fe_lagrange.C
Go to the documentation of this file.
00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 00003 00004 // This library is free software; you can redistribute it and/or 00005 // modify it under the terms of the GNU Lesser General Public 00006 // License as published by the Free Software Foundation; either 00007 // version 2.1 of the License, or (at your option) any later version. 00008 00009 // This library is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 // Lesser General Public License for more details. 00013 00014 // You should have received a copy of the GNU Lesser General Public 00015 // License along with this library; if not, write to the Free Software 00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00017 00018 00019 00020 // Local includes 00021 #include "libmesh/dof_map.h" 00022 #include "libmesh/fe.h" 00023 #include "libmesh/fe_interface.h" 00024 #include "libmesh/elem.h" 00025 #include "libmesh/threads.h" 00026 #include "libmesh/string_to_enum.h" 00027 00028 namespace libMesh 00029 { 00030 00031 // ------------------------------------------------------------ 00032 // Lagrange-specific implementations 00033 00034 00035 // Anonymous namespace for local helper functions 00036 namespace { 00037 void lagrange_nodal_soln(const Elem* elem, 00038 const Order order, 00039 const std::vector<Number>& elem_soln, 00040 std::vector<Number>& nodal_soln) 00041 { 00042 const unsigned int n_nodes = elem->n_nodes(); 00043 const ElemType type = elem->type(); 00044 00045 const Order totalorder = static_cast<Order>(order+elem->p_level()); 00046 00047 nodal_soln.resize(n_nodes); 00048 00049 00050 00051 switch (totalorder) 00052 { 00053 // linear Lagrange shape functions 00054 case FIRST: 00055 { 00056 switch (type) 00057 { 00058 case EDGE3: 00059 { 00060 libmesh_assert_equal_to (elem_soln.size(), 2); 00061 libmesh_assert_equal_to (nodal_soln.size(), 3); 00062 00063 nodal_soln[0] = elem_soln[0]; 00064 nodal_soln[1] = elem_soln[1]; 00065 nodal_soln[2] = .5*(elem_soln[0] + elem_soln[1]); 00066 00067 return; 00068 } 00069 00070 case EDGE4: 00071 { 00072 libmesh_assert_equal_to (elem_soln.size(), 2); 00073 libmesh_assert_equal_to (nodal_soln.size(), 4); 00074 00075 nodal_soln[0] = elem_soln[0]; 00076 nodal_soln[1] = elem_soln[1]; 00077 nodal_soln[2] = (2.*elem_soln[0] + elem_soln[1])/3.; 00078 nodal_soln[3] = (elem_soln[0] + 2.*elem_soln[1])/3.; 00079 00080 return; 00081 } 00082 00083 00084 case TRI6: 00085 { 00086 libmesh_assert_equal_to (elem_soln.size(), 3); 00087 libmesh_assert_equal_to (nodal_soln.size(), 6); 00088 00089 nodal_soln[0] = elem_soln[0]; 00090 nodal_soln[1] = elem_soln[1]; 00091 nodal_soln[2] = elem_soln[2]; 00092 nodal_soln[3] = .5*(elem_soln[0] + elem_soln[1]); 00093 nodal_soln[4] = .5*(elem_soln[1] + elem_soln[2]); 00094 nodal_soln[5] = .5*(elem_soln[2] + elem_soln[0]); 00095 00096 return; 00097 } 00098 00099 00100 case QUAD8: 00101 case QUAD9: 00102 { 00103 libmesh_assert_equal_to (elem_soln.size(), 4); 00104 00105 if (type == QUAD8) 00106 libmesh_assert_equal_to (nodal_soln.size(), 8); 00107 else 00108 libmesh_assert_equal_to (nodal_soln.size(), 9); 00109 00110 00111 nodal_soln[0] = elem_soln[0]; 00112 nodal_soln[1] = elem_soln[1]; 00113 nodal_soln[2] = elem_soln[2]; 00114 nodal_soln[3] = elem_soln[3]; 00115 nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]); 00116 nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]); 00117 nodal_soln[6] = .5*(elem_soln[2] + elem_soln[3]); 00118 nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]); 00119 00120 if (type == QUAD9) 00121 nodal_soln[8] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]); 00122 00123 return; 00124 } 00125 00126 00127 case TET10: 00128 { 00129 libmesh_assert_equal_to (elem_soln.size(), 4); 00130 libmesh_assert_equal_to (nodal_soln.size(), 10); 00131 00132 nodal_soln[0] = elem_soln[0]; 00133 nodal_soln[1] = elem_soln[1]; 00134 nodal_soln[2] = elem_soln[2]; 00135 nodal_soln[3] = elem_soln[3]; 00136 nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]); 00137 nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]); 00138 nodal_soln[6] = .5*(elem_soln[2] + elem_soln[0]); 00139 nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]); 00140 nodal_soln[8] = .5*(elem_soln[3] + elem_soln[1]); 00141 nodal_soln[9] = .5*(elem_soln[3] + elem_soln[2]); 00142 00143 return; 00144 } 00145 00146 00147 case HEX20: 00148 case HEX27: 00149 { 00150 libmesh_assert_equal_to (elem_soln.size(), 8); 00151 00152 if (type == HEX20) 00153 libmesh_assert_equal_to (nodal_soln.size(), 20); 00154 else 00155 libmesh_assert_equal_to (nodal_soln.size(), 27); 00156 00157 nodal_soln[0] = elem_soln[0]; 00158 nodal_soln[1] = elem_soln[1]; 00159 nodal_soln[2] = elem_soln[2]; 00160 nodal_soln[3] = elem_soln[3]; 00161 nodal_soln[4] = elem_soln[4]; 00162 nodal_soln[5] = elem_soln[5]; 00163 nodal_soln[6] = elem_soln[6]; 00164 nodal_soln[7] = elem_soln[7]; 00165 nodal_soln[8] = .5*(elem_soln[0] + elem_soln[1]); 00166 nodal_soln[9] = .5*(elem_soln[1] + elem_soln[2]); 00167 nodal_soln[10] = .5*(elem_soln[2] + elem_soln[3]); 00168 nodal_soln[11] = .5*(elem_soln[3] + elem_soln[0]); 00169 nodal_soln[12] = .5*(elem_soln[0] + elem_soln[4]); 00170 nodal_soln[13] = .5*(elem_soln[1] + elem_soln[5]); 00171 nodal_soln[14] = .5*(elem_soln[2] + elem_soln[6]); 00172 nodal_soln[15] = .5*(elem_soln[3] + elem_soln[7]); 00173 nodal_soln[16] = .5*(elem_soln[4] + elem_soln[5]); 00174 nodal_soln[17] = .5*(elem_soln[5] + elem_soln[6]); 00175 nodal_soln[18] = .5*(elem_soln[6] + elem_soln[7]); 00176 nodal_soln[19] = .5*(elem_soln[4] + elem_soln[7]); 00177 00178 if (type == HEX27) 00179 { 00180 nodal_soln[20] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]); 00181 nodal_soln[21] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[5]); 00182 nodal_soln[22] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[6]); 00183 nodal_soln[23] = .25*(elem_soln[2] + elem_soln[3] + elem_soln[6] + elem_soln[7]); 00184 nodal_soln[24] = .25*(elem_soln[3] + elem_soln[0] + elem_soln[7] + elem_soln[4]); 00185 nodal_soln[25] = .25*(elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]); 00186 00187 nodal_soln[26] = .125*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3] + 00188 elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]); 00189 } 00190 00191 return; 00192 } 00193 00194 00195 case PRISM15: 00196 case PRISM18: 00197 { 00198 libmesh_assert_equal_to (elem_soln.size(), 6); 00199 00200 if (type == PRISM15) 00201 libmesh_assert_equal_to (nodal_soln.size(), 15); 00202 else 00203 libmesh_assert_equal_to (nodal_soln.size(), 18); 00204 00205 nodal_soln[0] = elem_soln[0]; 00206 nodal_soln[1] = elem_soln[1]; 00207 nodal_soln[2] = elem_soln[2]; 00208 nodal_soln[3] = elem_soln[3]; 00209 nodal_soln[4] = elem_soln[4]; 00210 nodal_soln[5] = elem_soln[5]; 00211 nodal_soln[6] = .5*(elem_soln[0] + elem_soln[1]); 00212 nodal_soln[7] = .5*(elem_soln[1] + elem_soln[2]); 00213 nodal_soln[8] = .5*(elem_soln[0] + elem_soln[2]); 00214 nodal_soln[9] = .5*(elem_soln[0] + elem_soln[3]); 00215 nodal_soln[10] = .5*(elem_soln[1] + elem_soln[4]); 00216 nodal_soln[11] = .5*(elem_soln[2] + elem_soln[5]); 00217 nodal_soln[12] = .5*(elem_soln[3] + elem_soln[4]); 00218 nodal_soln[13] = .5*(elem_soln[4] + elem_soln[5]); 00219 nodal_soln[14] = .5*(elem_soln[3] + elem_soln[5]); 00220 00221 if (type == PRISM18) 00222 { 00223 nodal_soln[15] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[3]); 00224 nodal_soln[16] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[4]); 00225 nodal_soln[17] = .25*(elem_soln[2] + elem_soln[0] + elem_soln[3] + elem_soln[5]); 00226 } 00227 00228 return; 00229 } 00230 00231 case PYRAMID13: 00232 { 00233 libmesh_assert_equal_to (elem_soln.size(), 5); 00234 libmesh_assert_equal_to (nodal_soln.size(), 13); 00235 00236 nodal_soln[0] = elem_soln[0]; 00237 nodal_soln[1] = elem_soln[1]; 00238 nodal_soln[2] = elem_soln[2]; 00239 nodal_soln[3] = elem_soln[3]; 00240 nodal_soln[4] = elem_soln[4]; 00241 nodal_soln[5] = .5*(elem_soln[0] + elem_soln[1]); 00242 nodal_soln[6] = .5*(elem_soln[1] + elem_soln[2]); 00243 nodal_soln[7] = .5*(elem_soln[2] + elem_soln[3]); 00244 nodal_soln[8] = .5*(elem_soln[3] + elem_soln[0]); 00245 nodal_soln[9] = .5*(elem_soln[0] + elem_soln[4]); 00246 nodal_soln[10] = .5*(elem_soln[1] + elem_soln[4]); 00247 nodal_soln[11] = .5*(elem_soln[2] + elem_soln[4]); 00248 nodal_soln[12] = .5*(elem_soln[3] + elem_soln[4]); 00249 00250 return; 00251 } 00252 00253 case PYRAMID14: 00254 { 00255 libmesh_assert_equal_to (elem_soln.size(), 5); 00256 libmesh_assert_equal_to (nodal_soln.size(), 14); 00257 00258 nodal_soln[0] = elem_soln[0]; 00259 nodal_soln[1] = elem_soln[1]; 00260 nodal_soln[2] = elem_soln[2]; 00261 nodal_soln[3] = elem_soln[3]; 00262 nodal_soln[4] = elem_soln[4]; 00263 nodal_soln[5] = .5*(elem_soln[0] + elem_soln[1]); 00264 nodal_soln[6] = .5*(elem_soln[1] + elem_soln[2]); 00265 nodal_soln[7] = .5*(elem_soln[2] + elem_soln[3]); 00266 nodal_soln[8] = .5*(elem_soln[3] + elem_soln[0]); 00267 nodal_soln[9] = .5*(elem_soln[0] + elem_soln[4]); 00268 nodal_soln[10] = .5*(elem_soln[1] + elem_soln[4]); 00269 nodal_soln[11] = .5*(elem_soln[2] + elem_soln[4]); 00270 nodal_soln[12] = .5*(elem_soln[3] + elem_soln[4]); 00271 nodal_soln[13] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]); 00272 00273 return; 00274 } 00275 00276 default: 00277 { 00278 // By default the element solution _is_ nodal, 00279 // so just copy it. 00280 nodal_soln = elem_soln; 00281 00282 return; 00283 } 00284 } 00285 } 00286 00287 case SECOND: 00288 { 00289 switch (type) 00290 { 00291 case EDGE4: 00292 { 00293 libmesh_assert_equal_to (elem_soln.size(), 3); 00294 libmesh_assert_equal_to (nodal_soln.size(), 4); 00295 00296 // Project quadratic solution onto cubic element nodes 00297 nodal_soln[0] = elem_soln[0]; 00298 nodal_soln[1] = elem_soln[1]; 00299 nodal_soln[2] = (2.*elem_soln[0] - elem_soln[1] + 00300 8.*elem_soln[2])/9.; 00301 nodal_soln[3] = (-elem_soln[0] + 2.*elem_soln[1] + 00302 8.*elem_soln[2])/9.; 00303 return; 00304 } 00305 00306 default: 00307 { 00308 // By default the element solution _is_ nodal, 00309 // so just copy it. 00310 nodal_soln = elem_soln; 00311 00312 return; 00313 } 00314 } 00315 } 00316 00317 00318 00319 00320 default: 00321 { 00322 // By default the element solution _is_ nodal, 00323 // so just copy it. 00324 nodal_soln = elem_soln; 00325 00326 return; 00327 } 00328 } 00329 } 00330 00331 00332 00333 unsigned int lagrange_n_dofs(const ElemType t, const Order o) 00334 { 00335 switch (o) 00336 { 00337 00338 // linear Lagrange shape functions 00339 case FIRST: 00340 { 00341 switch (t) 00342 { 00343 case NODEELEM: 00344 return 1; 00345 00346 case EDGE2: 00347 case EDGE3: 00348 case EDGE4: 00349 return 2; 00350 00351 case TRI3: 00352 case TRI3SUBDIVISION: 00353 case TRI6: 00354 return 3; 00355 00356 case QUAD4: 00357 case QUAD8: 00358 case QUAD9: 00359 return 4; 00360 00361 case TET4: 00362 case TET10: 00363 return 4; 00364 00365 case HEX8: 00366 case HEX20: 00367 case HEX27: 00368 return 8; 00369 00370 case PRISM6: 00371 case PRISM15: 00372 case PRISM18: 00373 return 6; 00374 00375 case PYRAMID5: 00376 case PYRAMID13: 00377 case PYRAMID14: 00378 return 5; 00379 00380 case INVALID_ELEM: 00381 return 0; 00382 00383 default: 00384 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!"); 00385 } 00386 } 00387 00388 00389 // quadratic Lagrange shape functions 00390 case SECOND: 00391 { 00392 switch (t) 00393 { 00394 case NODEELEM: 00395 return 1; 00396 00397 case EDGE3: 00398 return 3; 00399 00400 case TRI6: 00401 return 6; 00402 00403 case QUAD8: 00404 return 8; 00405 00406 case QUAD9: 00407 return 9; 00408 00409 case TET10: 00410 return 10; 00411 00412 case HEX20: 00413 return 20; 00414 00415 case HEX27: 00416 return 27; 00417 00418 case PRISM15: 00419 return 15; 00420 00421 case PRISM18: 00422 return 18; 00423 00424 case PYRAMID13: 00425 return 13; 00426 00427 case PYRAMID14: 00428 return 14; 00429 00430 case INVALID_ELEM: 00431 return 0; 00432 00433 default: 00434 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!"); 00435 } 00436 } 00437 00438 case THIRD: 00439 { 00440 switch (t) 00441 { 00442 case NODEELEM: 00443 return 1; 00444 00445 case EDGE4: 00446 return 4; 00447 00448 case INVALID_ELEM: 00449 return 0; 00450 00451 default: 00452 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!"); 00453 } 00454 } 00455 00456 default: 00457 libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for LAGRANGE FE family!"); 00458 } 00459 00460 libmesh_error_msg("We'll never get here!"); 00461 return 0; 00462 } 00463 00464 00465 00466 00467 unsigned int lagrange_n_dofs_at_node(const ElemType t, 00468 const Order o, 00469 const unsigned int n) 00470 { 00471 switch (o) 00472 { 00473 00474 // linear Lagrange shape functions 00475 case FIRST: 00476 { 00477 switch (t) 00478 { 00479 case NODEELEM: 00480 return 1; 00481 00482 case EDGE2: 00483 case EDGE3: 00484 case EDGE4: 00485 { 00486 switch (n) 00487 { 00488 case 0: 00489 case 1: 00490 return 1; 00491 00492 default: 00493 return 0; 00494 } 00495 } 00496 00497 case TRI3: 00498 case TRI3SUBDIVISION: 00499 case TRI6: 00500 { 00501 switch (n) 00502 { 00503 case 0: 00504 case 1: 00505 case 2: 00506 return 1; 00507 00508 default: 00509 return 0; 00510 } 00511 } 00512 00513 case QUAD4: 00514 case QUAD8: 00515 case QUAD9: 00516 { 00517 switch (n) 00518 { 00519 case 0: 00520 case 1: 00521 case 2: 00522 case 3: 00523 return 1; 00524 00525 default: 00526 return 0; 00527 } 00528 } 00529 00530 00531 case TET4: 00532 case TET10: 00533 { 00534 switch (n) 00535 { 00536 case 0: 00537 case 1: 00538 case 2: 00539 case 3: 00540 return 1; 00541 00542 default: 00543 return 0; 00544 } 00545 } 00546 00547 case HEX8: 00548 case HEX20: 00549 case HEX27: 00550 { 00551 switch (n) 00552 { 00553 case 0: 00554 case 1: 00555 case 2: 00556 case 3: 00557 case 4: 00558 case 5: 00559 case 6: 00560 case 7: 00561 return 1; 00562 00563 default: 00564 return 0; 00565 } 00566 } 00567 00568 case PRISM6: 00569 case PRISM15: 00570 case PRISM18: 00571 { 00572 switch (n) 00573 { 00574 case 0: 00575 case 1: 00576 case 2: 00577 case 3: 00578 case 4: 00579 case 5: 00580 return 1; 00581 00582 default: 00583 return 0; 00584 } 00585 } 00586 00587 case PYRAMID5: 00588 case PYRAMID13: 00589 case PYRAMID14: 00590 { 00591 switch (n) 00592 { 00593 case 0: 00594 case 1: 00595 case 2: 00596 case 3: 00597 case 4: 00598 return 1; 00599 00600 default: 00601 return 0; 00602 } 00603 } 00604 00605 case INVALID_ELEM: 00606 return 0; 00607 00608 default: 00609 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!"); 00610 } 00611 } 00612 00613 // quadratic Lagrange shape functions 00614 case SECOND: 00615 { 00616 switch (t) 00617 { 00618 // quadratic lagrange has one dof at each node 00619 case NODEELEM: 00620 case EDGE3: 00621 case TRI6: 00622 case QUAD8: 00623 case QUAD9: 00624 case TET10: 00625 case HEX20: 00626 case HEX27: 00627 case PRISM15: 00628 case PRISM18: 00629 case PYRAMID13: 00630 case PYRAMID14: 00631 return 1; 00632 00633 case INVALID_ELEM: 00634 return 0; 00635 00636 default: 00637 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!"); 00638 } 00639 } 00640 00641 case THIRD: 00642 { 00643 switch (t) 00644 { 00645 case NODEELEM: 00646 case EDGE4: 00647 return 1; 00648 00649 case INVALID_ELEM: 00650 return 0; 00651 00652 default: 00653 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!"); 00654 } 00655 } 00656 00657 default: 00658 libmesh_error_msg("Unsupported order: " << o ); 00659 } 00660 00661 libmesh_error_msg("We'll never get here!"); 00662 return 0; 00663 00664 } 00665 00666 00667 00668 #ifdef LIBMESH_ENABLE_AMR 00669 void lagrange_compute_constraints (DofConstraints &constraints, 00670 DofMap &dof_map, 00671 const unsigned int variable_number, 00672 const Elem* elem, 00673 const unsigned Dim) 00674 { 00675 // Only constrain elements in 2,3D. 00676 if (Dim == 1) 00677 return; 00678 00679 libmesh_assert(elem); 00680 00681 // Only constrain active and ancestor elements 00682 if (elem->subactive()) 00683 return; 00684 00685 FEType fe_type = dof_map.variable_type(variable_number); 00686 fe_type.order = static_cast<Order>(fe_type.order + elem->p_level()); 00687 00688 std::vector<dof_id_type> my_dof_indices, parent_dof_indices; 00689 00690 // Look at the element faces. Check to see if we need to 00691 // build constraints. 00692 for (unsigned int s=0; s<elem->n_sides(); s++) 00693 if (elem->neighbor(s) != NULL) 00694 if (elem->neighbor(s)->level() < elem->level()) // constrain dofs shared between 00695 { // this element and ones coarser 00696 // than this element. 00697 // Get pointers to the elements of interest and its parent. 00698 const Elem* parent = elem->parent(); 00699 00700 // This can't happen... Only level-0 elements have NULL 00701 // parents, and no level-0 elements can be at a higher 00702 // level than their neighbors! 00703 libmesh_assert(parent); 00704 00705 const UniquePtr<Elem> my_side (elem->build_side(s)); 00706 const UniquePtr<Elem> parent_side (parent->build_side(s)); 00707 00708 // This function gets called element-by-element, so there 00709 // will be a lot of memory allocation going on. We can 00710 // at least minimize this for the case of the dof indices 00711 // by efficiently preallocating the requisite storage. 00712 my_dof_indices.reserve (my_side->n_nodes()); 00713 parent_dof_indices.reserve (parent_side->n_nodes()); 00714 00715 dof_map.dof_indices (my_side.get(), my_dof_indices, 00716 variable_number); 00717 dof_map.dof_indices (parent_side.get(), parent_dof_indices, 00718 variable_number); 00719 00720 for (unsigned int my_dof=0; 00721 my_dof<FEInterface::n_dofs(Dim-1, fe_type, my_side->type()); 00722 my_dof++) 00723 { 00724 libmesh_assert_less (my_dof, my_side->n_nodes()); 00725 00726 // My global dof index. 00727 const dof_id_type my_dof_g = my_dof_indices[my_dof]; 00728 00729 // Hunt for "constraining against myself" cases before 00730 // we bother creating a constraint row 00731 bool self_constraint = false; 00732 for (unsigned int their_dof=0; 00733 their_dof<FEInterface::n_dofs(Dim-1, fe_type, parent_side->type()); 00734 their_dof++) 00735 { 00736 libmesh_assert_less (their_dof, parent_side->n_nodes()); 00737 00738 // Their global dof index. 00739 const dof_id_type their_dof_g = 00740 parent_dof_indices[their_dof]; 00741 00742 if (their_dof_g == my_dof_g) 00743 { 00744 self_constraint = true; 00745 break; 00746 } 00747 } 00748 00749 if (self_constraint) 00750 continue; 00751 00752 DofConstraintRow* constraint_row; 00753 00754 // we may be running constraint methods concurrently 00755 // on multiple threads, so we need a lock to 00756 // ensure that this constraint is "ours" 00757 { 00758 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx); 00759 00760 if (dof_map.is_constrained_dof(my_dof_g)) 00761 continue; 00762 00763 constraint_row = &(constraints[my_dof_g]); 00764 libmesh_assert(constraint_row->empty()); 00765 } 00766 00767 // The support point of the DOF 00768 const Point& support_point = my_side->point(my_dof); 00769 00770 // Figure out where my node lies on their reference element. 00771 const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type, 00772 parent_side.get(), 00773 support_point); 00774 00775 // Compute the parent's side shape function values. 00776 for (unsigned int their_dof=0; 00777 their_dof<FEInterface::n_dofs(Dim-1, fe_type, parent_side->type()); 00778 their_dof++) 00779 { 00780 libmesh_assert_less (their_dof, parent_side->n_nodes()); 00781 00782 // Their global dof index. 00783 const dof_id_type their_dof_g = 00784 parent_dof_indices[their_dof]; 00785 00786 const Real their_dof_value = FEInterface::shape(Dim-1, 00787 fe_type, 00788 parent_side->type(), 00789 their_dof, 00790 mapped_point); 00791 00792 // Only add non-zero and non-identity values 00793 // for Lagrange basis functions. 00794 if ((std::abs(their_dof_value) > 1.e-5) && 00795 (std::abs(their_dof_value) < .999)) 00796 { 00797 constraint_row->insert(std::make_pair (their_dof_g, 00798 their_dof_value)); 00799 } 00800 #ifdef DEBUG 00801 // Protect for the case u_i = 0.999 u_j, 00802 // in which case i better equal j. 00803 else if (their_dof_value >= .999) 00804 libmesh_assert_equal_to (my_dof_g, their_dof_g); 00805 #endif 00806 } 00807 } 00808 } 00809 } // lagrange_compute_constrants() 00810 #endif // #ifdef LIBMESH_ENABLE_AMR 00811 00812 } // anonymous namespace 00813 00814 00815 // Do full-specialization for every dimension, instead 00816 // of explicit instantiation at the end of this file. 00817 // This could be macro-ified so that it fits on one line... 00818 template <> 00819 void FE<0,LAGRANGE>::nodal_soln(const Elem* elem, 00820 const Order order, 00821 const std::vector<Number>& elem_soln, 00822 std::vector<Number>& nodal_soln) 00823 { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); } 00824 00825 template <> 00826 void FE<1,LAGRANGE>::nodal_soln(const Elem* elem, 00827 const Order order, 00828 const std::vector<Number>& elem_soln, 00829 std::vector<Number>& nodal_soln) 00830 { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); } 00831 00832 template <> 00833 void FE<2,LAGRANGE>::nodal_soln(const Elem* elem, 00834 const Order order, 00835 const std::vector<Number>& elem_soln, 00836 std::vector<Number>& nodal_soln) 00837 { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); } 00838 00839 template <> 00840 void FE<3,LAGRANGE>::nodal_soln(const Elem* elem, 00841 const Order order, 00842 const std::vector<Number>& elem_soln, 00843 std::vector<Number>& nodal_soln) 00844 { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); } 00845 00846 00847 // Do full-specialization for every dimension, instead 00848 // of explicit instantiation at the end of this function. 00849 // This could be macro-ified. 00850 template <> unsigned int FE<0,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); } 00851 template <> unsigned int FE<1,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); } 00852 template <> unsigned int FE<2,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); } 00853 template <> unsigned int FE<3,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); } 00854 00855 00856 // Do full-specialization for every dimension, instead 00857 // of explicit instantiation at the end of this function. 00858 template <> unsigned int FE<0,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); } 00859 template <> unsigned int FE<1,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); } 00860 template <> unsigned int FE<2,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); } 00861 template <> unsigned int FE<3,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); } 00862 00863 00864 // Lagrange elements have no dofs per element 00865 // (just at the nodes) 00866 template <> unsigned int FE<0,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; } 00867 template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; } 00868 template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; } 00869 template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; } 00870 00871 // Lagrange FEMs are always C^0 continuous 00872 template <> FEContinuity FE<0,LAGRANGE>::get_continuity() const { return C_ZERO; } 00873 template <> FEContinuity FE<1,LAGRANGE>::get_continuity() const { return C_ZERO; } 00874 template <> FEContinuity FE<2,LAGRANGE>::get_continuity() const { return C_ZERO; } 00875 template <> FEContinuity FE<3,LAGRANGE>::get_continuity() const { return C_ZERO; } 00876 00877 // Lagrange FEMs are not hierarchic 00878 template <> bool FE<0,LAGRANGE>::is_hierarchic() const { return false; } 00879 template <> bool FE<1,LAGRANGE>::is_hierarchic() const { return false; } 00880 template <> bool FE<2,LAGRANGE>::is_hierarchic() const { return false; } 00881 template <> bool FE<3,LAGRANGE>::is_hierarchic() const { return false; } 00882 00883 // Lagrange FEM shapes do not need reinit (is this always true?) 00884 template <> bool FE<0,LAGRANGE>::shapes_need_reinit() const { return false; } 00885 template <> bool FE<1,LAGRANGE>::shapes_need_reinit() const { return false; } 00886 template <> bool FE<2,LAGRANGE>::shapes_need_reinit() const { return false; } 00887 template <> bool FE<3,LAGRANGE>::shapes_need_reinit() const { return false; } 00888 00889 // Methods for computing Lagrange constraints. Note: we pass the 00890 // dimension as the last argument to the anonymous helper function. 00891 // Also note: we only need instantiations of this function for 00892 // Dim==2 and 3. 00893 #ifdef LIBMESH_ENABLE_AMR 00894 template <> 00895 void FE<2,LAGRANGE>::compute_constraints (DofConstraints &constraints, 00896 DofMap &dof_map, 00897 const unsigned int variable_number, 00898 const Elem* elem) 00899 { lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/2); } 00900 00901 template <> 00902 void FE<3,LAGRANGE>::compute_constraints (DofConstraints &constraints, 00903 DofMap &dof_map, 00904 const unsigned int variable_number, 00905 const Elem* elem) 00906 { lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/3); } 00907 #endif // LIBMESH_ENABLE_AMR 00908 00909 } // namespace libMesh
generated by
