fe_lagrange_shape_2D.C
Go to the documentation of this file.00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2012 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 // C++ inlcludes 00020 00021 // Local includes 00022 #include "libmesh/fe.h" 00023 #include "libmesh/elem.h" 00024 00025 namespace libMesh 00026 { 00027 00028 00029 00030 00031 template <> 00032 Real FE<2,LAGRANGE>::shape(const ElemType type, 00033 const Order order, 00034 const unsigned int i, 00035 const Point& p) 00036 { 00037 #if LIBMESH_DIM > 1 00038 00039 switch (order) 00040 { 00041 // linear Lagrange shape functions 00042 case FIRST: 00043 { 00044 switch (type) 00045 { 00046 case QUAD4: 00047 case QUAD8: 00048 case QUAD9: 00049 { 00050 // Compute quad shape functions as a tensor-product 00051 const Real xi = p(0); 00052 const Real eta = p(1); 00053 00054 libmesh_assert_less (i, 4); 00055 00056 // 0 1 2 3 00057 static const unsigned int i0[] = {0, 1, 1, 0}; 00058 static const unsigned int i1[] = {0, 0, 1, 1}; 00059 00060 return (FE<1,LAGRANGE>::shape(EDGE2, FIRST, i0[i], xi)* 00061 FE<1,LAGRANGE>::shape(EDGE2, FIRST, i1[i], eta)); 00062 } 00063 00064 case TRI3: 00065 case TRI6: 00066 { 00067 const Real zeta1 = p(0); 00068 const Real zeta2 = p(1); 00069 const Real zeta0 = 1. - zeta1 - zeta2; 00070 00071 libmesh_assert_less (i, 3); 00072 00073 switch(i) 00074 { 00075 case 0: 00076 return zeta0; 00077 00078 case 1: 00079 return zeta1; 00080 00081 case 2: 00082 return zeta2; 00083 00084 default: 00085 libmesh_error(); 00086 00087 } 00088 } 00089 00090 default: 00091 { 00092 libMesh::err << "ERROR: Unsupported 2D element type!: " << type 00093 << std::endl; 00094 libmesh_error(); 00095 } 00096 } 00097 } 00098 00099 00100 // quadratic Lagrange shape functions 00101 case SECOND: 00102 { 00103 switch (type) 00104 { 00105 case QUAD8: 00106 { 00107 const Real xi = p(0); 00108 const Real eta = p(1); 00109 00110 libmesh_assert_less (i, 8); 00111 00112 switch (i) 00113 { 00114 case 0: 00115 return .25*(1. - xi)*(1. - eta)*(-1. - xi - eta); 00116 00117 case 1: 00118 return .25*(1. + xi)*(1. - eta)*(-1. + xi - eta); 00119 00120 case 2: 00121 return .25*(1. + xi)*(1. + eta)*(-1. + xi + eta); 00122 00123 case 3: 00124 return .25*(1. - xi)*(1. + eta)*(-1. - xi + eta); 00125 00126 case 4: 00127 return .5*(1. - xi*xi)*(1. - eta); 00128 00129 case 5: 00130 return .5*(1. + xi)*(1. - eta*eta); 00131 00132 case 6: 00133 return .5*(1. - xi*xi)*(1. + eta); 00134 00135 case 7: 00136 return .5*(1. - xi)*(1. - eta*eta); 00137 00138 default: 00139 libmesh_error(); 00140 } 00141 } 00142 00143 case QUAD9: 00144 { 00145 // Compute quad shape functions as a tensor-product 00146 const Real xi = p(0); 00147 const Real eta = p(1); 00148 00149 libmesh_assert_less (i, 9); 00150 00151 // 0 1 2 3 4 5 6 7 8 00152 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2}; 00153 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2}; 00154 00155 return (FE<1,LAGRANGE>::shape(EDGE3, SECOND, i0[i], xi)* 00156 FE<1,LAGRANGE>::shape(EDGE3, SECOND, i1[i], eta)); 00157 } 00158 00159 case TRI6: 00160 { 00161 const Real zeta1 = p(0); 00162 const Real zeta2 = p(1); 00163 const Real zeta0 = 1. - zeta1 - zeta2; 00164 00165 libmesh_assert_less (i, 6); 00166 00167 switch(i) 00168 { 00169 case 0: 00170 return 2.*zeta0*(zeta0-0.5); 00171 00172 case 1: 00173 return 2.*zeta1*(zeta1-0.5); 00174 00175 case 2: 00176 return 2.*zeta2*(zeta2-0.5); 00177 00178 case 3: 00179 return 4.*zeta0*zeta1; 00180 00181 case 4: 00182 return 4.*zeta1*zeta2; 00183 00184 case 5: 00185 return 4.*zeta2*zeta0; 00186 00187 default: 00188 libmesh_error(); 00189 00190 } 00191 } 00192 00193 default: 00194 { 00195 libMesh::err << "ERROR: Unsupported 2D element type!: " << type 00196 << std::endl; 00197 libmesh_error(); 00198 } 00199 } 00200 } 00201 00202 00203 00204 // unsupported order 00205 default: 00206 { 00207 libMesh::err << "ERROR: Unsupported 2D FE order!: " << order 00208 << std::endl; 00209 libmesh_error(); 00210 } 00211 } 00212 00213 libmesh_error(); 00214 #endif // LIBMESH_DIM > 1 00215 return 0.; 00216 } 00217 00218 00219 00220 template <> 00221 Real FE<2,LAGRANGE>::shape(const Elem* elem, 00222 const Order order, 00223 const unsigned int i, 00224 const Point& p) 00225 { 00226 libmesh_assert(elem); 00227 00228 // call the orientation-independent shape functions 00229 return FE<2,LAGRANGE>::shape(elem->type(), static_cast<Order>(order + elem->p_level()), i, p); 00230 } 00231 00232 00233 00234 template <> 00235 Real FE<2,LAGRANGE>::shape_deriv(const ElemType type, 00236 const Order order, 00237 const unsigned int i, 00238 const unsigned int j, 00239 const Point& p) 00240 { 00241 #if LIBMESH_DIM > 1 00242 00243 00244 libmesh_assert_less (j, 2); 00245 00246 switch (order) 00247 { 00248 // linear Lagrange shape functions 00249 case FIRST: 00250 { 00251 switch (type) 00252 { 00253 case QUAD4: 00254 case QUAD8: 00255 case QUAD9: 00256 { 00257 // Compute quad shape functions as a tensor-product 00258 const Real xi = p(0); 00259 const Real eta = p(1); 00260 00261 libmesh_assert_less (i, 4); 00262 00263 // 0 1 2 3 00264 static const unsigned int i0[] = {0, 1, 1, 0}; 00265 static const unsigned int i1[] = {0, 0, 1, 1}; 00266 00267 switch (j) 00268 { 00269 // d()/dxi 00270 case 0: 00271 return (FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, xi)* 00272 FE<1,LAGRANGE>::shape (EDGE2, FIRST, i1[i], eta)); 00273 00274 // d()/deta 00275 case 1: 00276 return (FE<1,LAGRANGE>::shape (EDGE2, FIRST, i0[i], xi)* 00277 FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i1[i], 0, eta)); 00278 00279 default: 00280 libmesh_error(); 00281 } 00282 } 00283 00284 case TRI3: 00285 case TRI6: 00286 { 00287 libmesh_assert_less (i, 3); 00288 00289 const Real dzeta0dxi = -1.; 00290 const Real dzeta1dxi = 1.; 00291 const Real dzeta2dxi = 0.; 00292 00293 const Real dzeta0deta = -1.; 00294 const Real dzeta1deta = 0.; 00295 const Real dzeta2deta = 1.; 00296 00297 switch (j) 00298 { 00299 // d()/dxi 00300 case 0: 00301 { 00302 switch(i) 00303 { 00304 case 0: 00305 return dzeta0dxi; 00306 00307 case 1: 00308 return dzeta1dxi; 00309 00310 case 2: 00311 return dzeta2dxi; 00312 00313 default: 00314 libmesh_error(); 00315 } 00316 } 00317 // d()/deta 00318 case 1: 00319 { 00320 switch(i) 00321 { 00322 case 0: 00323 return dzeta0deta; 00324 00325 case 1: 00326 return dzeta1deta; 00327 00328 case 2: 00329 return dzeta2deta; 00330 00331 default: 00332 libmesh_error(); 00333 00334 } 00335 } 00336 default: 00337 libmesh_error(); 00338 } 00339 } 00340 00341 default: 00342 { 00343 libMesh::err << "ERROR: Unsupported 2D element type!: " << type 00344 << std::endl; 00345 libmesh_error(); 00346 } 00347 } 00348 } 00349 00350 00351 // quadratic Lagrange shape functions 00352 case SECOND: 00353 { 00354 switch (type) 00355 { 00356 case QUAD8: 00357 { 00358 const Real xi = p(0); 00359 const Real eta = p(1); 00360 00361 libmesh_assert_less (i, 8); 00362 00363 switch (j) 00364 { 00365 // d/dxi 00366 case 0: 00367 switch (i) 00368 { 00369 case 0: 00370 return .25*(1. - eta)*((1. - xi)*(-1.) + 00371 (-1.)*(-1. - xi - eta)); 00372 00373 case 1: 00374 return .25*(1. - eta)*((1. + xi)*(1.) + 00375 (1.)*(-1. + xi - eta)); 00376 00377 case 2: 00378 return .25*(1. + eta)*((1. + xi)*(1.) + 00379 (1.)*(-1. + xi + eta)); 00380 00381 case 3: 00382 return .25*(1. + eta)*((1. - xi)*(-1.) + 00383 (-1.)*(-1. - xi + eta)); 00384 00385 case 4: 00386 return .5*(-2.*xi)*(1. - eta); 00387 00388 case 5: 00389 return .5*(1.)*(1. - eta*eta); 00390 00391 case 6: 00392 return .5*(-2.*xi)*(1. + eta); 00393 00394 case 7: 00395 return .5*(-1.)*(1. - eta*eta); 00396 00397 default: 00398 libmesh_error(); 00399 } 00400 00401 // d/deta 00402 case 1: 00403 switch (i) 00404 { 00405 case 0: 00406 return .25*(1. - xi)*((1. - eta)*(-1.) + 00407 (-1.)*(-1. - xi - eta)); 00408 00409 case 1: 00410 return .25*(1. + xi)*((1. - eta)*(-1.) + 00411 (-1.)*(-1. + xi - eta)); 00412 00413 case 2: 00414 return .25*(1. + xi)*((1. + eta)*(1.) + 00415 (1.)*(-1. + xi + eta)); 00416 00417 case 3: 00418 return .25*(1. - xi)*((1. + eta)*(1.) + 00419 (1.)*(-1. - xi + eta)); 00420 00421 case 4: 00422 return .5*(1. - xi*xi)*(-1.); 00423 00424 case 5: 00425 return .5*(1. + xi)*(-2.*eta); 00426 00427 case 6: 00428 return .5*(1. - xi*xi)*(1.); 00429 00430 case 7: 00431 return .5*(1. - xi)*(-2.*eta); 00432 00433 default: 00434 libmesh_error(); 00435 } 00436 00437 default: 00438 libmesh_error(); 00439 } 00440 } 00441 00442 case QUAD9: 00443 { 00444 // Compute quad shape functions as a tensor-product 00445 const Real xi = p(0); 00446 const Real eta = p(1); 00447 00448 libmesh_assert_less (i, 9); 00449 00450 // 0 1 2 3 4 5 6 7 8 00451 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2}; 00452 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2}; 00453 00454 switch (j) 00455 { 00456 // d()/dxi 00457 case 0: 00458 return (FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)* 00459 FE<1,LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)); 00460 00461 // d()/deta 00462 case 1: 00463 return (FE<1,LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)* 00464 FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)); 00465 00466 default: 00467 libmesh_error(); 00468 } 00469 } 00470 00471 case TRI6: 00472 { 00473 libmesh_assert_less (i, 6); 00474 00475 const Real zeta1 = p(0); 00476 const Real zeta2 = p(1); 00477 const Real zeta0 = 1. - zeta1 - zeta2; 00478 00479 const Real dzeta0dxi = -1.; 00480 const Real dzeta1dxi = 1.; 00481 const Real dzeta2dxi = 0.; 00482 00483 const Real dzeta0deta = -1.; 00484 const Real dzeta1deta = 0.; 00485 const Real dzeta2deta = 1.; 00486 00487 switch(j) 00488 { 00489 case 0: 00490 { 00491 switch(i) 00492 { 00493 case 0: 00494 return (4.*zeta0-1.)*dzeta0dxi; 00495 00496 case 1: 00497 return (4.*zeta1-1.)*dzeta1dxi; 00498 00499 case 2: 00500 return (4.*zeta2-1.)*dzeta2dxi; 00501 00502 case 3: 00503 return 4.*zeta1*dzeta0dxi + 4.*zeta0*dzeta1dxi; 00504 00505 case 4: 00506 return 4.*zeta2*dzeta1dxi + 4.*zeta1*dzeta2dxi; 00507 00508 case 5: 00509 return 4.*zeta2*dzeta0dxi + 4*zeta0*dzeta2dxi; 00510 00511 default: 00512 libmesh_error(); 00513 } 00514 } 00515 00516 case 1: 00517 { 00518 switch(i) 00519 { 00520 case 0: 00521 return (4.*zeta0-1.)*dzeta0deta; 00522 00523 case 1: 00524 return (4.*zeta1-1.)*dzeta1deta; 00525 00526 case 2: 00527 return (4.*zeta2-1.)*dzeta2deta; 00528 00529 case 3: 00530 return 4.*zeta1*dzeta0deta + 4.*zeta0*dzeta1deta; 00531 00532 case 4: 00533 return 4.*zeta2*dzeta1deta + 4.*zeta1*dzeta2deta; 00534 00535 case 5: 00536 return 4.*zeta2*dzeta0deta + 4*zeta0*dzeta2deta; 00537 00538 default: 00539 libmesh_error(); 00540 } 00541 } 00542 default: 00543 libmesh_error(); 00544 } 00545 } 00546 00547 default: 00548 { 00549 libMesh::err << "ERROR: Unsupported 2D element type!: " << type 00550 << std::endl; 00551 libmesh_error(); 00552 } 00553 } 00554 } 00555 00556 00557 00558 // unsupported order 00559 default: 00560 { 00561 libMesh::err << "ERROR: Unsupported 2D FE order!: " << order 00562 << std::endl; 00563 libmesh_error(); 00564 } 00565 } 00566 00567 00568 libmesh_error(); 00569 #endif // LIBMESH_DIM > 1 00570 return 0.; 00571 } 00572 00573 00574 00575 template <> 00576 Real FE<2,LAGRANGE>::shape_deriv(const Elem* elem, 00577 const Order order, 00578 const unsigned int i, 00579 const unsigned int j, 00580 const Point& p) 00581 { 00582 libmesh_assert(elem); 00583 00584 00585 // call the orientation-independent shape functions 00586 return FE<2,LAGRANGE>::shape_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p); 00587 } 00588 00589 00590 00591 00592 template <> 00593 Real FE<2,LAGRANGE>::shape_second_deriv(const ElemType type, 00594 const Order order, 00595 const unsigned int i, 00596 const unsigned int j, 00597 const Point& p) 00598 { 00599 #if LIBMESH_DIM > 1 00600 00601 // j = 0 ==> d^2 phi / dxi^2 00602 // j = 1 ==> d^2 phi / dxi deta 00603 // j = 2 ==> d^2 phi / deta^2 00604 libmesh_assert_less (j, 3); 00605 00606 switch (order) 00607 { 00608 // linear Lagrange shape functions 00609 case FIRST: 00610 { 00611 switch (type) 00612 { 00613 case QUAD4: 00614 case QUAD8: 00615 case QUAD9: 00616 { 00617 // Compute quad shape functions as a tensor-product 00618 const Real xi = p(0); 00619 const Real eta = p(1); 00620 00621 libmesh_assert_less (i, 4); 00622 00623 // 0 1 2 3 00624 static const unsigned int i0[] = {0, 1, 1, 0}; 00625 static const unsigned int i1[] = {0, 0, 1, 1}; 00626 00627 switch (j) 00628 { 00629 // d^2() / dxi^2 00630 case 0: 00631 return 0.; 00632 00633 // d^2() / dxi deta 00634 case 1: 00635 return (FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, xi)* 00636 FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i1[1], 0, eta)); 00637 00638 // d^2() / deta^2 00639 case 2: 00640 return 0.; 00641 00642 default: 00643 { 00644 libMesh::err << "ERROR: Invalid derivative requested! " 00645 << std::endl; 00646 libmesh_error(); 00647 } 00648 } 00649 } 00650 00651 case TRI3: 00652 case TRI6: 00653 { 00654 // All second derivatives for linear triangles are zero. 00655 return 0.; 00656 } 00657 00658 default: 00659 { 00660 libMesh::err << "ERROR: Unsupported 2D element type!: " << type 00661 << std::endl; 00662 libmesh_error(); 00663 } 00664 00665 } // end switch (type) 00666 } // end case FIRST 00667 00668 00669 // quadratic Lagrange shape functions 00670 case SECOND: 00671 { 00672 switch (type) 00673 { 00674 case QUAD8: 00675 { 00676 const Real xi = p(0); 00677 const Real eta = p(1); 00678 00679 libmesh_assert_less (j, 3); 00680 00681 switch (j) 00682 { 00683 // d^2() / dxi^2 00684 case 0: 00685 { 00686 switch (i) 00687 { 00688 case 0: 00689 case 1: 00690 return 0.5*(1.-eta); 00691 00692 case 2: 00693 case 3: 00694 return 0.5*(1.+eta); 00695 00696 case 4: 00697 return eta - 1.; 00698 00699 case 5: 00700 case 7: 00701 return 0.0; 00702 00703 case 6: 00704 return -1. - eta; 00705 00706 default: 00707 { 00708 libMesh::err << "Invalid shape function index requested!" 00709 << std::endl; 00710 libmesh_error(); 00711 } 00712 } 00713 } 00714 00715 // d^2() / dxi deta 00716 case 1: 00717 { 00718 switch (i) 00719 { 00720 case 0: 00721 return 0.25*( 1. - 2.*xi - 2.*eta); 00722 00723 case 1: 00724 return 0.25*(-1. - 2.*xi + 2.*eta); 00725 00726 case 2: 00727 return 0.25*( 1. + 2.*xi + 2.*eta); 00728 00729 case 3: 00730 return 0.25*(-1. + 2.*xi - 2.*eta); 00731 00732 case 4: 00733 return xi; 00734 00735 case 5: 00736 return -eta; 00737 00738 case 6: 00739 return -xi; 00740 00741 case 7: 00742 return eta; 00743 00744 default: 00745 { 00746 libMesh::err << "Invalid shape function index requested!" 00747 << std::endl; 00748 libmesh_error(); 00749 } 00750 } 00751 } 00752 00753 // d^2() / deta^2 00754 case 2: 00755 { 00756 switch (i) 00757 { 00758 case 0: 00759 case 3: 00760 return 0.5*(1.-xi); 00761 00762 case 1: 00763 case 2: 00764 return 0.5*(1.+xi); 00765 00766 case 4: 00767 case 6: 00768 return 0.0; 00769 00770 case 5: 00771 return -1.0 - xi; 00772 00773 case 7: 00774 return xi - 1.0; 00775 00776 default: 00777 { 00778 libMesh::err << "Invalid shape function index requested!" 00779 << std::endl; 00780 libmesh_error(); 00781 } 00782 } 00783 } 00784 00785 00786 default: 00787 { 00788 libMesh::err << "ERROR: Invalid derivative requested! " 00789 << std::endl; 00790 libmesh_error(); 00791 } 00792 } // end switch (j) 00793 } // end case QUAD8 00794 00795 case QUAD9: 00796 { 00797 // Compute QUAD9 second derivatives as tensor product 00798 const Real xi = p(0); 00799 const Real eta = p(1); 00800 00801 libmesh_assert_less (i, 9); 00802 00803 // 0 1 2 3 4 5 6 7 8 00804 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2}; 00805 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2}; 00806 00807 switch (j) 00808 { 00809 // d^2() / dxi^2 00810 case 0: 00811 return (FE<1,LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i0[i], 0, xi)* 00812 FE<1,LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)); 00813 00814 // d^2() / dxi deta 00815 case 1: 00816 return (FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)* 00817 FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)); 00818 00819 // d^2() / deta^2 00820 case 2: 00821 return (FE<1,LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)* 00822 FE<1,LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i1[i], 0, eta)); 00823 00824 default: 00825 { 00826 libMesh::err << "ERROR: Invalid derivative requested! " 00827 << std::endl; 00828 libmesh_error(); 00829 } 00830 } // end switch (j) 00831 } // end case QUAD9 00832 00833 case TRI6: 00834 { 00835 const Real dzeta0dxi = -1.; 00836 const Real dzeta1dxi = 1.; 00837 const Real dzeta2dxi = 0.; 00838 00839 const Real dzeta0deta = -1.; 00840 const Real dzeta1deta = 0.; 00841 const Real dzeta2deta = 1.; 00842 00843 libmesh_assert_less (j, 3); 00844 00845 switch (j) 00846 { 00847 // d^2() / dxi^2 00848 case 0: 00849 { 00850 switch (i) 00851 { 00852 case 0: 00853 return 4.*dzeta0dxi*dzeta0dxi; 00854 00855 case 1: 00856 return 4.*dzeta1dxi*dzeta1dxi; 00857 00858 case 2: 00859 return 4.*dzeta2dxi*dzeta2dxi; 00860 00861 case 3: 00862 return 8.*dzeta0dxi*dzeta1dxi; 00863 00864 case 4: 00865 return 8.*dzeta1dxi*dzeta2dxi; 00866 00867 case 5: 00868 return 8.*dzeta0dxi*dzeta2dxi; 00869 00870 default: 00871 { 00872 libMesh::err << "Invalid shape function index requested!" 00873 << std::endl; 00874 libmesh_error(); 00875 } 00876 } 00877 } 00878 00879 // d^2() / dxi deta 00880 case 1: 00881 { 00882 switch (i) 00883 { 00884 case 0: 00885 return 4.*dzeta0dxi*dzeta0deta; 00886 00887 case 1: 00888 return 4.*dzeta1dxi*dzeta1deta; 00889 00890 case 2: 00891 return 4.*dzeta2dxi*dzeta2deta; 00892 00893 case 3: 00894 return 4.*dzeta1deta*dzeta0dxi + 4.*dzeta0deta*dzeta1dxi; 00895 00896 case 4: 00897 return 4.*dzeta2deta*dzeta1dxi + 4.*dzeta1deta*dzeta2dxi; 00898 00899 case 5: 00900 return 4.*dzeta2deta*dzeta0dxi + 4.*dzeta0deta*dzeta2dxi; 00901 00902 default: 00903 { 00904 libMesh::err << "Invalid shape function index requested!" 00905 << std::endl; 00906 libmesh_error(); 00907 } 00908 } 00909 } 00910 00911 // d^2() / deta^2 00912 case 2: 00913 { 00914 switch (i) 00915 { 00916 case 0: 00917 return 4.*dzeta0deta*dzeta0deta; 00918 00919 case 1: 00920 return 4.*dzeta1deta*dzeta1deta; 00921 00922 case 2: 00923 return 4.*dzeta2deta*dzeta2deta; 00924 00925 case 3: 00926 return 8.*dzeta0deta*dzeta1deta; 00927 00928 case 4: 00929 return 8.*dzeta1deta*dzeta2deta; 00930 00931 case 5: 00932 return 8.*dzeta0deta*dzeta2deta; 00933 00934 default: 00935 { 00936 libMesh::err << "Invalid shape function index requested!" 00937 << std::endl; 00938 libmesh_error(); 00939 } 00940 } 00941 } 00942 00943 default: 00944 { 00945 libMesh::err << "ERROR: Invalid derivative requested! " 00946 << std::endl; 00947 libmesh_error(); 00948 } 00949 } // end switch (j) 00950 } // end case TRI6 00951 00952 default: 00953 { 00954 libMesh::err << "ERROR: Unsupported 2D element type!: " << type 00955 << std::endl; 00956 libmesh_error(); 00957 } 00958 } 00959 } // end case SECOND 00960 00961 00962 00963 // unsupported order 00964 default: 00965 { 00966 libMesh::err << "ERROR: Unsupported 2D FE order!: " << order 00967 << std::endl; 00968 libmesh_error(); 00969 } 00970 00971 } // end switch (order) 00972 00973 00974 libmesh_error(); 00975 #endif // LIBMESH_DIM > 1 00976 return 0.; 00977 } 00978 00979 00980 00981 template <> 00982 Real FE<2,LAGRANGE>::shape_second_deriv(const Elem* elem, 00983 const Order order, 00984 const unsigned int i, 00985 const unsigned int j, 00986 const Point& p) 00987 { 00988 libmesh_assert(elem); 00989 00990 // call the orientation-independent shape functions 00991 return FE<2,LAGRANGE>::shape_second_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p); 00992 } 00993 00994 } // namespace libMesh
Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:46 UTC
Hosted By: