fe_l2_lagrange_shape_3D.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 // 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<3,L2_LAGRANGE>::shape(const ElemType type, 00033 const Order order, 00034 const unsigned int i, 00035 const Point& p) 00036 { 00037 #if LIBMESH_DIM == 3 00038 00039 00040 switch (order) 00041 { 00042 // linear Lagrange shape functions 00043 case FIRST: 00044 { 00045 switch (type) 00046 { 00047 // trilinear hexahedral shape functions 00048 case HEX8: 00049 case HEX20: 00050 case HEX27: 00051 { 00052 libmesh_assert_less (i, 8); 00053 00054 // Compute hex shape functions as a tensor-product 00055 const Real xi = p(0); 00056 const Real eta = p(1); 00057 const Real zeta = p(2); 00058 00059 // 0 1 2 3 4 5 6 7 00060 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0}; 00061 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1}; 00062 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1}; 00063 00064 return (FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i0[i], xi)* 00065 FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i1[i], eta)* 00066 FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i2[i], zeta)); 00067 } 00068 00069 // linear tetrahedral shape functions 00070 case TET4: 00071 case TET10: 00072 { 00073 libmesh_assert_less (i, 4); 00074 00075 // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM 00076 const Real zeta1 = p(0); 00077 const Real zeta2 = p(1); 00078 const Real zeta3 = p(2); 00079 const Real zeta0 = 1. - zeta1 - zeta2 - zeta3; 00080 00081 switch(i) 00082 { 00083 case 0: 00084 return zeta0; 00085 00086 case 1: 00087 return zeta1; 00088 00089 case 2: 00090 return zeta2; 00091 00092 case 3: 00093 return zeta3; 00094 00095 default: 00096 libmesh_error_msg("Invalid i = " << i); 00097 } 00098 } 00099 00100 // linear prism shape functions 00101 case PRISM6: 00102 case PRISM15: 00103 case PRISM18: 00104 { 00105 libmesh_assert_less (i, 6); 00106 00107 // Compute prism shape functions as a tensor-product 00108 // of a triangle and an edge 00109 00110 Point p2d(p(0),p(1)); 00111 Point p1d(p(2)); 00112 00113 // 0 1 2 3 4 5 00114 static const unsigned int i0[] = {0, 0, 0, 1, 1, 1}; 00115 static const unsigned int i1[] = {0, 1, 2, 0, 1, 2}; 00116 00117 return (FE<2,L2_LAGRANGE>::shape(TRI3, FIRST, i1[i], p2d)* 00118 FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d)); 00119 } 00120 00121 // linear pyramid shape functions 00122 case PYRAMID5: 00123 { 00124 libmesh_assert_less (i, 5); 00125 00126 const Real xi = p(0); 00127 const Real eta = p(1); 00128 const Real zeta = p(2); 00129 const Real eps = 1.e-35; 00130 00131 switch(i) 00132 { 00133 case 0: 00134 return .25*(zeta + xi - 1.)*(zeta + eta - 1.)/((1. - zeta) + eps); 00135 00136 case 1: 00137 return .25*(zeta - xi - 1.)*(zeta + eta - 1.)/((1. - zeta) + eps); 00138 00139 case 2: 00140 return .25*(zeta - xi - 1.)*(zeta - eta - 1.)/((1. - zeta) + eps); 00141 00142 case 3: 00143 return .25*(zeta + xi - 1.)*(zeta - eta - 1.)/((1. - zeta) + eps); 00144 00145 case 4: 00146 return zeta; 00147 00148 default: 00149 libmesh_error_msg("Invalid i = " << i); 00150 } 00151 } 00152 00153 00154 default: 00155 libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type); 00156 } 00157 } 00158 00159 00160 // quadratic Lagrange shape functions 00161 case SECOND: 00162 { 00163 switch (type) 00164 { 00165 00166 // serendipity hexahedral quadratic shape functions 00167 case HEX20: 00168 { 00169 libmesh_assert_less (i, 20); 00170 00171 const Real xi = p(0); 00172 const Real eta = p(1); 00173 const Real zeta = p(2); 00174 00175 // these functions are defined for (x,y,z) in [0,1]^3 00176 // so transform the locations 00177 const Real x = .5*(xi + 1.); 00178 const Real y = .5*(eta + 1.); 00179 const Real z = .5*(zeta + 1.); 00180 00181 switch (i) 00182 { 00183 case 0: 00184 return (1. - x)*(1. - y)*(1. - z)*(1. - 2.*x - 2.*y - 2.*z); 00185 00186 case 1: 00187 return x*(1. - y)*(1. - z)*(2.*x - 2.*y - 2.*z - 1.); 00188 00189 case 2: 00190 return x*y*(1. - z)*(2.*x + 2.*y - 2.*z - 3.); 00191 00192 case 3: 00193 return (1. - x)*y*(1. - z)*(2.*y - 2.*x - 2.*z - 1.); 00194 00195 case 4: 00196 return (1. - x)*(1. - y)*z*(2.*z - 2.*x - 2.*y - 1.); 00197 00198 case 5: 00199 return x*(1. - y)*z*(2.*x - 2.*y + 2.*z - 3.); 00200 00201 case 6: 00202 return x*y*z*(2.*x + 2.*y + 2.*z - 5.); 00203 00204 case 7: 00205 return (1. - x)*y*z*(2.*y - 2.*x + 2.*z - 3.); 00206 00207 case 8: 00208 return 4.*x*(1. - x)*(1. - y)*(1. - z); 00209 00210 case 9: 00211 return 4.*x*y*(1. - y)*(1. - z); 00212 00213 case 10: 00214 return 4.*x*(1. - x)*y*(1. - z); 00215 00216 case 11: 00217 return 4.*(1. - x)*y*(1. - y)*(1. - z); 00218 00219 case 12: 00220 return 4.*(1. - x)*(1. - y)*z*(1. - z); 00221 00222 case 13: 00223 return 4.*x*(1. - y)*z*(1. - z); 00224 00225 case 14: 00226 return 4.*x*y*z*(1. - z); 00227 00228 case 15: 00229 return 4.*(1. - x)*y*z*(1. - z); 00230 00231 case 16: 00232 return 4.*x*(1. - x)*(1. - y)*z; 00233 00234 case 17: 00235 return 4.*x*y*(1. - y)*z; 00236 00237 case 18: 00238 return 4.*x*(1. - x)*y*z; 00239 00240 case 19: 00241 return 4.*(1. - x)*y*(1. - y)*z; 00242 00243 default: 00244 libmesh_error_msg("Invalid i = " << i); 00245 } 00246 } 00247 00248 // triquadraic hexahedral shape funcions 00249 case HEX27: 00250 { 00251 libmesh_assert_less (i, 27); 00252 00253 // Compute hex shape functions as a tensor-product 00254 const Real xi = p(0); 00255 const Real eta = p(1); 00256 const Real zeta = p(2); 00257 00258 // The only way to make any sense of this 00259 // is to look at the mgflo/mg2/mgf documentation 00260 // and make the cut-out cube! 00261 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 00262 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2}; 00263 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2}; 00264 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2}; 00265 00266 return (FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], xi)* 00267 FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i1[i], eta)* 00268 FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i2[i], zeta)); 00269 } 00270 00271 // quadratic tetrahedral shape functions 00272 case TET10: 00273 { 00274 libmesh_assert_less (i, 10); 00275 00276 // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM 00277 const Real zeta1 = p(0); 00278 const Real zeta2 = p(1); 00279 const Real zeta3 = p(2); 00280 const Real zeta0 = 1. - zeta1 - zeta2 - zeta3; 00281 00282 switch(i) 00283 { 00284 case 0: 00285 return zeta0*(2.*zeta0 - 1.); 00286 00287 case 1: 00288 return zeta1*(2.*zeta1 - 1.); 00289 00290 case 2: 00291 return zeta2*(2.*zeta2 - 1.); 00292 00293 case 3: 00294 return zeta3*(2.*zeta3 - 1.); 00295 00296 case 4: 00297 return 4.*zeta0*zeta1; 00298 00299 case 5: 00300 return 4.*zeta1*zeta2; 00301 00302 case 6: 00303 return 4.*zeta2*zeta0; 00304 00305 case 7: 00306 return 4.*zeta0*zeta3; 00307 00308 case 8: 00309 return 4.*zeta1*zeta3; 00310 00311 case 9: 00312 return 4.*zeta2*zeta3; 00313 00314 default: 00315 libmesh_error_msg("Invalid i = " << i); 00316 } 00317 } 00318 00319 // quadradic prism shape functions 00320 case PRISM18: 00321 { 00322 libmesh_assert_less (i, 18); 00323 00324 // Compute prism shape functions as a tensor-product 00325 // of a triangle and an edge 00326 00327 Point p2d(p(0),p(1)); 00328 Point p1d(p(2)); 00329 00330 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 00331 static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2}; 00332 static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5}; 00333 00334 return (FE<2,L2_LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)* 00335 FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d)); 00336 } 00337 00338 00339 default: 00340 libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type); 00341 } 00342 } 00343 00344 00345 // unsupported order 00346 default: 00347 libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order); 00348 } 00349 00350 #endif 00351 00352 libmesh_error_msg("We'll never get here!"); 00353 return 0.; 00354 } 00355 00356 00357 00358 template <> 00359 Real FE<3,L2_LAGRANGE>::shape(const Elem* elem, 00360 const Order order, 00361 const unsigned int i, 00362 const Point& p) 00363 { 00364 libmesh_assert(elem); 00365 00366 // call the orientation-independent shape functions 00367 return FE<3,L2_LAGRANGE>::shape(elem->type(), static_cast<Order>(order + elem->p_level()), i, p); 00368 } 00369 00370 00371 00372 00373 template <> 00374 Real FE<3,L2_LAGRANGE>::shape_deriv(const ElemType type, 00375 const Order order, 00376 const unsigned int i, 00377 const unsigned int j, 00378 const Point& p) 00379 { 00380 #if LIBMESH_DIM == 3 00381 00382 libmesh_assert_less (j, 3); 00383 00384 switch (order) 00385 { 00386 // linear Lagrange shape functions 00387 case FIRST: 00388 { 00389 switch (type) 00390 { 00391 // trilinear hexahedral shape functions 00392 case HEX8: 00393 case HEX20: 00394 case HEX27: 00395 { 00396 libmesh_assert_less (i, 8); 00397 00398 // Compute hex shape functions as a tensor-product 00399 const Real xi = p(0); 00400 const Real eta = p(1); 00401 const Real zeta = p(2); 00402 00403 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0}; 00404 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1}; 00405 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1}; 00406 00407 switch(j) 00408 { 00409 case 0: 00410 return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, xi)* 00411 FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i1[i], eta)* 00412 FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i2[i], zeta)); 00413 00414 case 1: 00415 return (FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i0[i], xi)* 00416 FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i1[i], 0, eta)* 00417 FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i2[i], zeta)); 00418 00419 case 2: 00420 return (FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i0[i], xi)* 00421 FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i1[i], eta)* 00422 FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i2[i], 0, zeta)); 00423 00424 default: 00425 libmesh_error_msg("Invalid j = " << j); 00426 } 00427 } 00428 00429 // linear tetrahedral shape functions 00430 case TET4: 00431 case TET10: 00432 { 00433 libmesh_assert_less (i, 4); 00434 00435 // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM 00436 const Real dzeta0dxi = -1.; 00437 const Real dzeta1dxi = 1.; 00438 const Real dzeta2dxi = 0.; 00439 const Real dzeta3dxi = 0.; 00440 00441 const Real dzeta0deta = -1.; 00442 const Real dzeta1deta = 0.; 00443 const Real dzeta2deta = 1.; 00444 const Real dzeta3deta = 0.; 00445 00446 const Real dzeta0dzeta = -1.; 00447 const Real dzeta1dzeta = 0.; 00448 const Real dzeta2dzeta = 0.; 00449 const Real dzeta3dzeta = 1.; 00450 00451 switch (j) 00452 { 00453 // d()/dxi 00454 case 0: 00455 { 00456 switch(i) 00457 { 00458 case 0: 00459 return dzeta0dxi; 00460 00461 case 1: 00462 return dzeta1dxi; 00463 00464 case 2: 00465 return dzeta2dxi; 00466 00467 case 3: 00468 return dzeta3dxi; 00469 00470 default: 00471 libmesh_error_msg("Invalid i = " << i); 00472 } 00473 } 00474 00475 // d()/deta 00476 case 1: 00477 { 00478 switch(i) 00479 { 00480 case 0: 00481 return dzeta0deta; 00482 00483 case 1: 00484 return dzeta1deta; 00485 00486 case 2: 00487 return dzeta2deta; 00488 00489 case 3: 00490 return dzeta3deta; 00491 00492 default: 00493 libmesh_error_msg("Invalid i = " << i); 00494 } 00495 } 00496 00497 // d()/dzeta 00498 case 2: 00499 { 00500 switch(i) 00501 { 00502 case 0: 00503 return dzeta0dzeta; 00504 00505 case 1: 00506 return dzeta1dzeta; 00507 00508 case 2: 00509 return dzeta2dzeta; 00510 00511 case 3: 00512 return dzeta3dzeta; 00513 00514 default: 00515 libmesh_error_msg("Invalid i = " << i); 00516 } 00517 } 00518 00519 default: 00520 libmesh_error_msg("Invalid shape function derivative j = " << j); 00521 } 00522 } 00523 00524 // linear prism shape functions 00525 case PRISM6: 00526 case PRISM15: 00527 case PRISM18: 00528 { 00529 libmesh_assert_less (i, 6); 00530 00531 // Compute prism shape functions as a tensor-product 00532 // of a triangle and an edge 00533 00534 Point p2d(p(0),p(1)); 00535 Point p1d(p(2)); 00536 00537 // 0 1 2 3 4 5 00538 static const unsigned int i0[] = {0, 0, 0, 1, 1, 1}; 00539 static const unsigned int i1[] = {0, 1, 2, 0, 1, 2}; 00540 00541 switch (j) 00542 { 00543 // d()/dxi 00544 case 0: 00545 return (FE<2,L2_LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 0, p2d)* 00546 FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d)); 00547 00548 // d()/deta 00549 case 1: 00550 return (FE<2,L2_LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 1, p2d)* 00551 FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d)); 00552 00553 // d()/dzeta 00554 case 2: 00555 return (FE<2,L2_LAGRANGE>::shape(TRI3, FIRST, i1[i], p2d)* 00556 FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, p1d)); 00557 00558 default: 00559 libmesh_error_msg("Invalid shape function derivative j = " << j); 00560 } 00561 } 00562 00563 // linear pyramid shape functions 00564 case PYRAMID5: 00565 { 00566 libmesh_assert_less (i, 5); 00567 00568 const Real xi = p(0); 00569 const Real eta = p(1); 00570 const Real zeta = p(2); 00571 const Real eps = 1.e-35; 00572 00573 switch (j) 00574 { 00575 // d/dxi 00576 case 0: 00577 switch(i) 00578 { 00579 case 0: 00580 return .25*(zeta + eta - 1.)/((1. - zeta) + eps); 00581 00582 case 1: 00583 return -.25*(zeta + eta - 1.)/((1. - zeta) + eps); 00584 00585 case 2: 00586 return -.25*(zeta - eta - 1.)/((1. - zeta) + eps); 00587 00588 case 3: 00589 return .25*(zeta - eta - 1.)/((1. - zeta) + eps); 00590 00591 case 4: 00592 return 0; 00593 00594 default: 00595 libmesh_error_msg("Invalid i = " << i); 00596 } 00597 00598 00599 // d/deta 00600 case 1: 00601 switch(i) 00602 { 00603 case 0: 00604 return .25*(zeta + xi - 1.)/((1. - zeta) + eps); 00605 00606 case 1: 00607 return .25*(zeta - xi - 1.)/((1. - zeta) + eps); 00608 00609 case 2: 00610 return -.25*(zeta - xi - 1.)/((1. - zeta) + eps); 00611 00612 case 3: 00613 return -.25*(zeta + xi - 1.)/((1. - zeta) + eps); 00614 00615 case 4: 00616 return 0; 00617 00618 default: 00619 libmesh_error_msg("Invalid i = " << i); 00620 } 00621 00622 00623 // d/dzeta 00624 case 2: 00625 switch(i) 00626 { 00627 case 0: 00628 { 00629 const Real a=1.; 00630 const Real b=1.; 00631 00632 return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) + 00633 (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/ 00634 ((1. - zeta)*(1. - zeta) + eps)); 00635 } 00636 00637 case 1: 00638 { 00639 const Real a=-1.; 00640 const Real b=1.; 00641 00642 return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) + 00643 (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/ 00644 ((1. - zeta)*(1. - zeta) + eps)); 00645 } 00646 00647 case 2: 00648 { 00649 const Real a=-1.; 00650 const Real b=-1.; 00651 00652 return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) + 00653 (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/ 00654 ((1. - zeta)*(1. - zeta) + eps)); 00655 } 00656 00657 case 3: 00658 { 00659 const Real a=1.; 00660 const Real b=-1.; 00661 00662 return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) + 00663 (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/ 00664 ((1. - zeta)*(1. - zeta) + eps)); 00665 } 00666 00667 case 4: 00668 return 1.; 00669 00670 default: 00671 libmesh_error_msg("Invalid i = " << i); 00672 } 00673 00674 00675 default: 00676 libmesh_error_msg("Invalid j = " << j); 00677 } 00678 } 00679 00680 00681 default: 00682 libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type); 00683 } 00684 } 00685 00686 00687 // quadratic Lagrange shape functions 00688 case SECOND: 00689 { 00690 switch (type) 00691 { 00692 00693 // serendipity hexahedral quadratic shape functions 00694 case HEX20: 00695 { 00696 libmesh_assert_less (i, 20); 00697 00698 const Real xi = p(0); 00699 const Real eta = p(1); 00700 const Real zeta = p(2); 00701 00702 // these functions are defined for (x,y,z) in [0,1]^3 00703 // so transform the locations 00704 const Real x = .5*(xi + 1.); 00705 const Real y = .5*(eta + 1.); 00706 const Real z = .5*(zeta + 1.); 00707 00708 // and don't forget the chain rule! 00709 00710 switch (j) 00711 { 00712 00713 // d/dx*dx/dxi 00714 case 0: 00715 switch (i) 00716 { 00717 case 0: 00718 return .5*(1. - y)*(1. - z)*((1. - x)*(-2.) + 00719 (-1.)*(1. - 2.*x - 2.*y - 2.*z)); 00720 00721 case 1: 00722 return .5*(1. - y)*(1. - z)*(x*(2.) + 00723 (1.)*(2.*x - 2.*y - 2.*z - 1.)); 00724 00725 case 2: 00726 return .5*y*(1. - z)*(x*(2.) + 00727 (1.)*(2.*x + 2.*y - 2.*z - 3.)); 00728 00729 case 3: 00730 return .5*y*(1. - z)*((1. - x)*(-2.) + 00731 (-1.)*(2.*y - 2.*x - 2.*z - 1.)); 00732 00733 case 4: 00734 return .5*(1. - y)*z*((1. - x)*(-2.) + 00735 (-1.)*(2.*z - 2.*x - 2.*y - 1.)); 00736 00737 case 5: 00738 return .5*(1. - y)*z*(x*(2.) + 00739 (1.)*(2.*x - 2.*y + 2.*z - 3.)); 00740 00741 case 6: 00742 return .5*y*z*(x*(2.) + 00743 (1.)*(2.*x + 2.*y + 2.*z - 5.)); 00744 00745 case 7: 00746 return .5*y*z*((1. - x)*(-2.) + 00747 (-1.)*(2.*y - 2.*x + 2.*z - 3.)); 00748 00749 case 8: 00750 return 2.*(1. - y)*(1. - z)*(1. - 2.*x); 00751 00752 case 9: 00753 return 2.*y*(1. - y)*(1. - z); 00754 00755 case 10: 00756 return 2.*y*(1. - z)*(1. - 2.*x); 00757 00758 case 11: 00759 return 2.*y*(1. - y)*(1. - z)*(-1.); 00760 00761 case 12: 00762 return 2.*(1. - y)*z*(1. - z)*(-1.); 00763 00764 case 13: 00765 return 2.*(1. - y)*z*(1. - z); 00766 00767 case 14: 00768 return 2.*y*z*(1. - z); 00769 00770 case 15: 00771 return 2.*y*z*(1. - z)*(-1.); 00772 00773 case 16: 00774 return 2.*(1. - y)*z*(1. - 2.*x); 00775 00776 case 17: 00777 return 2.*y*(1. - y)*z; 00778 00779 case 18: 00780 return 2.*y*z*(1. - 2.*x); 00781 00782 case 19: 00783 return 2.*y*(1. - y)*z*(-1.); 00784 00785 default: 00786 libmesh_error_msg("Invalid i = " << i); 00787 } 00788 00789 00790 // d/dy*dy/deta 00791 case 1: 00792 switch (i) 00793 { 00794 case 0: 00795 return .5*(1. - x)*(1. - z)*((1. - y)*(-2.) + 00796 (-1.)*(1. - 2.*x - 2.*y - 2.*z)); 00797 00798 case 1: 00799 return .5*x*(1. - z)*((1. - y)*(-2.) + 00800 (-1.)*(2.*x - 2.*y - 2.*z - 1.)); 00801 00802 case 2: 00803 return .5*x*(1. - z)*(y*(2.) + 00804 (1.)*(2.*x + 2.*y - 2.*z - 3.)); 00805 00806 case 3: 00807 return .5*(1. - x)*(1. - z)*(y*(2.) + 00808 (1.)*(2.*y - 2.*x - 2.*z - 1.)); 00809 00810 case 4: 00811 return .5*(1. - x)*z*((1. - y)*(-2.) + 00812 (-1.)*(2.*z - 2.*x - 2.*y - 1.)); 00813 00814 case 5: 00815 return .5*x*z*((1. - y)*(-2.) + 00816 (-1.)*(2.*x - 2.*y + 2.*z - 3.)); 00817 00818 case 6: 00819 return .5*x*z*(y*(2.) + 00820 (1.)*(2.*x + 2.*y + 2.*z - 5.)); 00821 00822 case 7: 00823 return .5*(1. - x)*z*(y*(2.) + 00824 (1.)*(2.*y - 2.*x + 2.*z - 3.)); 00825 00826 case 8: 00827 return 2.*x*(1. - x)*(1. - z)*(-1.); 00828 00829 case 9: 00830 return 2.*x*(1. - z)*(1. - 2.*y); 00831 00832 case 10: 00833 return 2.*x*(1. - x)*(1. - z); 00834 00835 case 11: 00836 return 2.*(1. - x)*(1. - z)*(1. - 2.*y); 00837 00838 case 12: 00839 return 2.*(1. - x)*z*(1. - z)*(-1.); 00840 00841 case 13: 00842 return 2.*x*z*(1. - z)*(-1.); 00843 00844 case 14: 00845 return 2.*x*z*(1. - z); 00846 00847 case 15: 00848 return 2.*(1. - x)*z*(1. - z); 00849 00850 case 16: 00851 return 2.*x*(1. - x)*z*(-1.); 00852 00853 case 17: 00854 return 2.*x*z*(1. - 2.*y); 00855 00856 case 18: 00857 return 2.*x*(1. - x)*z; 00858 00859 case 19: 00860 return 2.*(1. - x)*z*(1. - 2.*y); 00861 00862 default: 00863 libmesh_error_msg("Invalid i = " << i); 00864 } 00865 00866 00867 // d/dz*dz/dzeta 00868 case 2: 00869 switch (i) 00870 { 00871 case 0: 00872 return .5*(1. - x)*(1. - y)*((1. - z)*(-2.) + 00873 (-1.)*(1. - 2.*x - 2.*y - 2.*z)); 00874 00875 case 1: 00876 return .5*x*(1. - y)*((1. - z)*(-2.) + 00877 (-1.)*(2.*x - 2.*y - 2.*z - 1.)); 00878 00879 case 2: 00880 return .5*x*y*((1. - z)*(-2.) + 00881 (-1.)*(2.*x + 2.*y - 2.*z - 3.)); 00882 00883 case 3: 00884 return .5*(1. - x)*y*((1. - z)*(-2.) + 00885 (-1.)*(2.*y - 2.*x - 2.*z - 1.)); 00886 00887 case 4: 00888 return .5*(1. - x)*(1. - y)*(z*(2.) + 00889 (1.)*(2.*z - 2.*x - 2.*y - 1.)); 00890 00891 case 5: 00892 return .5*x*(1. - y)*(z*(2.) + 00893 (1.)*(2.*x - 2.*y + 2.*z - 3.)); 00894 00895 case 6: 00896 return .5*x*y*(z*(2.) + 00897 (1.)*(2.*x + 2.*y + 2.*z - 5.)); 00898 00899 case 7: 00900 return .5*(1. - x)*y*(z*(2.) + 00901 (1.)*(2.*y - 2.*x + 2.*z - 3.)); 00902 00903 case 8: 00904 return 2.*x*(1. - x)*(1. - y)*(-1.); 00905 00906 case 9: 00907 return 2.*x*y*(1. - y)*(-1.); 00908 00909 case 10: 00910 return 2.*x*(1. - x)*y*(-1.); 00911 00912 case 11: 00913 return 2.*(1. - x)*y*(1. - y)*(-1.); 00914 00915 case 12: 00916 return 2.*(1. - x)*(1. - y)*(1. - 2.*z); 00917 00918 case 13: 00919 return 2.*x*(1. - y)*(1. - 2.*z); 00920 00921 case 14: 00922 return 2.*x*y*(1. - 2.*z); 00923 00924 case 15: 00925 return 2.*(1. - x)*y*(1. - 2.*z); 00926 00927 case 16: 00928 return 2.*x*(1. - x)*(1. - y); 00929 00930 case 17: 00931 return 2.*x*y*(1. - y); 00932 00933 case 18: 00934 return 2.*x*(1. - x)*y; 00935 00936 case 19: 00937 return 2.*(1. - x)*y*(1. - y); 00938 00939 default: 00940 libmesh_error_msg("Invalid i = " << i); 00941 } 00942 00943 default: 00944 libmesh_error_msg("Invalid shape function derivative j = " << j); 00945 } 00946 } 00947 00948 // triquadraic hexahedral shape funcions 00949 case HEX27: 00950 { 00951 libmesh_assert_less (i, 27); 00952 00953 // Compute hex shape functions as a tensor-product 00954 const Real xi = p(0); 00955 const Real eta = p(1); 00956 const Real zeta = p(2); 00957 00958 // The only way to make any sense of this 00959 // is to look at the mgflo/mg2/mgf documentation 00960 // and make the cut-out cube! 00961 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 00962 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2}; 00963 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2}; 00964 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2}; 00965 00966 switch(j) 00967 { 00968 case 0: 00969 return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)* 00970 FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)* 00971 FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta)); 00972 00973 case 1: 00974 return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)* 00975 FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)* 00976 FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta)); 00977 00978 case 2: 00979 return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)* 00980 FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)* 00981 FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta)); 00982 00983 default: 00984 libmesh_error_msg("Invalid j = " << j); 00985 } 00986 } 00987 00988 // quadratic tetrahedral shape functions 00989 case TET10: 00990 { 00991 libmesh_assert_less (i, 10); 00992 00993 // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM 00994 const Real zeta1 = p(0); 00995 const Real zeta2 = p(1); 00996 const Real zeta3 = p(2); 00997 const Real zeta0 = 1. - zeta1 - zeta2 - zeta3; 00998 00999 const Real dzeta0dxi = -1.; 01000 const Real dzeta1dxi = 1.; 01001 const Real dzeta2dxi = 0.; 01002 const Real dzeta3dxi = 0.; 01003 01004 const Real dzeta0deta = -1.; 01005 const Real dzeta1deta = 0.; 01006 const Real dzeta2deta = 1.; 01007 const Real dzeta3deta = 0.; 01008 01009 const Real dzeta0dzeta = -1.; 01010 const Real dzeta1dzeta = 0.; 01011 const Real dzeta2dzeta = 0.; 01012 const Real dzeta3dzeta = 1.; 01013 01014 switch (j) 01015 { 01016 // d()/dxi 01017 case 0: 01018 { 01019 switch(i) 01020 { 01021 case 0: 01022 return (4.*zeta0 - 1.)*dzeta0dxi; 01023 01024 case 1: 01025 return (4.*zeta1 - 1.)*dzeta1dxi; 01026 01027 case 2: 01028 return (4.*zeta2 - 1.)*dzeta2dxi; 01029 01030 case 3: 01031 return (4.*zeta3 - 1.)*dzeta3dxi; 01032 01033 case 4: 01034 return 4.*(zeta0*dzeta1dxi + dzeta0dxi*zeta1); 01035 01036 case 5: 01037 return 4.*(zeta1*dzeta2dxi + dzeta1dxi*zeta2); 01038 01039 case 6: 01040 return 4.*(zeta0*dzeta2dxi + dzeta0dxi*zeta2); 01041 01042 case 7: 01043 return 4.*(zeta0*dzeta3dxi + dzeta0dxi*zeta3); 01044 01045 case 8: 01046 return 4.*(zeta1*dzeta3dxi + dzeta1dxi*zeta3); 01047 01048 case 9: 01049 return 4.*(zeta2*dzeta3dxi + dzeta2dxi*zeta3); 01050 01051 default: 01052 libmesh_error_msg("Invalid i = " << i); 01053 } 01054 } 01055 01056 // d()/deta 01057 case 1: 01058 { 01059 switch(i) 01060 { 01061 case 0: 01062 return (4.*zeta0 - 1.)*dzeta0deta; 01063 01064 case 1: 01065 return (4.*zeta1 - 1.)*dzeta1deta; 01066 01067 case 2: 01068 return (4.*zeta2 - 1.)*dzeta2deta; 01069 01070 case 3: 01071 return (4.*zeta3 - 1.)*dzeta3deta; 01072 01073 case 4: 01074 return 4.*(zeta0*dzeta1deta + dzeta0deta*zeta1); 01075 01076 case 5: 01077 return 4.*(zeta1*dzeta2deta + dzeta1deta*zeta2); 01078 01079 case 6: 01080 return 4.*(zeta0*dzeta2deta + dzeta0deta*zeta2); 01081 01082 case 7: 01083 return 4.*(zeta0*dzeta3deta + dzeta0deta*zeta3); 01084 01085 case 8: 01086 return 4.*(zeta1*dzeta3deta + dzeta1deta*zeta3); 01087 01088 case 9: 01089 return 4.*(zeta2*dzeta3deta + dzeta2deta*zeta3); 01090 01091 default: 01092 libmesh_error_msg("Invalid i = " << i); 01093 } 01094 } 01095 01096 // d()/dzeta 01097 case 2: 01098 { 01099 switch(i) 01100 { 01101 case 0: 01102 return (4.*zeta0 - 1.)*dzeta0dzeta; 01103 01104 case 1: 01105 return (4.*zeta1 - 1.)*dzeta1dzeta; 01106 01107 case 2: 01108 return (4.*zeta2 - 1.)*dzeta2dzeta; 01109 01110 case 3: 01111 return (4.*zeta3 - 1.)*dzeta3dzeta; 01112 01113 case 4: 01114 return 4.*(zeta0*dzeta1dzeta + dzeta0dzeta*zeta1); 01115 01116 case 5: 01117 return 4.*(zeta1*dzeta2dzeta + dzeta1dzeta*zeta2); 01118 01119 case 6: 01120 return 4.*(zeta0*dzeta2dzeta + dzeta0dzeta*zeta2); 01121 01122 case 7: 01123 return 4.*(zeta0*dzeta3dzeta + dzeta0dzeta*zeta3); 01124 01125 case 8: 01126 return 4.*(zeta1*dzeta3dzeta + dzeta1dzeta*zeta3); 01127 01128 case 9: 01129 return 4.*(zeta2*dzeta3dzeta + dzeta2dzeta*zeta3); 01130 01131 default: 01132 libmesh_error_msg("Invalid i = " << i); 01133 } 01134 } 01135 01136 default: 01137 libmesh_error_msg("Invalid j = " << j); 01138 } 01139 } 01140 01141 01142 01143 // quadradic prism shape functions 01144 case PRISM18: 01145 { 01146 libmesh_assert_less (i, 18); 01147 01148 // Compute prism shape functions as a tensor-product 01149 // of a triangle and an edge 01150 01151 Point p2d(p(0),p(1)); 01152 Point p1d(p(2)); 01153 01154 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 01155 static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2}; 01156 static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5}; 01157 01158 switch (j) 01159 { 01160 // d()/dxi 01161 case 0: 01162 return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)* 01163 FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d)); 01164 01165 // d()/deta 01166 case 1: 01167 return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)* 01168 FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d)); 01169 01170 // d()/dzeta 01171 case 2: 01172 return (FE<2,L2_LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)* 01173 FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d)); 01174 01175 default: 01176 libmesh_error_msg("Invalid shape function derivative j = " << j); 01177 } 01178 } 01179 01180 01181 default: 01182 libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type); 01183 } 01184 } 01185 01186 01187 // unsupported order 01188 default: 01189 libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order); 01190 } 01191 01192 #endif 01193 01194 libmesh_error_msg("We'll never get here!"); 01195 return 0.; 01196 } 01197 01198 01199 01200 template <> 01201 Real FE<3,L2_LAGRANGE>::shape_deriv(const Elem* elem, 01202 const Order order, 01203 const unsigned int i, 01204 const unsigned int j, 01205 const Point& p) 01206 { 01207 libmesh_assert(elem); 01208 01209 // call the orientation-independent shape function derivatives 01210 return FE<3,L2_LAGRANGE>::shape_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p); 01211 } 01212 01213 01214 01215 template <> 01216 Real FE<3,L2_LAGRANGE>::shape_second_deriv(const ElemType type, 01217 const Order order, 01218 const unsigned int i, 01219 const unsigned int j, 01220 const Point& p) 01221 { 01222 #if LIBMESH_DIM == 3 01223 01224 libmesh_assert_less (j, 6); 01225 01226 switch (order) 01227 { 01228 // linear Lagrange shape functions 01229 case FIRST: 01230 { 01231 return 0.; 01232 } 01233 01234 // quadratic Lagrange shape functions 01235 case SECOND: 01236 { 01237 switch (type) 01238 { 01239 01240 // serendipity hexahedral quadratic shape functions 01241 case HEX20: 01242 { 01243 static bool warning_given_HEX20 = false; 01244 01245 if (!warning_given_HEX20) 01246 libMesh::err << "Second derivatives for 3D Lagrangian HEX20" 01247 << " elements are not yet implemented!" 01248 << std::endl; 01249 warning_given_HEX20 = true; 01250 } 01251 01252 // triquadraic hexahedral shape funcions 01253 case HEX27: 01254 { 01255 libmesh_assert_less (i, 27); 01256 01257 // Compute hex shape functions as a tensor-product 01258 const Real xi = p(0); 01259 const Real eta = p(1); 01260 const Real zeta = p(2); 01261 01262 // The only way to make any sense of this 01263 // is to look at the mgflo/mg2/mgf documentation 01264 // and make the cut-out cube! 01265 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 01266 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2}; 01267 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2}; 01268 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2}; 01269 01270 switch(j) 01271 { 01272 // d^2()/dxi^2 01273 case 0: 01274 return (FE<1,L2_LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i0[i], 0, xi)* 01275 FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)* 01276 FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta)); 01277 01278 // d^2()/dxideta 01279 case 1: 01280 return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)* 01281 FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)* 01282 FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta)); 01283 01284 // d^2()/deta^2 01285 case 2: 01286 return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)* 01287 FE<1,L2_LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i1[i], 0, eta)* 01288 FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta)); 01289 01290 // d^2()/dxidzeta 01291 case 3: 01292 return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)* 01293 FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)* 01294 FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta)); 01295 01296 // d^2()/detadzeta 01297 case 4: 01298 return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)* 01299 FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)* 01300 FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta)); 01301 01302 // d^2()/dzeta^2 01303 case 5: 01304 return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)* 01305 FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)* 01306 FE<1,L2_LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i2[i], 0, zeta)); 01307 01308 default: 01309 libmesh_error_msg("Invalid j = " << j); 01310 } 01311 } 01312 01313 // quadratic tetrahedral shape functions 01314 case TET10: 01315 { 01316 // The area coordinates are the same as used for the 01317 // shape() and shape_deriv() functions. 01318 // const Real zeta0 = 1. - zeta1 - zeta2 - zeta3; 01319 // const Real zeta1 = p(0); 01320 // const Real zeta2 = p(1); 01321 // const Real zeta3 = p(2); 01322 static const Real dzetadxi[4][3] = 01323 { 01324 {-1., -1., -1.}, 01325 {1., 0., 0.}, 01326 {0., 1., 0.}, 01327 {0., 0., 1.} 01328 }; 01329 01330 // Convert from j -> (j,k) indices for independent variable 01331 // (0=xi, 1=eta, 2=zeta) 01332 static const unsigned short int independent_var_indices[6][2] = 01333 { 01334 {0, 0}, // d^2 phi / dxi^2 01335 {0, 1}, // d^2 phi / dxi deta 01336 {1, 1}, // d^2 phi / deta^2 01337 {0, 2}, // d^2 phi / dxi dzeta 01338 {1, 2}, // d^2 phi / deta dzeta 01339 {2, 2} // d^2 phi / dzeta^2 01340 }; 01341 01342 // Convert from i -> zeta indices. Each quadratic shape 01343 // function for the Tet10 depends on up to two of the zeta 01344 // area coordinate functions (see the shape() function above). 01345 // This table just tells which two area coords it uses. 01346 static const unsigned short int zeta_indices[10][2] = 01347 { 01348 {0, 0}, 01349 {1, 1}, 01350 {2, 2}, 01351 {3, 3}, 01352 {0, 1}, 01353 {1, 2}, 01354 {2, 0}, 01355 {0, 3}, 01356 {1, 3}, 01357 {2, 3}, 01358 }; 01359 01360 // Look up the independent variable indices for this value of j. 01361 const unsigned int my_j = independent_var_indices[j][0]; 01362 const unsigned int my_k = independent_var_indices[j][1]; 01363 01364 if (i<4) 01365 { 01366 return 4.*dzetadxi[i][my_j]*dzetadxi[i][my_k]; 01367 } 01368 01369 else if (i<10) 01370 { 01371 const unsigned short int my_m = zeta_indices[i][0]; 01372 const unsigned short int my_n = zeta_indices[i][1]; 01373 01374 return 4.*(dzetadxi[my_n][my_j]*dzetadxi[my_m][my_k] + 01375 dzetadxi[my_m][my_j]*dzetadxi[my_n][my_k] ); 01376 } 01377 else 01378 libmesh_error_msg("Invalid shape function index " << i); 01379 } 01380 01381 01382 // quadradic prism shape functions 01383 case PRISM18: 01384 { 01385 libmesh_assert_less (i, 18); 01386 01387 // Compute prism shape functions as a tensor-product 01388 // of a triangle and an edge 01389 01390 Point p2d(p(0),p(1)); 01391 Point p1d(p(2)); 01392 01393 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 01394 static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2}; 01395 static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5}; 01396 01397 switch (j) 01398 { 01399 // d^2()/dxi^2 01400 case 0: 01401 return (FE<2,L2_LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 0, p2d)* 01402 FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d)); 01403 01404 // d^2()/dxideta 01405 case 1: 01406 return (FE<2,L2_LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 1, p2d)* 01407 FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d)); 01408 01409 // d^2()/deta^2 01410 case 2: 01411 return (FE<2,L2_LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 2, p2d)* 01412 FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d)); 01413 01414 // d^2()/dxidzeta 01415 case 3: 01416 return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)* 01417 FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d)); 01418 01419 // d^2()/detadzeta 01420 case 4: 01421 return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)* 01422 FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d)); 01423 01424 // d^2()/dzeta^2 01425 case 5: 01426 return (FE<2,L2_LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)* 01427 FE<1,L2_LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i0[i], 0, p1d)); 01428 01429 default: 01430 libmesh_error_msg("Invalid shape function derivative j = " << j); 01431 } 01432 } 01433 01434 01435 01436 default: 01437 libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type); 01438 } 01439 } 01440 01441 01442 // unsupported order 01443 default: 01444 libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order); 01445 } 01446 01447 #endif 01448 01449 libmesh_error_msg("We'll never get here!"); 01450 return 0.; 01451 } 01452 01453 01454 01455 template <> 01456 Real FE<3,L2_LAGRANGE>::shape_second_deriv(const Elem* elem, 01457 const Order order, 01458 const unsigned int i, 01459 const unsigned int j, 01460 const Point& p) 01461 { 01462 libmesh_assert(elem); 01463 01464 // call the orientation-independent shape function derivatives 01465 return FE<3,L2_LAGRANGE>::shape_second_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p); 01466 } 01467 01468 } // namespace libMesh
generated by
