fe_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,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,LAGRANGE>::shape(EDGE2, FIRST, i0[i], xi)* 00065 FE<1,LAGRANGE>::shape(EDGE2, FIRST, i1[i], eta)* 00066 FE<1,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,LAGRANGE>::shape(TRI3, FIRST, i1[i], p2d)* 00118 FE<1,LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d)); 00119 } 00120 00121 // linear pyramid shape functions 00122 case PYRAMID5: 00123 case PYRAMID13: 00124 case PYRAMID14: 00125 { 00126 libmesh_assert_less (i, 5); 00127 00128 const Real xi = p(0); 00129 const Real eta = p(1); 00130 const Real zeta = p(2); 00131 const Real eps = 1.e-35; 00132 00133 switch(i) 00134 { 00135 case 0: 00136 return .25*(zeta + xi - 1.)*(zeta + eta - 1.)/((1. - zeta) + eps); 00137 00138 case 1: 00139 return .25*(zeta - xi - 1.)*(zeta + eta - 1.)/((1. - zeta) + eps); 00140 00141 case 2: 00142 return .25*(zeta - xi - 1.)*(zeta - eta - 1.)/((1. - zeta) + eps); 00143 00144 case 3: 00145 return .25*(zeta + xi - 1.)*(zeta - eta - 1.)/((1. - zeta) + eps); 00146 00147 case 4: 00148 return zeta; 00149 00150 default: 00151 libmesh_error_msg("Invalid i = " << i); 00152 } 00153 } 00154 00155 00156 default: 00157 libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type); 00158 } 00159 } 00160 00161 00162 // quadratic Lagrange shape functions 00163 case SECOND: 00164 { 00165 switch (type) 00166 { 00167 00168 // serendipity hexahedral quadratic shape functions 00169 case HEX20: 00170 { 00171 libmesh_assert_less (i, 20); 00172 00173 const Real xi = p(0); 00174 const Real eta = p(1); 00175 const Real zeta = p(2); 00176 00177 // these functions are defined for (x,y,z) in [0,1]^3 00178 // so transform the locations 00179 const Real x = .5*(xi + 1.); 00180 const Real y = .5*(eta + 1.); 00181 const Real z = .5*(zeta + 1.); 00182 00183 switch (i) 00184 { 00185 case 0: 00186 return (1. - x)*(1. - y)*(1. - z)*(1. - 2.*x - 2.*y - 2.*z); 00187 00188 case 1: 00189 return x*(1. - y)*(1. - z)*(2.*x - 2.*y - 2.*z - 1.); 00190 00191 case 2: 00192 return x*y*(1. - z)*(2.*x + 2.*y - 2.*z - 3.); 00193 00194 case 3: 00195 return (1. - x)*y*(1. - z)*(2.*y - 2.*x - 2.*z - 1.); 00196 00197 case 4: 00198 return (1. - x)*(1. - y)*z*(2.*z - 2.*x - 2.*y - 1.); 00199 00200 case 5: 00201 return x*(1. - y)*z*(2.*x - 2.*y + 2.*z - 3.); 00202 00203 case 6: 00204 return x*y*z*(2.*x + 2.*y + 2.*z - 5.); 00205 00206 case 7: 00207 return (1. - x)*y*z*(2.*y - 2.*x + 2.*z - 3.); 00208 00209 case 8: 00210 return 4.*x*(1. - x)*(1. - y)*(1. - z); 00211 00212 case 9: 00213 return 4.*x*y*(1. - y)*(1. - z); 00214 00215 case 10: 00216 return 4.*x*(1. - x)*y*(1. - z); 00217 00218 case 11: 00219 return 4.*(1. - x)*y*(1. - y)*(1. - z); 00220 00221 case 12: 00222 return 4.*(1. - x)*(1. - y)*z*(1. - z); 00223 00224 case 13: 00225 return 4.*x*(1. - y)*z*(1. - z); 00226 00227 case 14: 00228 return 4.*x*y*z*(1. - z); 00229 00230 case 15: 00231 return 4.*(1. - x)*y*z*(1. - z); 00232 00233 case 16: 00234 return 4.*x*(1. - x)*(1. - y)*z; 00235 00236 case 17: 00237 return 4.*x*y*(1. - y)*z; 00238 00239 case 18: 00240 return 4.*x*(1. - x)*y*z; 00241 00242 case 19: 00243 return 4.*(1. - x)*y*(1. - y)*z; 00244 00245 default: 00246 libmesh_error_msg("Invalid i = " << i); 00247 } 00248 } 00249 00250 // triquadraic hexahedral shape funcions 00251 case HEX27: 00252 { 00253 libmesh_assert_less (i, 27); 00254 00255 // Compute hex shape functions as a tensor-product 00256 const Real xi = p(0); 00257 const Real eta = p(1); 00258 const Real zeta = p(2); 00259 00260 // The only way to make any sense of this 00261 // is to look at the mgflo/mg2/mgf documentation 00262 // and make the cut-out cube! 00263 // 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 00264 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}; 00265 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}; 00266 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}; 00267 00268 return (FE<1,LAGRANGE>::shape(EDGE3, SECOND, i0[i], xi)* 00269 FE<1,LAGRANGE>::shape(EDGE3, SECOND, i1[i], eta)* 00270 FE<1,LAGRANGE>::shape(EDGE3, SECOND, i2[i], zeta)); 00271 } 00272 00273 // quadratic tetrahedral shape functions 00274 case TET10: 00275 { 00276 libmesh_assert_less (i, 10); 00277 00278 // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM 00279 const Real zeta1 = p(0); 00280 const Real zeta2 = p(1); 00281 const Real zeta3 = p(2); 00282 const Real zeta0 = 1. - zeta1 - zeta2 - zeta3; 00283 00284 switch(i) 00285 { 00286 case 0: 00287 return zeta0*(2.*zeta0 - 1.); 00288 00289 case 1: 00290 return zeta1*(2.*zeta1 - 1.); 00291 00292 case 2: 00293 return zeta2*(2.*zeta2 - 1.); 00294 00295 case 3: 00296 return zeta3*(2.*zeta3 - 1.); 00297 00298 case 4: 00299 return 4.*zeta0*zeta1; 00300 00301 case 5: 00302 return 4.*zeta1*zeta2; 00303 00304 case 6: 00305 return 4.*zeta2*zeta0; 00306 00307 case 7: 00308 return 4.*zeta0*zeta3; 00309 00310 case 8: 00311 return 4.*zeta1*zeta3; 00312 00313 case 9: 00314 return 4.*zeta2*zeta3; 00315 00316 default: 00317 libmesh_error_msg("Invalid i = " << i); 00318 } 00319 } 00320 00321 // "serendipity" prism 00322 case PRISM15: 00323 { 00324 libmesh_assert_less (i, 15); 00325 00326 const Real xi = p(0); 00327 const Real eta = p(1); 00328 const Real zeta = p(2); 00329 00330 switch(i) 00331 { 00332 case 0: 00333 return (1. - zeta)*(xi + eta - 1.)*(xi + eta + 0.5*zeta); 00334 00335 case 1: 00336 return (1. - zeta)*xi*(xi - 1. - 0.5*zeta); 00337 00338 case 2: // phi1 with xi <- eta 00339 return (1. - zeta)*eta*(eta - 1. - 0.5*zeta); 00340 00341 case 3: // phi0 with zeta <- (-zeta) 00342 return (1. + zeta)*(xi + eta - 1.)*(xi + eta - 0.5*zeta); 00343 00344 case 4: // phi1 with zeta <- (-zeta) 00345 return (1. + zeta)*xi*(xi - 1. + 0.5*zeta); 00346 00347 case 5: // phi4 with xi <- eta 00348 return (1. + zeta)*eta*(eta - 1. + 0.5*zeta); 00349 00350 case 6: 00351 return 2.*(1. - zeta)*xi*(1. - xi - eta); 00352 00353 case 7: 00354 return 2.*(1. - zeta)*xi*eta; 00355 00356 case 8: 00357 return 2.*(1. - zeta)*eta*(1. - xi - eta); 00358 00359 case 9: 00360 return (1. - zeta)*(1. + zeta)*(1. - xi - eta); 00361 00362 case 10: 00363 return (1. - zeta)*(1. + zeta)*xi; 00364 00365 case 11: // phi10 with xi <-> eta 00366 return (1. - zeta)*(1. + zeta)*eta; 00367 00368 case 12: // phi6 with zeta <- (-zeta) 00369 return 2.*(1. + zeta)*xi*(1. - xi - eta); 00370 00371 case 13: // phi7 with zeta <- (-zeta) 00372 return 2.*(1. + zeta)*xi*eta; 00373 00374 case 14: // phi8 with zeta <- (-zeta) 00375 return 2.*(1. + zeta)*eta*(1. - xi - eta); 00376 00377 default: 00378 libmesh_error_msg("Invalid i = " << i); 00379 } 00380 } 00381 00382 // quadradic prism shape functions 00383 case PRISM18: 00384 { 00385 libmesh_assert_less (i, 18); 00386 00387 // Compute prism shape functions as a tensor-product 00388 // of a triangle and an edge 00389 00390 Point p2d(p(0),p(1)); 00391 Point p1d(p(2)); 00392 00393 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 00394 static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2}; 00395 static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5}; 00396 00397 return (FE<2,LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)* 00398 FE<1,LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d)); 00399 } 00400 00401 // G. Bedrosian, "Shape functions and integration formulas for 00402 // three-dimensional finite element analysis", Int. J. Numerical 00403 // Methods Engineering, vol 35, p. 95-108, 1992. 00404 case PYRAMID13: 00405 { 00406 libmesh_assert_less (i, 13); 00407 00408 const Real xi = p(0); 00409 const Real eta = p(1); 00410 const Real zeta = p(2); 00411 const Real eps = 1.e-35; 00412 00413 // Denominators are perturbed by epsilon to avoid 00414 // divide-by-zero issues. 00415 Real den = (1. - zeta + eps); 00416 00417 switch(i) 00418 { 00419 case 0: 00420 return 0.25*(-xi - eta - 1.)*((1. - xi)*(1. - eta) - zeta + xi*eta*zeta/den); 00421 00422 case 1: 00423 return 0.25*(-eta + xi - 1.)*((1. + xi)*(1. - eta) - zeta - xi*eta*zeta/den); 00424 00425 case 2: 00426 return 0.25*(xi + eta - 1.)*((1. + xi)*(1. + eta) - zeta + xi*eta*zeta/den); 00427 00428 case 3: 00429 return 0.25*(eta - xi - 1.)*((1. - xi)*(1. + eta) - zeta - xi*eta*zeta/den); 00430 00431 case 4: 00432 return zeta*(2.*zeta - 1.); 00433 00434 case 5: 00435 return 0.5*(1. + xi - zeta)*(1. - xi - zeta)*(1. - eta - zeta)/den; 00436 00437 case 6: 00438 return 0.5*(1. + eta - zeta)*(1. - eta - zeta)*(1. + xi - zeta)/den; 00439 00440 case 7: 00441 return 0.5*(1. + xi - zeta)*(1. - xi - zeta)*(1. + eta - zeta)/den; 00442 00443 case 8: 00444 return 0.5*(1. + eta - zeta)*(1. - eta - zeta)*(1. - xi - zeta)/den; 00445 00446 case 9: 00447 return zeta*(1. - xi - zeta)*(1. - eta - zeta)/den; 00448 00449 case 10: 00450 return zeta*(1. + xi - zeta)*(1. - eta - zeta)/den; 00451 00452 case 11: 00453 return zeta*(1. + eta - zeta)*(1. + xi - zeta)/den; 00454 00455 case 12: 00456 return zeta*(1. - xi - zeta)*(1. + eta - zeta)/den; 00457 00458 default: 00459 libmesh_error_msg("Invalid i = " << i); 00460 } 00461 } 00462 00463 // Quadratic shape functions, as defined in R. Graglia, "Higher order 00464 // bases on pyramidal elements", IEEE Trans Antennas and Propagation, 00465 // vol 47, no 5, May 1999. 00466 case PYRAMID14: 00467 { 00468 libmesh_assert_less (i, 14); 00469 00470 const Real xi = p(0); 00471 const Real eta = p(1); 00472 const Real zeta = p(2); 00473 const Real eps = 1.e-35; 00474 00475 // The "normalized coordinates" defined by Graglia. These are 00476 // the planes which define the faces of the pyramid. 00477 Real 00478 p1 = 0.5*(1. - eta - zeta), // back 00479 p2 = 0.5*(1. + xi - zeta), // left 00480 p3 = 0.5*(1. + eta - zeta), // front 00481 p4 = 0.5*(1. - xi - zeta); // right 00482 00483 // Denominators are perturbed by epsilon to avoid 00484 // divide-by-zero issues. 00485 Real 00486 den = (-1. + zeta + eps), 00487 den2 = den*den; 00488 00489 switch(i) 00490 { 00491 case 0: 00492 return p4*p1*(xi*eta - zeta + zeta*zeta)/den2; 00493 00494 case 1: 00495 return -p1*p2*(xi*eta + zeta - zeta*zeta)/den2; 00496 00497 case 2: 00498 return p2*p3*(xi*eta - zeta + zeta*zeta)/den2; 00499 00500 case 3: 00501 return -p3*p4*(xi*eta + zeta - zeta*zeta)/den2; 00502 00503 case 4: 00504 return zeta*(2.*zeta - 1.); 00505 00506 case 5: 00507 return -4.*p2*p1*p4*eta/den2; 00508 00509 case 6: 00510 return 4.*p1*p2*p3*xi/den2; 00511 00512 case 7: 00513 return 4.*p2*p3*p4*eta/den2; 00514 00515 case 8: 00516 return -4.*p3*p4*p1*xi/den2; 00517 00518 case 9: 00519 return -4.*p1*p4*zeta/den; 00520 00521 case 10: 00522 return -4.*p2*p1*zeta/den; 00523 00524 case 11: 00525 return -4.*p3*p2*zeta/den; 00526 00527 case 12: 00528 return -4.*p4*p3*zeta/den; 00529 00530 case 13: 00531 return 16.*p1*p2*p3*p4/den2; 00532 00533 default: 00534 libmesh_error_msg("Invalid i = " << i); 00535 } 00536 } 00537 00538 00539 default: 00540 libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type); 00541 } 00542 } 00543 00544 00545 // unsupported order 00546 default: 00547 libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order); 00548 } 00549 00550 #endif 00551 00552 libmesh_error_msg("We'll never get here!"); 00553 return 0.; 00554 } 00555 00556 00557 00558 template <> 00559 Real FE<3,LAGRANGE>::shape(const Elem* elem, 00560 const Order order, 00561 const unsigned int i, 00562 const Point& p) 00563 { 00564 libmesh_assert(elem); 00565 00566 // call the orientation-independent shape functions 00567 return FE<3,LAGRANGE>::shape(elem->type(), static_cast<Order>(order + elem->p_level()), i, p); 00568 } 00569 00570 00571 00572 00573 template <> 00574 Real FE<3,LAGRANGE>::shape_deriv(const ElemType type, 00575 const Order order, 00576 const unsigned int i, 00577 const unsigned int j, 00578 const Point& p) 00579 { 00580 #if LIBMESH_DIM == 3 00581 00582 libmesh_assert_less (j, 3); 00583 00584 switch (order) 00585 { 00586 // linear Lagrange shape functions 00587 case FIRST: 00588 { 00589 switch (type) 00590 { 00591 // trilinear hexahedral shape functions 00592 case HEX8: 00593 case HEX20: 00594 case HEX27: 00595 { 00596 libmesh_assert_less (i, 8); 00597 00598 // Compute hex shape functions as a tensor-product 00599 const Real xi = p(0); 00600 const Real eta = p(1); 00601 const Real zeta = p(2); 00602 00603 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0}; 00604 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1}; 00605 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1}; 00606 00607 switch(j) 00608 { 00609 case 0: 00610 return (FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, xi)* 00611 FE<1,LAGRANGE>::shape (EDGE2, FIRST, i1[i], eta)* 00612 FE<1,LAGRANGE>::shape (EDGE2, FIRST, i2[i], zeta)); 00613 00614 case 1: 00615 return (FE<1,LAGRANGE>::shape (EDGE2, FIRST, i0[i], xi)* 00616 FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i1[i], 0, eta)* 00617 FE<1,LAGRANGE>::shape (EDGE2, FIRST, i2[i], zeta)); 00618 00619 case 2: 00620 return (FE<1,LAGRANGE>::shape (EDGE2, FIRST, i0[i], xi)* 00621 FE<1,LAGRANGE>::shape (EDGE2, FIRST, i1[i], eta)* 00622 FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i2[i], 0, zeta)); 00623 00624 default: 00625 libmesh_error_msg("Invalid j = " << j); 00626 } 00627 } 00628 00629 // linear tetrahedral shape functions 00630 case TET4: 00631 case TET10: 00632 { 00633 libmesh_assert_less (i, 4); 00634 00635 // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM 00636 const Real dzeta0dxi = -1.; 00637 const Real dzeta1dxi = 1.; 00638 const Real dzeta2dxi = 0.; 00639 const Real dzeta3dxi = 0.; 00640 00641 const Real dzeta0deta = -1.; 00642 const Real dzeta1deta = 0.; 00643 const Real dzeta2deta = 1.; 00644 const Real dzeta3deta = 0.; 00645 00646 const Real dzeta0dzeta = -1.; 00647 const Real dzeta1dzeta = 0.; 00648 const Real dzeta2dzeta = 0.; 00649 const Real dzeta3dzeta = 1.; 00650 00651 switch (j) 00652 { 00653 // d()/dxi 00654 case 0: 00655 { 00656 switch(i) 00657 { 00658 case 0: 00659 return dzeta0dxi; 00660 00661 case 1: 00662 return dzeta1dxi; 00663 00664 case 2: 00665 return dzeta2dxi; 00666 00667 case 3: 00668 return dzeta3dxi; 00669 00670 default: 00671 libmesh_error_msg("Invalid i = " << i); 00672 } 00673 } 00674 00675 // d()/deta 00676 case 1: 00677 { 00678 switch(i) 00679 { 00680 case 0: 00681 return dzeta0deta; 00682 00683 case 1: 00684 return dzeta1deta; 00685 00686 case 2: 00687 return dzeta2deta; 00688 00689 case 3: 00690 return dzeta3deta; 00691 00692 default: 00693 libmesh_error_msg("Invalid i = " << i); 00694 } 00695 } 00696 00697 // d()/dzeta 00698 case 2: 00699 { 00700 switch(i) 00701 { 00702 case 0: 00703 return dzeta0dzeta; 00704 00705 case 1: 00706 return dzeta1dzeta; 00707 00708 case 2: 00709 return dzeta2dzeta; 00710 00711 case 3: 00712 return dzeta3dzeta; 00713 00714 default: 00715 libmesh_error_msg("Invalid i = " << i); 00716 } 00717 } 00718 00719 default: 00720 libmesh_error_msg("Invalid shape function derivative j = " << j); 00721 } 00722 } 00723 00724 // linear prism shape functions 00725 case PRISM6: 00726 case PRISM15: 00727 case PRISM18: 00728 { 00729 libmesh_assert_less (i, 6); 00730 00731 // Compute prism shape functions as a tensor-product 00732 // of a triangle and an edge 00733 00734 Point p2d(p(0),p(1)); 00735 Point p1d(p(2)); 00736 00737 // 0 1 2 3 4 5 00738 static const unsigned int i0[] = {0, 0, 0, 1, 1, 1}; 00739 static const unsigned int i1[] = {0, 1, 2, 0, 1, 2}; 00740 00741 switch (j) 00742 { 00743 // d()/dxi 00744 case 0: 00745 return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 0, p2d)* 00746 FE<1,LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d)); 00747 00748 // d()/deta 00749 case 1: 00750 return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 1, p2d)* 00751 FE<1,LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d)); 00752 00753 // d()/dzeta 00754 case 2: 00755 return (FE<2,LAGRANGE>::shape(TRI3, FIRST, i1[i], p2d)* 00756 FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, p1d)); 00757 00758 default: 00759 libmesh_error_msg("Invalid shape function derivative j = " << j); 00760 } 00761 } 00762 00763 // linear pyramid shape functions 00764 case PYRAMID5: 00765 case PYRAMID13: 00766 case PYRAMID14: 00767 { 00768 libmesh_assert_less (i, 5); 00769 00770 const Real xi = p(0); 00771 const Real eta = p(1); 00772 const Real zeta = p(2); 00773 const Real eps = 1.e-35; 00774 00775 switch (j) 00776 { 00777 // d/dxi 00778 case 0: 00779 switch(i) 00780 { 00781 case 0: 00782 return .25*(zeta + eta - 1.)/((1. - zeta) + eps); 00783 00784 case 1: 00785 return -.25*(zeta + eta - 1.)/((1. - zeta) + eps); 00786 00787 case 2: 00788 return -.25*(zeta - eta - 1.)/((1. - zeta) + eps); 00789 00790 case 3: 00791 return .25*(zeta - eta - 1.)/((1. - zeta) + eps); 00792 00793 case 4: 00794 return 0; 00795 00796 default: 00797 libmesh_error_msg("Invalid i = " << i); 00798 } 00799 00800 00801 // d/deta 00802 case 1: 00803 switch(i) 00804 { 00805 case 0: 00806 return .25*(zeta + xi - 1.)/((1. - zeta) + eps); 00807 00808 case 1: 00809 return .25*(zeta - xi - 1.)/((1. - zeta) + eps); 00810 00811 case 2: 00812 return -.25*(zeta - xi - 1.)/((1. - zeta) + eps); 00813 00814 case 3: 00815 return -.25*(zeta + xi - 1.)/((1. - zeta) + eps); 00816 00817 case 4: 00818 return 0; 00819 00820 default: 00821 libmesh_error_msg("Invalid i = " << i); 00822 } 00823 00824 00825 // d/dzeta 00826 case 2: 00827 { 00828 // We computed the derivatives with general eps and 00829 // then let eps tend to zero in the numerators... 00830 Real 00831 num = zeta*(2. - zeta) - 1., 00832 den = (1. - zeta + eps)*(1. - zeta + eps); 00833 00834 switch(i) 00835 { 00836 case 0: 00837 case 2: 00838 return .25*(num + xi*eta)/den; 00839 00840 case 1: 00841 case 3: 00842 return .25*(num - xi*eta)/den; 00843 00844 case 4: 00845 return 1.; 00846 00847 default: 00848 libmesh_error_msg("Invalid i = " << i); 00849 } 00850 } 00851 00852 default: 00853 libmesh_error_msg("Invalid j = " << j); 00854 } 00855 } 00856 00857 00858 default: 00859 libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type); 00860 } 00861 } 00862 00863 00864 // quadratic Lagrange shape functions 00865 case SECOND: 00866 { 00867 switch (type) 00868 { 00869 00870 // serendipity hexahedral quadratic shape functions 00871 case HEX20: 00872 { 00873 libmesh_assert_less (i, 20); 00874 00875 const Real xi = p(0); 00876 const Real eta = p(1); 00877 const Real zeta = p(2); 00878 00879 // these functions are defined for (x,y,z) in [0,1]^3 00880 // so transform the locations 00881 const Real x = .5*(xi + 1.); 00882 const Real y = .5*(eta + 1.); 00883 const Real z = .5*(zeta + 1.); 00884 00885 // and don't forget the chain rule! 00886 00887 switch (j) 00888 { 00889 00890 // d/dx*dx/dxi 00891 case 0: 00892 switch (i) 00893 { 00894 case 0: 00895 return .5*(1. - y)*(1. - z)*((1. - x)*(-2.) + 00896 (-1.)*(1. - 2.*x - 2.*y - 2.*z)); 00897 00898 case 1: 00899 return .5*(1. - y)*(1. - z)*(x*(2.) + 00900 (1.)*(2.*x - 2.*y - 2.*z - 1.)); 00901 00902 case 2: 00903 return .5*y*(1. - z)*(x*(2.) + 00904 (1.)*(2.*x + 2.*y - 2.*z - 3.)); 00905 00906 case 3: 00907 return .5*y*(1. - z)*((1. - x)*(-2.) + 00908 (-1.)*(2.*y - 2.*x - 2.*z - 1.)); 00909 00910 case 4: 00911 return .5*(1. - y)*z*((1. - x)*(-2.) + 00912 (-1.)*(2.*z - 2.*x - 2.*y - 1.)); 00913 00914 case 5: 00915 return .5*(1. - y)*z*(x*(2.) + 00916 (1.)*(2.*x - 2.*y + 2.*z - 3.)); 00917 00918 case 6: 00919 return .5*y*z*(x*(2.) + 00920 (1.)*(2.*x + 2.*y + 2.*z - 5.)); 00921 00922 case 7: 00923 return .5*y*z*((1. - x)*(-2.) + 00924 (-1.)*(2.*y - 2.*x + 2.*z - 3.)); 00925 00926 case 8: 00927 return 2.*(1. - y)*(1. - z)*(1. - 2.*x); 00928 00929 case 9: 00930 return 2.*y*(1. - y)*(1. - z); 00931 00932 case 10: 00933 return 2.*y*(1. - z)*(1. - 2.*x); 00934 00935 case 11: 00936 return 2.*y*(1. - y)*(1. - z)*(-1.); 00937 00938 case 12: 00939 return 2.*(1. - y)*z*(1. - z)*(-1.); 00940 00941 case 13: 00942 return 2.*(1. - y)*z*(1. - z); 00943 00944 case 14: 00945 return 2.*y*z*(1. - z); 00946 00947 case 15: 00948 return 2.*y*z*(1. - z)*(-1.); 00949 00950 case 16: 00951 return 2.*(1. - y)*z*(1. - 2.*x); 00952 00953 case 17: 00954 return 2.*y*(1. - y)*z; 00955 00956 case 18: 00957 return 2.*y*z*(1. - 2.*x); 00958 00959 case 19: 00960 return 2.*y*(1. - y)*z*(-1.); 00961 00962 default: 00963 libmesh_error_msg("Invalid i = " << i); 00964 } 00965 00966 00967 // d/dy*dy/deta 00968 case 1: 00969 switch (i) 00970 { 00971 case 0: 00972 return .5*(1. - x)*(1. - z)*((1. - y)*(-2.) + 00973 (-1.)*(1. - 2.*x - 2.*y - 2.*z)); 00974 00975 case 1: 00976 return .5*x*(1. - z)*((1. - y)*(-2.) + 00977 (-1.)*(2.*x - 2.*y - 2.*z - 1.)); 00978 00979 case 2: 00980 return .5*x*(1. - z)*(y*(2.) + 00981 (1.)*(2.*x + 2.*y - 2.*z - 3.)); 00982 00983 case 3: 00984 return .5*(1. - x)*(1. - z)*(y*(2.) + 00985 (1.)*(2.*y - 2.*x - 2.*z - 1.)); 00986 00987 case 4: 00988 return .5*(1. - x)*z*((1. - y)*(-2.) + 00989 (-1.)*(2.*z - 2.*x - 2.*y - 1.)); 00990 00991 case 5: 00992 return .5*x*z*((1. - y)*(-2.) + 00993 (-1.)*(2.*x - 2.*y + 2.*z - 3.)); 00994 00995 case 6: 00996 return .5*x*z*(y*(2.) + 00997 (1.)*(2.*x + 2.*y + 2.*z - 5.)); 00998 00999 case 7: 01000 return .5*(1. - x)*z*(y*(2.) + 01001 (1.)*(2.*y - 2.*x + 2.*z - 3.)); 01002 01003 case 8: 01004 return 2.*x*(1. - x)*(1. - z)*(-1.); 01005 01006 case 9: 01007 return 2.*x*(1. - z)*(1. - 2.*y); 01008 01009 case 10: 01010 return 2.*x*(1. - x)*(1. - z); 01011 01012 case 11: 01013 return 2.*(1. - x)*(1. - z)*(1. - 2.*y); 01014 01015 case 12: 01016 return 2.*(1. - x)*z*(1. - z)*(-1.); 01017 01018 case 13: 01019 return 2.*x*z*(1. - z)*(-1.); 01020 01021 case 14: 01022 return 2.*x*z*(1. - z); 01023 01024 case 15: 01025 return 2.*(1. - x)*z*(1. - z); 01026 01027 case 16: 01028 return 2.*x*(1. - x)*z*(-1.); 01029 01030 case 17: 01031 return 2.*x*z*(1. - 2.*y); 01032 01033 case 18: 01034 return 2.*x*(1. - x)*z; 01035 01036 case 19: 01037 return 2.*(1. - x)*z*(1. - 2.*y); 01038 01039 default: 01040 libmesh_error_msg("Invalid i = " << i); 01041 } 01042 01043 01044 // d/dz*dz/dzeta 01045 case 2: 01046 switch (i) 01047 { 01048 case 0: 01049 return .5*(1. - x)*(1. - y)*((1. - z)*(-2.) + 01050 (-1.)*(1. - 2.*x - 2.*y - 2.*z)); 01051 01052 case 1: 01053 return .5*x*(1. - y)*((1. - z)*(-2.) + 01054 (-1.)*(2.*x - 2.*y - 2.*z - 1.)); 01055 01056 case 2: 01057 return .5*x*y*((1. - z)*(-2.) + 01058 (-1.)*(2.*x + 2.*y - 2.*z - 3.)); 01059 01060 case 3: 01061 return .5*(1. - x)*y*((1. - z)*(-2.) + 01062 (-1.)*(2.*y - 2.*x - 2.*z - 1.)); 01063 01064 case 4: 01065 return .5*(1. - x)*(1. - y)*(z*(2.) + 01066 (1.)*(2.*z - 2.*x - 2.*y - 1.)); 01067 01068 case 5: 01069 return .5*x*(1. - y)*(z*(2.) + 01070 (1.)*(2.*x - 2.*y + 2.*z - 3.)); 01071 01072 case 6: 01073 return .5*x*y*(z*(2.) + 01074 (1.)*(2.*x + 2.*y + 2.*z - 5.)); 01075 01076 case 7: 01077 return .5*(1. - x)*y*(z*(2.) + 01078 (1.)*(2.*y - 2.*x + 2.*z - 3.)); 01079 01080 case 8: 01081 return 2.*x*(1. - x)*(1. - y)*(-1.); 01082 01083 case 9: 01084 return 2.*x*y*(1. - y)*(-1.); 01085 01086 case 10: 01087 return 2.*x*(1. - x)*y*(-1.); 01088 01089 case 11: 01090 return 2.*(1. - x)*y*(1. - y)*(-1.); 01091 01092 case 12: 01093 return 2.*(1. - x)*(1. - y)*(1. - 2.*z); 01094 01095 case 13: 01096 return 2.*x*(1. - y)*(1. - 2.*z); 01097 01098 case 14: 01099 return 2.*x*y*(1. - 2.*z); 01100 01101 case 15: 01102 return 2.*(1. - x)*y*(1. - 2.*z); 01103 01104 case 16: 01105 return 2.*x*(1. - x)*(1. - y); 01106 01107 case 17: 01108 return 2.*x*y*(1. - y); 01109 01110 case 18: 01111 return 2.*x*(1. - x)*y; 01112 01113 case 19: 01114 return 2.*(1. - x)*y*(1. - y); 01115 01116 default: 01117 libmesh_error_msg("Invalid i = " << i); 01118 } 01119 01120 default: 01121 libmesh_error_msg("Invalid shape function derivative j = " << j); 01122 } 01123 } 01124 01125 // triquadraic hexahedral shape funcions 01126 case HEX27: 01127 { 01128 libmesh_assert_less (i, 27); 01129 01130 // Compute hex shape functions as a tensor-product 01131 const Real xi = p(0); 01132 const Real eta = p(1); 01133 const Real zeta = p(2); 01134 01135 // The only way to make any sense of this 01136 // is to look at the mgflo/mg2/mgf documentation 01137 // and make the cut-out cube! 01138 // 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 01139 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}; 01140 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}; 01141 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}; 01142 01143 switch(j) 01144 { 01145 case 0: 01146 return (FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)* 01147 FE<1,LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)* 01148 FE<1,LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta)); 01149 01150 case 1: 01151 return (FE<1,LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)* 01152 FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)* 01153 FE<1,LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta)); 01154 01155 case 2: 01156 return (FE<1,LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)* 01157 FE<1,LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)* 01158 FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta)); 01159 01160 default: 01161 libmesh_error_msg("Invalid j = " << j); 01162 } 01163 } 01164 01165 // quadratic tetrahedral shape functions 01166 case TET10: 01167 { 01168 libmesh_assert_less (i, 10); 01169 01170 // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM 01171 const Real zeta1 = p(0); 01172 const Real zeta2 = p(1); 01173 const Real zeta3 = p(2); 01174 const Real zeta0 = 1. - zeta1 - zeta2 - zeta3; 01175 01176 const Real dzeta0dxi = -1.; 01177 const Real dzeta1dxi = 1.; 01178 const Real dzeta2dxi = 0.; 01179 const Real dzeta3dxi = 0.; 01180 01181 const Real dzeta0deta = -1.; 01182 const Real dzeta1deta = 0.; 01183 const Real dzeta2deta = 1.; 01184 const Real dzeta3deta = 0.; 01185 01186 const Real dzeta0dzeta = -1.; 01187 const Real dzeta1dzeta = 0.; 01188 const Real dzeta2dzeta = 0.; 01189 const Real dzeta3dzeta = 1.; 01190 01191 switch (j) 01192 { 01193 // d()/dxi 01194 case 0: 01195 { 01196 switch(i) 01197 { 01198 case 0: 01199 return (4.*zeta0 - 1.)*dzeta0dxi; 01200 01201 case 1: 01202 return (4.*zeta1 - 1.)*dzeta1dxi; 01203 01204 case 2: 01205 return (4.*zeta2 - 1.)*dzeta2dxi; 01206 01207 case 3: 01208 return (4.*zeta3 - 1.)*dzeta3dxi; 01209 01210 case 4: 01211 return 4.*(zeta0*dzeta1dxi + dzeta0dxi*zeta1); 01212 01213 case 5: 01214 return 4.*(zeta1*dzeta2dxi + dzeta1dxi*zeta2); 01215 01216 case 6: 01217 return 4.*(zeta0*dzeta2dxi + dzeta0dxi*zeta2); 01218 01219 case 7: 01220 return 4.*(zeta0*dzeta3dxi + dzeta0dxi*zeta3); 01221 01222 case 8: 01223 return 4.*(zeta1*dzeta3dxi + dzeta1dxi*zeta3); 01224 01225 case 9: 01226 return 4.*(zeta2*dzeta3dxi + dzeta2dxi*zeta3); 01227 01228 default: 01229 libmesh_error_msg("Invalid i = " << i); 01230 } 01231 } 01232 01233 // d()/deta 01234 case 1: 01235 { 01236 switch(i) 01237 { 01238 case 0: 01239 return (4.*zeta0 - 1.)*dzeta0deta; 01240 01241 case 1: 01242 return (4.*zeta1 - 1.)*dzeta1deta; 01243 01244 case 2: 01245 return (4.*zeta2 - 1.)*dzeta2deta; 01246 01247 case 3: 01248 return (4.*zeta3 - 1.)*dzeta3deta; 01249 01250 case 4: 01251 return 4.*(zeta0*dzeta1deta + dzeta0deta*zeta1); 01252 01253 case 5: 01254 return 4.*(zeta1*dzeta2deta + dzeta1deta*zeta2); 01255 01256 case 6: 01257 return 4.*(zeta0*dzeta2deta + dzeta0deta*zeta2); 01258 01259 case 7: 01260 return 4.*(zeta0*dzeta3deta + dzeta0deta*zeta3); 01261 01262 case 8: 01263 return 4.*(zeta1*dzeta3deta + dzeta1deta*zeta3); 01264 01265 case 9: 01266 return 4.*(zeta2*dzeta3deta + dzeta2deta*zeta3); 01267 01268 default: 01269 libmesh_error_msg("Invalid i = " << i); 01270 } 01271 } 01272 01273 // d()/dzeta 01274 case 2: 01275 { 01276 switch(i) 01277 { 01278 case 0: 01279 return (4.*zeta0 - 1.)*dzeta0dzeta; 01280 01281 case 1: 01282 return (4.*zeta1 - 1.)*dzeta1dzeta; 01283 01284 case 2: 01285 return (4.*zeta2 - 1.)*dzeta2dzeta; 01286 01287 case 3: 01288 return (4.*zeta3 - 1.)*dzeta3dzeta; 01289 01290 case 4: 01291 return 4.*(zeta0*dzeta1dzeta + dzeta0dzeta*zeta1); 01292 01293 case 5: 01294 return 4.*(zeta1*dzeta2dzeta + dzeta1dzeta*zeta2); 01295 01296 case 6: 01297 return 4.*(zeta0*dzeta2dzeta + dzeta0dzeta*zeta2); 01298 01299 case 7: 01300 return 4.*(zeta0*dzeta3dzeta + dzeta0dzeta*zeta3); 01301 01302 case 8: 01303 return 4.*(zeta1*dzeta3dzeta + dzeta1dzeta*zeta3); 01304 01305 case 9: 01306 return 4.*(zeta2*dzeta3dzeta + dzeta2dzeta*zeta3); 01307 01308 default: 01309 libmesh_error_msg("Invalid i = " << i); 01310 } 01311 } 01312 01313 default: 01314 libmesh_error_msg("Invalid j = " << j); 01315 } 01316 } 01317 01318 01319 // "serendipity" prism 01320 case PRISM15: 01321 { 01322 libmesh_assert_less (i, 15); 01323 01324 const Real xi = p(0); 01325 const Real eta = p(1); 01326 const Real zeta = p(2); 01327 01328 switch (j) 01329 { 01330 // d()/dxi 01331 case 0: 01332 { 01333 switch(i) 01334 { 01335 case 0: 01336 return (2.*xi + 2.*eta + 0.5*zeta - 1.)*(1. - zeta); 01337 case 1: 01338 return (2.*xi - 1. - 0.5*zeta)*(1. - zeta); 01339 case 2: 01340 return 0.; 01341 case 3: 01342 return (2.*xi + 2.*eta - 0.5*zeta - 1.)*(1. + zeta); 01343 case 4: 01344 return (2.*xi - 1. + 0.5*zeta)*(1. + zeta); 01345 case 5: 01346 return 0.; 01347 case 6: 01348 return (4.*xi + 2.*eta - 2.)*(zeta - 1.); 01349 case 7: 01350 return -2.*(zeta - 1.)*eta; 01351 case 8: 01352 return 2.*(zeta - 1.)*eta; 01353 case 9: 01354 return (zeta - 1.)*(1. + zeta); 01355 case 10: 01356 return (1. - zeta)*(1. + zeta); 01357 case 11: 01358 return 0.; 01359 case 12: 01360 return (-4.*xi - 2.*eta + 2.)*(1. + zeta); 01361 case 13: 01362 return 2.*(1. + zeta)*eta; 01363 case 14: 01364 return -2.*(1. + zeta)*eta; 01365 default: 01366 libmesh_error_msg("Invalid i = " << i); 01367 } 01368 } 01369 01370 // d()/deta 01371 case 1: 01372 { 01373 switch(i) 01374 { 01375 case 0: 01376 return (2.*xi + 2.*eta + 0.5*zeta - 1.)*(1. - zeta); 01377 case 1: 01378 return 0.; 01379 case 2: 01380 return (2.*eta - 1. - 0.5*zeta)*(1. - zeta); 01381 case 3: 01382 return (2.*xi + 2.*eta - 0.5*zeta - 1.)*(1. + zeta); 01383 case 4: 01384 return 0.; 01385 case 5: 01386 return (2.*eta - 1. + 0.5*zeta)*(1. + zeta); 01387 case 6: 01388 return 2.*(zeta - 1.)*xi; 01389 case 7: 01390 return 2.*(1. - zeta)*xi; 01391 case 8: 01392 return (2.*xi + 4.*eta - 2.)*(zeta - 1.); 01393 case 9: 01394 return (zeta - 1.)*(1. + zeta); 01395 case 10: 01396 return 0.; 01397 case 11: 01398 return (1. - zeta)*(1. + zeta); 01399 case 12: 01400 return -2.*(1. + zeta)*xi; 01401 case 13: 01402 return 2.*(1. + zeta)*xi; 01403 case 14: 01404 return (-2.*xi - 4.*eta + 2.)*(1. + zeta); 01405 01406 default: 01407 libmesh_error_msg("Invalid i = " << i); 01408 } 01409 } 01410 01411 // d()/dzeta 01412 case 2: 01413 { 01414 switch(i) 01415 { 01416 case 0: 01417 return (-xi - eta - zeta + 0.5)*(xi + eta - 1.); 01418 case 1: 01419 return -0.5*xi*(2.*xi - 1. - 2.*zeta); 01420 case 2: 01421 return -0.5*eta*(2.*eta - 1. - 2.*zeta); 01422 case 3: 01423 return (xi + eta - zeta - 0.5)*(xi + eta - 1.); 01424 case 4: 01425 return 0.5*xi*(2.*xi - 1. + 2.*zeta); 01426 case 5: 01427 return 0.5*eta*(2.*eta - 1. + 2.*zeta); 01428 case 6: 01429 return 2.*xi*(xi + eta - 1.); 01430 case 7: 01431 return -2.*xi*eta; 01432 case 8: 01433 return 2.*eta*(xi + eta - 1.); 01434 case 9: 01435 return 2.*zeta*(xi + eta - 1.); 01436 case 10: 01437 return -2.*xi*zeta; 01438 case 11: 01439 return -2.*eta*zeta; 01440 case 12: 01441 return 2.*xi*(1. - xi - eta); 01442 case 13: 01443 return 2.*xi*eta; 01444 case 14: 01445 return 2.*eta*(1. - xi - eta); 01446 01447 default: 01448 libmesh_error_msg("Invalid i = " << i); 01449 } 01450 } 01451 01452 default: 01453 libmesh_error_msg("Invalid j = " << j); 01454 } 01455 } 01456 01457 01458 01459 // quadradic prism shape functions 01460 case PRISM18: 01461 { 01462 libmesh_assert_less (i, 18); 01463 01464 // Compute prism shape functions as a tensor-product 01465 // of a triangle and an edge 01466 01467 Point p2d(p(0),p(1)); 01468 Point p1d(p(2)); 01469 01470 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 01471 static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2}; 01472 static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5}; 01473 01474 switch (j) 01475 { 01476 // d()/dxi 01477 case 0: 01478 return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)* 01479 FE<1,LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d)); 01480 01481 // d()/deta 01482 case 1: 01483 return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)* 01484 FE<1,LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d)); 01485 01486 // d()/dzeta 01487 case 2: 01488 return (FE<2,LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)* 01489 FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d)); 01490 01491 default: 01492 libmesh_error_msg("Invalid shape function derivative j = " << j); 01493 } 01494 } 01495 01496 // G. Bedrosian, "Shape functions and integration formulas for 01497 // three-dimensional finite element analysis", Int. J. Numerical 01498 // Methods Engineering, vol 35, p. 95-108, 1992. 01499 case PYRAMID13: 01500 { 01501 libmesh_assert_less (i, 13); 01502 01503 const Real xi = p(0); 01504 const Real eta = p(1); 01505 const Real zeta = p(2); 01506 const Real eps = 1.e-35; 01507 01508 // Denominators are perturbed by epsilon to avoid 01509 // divide-by-zero issues. 01510 Real 01511 den = (-1. + zeta + eps), 01512 den2 = den*den, 01513 xi2 = xi*xi, 01514 eta2 = eta*eta, 01515 zeta2 = zeta*zeta, 01516 zeta3 = zeta2*zeta; 01517 01518 switch (j) 01519 { 01520 // d/dxi 01521 case 0: 01522 switch(i) 01523 { 01524 case 0: 01525 return 0.25*(-zeta - eta + 2.*eta*zeta - 2.*xi + 2.*zeta*xi + 2.*eta*xi + zeta2 + eta2)/den; 01526 01527 case 1: 01528 return -0.25*(-zeta - eta + 2.*eta*zeta + 2.*xi - 2.*zeta*xi - 2.*eta*xi + zeta2 + eta2)/den; 01529 01530 case 2: 01531 return -0.25*(-zeta + eta - 2.*eta*zeta + 2.*xi - 2.*zeta*xi + 2.*eta*xi + zeta2 + eta2)/den; 01532 01533 case 3: 01534 return 0.25*(-zeta + eta - 2.*eta*zeta - 2.*xi + 2.*zeta*xi - 2.*eta*xi + zeta2 + eta2)/den; 01535 01536 case 4: 01537 return 0.; 01538 01539 case 5: 01540 return -(-1. + eta + zeta)*xi/den; 01541 01542 case 6: 01543 return 0.5*(-1. + eta + zeta)*(1. + eta - zeta)/den; 01544 01545 case 7: 01546 return (1. + eta - zeta)*xi/den; 01547 01548 case 8: 01549 return -0.5*(-1. + eta + zeta)*(1. + eta - zeta)/den; 01550 01551 case 9: 01552 return -(-1. + eta + zeta)*zeta/den; 01553 01554 case 10: 01555 return (-1. + eta + zeta)*zeta/den; 01556 01557 case 11: 01558 return -(1. + eta - zeta)*zeta/den; 01559 01560 case 12: 01561 return (1. + eta - zeta)*zeta/den; 01562 01563 default: 01564 libmesh_error_msg("Invalid i = " << i); 01565 } 01566 01567 // d/deta 01568 case 1: 01569 switch(i) 01570 { 01571 case 0: 01572 return 0.25*(-zeta - 2.*eta + 2.*eta*zeta - xi + 2.*zeta*xi + 2.*eta*xi + zeta2 + xi2)/den; 01573 01574 case 1: 01575 return -0.25*(zeta + 2.*eta - 2.*eta*zeta - xi + 2.*zeta*xi + 2.*eta*xi - zeta2 - xi2)/den; 01576 01577 case 2: 01578 return -0.25*(-zeta + 2.*eta - 2.*eta*zeta + xi - 2.*zeta*xi + 2.*eta*xi + zeta2 + xi2)/den; 01579 01580 case 3: 01581 return 0.25*(zeta - 2.*eta + 2.*eta*zeta + xi - 2.*zeta*xi + 2.*eta*xi - zeta2 - xi2)/den; 01582 01583 case 4: 01584 return 0.; 01585 01586 case 5: 01587 return -0.5*(-1. + xi + zeta)*(1. + xi - zeta)/den; 01588 01589 case 6: 01590 return (1. + xi - zeta)*eta/den; 01591 01592 case 7: 01593 return 0.5*(-1. + xi + zeta)*(1. + xi - zeta)/den; 01594 01595 case 8: 01596 return -(-1. + xi + zeta)*eta/den; 01597 01598 case 9: 01599 return -(-1. + xi + zeta)*zeta/den; 01600 01601 case 10: 01602 return (1. + xi - zeta)*zeta/den; 01603 01604 case 11: 01605 return -(1. + xi - zeta)*zeta/den; 01606 01607 case 12: 01608 return (-1. + xi + zeta)*zeta/den; 01609 01610 default: 01611 libmesh_error_msg("Invalid i = " << i); 01612 } 01613 01614 // d/dzeta 01615 case 2: 01616 { 01617 switch(i) 01618 { 01619 case 0: 01620 return -0.25*(xi + eta + 1.)*(-1. + 2.*zeta - zeta2 + eta*xi)/den2; 01621 01622 case 1: 01623 return 0.25*(eta - xi + 1.)*(1. - 2.*zeta + zeta2 + eta*xi)/den2; 01624 01625 case 2: 01626 return 0.25*(xi + eta - 1.)*(-1. + 2.*zeta - zeta2 + eta*xi)/den2; 01627 01628 case 3: 01629 return -0.25*(eta - xi - 1.)*(1. - 2.*zeta + zeta2 + eta*xi)/den2; 01630 01631 case 4: 01632 return 4.*zeta - 1.; 01633 01634 case 5: 01635 return 0.5*(-2 + eta + 6.*zeta + eta*xi2 + eta*zeta2 - 6.*zeta2 + 2.*zeta3 - 2.*eta*zeta)/den2; 01636 01637 case 6: 01638 return -0.5*(2 - 6.*zeta + xi + xi*zeta2 + eta2*xi + 6.*zeta2 - 2.*zeta3 - 2.*zeta*xi)/den2; 01639 01640 case 7: 01641 return -0.5*(2 + eta - 6.*zeta + eta*xi2 + eta*zeta2 + 6.*zeta2 - 2.*zeta3 - 2.*eta*zeta)/den2; 01642 01643 case 8: 01644 return 0.5*(-2 + 6.*zeta + xi + xi*zeta2 + eta2*xi - 6.*zeta2 + 2.*zeta3 - 2.*zeta*xi)/den2; 01645 01646 case 9: 01647 return (1. - eta - 4.*zeta - xi - xi*zeta2 - eta*zeta2 + eta*xi + 5.*zeta2 - 2.*zeta3 + 2.*eta*zeta + 2.*zeta*xi)/den2; 01648 01649 case 10: 01650 return -(-1. + eta + 4.*zeta - xi - xi*zeta2 + eta*zeta2 + eta*xi - 5.*zeta2 + 2.*zeta3 - 2.*eta*zeta + 2.*zeta*xi)/den2; 01651 01652 case 11: 01653 return (1. + eta - 4.*zeta + xi + xi*zeta2 + eta*zeta2 + eta*xi + 5.*zeta2 - 2.*zeta3 - 2.*eta*zeta - 2.*zeta*xi)/den2; 01654 01655 case 12: 01656 return -(-1. - eta + 4.*zeta + xi + xi*zeta2 - eta*zeta2 + eta*xi - 5.*zeta2 + 2.*zeta3 + 2.*eta*zeta - 2.*zeta*xi)/den2; 01657 01658 default: 01659 libmesh_error_msg("Invalid i = " << i); 01660 } 01661 } 01662 01663 default: 01664 libmesh_error_msg("Invalid j = " << j); 01665 } 01666 } 01667 01668 // Quadratic shape functions, as defined in R. Graglia, "Higher order 01669 // bases on pyramidal elements", IEEE Trans Antennas and Propagation, 01670 // vol 47, no 5, May 1999. 01671 case PYRAMID14: 01672 { 01673 libmesh_assert_less (i, 14); 01674 01675 const Real xi = p(0); 01676 const Real eta = p(1); 01677 const Real zeta = p(2); 01678 const Real eps = 1.e-35; 01679 01680 // The "normalized coordinates" defined by Graglia. These are 01681 // the planes which define the faces of the pyramid. 01682 Real 01683 p1 = 0.5*(1. - eta - zeta), // back 01684 p2 = 0.5*(1. + xi - zeta), // left 01685 p3 = 0.5*(1. + eta - zeta), // front 01686 p4 = 0.5*(1. - xi - zeta); // right 01687 01688 // Denominators are perturbed by epsilon to avoid 01689 // divide-by-zero issues. 01690 Real 01691 den = (-1. + zeta + eps), 01692 den2 = den*den, 01693 den3 = den2*den; 01694 01695 switch (j) 01696 { 01697 // d/dxi 01698 case 0: 01699 switch(i) 01700 { 01701 case 0: 01702 return 0.5*p1*(-xi*eta + zeta - zeta*zeta + 2.*p4*eta)/den2; 01703 01704 case 1: 01705 return -0.5*p1*(xi*eta + zeta - zeta*zeta + 2.*p2*eta)/den2; 01706 01707 case 2: 01708 return 0.5*p3*(xi*eta - zeta + zeta*zeta + 2.*p2*eta)/den2; 01709 01710 case 3: 01711 return -0.5*p3*(-xi*eta - zeta + zeta*zeta + 2.*p4*eta)/den2; 01712 01713 case 4: 01714 return 0.; 01715 01716 case 5: 01717 return 2.*p1*eta*xi/den2; 01718 01719 case 6: 01720 return 2.*p1*p3*(xi + 2.*p2)/den2; 01721 01722 case 7: 01723 return -2.*p3*eta*xi/den2; 01724 01725 case 8: 01726 return -2.*p1*p3*(-xi + 2.*p4)/den2; 01727 01728 case 9: 01729 return 2.*p1*zeta/den; 01730 01731 case 10: 01732 return -2.*p1*zeta/den; 01733 01734 case 11: 01735 return -2.*p3*zeta/den; 01736 01737 case 12: 01738 return 2.*p3*zeta/den; 01739 01740 case 13: 01741 return -8.*p1*p3*xi/den2; 01742 01743 default: 01744 libmesh_error_msg("Invalid i = " << i); 01745 } 01746 01747 // d/deta 01748 case 1: 01749 switch(i) 01750 { 01751 case 0: 01752 return -0.5*p4*(xi*eta - zeta + zeta*zeta - 2.*p1*xi)/den2; 01753 01754 case 1: 01755 return 0.5*p2*(xi*eta + zeta - zeta*zeta - 2.*p1*xi)/den2; 01756 01757 case 2: 01758 return 0.5*p2*(xi*eta - zeta + zeta*zeta + 2.*p3*xi)/den2; 01759 01760 case 3: 01761 return -0.5*p4*(xi*eta + zeta - zeta*zeta + 2.*p3*xi)/den2; 01762 01763 case 4: 01764 return 0.; 01765 01766 case 5: 01767 return 2.*p2*p4*(eta - 2.*p1)/den2; 01768 01769 case 6: 01770 return -2.*p2*xi*eta/den2; 01771 01772 case 7: 01773 return 2.*p2*p4*(eta + 2.*p3)/den2; 01774 01775 case 8: 01776 return 2.*p4*xi*eta/den2; 01777 01778 case 9: 01779 return 2.*p4*zeta/den; 01780 01781 case 10: 01782 return 2.*p2*zeta/den; 01783 01784 case 11: 01785 return -2.*p2*zeta/den; 01786 01787 case 12: 01788 return -2.*p4*zeta/den; 01789 01790 case 13: 01791 return -8.*p2*p4*eta/den2; 01792 01793 default: 01794 libmesh_error_msg("Invalid i = " << i); 01795 } 01796 01797 01798 // d/dzeta 01799 case 2: 01800 { 01801 switch(i) 01802 { 01803 case 0: 01804 return -0.5*p1*(xi*eta - zeta + zeta*zeta)/den2 01805 - 0.5*p4*(xi*eta - zeta + zeta*zeta)/den2 01806 + p4*p1*(2.*zeta - 1)/den2 01807 - 2.*p4*p1*(xi*eta - zeta + zeta*zeta)/den3; 01808 01809 case 1: 01810 return 0.5*p2*(xi*eta + zeta - zeta*zeta)/den2 01811 + 0.5*p1*(xi*eta + zeta - zeta*zeta)/den2 01812 - p1*p2*(1 - 2.*zeta)/den2 01813 + 2.*p1*p2*(xi*eta + zeta - zeta*zeta)/den3; 01814 01815 case 2: 01816 return -0.5*p3*(xi*eta - zeta + zeta*zeta)/den2 01817 - 0.5*p2*(xi*eta - zeta + zeta*zeta)/den2 01818 + p2*p3*(2.*zeta - 1)/den2 01819 - 2.*p2*p3*(xi*eta - zeta + zeta*zeta)/den3; 01820 01821 case 3: 01822 return 0.5*p4*(xi*eta + zeta - zeta*zeta)/den2 01823 + 0.5*p3*(xi*eta + zeta - zeta*zeta)/den2 01824 - p3*p4*(1 - 2.*zeta)/den2 01825 + 2.*p3*p4*(xi*eta + zeta - zeta*zeta)/den3; 01826 01827 case 4: 01828 return 4.*zeta - 1.; 01829 01830 case 5: 01831 return 2.*p4*p1*eta/den2 01832 + 2.*p2*p4*eta/den2 01833 + 2.*p1*p2*eta/den2 01834 + 8.*p2*p1*p4*eta/den3; 01835 01836 case 6: 01837 return -2.*p2*p3*xi/den2 01838 - 2.*p1*p3*xi/den2 01839 - 2.*p1*p2*xi/den2 01840 - 8.*p1*p2*p3*xi/den3; 01841 01842 case 7: 01843 return -2.*p3*p4*eta/den2 01844 - 2.*p2*p4*eta/den2 01845 - 2.*p2*p3*eta/den2 01846 - 8.*p2*p3*p4*eta/den3; 01847 01848 case 8: 01849 return 2.*p4*p1*xi/den2 01850 + 2.*p1*p3*xi/den2 01851 + 2.*p3*p4*xi/den2 01852 + 8.*p3*p4*p1*xi/den3; 01853 01854 case 9: 01855 return 2.*p4*zeta/den 01856 + 2.*p1*zeta/den 01857 - 4.*p1*p4/den 01858 + 4.*p1*p4*zeta/den2; 01859 01860 case 10: 01861 return 2.*p1*zeta/den 01862 + 2.*p2*zeta/den 01863 - 4.*p2*p1/den 01864 + 4.*p2*p1*zeta/den2; 01865 01866 case 11: 01867 return 2.*p2*zeta/den 01868 + 2.*p3*zeta/den 01869 - 4.*p3*p2/den 01870 + 4.*p3*p2*zeta/den2; 01871 01872 case 12: 01873 return 2.*p3*zeta/den 01874 + 2.*p4*zeta/den 01875 - 4.*p4*p3/den 01876 + 4.*p4*p3*zeta/den2; 01877 01878 case 13: 01879 return -8.*p2*p3*p4/den2 01880 - 8.*p3*p4*p1/den2 01881 - 8.*p2*p1*p4/den2 01882 - 8.*p1*p2*p3/den2 01883 - 32.*p1*p2*p3*p4/den3; 01884 01885 default: 01886 libmesh_error_msg("Invalid i = " << i); 01887 } 01888 } 01889 01890 default: 01891 libmesh_error_msg("Invalid j = " << j); 01892 } 01893 } 01894 01895 01896 default: 01897 libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type); 01898 } 01899 } 01900 01901 01902 // unsupported order 01903 default: 01904 libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order); 01905 } 01906 01907 #endif 01908 01909 libmesh_error_msg("We'll never get here!"); 01910 return 0.; 01911 } 01912 01913 01914 01915 template <> 01916 Real FE<3,LAGRANGE>::shape_deriv(const Elem* elem, 01917 const Order order, 01918 const unsigned int i, 01919 const unsigned int j, 01920 const Point& p) 01921 { 01922 libmesh_assert(elem); 01923 01924 // call the orientation-independent shape function derivatives 01925 return FE<3,LAGRANGE>::shape_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p); 01926 } 01927 01928 01929 01930 template <> 01931 Real FE<3,LAGRANGE>::shape_second_deriv(const ElemType type, 01932 const Order order, 01933 const unsigned int i, 01934 const unsigned int j, 01935 const Point& p) 01936 { 01937 #if LIBMESH_DIM == 3 01938 01939 libmesh_assert_less (j, 6); 01940 01941 switch (order) 01942 { 01943 // linear Lagrange shape functions 01944 case FIRST: 01945 { 01946 switch (type) 01947 { 01948 // Linear tets have all second derivatives = 0 01949 case TET4: 01950 case TET10: 01951 { 01952 return 0.; 01953 } 01954 01955 // The following elements use either tensor product or 01956 // rational basis functions, and therefore probably have 01957 // second derivatives, but we have not implemented them 01958 // yet... 01959 case PRISM6: 01960 case PRISM15: 01961 case PRISM18: 01962 { 01963 libmesh_assert_less (i, 6); 01964 01965 // Compute prism shape functions as a tensor-product 01966 // of a triangle and an edge 01967 01968 Point p2d(p(0),p(1)); 01969 Point p1d(p(2)); 01970 01971 // 0 1 2 3 4 5 01972 static const unsigned int i0[] = {0, 0, 0, 1, 1, 1}; 01973 static const unsigned int i1[] = {0, 1, 2, 0, 1, 2}; 01974 01975 switch (j) 01976 { 01977 // All repeated second derivatives and the xi-eta derivative are zero on PRISMs 01978 case 0: // d^2()/dxi^2 01979 case 1: // d^2()/dxideta 01980 case 2: // d^2()/deta^2 01981 case 5: // d^2()/dzeta^2 01982 { 01983 return 0.; 01984 } 01985 01986 case 3: // d^2()/dxidzeta 01987 return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 0, p2d)* 01988 FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, p1d)); 01989 01990 case 4: // d^2()/detadzeta 01991 return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 1, p2d)* 01992 FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, p1d)); 01993 01994 default: 01995 libmesh_error_msg("Invalid j = " << j); 01996 } 01997 } 01998 01999 case PYRAMID5: 02000 case PYRAMID13: 02001 case PYRAMID14: 02002 { 02003 libmesh_assert_less (i, 5); 02004 02005 const Real xi = p(0); 02006 const Real eta = p(1); 02007 const Real zeta = p(2); 02008 const Real eps = 1.e-35; 02009 02010 switch (j) 02011 { 02012 // xi-xi and eta-eta derivatives are all zero for PYRAMID5. 02013 case 0: // d^2()/dxi^2 02014 case 2: // d^2()/deta^2 02015 return 0.; 02016 02017 case 1: // d^2()/dxideta 02018 { 02019 switch (i) 02020 { 02021 case 0: 02022 case 2: 02023 return 0.25/(1. - zeta + eps); 02024 case 1: 02025 case 3: 02026 return -0.25/(1. - zeta + eps); 02027 case 4: 02028 return 0.; 02029 default: 02030 libmesh_error_msg("Invalid i = " << i); 02031 } 02032 } 02033 02034 case 3: // d^2()/dxidzeta 02035 { 02036 Real den = (1. - zeta + eps)*(1. - zeta + eps); 02037 02038 switch (i) 02039 { 02040 case 0: 02041 case 2: 02042 return 0.25*eta/den; 02043 case 1: 02044 case 3: 02045 return -0.25*eta/den; 02046 case 4: 02047 return 0.; 02048 default: 02049 libmesh_error_msg("Invalid i = " << i); 02050 } 02051 } 02052 02053 case 4: // d^2()/detadzeta 02054 { 02055 Real den = (1. - zeta + eps)*(1. - zeta + eps); 02056 02057 switch (i) 02058 { 02059 case 0: 02060 case 2: 02061 return 0.25*xi/den; 02062 case 1: 02063 case 3: 02064 return -0.25*xi/den; 02065 case 4: 02066 return 0.; 02067 default: 02068 libmesh_error_msg("Invalid i = " << i); 02069 } 02070 } 02071 02072 case 5: // d^2()/dzeta^2 02073 { 02074 Real den = (1. - zeta + eps)*(1. - zeta + eps)*(1. - zeta + eps); 02075 02076 switch (i) 02077 { 02078 case 0: 02079 case 2: 02080 return 0.5*xi*eta/den; 02081 case 1: 02082 case 3: 02083 return -0.5*xi*eta/den; 02084 case 4: 02085 return 0.; 02086 default: 02087 libmesh_error_msg("Invalid i = " << i); 02088 } 02089 } 02090 02091 default: 02092 libmesh_error_msg("Invalid j = " << j); 02093 } 02094 } 02095 02096 // Trilinear shape functions on HEX8s have nonzero mixed second derivatives 02097 case HEX8: 02098 case HEX20: 02099 case HEX27: 02100 { 02101 libmesh_assert_less (i, 8); 02102 02103 // Compute hex shape functions as a tensor-product 02104 const Real xi = p(0); 02105 const Real eta = p(1); 02106 const Real zeta = p(2); 02107 02108 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0}; 02109 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1}; 02110 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1}; 02111 02112 switch (j) 02113 { 02114 // All repeated second derivatives are zero on HEX8 02115 case 0: // d^2()/dxi^2 02116 case 2: // d^2()/deta^2 02117 case 5: // d^2()/dzeta^2 02118 { 02119 return 0.; 02120 } 02121 02122 case 1: // d^2()/dxideta 02123 return (FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, xi)* 02124 FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i1[i], 0, eta)* 02125 FE<1,LAGRANGE>::shape (EDGE2, FIRST, i2[i], zeta)); 02126 02127 case 3: // d^2()/dxidzeta 02128 return (FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, xi)* 02129 FE<1,LAGRANGE>::shape (EDGE2, FIRST, i1[i], eta)* 02130 FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i2[i], 0, zeta)); 02131 02132 case 4: // d^2()/detadzeta 02133 return (FE<1,LAGRANGE>::shape (EDGE2, FIRST, i0[i], xi)* 02134 FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i1[i], 0, eta)* 02135 FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i2[i], 0, zeta)); 02136 02137 default: 02138 libmesh_error_msg("Invalid j = " << j); 02139 } 02140 } 02141 02142 default: 02143 libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type); 02144 } 02145 02146 } 02147 02148 // quadratic Lagrange shape functions 02149 case SECOND: 02150 { 02151 switch (type) 02152 { 02153 02154 // serendipity hexahedral quadratic shape functions 02155 case HEX20: 02156 { 02157 libmesh_assert_less (i, 20); 02158 02159 const Real xi = p(0); 02160 const Real eta = p(1); 02161 const Real zeta = p(2); 02162 02163 // these functions are defined for (x,y,z) in [0,1]^3 02164 // so transform the locations 02165 const Real x = .5*(xi + 1.); 02166 const Real y = .5*(eta + 1.); 02167 const Real z = .5*(zeta + 1.); 02168 02169 switch(j) 02170 { 02171 case 0: // d^2()/dxi^2 02172 { 02173 switch(i) 02174 { 02175 case 0: 02176 case 1: 02177 return (1. - y) * (1. - z); 02178 case 2: 02179 case 3: 02180 return y * (1. - z); 02181 case 4: 02182 case 5: 02183 return (1. - y) * z; 02184 case 6: 02185 case 7: 02186 return y * z; 02187 case 8: 02188 return -2. * (1. - y) * (1. - z); 02189 case 10: 02190 return -2. * y * (1. - z); 02191 case 16: 02192 return -2. * (1. - y) * z; 02193 case 18: 02194 return -2. * y * z; 02195 case 9: 02196 case 11: 02197 case 12: 02198 case 13: 02199 case 14: 02200 case 15: 02201 case 17: 02202 case 19: 02203 return 0; 02204 default: 02205 libmesh_error_msg("Invalid i = " << i); 02206 } 02207 } 02208 case 1: // d^2()/dxideta 02209 { 02210 switch(i) 02211 { 02212 case 0: 02213 return (1.25 - x - y - .5*z) * (1. - z); 02214 case 1: 02215 return (-x + y + .5*z - .25) * (1. - z); 02216 case 2: 02217 return (x + y - .5*z - .75) * (1. - z); 02218 case 3: 02219 return (-y + x + .5*z - .25) * (1. - z); 02220 case 4: 02221 return -.25*z * (4.*x + 4.*y - 2.*z - 3); 02222 case 5: 02223 return -.25*z * (-4.*y + 4.*x + 2.*z - 1.); 02224 case 6: 02225 return .25*z * (-5 + 4.*x + 4.*y + 2.*z); 02226 case 7: 02227 return .25*z * (4.*x - 4.*y - 2.*z + 1.); 02228 case 8: 02229 return (-1. + 2.*x) * (1. - z); 02230 case 9: 02231 return (1. - 2.*y) * (1. - z); 02232 case 10: 02233 return (1. - 2.*x) * (1. - z); 02234 case 11: 02235 return (-1. + 2.*y) * (1. - z); 02236 case 12: 02237 return z * (1. - z); 02238 case 13: 02239 return -z * (1. - z); 02240 case 14: 02241 return z * (1. - z); 02242 case 15: 02243 return -z * (1. - z); 02244 case 16: 02245 return (-1. + 2.*x) * z; 02246 case 17: 02247 return (1. - 2.*y) * z; 02248 case 18: 02249 return (1. - 2.*x) * z; 02250 case 19: 02251 return (-1. + 2.*y) * z; 02252 default: 02253 libmesh_error_msg("Invalid i = " << i); 02254 } 02255 } 02256 case 2: // d^2()/deta^2 02257 switch(i) 02258 { 02259 case 0: 02260 case 3: 02261 return (1. - x) * (1. - z); 02262 case 1: 02263 case 2: 02264 return x * (1. - z); 02265 case 4: 02266 case 7: 02267 return (1. - x) * z; 02268 case 5: 02269 case 6: 02270 return x * z; 02271 case 9: 02272 return -2. * x * (1. - z); 02273 case 11: 02274 return -2. * (1. - x) * (1. - z); 02275 case 17: 02276 return -2. * x * z; 02277 case 19: 02278 return -2. * (1. - x) * z; 02279 case 8: 02280 case 10: 02281 case 12: 02282 case 13: 02283 case 14: 02284 case 15: 02285 case 16: 02286 case 18: 02287 return 0.; 02288 default: 02289 libmesh_error_msg("Invalid i = " << i); 02290 } 02291 case 3: // d^2()/dxidzeta 02292 switch(i) 02293 { 02294 case 0: 02295 return (1.25 - x - .5*y - z) * (1. - y); 02296 case 1: 02297 return (-x + .5*y + z - .25) * (1. - y); 02298 case 2: 02299 return -.25*y * (2.*y + 4.*x - 4.*z - 1.); 02300 case 3: 02301 return -.25*y * (-2.*y + 4.*x + 4.*z - 3); 02302 case 4: 02303 return (-z + x + .5*y - .25) * (1. - y); 02304 case 5: 02305 return (x - .5*y + z - .75) * (1. - y); 02306 case 6: 02307 return .25*y * (2.*y + 4.*x + 4.*z - 5); 02308 case 7: 02309 return .25*y * (-2.*y + 4.*x - 4.*z + 1.); 02310 case 8: 02311 return (-1. + 2.*x) * (1. - y); 02312 case 9: 02313 return -y * (1. - y); 02314 case 10: 02315 return (-1. + 2.*x) * y; 02316 case 11: 02317 return y * (1. - y); 02318 case 12: 02319 return (-1. + 2.*z) * (1. - y); 02320 case 13: 02321 return (1. - 2.*z) * (1. - y); 02322 case 14: 02323 return (1. - 2.*z) * y; 02324 case 15: 02325 return (-1. + 2.*z) * y; 02326 case 16: 02327 return (1. - 2.*x) * (1. - y); 02328 case 17: 02329 return y * (1. - y); 02330 case 18: 02331 return (1. - 2.*x) * y; 02332 case 19: 02333 return -y * (1. - y); 02334 default: 02335 libmesh_error_msg("Invalid i = " << i); 02336 } 02337 case 4: // d^2()/detadzeta 02338 switch(i) 02339 { 02340 case 0: 02341 return (1.25 - .5*x - y - z) * (1. - x); 02342 case 1: 02343 return .25*x * (2.*x - 4.*y - 4.*z + 3.); 02344 case 2: 02345 return -.25*x * (2.*x + 4.*y - 4.*z - 1.); 02346 case 3: 02347 return (-y + .5*x + z - .25) * (1. - x); 02348 case 4: 02349 return (-z + .5*x + y - .25) * (1. - x); 02350 case 5: 02351 return -.25*x * (2.*x - 4.*y + 4.*z - 1.); 02352 case 6: 02353 return .25*x * (2.*x + 4.*y + 4.*z - 5); 02354 case 7: 02355 return (y - .5*x + z - .75) * (1. - x); 02356 case 8: 02357 return x * (1. - x); 02358 case 9: 02359 return (-1. + 2.*y) * x; 02360 case 10: 02361 return -x * (1. - x); 02362 case 11: 02363 return (-1. + 2.*y) * (1. - x); 02364 case 12: 02365 return (-1. + 2.*z) * (1. - x); 02366 case 13: 02367 return (-1. + 2.*z) * x; 02368 case 14: 02369 return (1. - 2.*z) * x; 02370 case 15: 02371 return (1. - 2.*z) * (1. - x); 02372 case 16: 02373 return -x * (1. - x); 02374 case 17: 02375 return (1. - 2.*y) * x; 02376 case 18: 02377 return x * (1. - x); 02378 case 19: 02379 return (1. - 2.*y) * (1. - x); 02380 default: 02381 libmesh_error_msg("Invalid i = " << i); 02382 } 02383 case 5: // d^2()/dzeta^2 02384 switch(i) 02385 { 02386 case 0: 02387 case 4: 02388 return (1. - x) * (1. - y); 02389 case 1: 02390 case 5: 02391 return x * (1. - y); 02392 case 2: 02393 case 6: 02394 return x * y; 02395 case 3: 02396 case 7: 02397 return (1. - x) * y; 02398 case 12: 02399 return -2. * (1. - x) * (1. - y); 02400 case 13: 02401 return -2. * x * (1. - y); 02402 case 14: 02403 return -2. * x * y; 02404 case 15: 02405 return -2. * (1. - x) * y; 02406 case 8: 02407 case 9: 02408 case 10: 02409 case 11: 02410 case 16: 02411 case 17: 02412 case 18: 02413 case 19: 02414 return 0.; 02415 default: 02416 libmesh_error_msg("Invalid i = " << i); 02417 } 02418 default: 02419 libmesh_error_msg("Invalid j = " << j); 02420 } 02421 } 02422 02423 // triquadraic hexahedral shape funcions 02424 case HEX27: 02425 { 02426 libmesh_assert_less (i, 27); 02427 02428 // Compute hex shape functions as a tensor-product 02429 const Real xi = p(0); 02430 const Real eta = p(1); 02431 const Real zeta = p(2); 02432 02433 // The only way to make any sense of this 02434 // is to look at the mgflo/mg2/mgf documentation 02435 // and make the cut-out cube! 02436 // 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 02437 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}; 02438 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}; 02439 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}; 02440 02441 switch(j) 02442 { 02443 // d^2()/dxi^2 02444 case 0: 02445 return (FE<1,LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i0[i], 0, xi)* 02446 FE<1,LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)* 02447 FE<1,LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta)); 02448 02449 // d^2()/dxideta 02450 case 1: 02451 return (FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)* 02452 FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)* 02453 FE<1,LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta)); 02454 02455 // d^2()/deta^2 02456 case 2: 02457 return (FE<1,LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)* 02458 FE<1,LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i1[i], 0, eta)* 02459 FE<1,LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta)); 02460 02461 // d^2()/dxidzeta 02462 case 3: 02463 return (FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)* 02464 FE<1,LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)* 02465 FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta)); 02466 02467 // d^2()/detadzeta 02468 case 4: 02469 return (FE<1,LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)* 02470 FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)* 02471 FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta)); 02472 02473 // d^2()/dzeta^2 02474 case 5: 02475 return (FE<1,LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)* 02476 FE<1,LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)* 02477 FE<1,LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i2[i], 0, zeta)); 02478 02479 default: 02480 libmesh_error_msg("Invalid j = " << j); 02481 } 02482 } 02483 02484 // quadratic tetrahedral shape functions 02485 case TET10: 02486 { 02487 // The area coordinates are the same as used for the 02488 // shape() and shape_deriv() functions. 02489 // const Real zeta0 = 1. - zeta1 - zeta2 - zeta3; 02490 // const Real zeta1 = p(0); 02491 // const Real zeta2 = p(1); 02492 // const Real zeta3 = p(2); 02493 static const Real dzetadxi[4][3] = 02494 { 02495 {-1., -1., -1.}, 02496 {1., 0., 0.}, 02497 {0., 1., 0.}, 02498 {0., 0., 1.} 02499 }; 02500 02501 // Convert from j -> (j,k) indices for independent variable 02502 // (0=xi, 1=eta, 2=zeta) 02503 static const unsigned short int independent_var_indices[6][2] = 02504 { 02505 {0, 0}, // d^2 phi / dxi^2 02506 {0, 1}, // d^2 phi / dxi deta 02507 {1, 1}, // d^2 phi / deta^2 02508 {0, 2}, // d^2 phi / dxi dzeta 02509 {1, 2}, // d^2 phi / deta dzeta 02510 {2, 2} // d^2 phi / dzeta^2 02511 }; 02512 02513 // Convert from i -> zeta indices. Each quadratic shape 02514 // function for the Tet10 depends on up to two of the zeta 02515 // area coordinate functions (see the shape() function above). 02516 // This table just tells which two area coords it uses. 02517 static const unsigned short int zeta_indices[10][2] = 02518 { 02519 {0, 0}, 02520 {1, 1}, 02521 {2, 2}, 02522 {3, 3}, 02523 {0, 1}, 02524 {1, 2}, 02525 {2, 0}, 02526 {0, 3}, 02527 {1, 3}, 02528 {2, 3}, 02529 }; 02530 02531 // Look up the independent variable indices for this value of j. 02532 const unsigned int my_j = independent_var_indices[j][0]; 02533 const unsigned int my_k = independent_var_indices[j][1]; 02534 02535 if (i<4) 02536 { 02537 return 4.*dzetadxi[i][my_j]*dzetadxi[i][my_k]; 02538 } 02539 02540 else if (i<10) 02541 { 02542 const unsigned short int my_m = zeta_indices[i][0]; 02543 const unsigned short int my_n = zeta_indices[i][1]; 02544 02545 return 4.*(dzetadxi[my_n][my_j]*dzetadxi[my_m][my_k] + 02546 dzetadxi[my_m][my_j]*dzetadxi[my_n][my_k] ); 02547 } 02548 else 02549 libmesh_error_msg("Invalid shape function index " << i); 02550 } 02551 02552 02553 02554 // "serendipity" prism 02555 case PRISM15: 02556 { 02557 libmesh_assert_less (i, 15); 02558 02559 const Real xi = p(0); 02560 const Real eta = p(1); 02561 const Real zeta = p(2); 02562 02563 switch (j) 02564 { 02565 // d^2()/dxi^2 02566 case 0: 02567 { 02568 switch(i) 02569 { 02570 case 0: 02571 case 1: 02572 return 2.*(1. - zeta); 02573 case 2: 02574 case 5: 02575 case 7: 02576 case 8: 02577 case 9: 02578 case 10: 02579 case 11: 02580 case 13: 02581 case 14: 02582 return 0.; 02583 case 3: 02584 case 4: 02585 return 2.*(1. + zeta); 02586 case 6: 02587 return 4.*(zeta - 1); 02588 case 12: 02589 return -4.*(1. + zeta); 02590 default: 02591 libmesh_error_msg("Invalid i = " << i); 02592 } 02593 } 02594 02595 // d^2()/dxideta 02596 case 1: 02597 { 02598 switch(i) 02599 { 02600 case 0: 02601 case 7: 02602 return 2.*(1. - zeta); 02603 case 1: 02604 case 2: 02605 case 4: 02606 case 5: 02607 case 9: 02608 case 10: 02609 case 11: 02610 return 0.; 02611 case 3: 02612 case 13: 02613 return 2.*(1. + zeta); 02614 case 6: 02615 case 8: 02616 return 2.*(zeta - 1.); 02617 case 12: 02618 case 14: 02619 return -2.*(1. + zeta); 02620 default: 02621 libmesh_error_msg("Invalid i = " << i); 02622 } 02623 } 02624 02625 // d^2()/deta^2 02626 case 2: 02627 { 02628 switch(i) 02629 { 02630 case 0: 02631 case 2: 02632 return 2.*(1. - zeta); 02633 case 1: 02634 case 4: 02635 case 6: 02636 case 7: 02637 case 9: 02638 case 10: 02639 case 11: 02640 case 12: 02641 case 13: 02642 return 0.; 02643 case 3: 02644 case 5: 02645 return 2.*(1. + zeta); 02646 case 8: 02647 return 4.*(zeta - 1.); 02648 case 14: 02649 return -4.*(1. + zeta); 02650 default: 02651 libmesh_error_msg("Invalid i = " << i); 02652 } 02653 } 02654 02655 // d^2()/dxidzeta 02656 case 3: 02657 { 02658 switch(i) 02659 { 02660 case 0: 02661 return 1.5 - zeta - 2.*xi - 2.*eta; 02662 case 1: 02663 return 0.5 + zeta - 2.*xi; 02664 case 2: 02665 case 5: 02666 case 11: 02667 return 0.; 02668 case 3: 02669 return -1.5 - zeta + 2.*xi + 2.*eta; 02670 case 4: 02671 return -0.5 + zeta + 2.*xi; 02672 case 6: 02673 return 4.*xi + 2.*eta - 2.; 02674 case 7: 02675 return -2.*eta; 02676 case 8: 02677 return 2.*eta; 02678 case 9: 02679 return 2.*zeta; 02680 case 10: 02681 return -2.*zeta; 02682 case 12: 02683 return -4.*xi - 2.*eta + 2.; 02684 case 13: 02685 return 2.*eta; 02686 case 14: 02687 return -2.*eta; 02688 default: 02689 libmesh_error_msg("Invalid i = " << i); 02690 } 02691 } 02692 02693 // d^2()/detadzeta 02694 case 4: 02695 { 02696 switch(i) 02697 { 02698 case 0: 02699 return 1.5 - zeta - 2.*xi - 2.*eta; 02700 case 1: 02701 case 4: 02702 case 10: 02703 return 0.; 02704 case 2: 02705 return .5 + zeta - 2.*eta; 02706 case 3: 02707 return -1.5 - zeta + 2.*xi + 2.*eta; 02708 case 5: 02709 return -.5 + zeta + 2.*eta; 02710 case 6: 02711 return 2.*xi; 02712 case 7: 02713 return -2.*xi; 02714 case 8: 02715 return 2.*xi + 4.*eta - 2.; 02716 case 9: 02717 return 2.*zeta; 02718 case 11: 02719 return -2.*zeta; 02720 case 12: 02721 return -2.*xi; 02722 case 13: 02723 return 2.*xi; 02724 case 14: 02725 return -2.*xi - 4.*eta + 2.; 02726 default: 02727 libmesh_error_msg("Invalid i = " << i); 02728 } 02729 } 02730 02731 // d^2()/dzeta^2 02732 case 5: 02733 { 02734 switch(i) 02735 { 02736 case 0: 02737 case 3: 02738 return 1. - xi - eta; 02739 case 1: 02740 case 4: 02741 return xi; 02742 case 2: 02743 case 5: 02744 return eta; 02745 case 6: 02746 case 7: 02747 case 8: 02748 case 12: 02749 case 13: 02750 case 14: 02751 return 0.; 02752 case 9: 02753 return 2.*xi + 2.*eta - 2.; 02754 case 10: 02755 return -2.*xi; 02756 case 11: 02757 return -2.*eta; 02758 default: 02759 libmesh_error_msg("Invalid i = " << i); 02760 } 02761 } 02762 02763 default: 02764 libmesh_error_msg("Invalid j = " << j); 02765 } 02766 } 02767 02768 02769 02770 // quadradic prism shape functions 02771 case PRISM18: 02772 { 02773 libmesh_assert_less (i, 18); 02774 02775 // Compute prism shape functions as a tensor-product 02776 // of a triangle and an edge 02777 02778 Point p2d(p(0),p(1)); 02779 Point p1d(p(2)); 02780 02781 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 02782 static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2}; 02783 static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5}; 02784 02785 switch (j) 02786 { 02787 // d^2()/dxi^2 02788 case 0: 02789 return (FE<2,LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 0, p2d)* 02790 FE<1,LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d)); 02791 02792 // d^2()/dxideta 02793 case 1: 02794 return (FE<2,LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 1, p2d)* 02795 FE<1,LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d)); 02796 02797 // d^2()/deta^2 02798 case 2: 02799 return (FE<2,LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 2, p2d)* 02800 FE<1,LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d)); 02801 02802 // d^2()/dxidzeta 02803 case 3: 02804 return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)* 02805 FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d)); 02806 02807 // d^2()/detadzeta 02808 case 4: 02809 return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)* 02810 FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d)); 02811 02812 // d^2()/dzeta^2 02813 case 5: 02814 return (FE<2,LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)* 02815 FE<1,LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i0[i], 0, p1d)); 02816 02817 default: 02818 libmesh_error_msg("Invalid shape function derivative j = " << j); 02819 } 02820 } 02821 02822 02823 // Quadratic shape functions, as defined in R. Graglia, "Higher order 02824 // bases on pyramidal elements", IEEE Trans Antennas and Propagation, 02825 // vol 47, no 5, May 1999. 02826 case PYRAMID14: 02827 { 02828 libmesh_assert_less (i, 14); 02829 02830 const Real xi = p(0); 02831 const Real eta = p(1); 02832 const Real zeta = p(2); 02833 const Real eps = 1.e-35; 02834 02835 // The "normalized coordinates" defined by Graglia. These are 02836 // the planes which define the faces of the pyramid. 02837 Real 02838 p1 = 0.5*(1. - eta - zeta), // back 02839 p2 = 0.5*(1. + xi - zeta), // left 02840 p3 = 0.5*(1. + eta - zeta), // front 02841 p4 = 0.5*(1. - xi - zeta); // right 02842 02843 // Denominators are perturbed by epsilon to avoid 02844 // divide-by-zero issues. 02845 Real 02846 den = (-1. + zeta + eps), 02847 den2 = den*den, 02848 den3 = den2*den, 02849 den4 = den2*den2; 02850 02851 // These terms are used in several of the derivatives 02852 Real 02853 numer_mp = xi*eta - zeta + zeta*zeta, 02854 numer_pm = xi*eta + zeta - zeta*zeta; 02855 02856 switch (j) 02857 { 02858 case 0: // d^2()/dxi^2 02859 { 02860 switch(i) 02861 { 02862 case 0: 02863 case 1: 02864 return -p1*eta/den2; 02865 case 2: 02866 case 3: 02867 return p3*eta/den2; 02868 case 4: 02869 case 9: 02870 case 10: 02871 case 11: 02872 case 12: 02873 return 0.; 02874 case 5: 02875 return 2.*p1*eta/den2; 02876 case 6: 02877 case 8: 02878 return 4.*p1*p3/den2; 02879 case 7: 02880 return -2.*p3*eta/den2; 02881 case 13: 02882 return -8.*p1*p3/den2; 02883 default: 02884 libmesh_error_msg("Invalid i = " << i); 02885 } 02886 } 02887 02888 case 1: // d^2()/dxideta 02889 { 02890 switch(i) 02891 { 02892 case 0: 02893 return 0.25*numer_mp/den2 02894 - 0.5*p1*xi/den2 02895 - 0.5*p4*eta/den2 02896 + p4*p1/den2; 02897 02898 case 1: 02899 return 0.25*numer_pm/den2 02900 - 0.5*p1*xi/den2 02901 + 0.5*p2*eta/den2 02902 - p1*p2/den2; 02903 02904 case 2: 02905 return 0.25*numer_mp/den2 02906 + 0.5*p3*xi/den2 02907 + 0.5*p2*eta/den2 02908 + p2*p3/den2; 02909 02910 case 3: 02911 return 0.25*numer_pm/den2 02912 + 0.5*p3*xi/den2 02913 - 0.5*p4*eta/den2 02914 - p3*p4/den2; 02915 02916 case 4: 02917 return 0.; 02918 02919 case 5: 02920 return p4*eta/den2 02921 - 2.*p4*p1/den2 02922 - p2*eta/den2 02923 + 2.*p1*p2/den2; 02924 02925 case 6: 02926 return -p3*xi/den2 02927 + p1*xi/den2 02928 - 2.*p2*p3/den2 02929 + 2.*p1*p2/den2; 02930 02931 case 7: 02932 return p4*eta/den2 02933 + 2.*p3*p4/den2 02934 - p2*eta/den2 02935 - 2.*p2*p3/den2; 02936 02937 case 8: 02938 return -p3*xi/den2 02939 + p1*xi/den2 02940 - 2.*p4*p1/den2 02941 + 2.*p3*p4/den2; 02942 02943 case 9: 02944 case 11: 02945 return -zeta/den; 02946 02947 case 10: 02948 case 12: 02949 return zeta/den; 02950 02951 case 13: 02952 return 4.*p4*p1/den2 02953 - 4.*p3*p4/den2 02954 + 4.*p2*p3/den2 02955 - 4.*p1*p2/den2; 02956 02957 default: 02958 libmesh_error_msg("Invalid i = " << i); 02959 } 02960 } 02961 02962 02963 case 2: // d^2()/deta^2 02964 { 02965 switch(i) 02966 { 02967 case 0: 02968 case 3: 02969 return -p4*xi/den2; 02970 case 1: 02971 case 2: 02972 return p2*xi/den2; 02973 case 4: 02974 case 9: 02975 case 10: 02976 case 11: 02977 case 12: 02978 return 0.; 02979 case 5: 02980 case 7: 02981 return 4.*p2*p4/den2; 02982 case 6: 02983 return -2.*p2*xi/den2; 02984 case 8: 02985 return 2.*p4*xi/den2; 02986 case 13: 02987 return -8.*p2*p4/den2; 02988 default: 02989 libmesh_error_msg("Invalid i = " << i); 02990 } 02991 } 02992 02993 02994 case 3: // d^2()/dxidzeta 02995 { 02996 switch(i) 02997 { 02998 case 0: 02999 return 0.25*numer_mp/den2 03000 - 0.5*p1*(2.*zeta - 1.)/den2 03001 + p1*numer_mp/den3 03002 - 0.5*p1*eta/den2 03003 - 0.5*p4*eta/den2 03004 - 2.*p4*p1*eta/den3; 03005 03006 case 1: 03007 return 0.25*numer_pm/den2 03008 - 0.5*p1*(1 - 2.*zeta)/den2 03009 + p1*numer_pm/den3 03010 + 0.5*p2*eta/den2 03011 + 0.5*p1*eta/den2 03012 + 2.*p1*p2*eta/den3; 03013 03014 case 2: 03015 return -0.25*numer_mp/den2 03016 + 0.5*p3*(2.*zeta - 1.)/den2 03017 - p3*numer_mp/den3 03018 - 0.5*p3*eta/den2 03019 - 0.5*p2*eta/den2 03020 - 2.*p2*p3*eta/den3; 03021 03022 case 3: 03023 return -0.25*numer_pm/den2 03024 + 0.5*p3*(1 - 2.*zeta)/den2 03025 - p3*numer_pm/den3 03026 + 0.5*p4*eta/den2 03027 + 0.5*p3*eta/den2 03028 + 2.*p3*p4*eta/den3; 03029 03030 case 4: 03031 return 0.; 03032 03033 case 5: 03034 return p4*eta/den2 03035 + 4.*p4*p1*eta/den3 03036 - p2*eta/den2 03037 - 4.*p1*p2*eta/den3; 03038 03039 case 6: 03040 return -p3*xi/den2 03041 - p1*xi/den2 03042 - 4.*p1*p3*xi/den3 03043 - 2.*p2*p3/den2 03044 - 2.*p1*p3/den2 03045 - 2.*p1*p2/den2 03046 - 8.*p1*p2*p3/den3; 03047 03048 case 7: 03049 return -p4*eta/den2 03050 - 4.*p3*p4*eta/den3 03051 + p2*eta/den2 03052 + 4.*p2*p3*eta/den3; 03053 03054 case 8: 03055 return -p3*xi/den2 03056 - p1*xi/den2 03057 - 4.*p1*p3*xi/den3 03058 + 2.*p4*p1/den2 03059 + 2.*p1*p3/den2 03060 + 2.*p3*p4/den2 03061 + 8.*p3*p4*p1/den3; 03062 03063 case 9: 03064 return -zeta/den 03065 + 2.*p1/den 03066 - 2.*p1*zeta/den2; 03067 03068 case 10: 03069 return zeta/den 03070 - 2.*p1/den 03071 + 2.*p1*zeta/den2; 03072 03073 case 11: 03074 return zeta/den 03075 - 2.*p3/den 03076 + 2.*p3*zeta/den2; 03077 03078 case 12: 03079 return -zeta/den 03080 + 2.*p3/den 03081 - 2.*p3*zeta/den2; 03082 03083 case 13: 03084 return -4.*p4*p1/den2 03085 - 4.*p3*p4/den2 03086 - 16.*p3*p4*p1/den3 03087 + 4.*p2*p3/den2 03088 + 4.*p1*p2/den2 03089 + 16.*p1*p2*p3/den3; 03090 03091 default: 03092 libmesh_error_msg("Invalid i = " << i); 03093 } 03094 } 03095 03096 case 4: // d^2()/detadzeta 03097 { 03098 switch(i) 03099 { 03100 case 0: 03101 return 0.25*numer_mp/den2 03102 - 0.5*p4*(2.*zeta - 1.)/den2 03103 + p4*numer_mp/den3 03104 - 0.5*p1*xi/den2 03105 - 0.5*p4*xi/den2 03106 - 2.*p4*p1*xi/den3; 03107 03108 case 1: 03109 return -0.25*numer_pm/den2 03110 + 0.5*p2*(1. - 2.*zeta)/den2 03111 - p2*numer_pm/den3 03112 + 0.5*p2*xi/den2 03113 + 0.5*p1*xi/den2 03114 + 2.*p1*p2*xi/den3; 03115 03116 case 2: 03117 return -0.25*numer_mp/den2 03118 + 0.5*p2*(2.*zeta - 1.)/den2 03119 - p2*numer_mp/den3 03120 - 0.5*p3*xi/den2 03121 - 0.5*p2*xi/den2 03122 - 2.*p2*p3*xi/den3; 03123 03124 case 3: 03125 return 0.25*numer_pm/den2 03126 - 0.5*p4*(1. - 2.*zeta)/den2 03127 + p4*numer_pm/den3 03128 + 0.5*p4*xi/den2 03129 + 0.5*p3*xi/den2 03130 + 2.*p3*p4*xi/den3; 03131 03132 case 4: 03133 return 0.; 03134 03135 case 5: 03136 return -p4*eta/den2 03137 - p2*eta/den2 03138 - 4.*p2*p4*eta/den3 03139 + 2.*p4*p1/den2 03140 + 2.*p2*p4/den2 03141 + 2.*p1*p2/den2 03142 + 8.*p2*p1*p4/den3; 03143 03144 case 6: 03145 return p3*xi/den2 03146 + 4.*p2*p3*xi/den3 03147 - p1*xi/den2 03148 - 4.*p1*p2*xi/den3; 03149 03150 case 7: 03151 return -p4*eta/den2 03152 - p2*eta/den2 03153 - 4.*p2*p4*eta/den3 03154 - 2.*p3*p4/den2 03155 - 2.*p2*p4/den2 03156 - 2.*p2*p3/den2 03157 - 8.*p2*p3*p4/den3; 03158 03159 case 8: 03160 return p1*xi/den2 03161 + 4.*p4*p1*xi/den3 03162 - p3*xi/den2 03163 - 4.*p3*p4*xi/den3; 03164 03165 case 9: 03166 return -zeta/den 03167 + 2.*p4/den 03168 - 2.*p4*zeta/den2; 03169 03170 case 10: 03171 return -zeta/den 03172 + 2.*p2/den 03173 - 2.*p2*zeta/den2; 03174 03175 case 11: 03176 return zeta/den 03177 - 2.*p2/den 03178 + 2.*p2*zeta/den2; 03179 03180 case 12: 03181 return zeta/den 03182 - 2.*p4/den 03183 + 2.*p4*zeta/den2; 03184 03185 case 13: 03186 return 4.*p3*p4/den2 03187 + 4.*p2*p3/den2 03188 + 16.*p2*p3*p4/den3 03189 - 4.*p4*p1/den2 03190 - 4.*p1*p2/den2 03191 - 16.*p2*p1*p4/den3; 03192 03193 default: 03194 libmesh_error_msg("Invalid i = " << i); 03195 } 03196 } 03197 03198 case 5: // d^2()/dzeta^2 03199 { 03200 switch(i) 03201 { 03202 case 0: 03203 return 0.5*numer_mp/den2 03204 - p1*(2.*zeta - 1.)/den2 03205 + 2.*p1*numer_mp/den3 03206 - p4*(2.*zeta - 1.)/den2 03207 + 2.*p4*numer_mp/den3 03208 + 2.*p4*p1/den2 03209 - 4.*p4*p1*(2.*zeta - 1.)/den3 03210 + 6.*p4*p1*numer_mp/den4; 03211 03212 case 1: 03213 return -0.5*numer_pm/den2 03214 + p2*(1 - 2.*zeta)/den2 03215 - 2.*p2*numer_pm/den3 03216 + p1*(1 - 2.*zeta)/den2 03217 - 2.*p1*numer_pm/den3 03218 + 2.*p1*p2/den2 03219 + 4.*p1*p2*(1 - 2.*zeta)/den3 03220 - 6.*p1*p2*numer_pm/den4; 03221 03222 case 2: 03223 return 0.5*numer_mp/den2 03224 - p3*(2.*zeta - 1.)/den2 03225 + 2.*p3*numer_mp/den3 03226 - p2*(2.*zeta - 1.)/den2 03227 + 2.*p2*numer_mp/den3 03228 + 2.*p2*p3/den2 03229 - 4.*p2*p3*(2.*zeta - 1.)/den3 03230 + 6.*p2*p3*numer_mp/den4; 03231 03232 case 3: 03233 return -0.5*numer_pm/den2 03234 + p4*(1 - 2.*zeta)/den2 03235 - 2.*p4*numer_pm/den3 03236 + p3*(1 - 2.*zeta)/den2 03237 - 2.*p3*numer_pm/den3 03238 + 2.*p3*p4/den2 03239 + 4.*p3*p4*(1 - 2.*zeta)/den3 03240 - 6.*p3*p4*numer_pm/den4; 03241 03242 case 4: 03243 return 4.; 03244 03245 case 5: 03246 return -2.*p1*eta/den2 03247 - 2.*p4*eta/den2 03248 - 8.*p4*p1*eta/den3 03249 - 2.*p2*eta/den2 03250 - 8.*p2*p4*eta/den3 03251 - 8.*p1*p2*eta/den3 03252 - 24.*p2*p1*p4*eta/den4; 03253 03254 case 6: 03255 return 2.*p3*xi/den2 03256 + 2.*p2*xi/den2 03257 + 8.*p2*p3*xi/den3 03258 + 2.*p1*xi/den2 03259 + 8.*p1*p3*xi/den3 03260 + 8.*p1*p2*xi/den3 03261 + 24.*p1*p2*p3*xi/den4; 03262 03263 case 7: 03264 return 2.*p4*eta/den2 03265 + 2.*p3*eta/den2 03266 + 8.*p3*p4*eta/den3 03267 + 2.*p2*eta/den2 03268 + 8.*p2*p4*eta/den3 03269 + 8.*p2*p3*eta/den3 03270 + 24.*p2*p3*p4*eta/den4; 03271 03272 case 8: 03273 return -2.*p1*xi/den2 03274 - 2.*p4*xi/den2 03275 - 8.*p4*p1*xi/den3 03276 - 2.*p3*xi/den2 03277 - 8.*p1*p3*xi/den3 03278 - 8.*p3*p4*xi/den3 03279 - 24.*p3*p4*p1*xi/den4; 03280 03281 case 9: 03282 return -2.*zeta/den 03283 + 4.*p4/den 03284 - 4.*p4*zeta/den2 03285 + 4.*p1/den 03286 - 4.*p1*zeta/den2 03287 + 8.*p4*p1/den2 03288 - 8.*p1*p4*zeta/den3; 03289 03290 case 10: 03291 return -2.*zeta/den 03292 + 4.*p1/den 03293 - 4.*p1*zeta/den2 03294 + 4.*p2/den 03295 - 4.*p2*zeta/den2 03296 + 8.*p1*p2/den2 03297 - 8.*p2*p1*zeta/den3; 03298 03299 case 11: 03300 return -2.*zeta/den 03301 + 4.*p2/den 03302 - 4.*p2*zeta/den2 03303 + 4.*p3/den 03304 - 4.*p3*zeta/den2 03305 + 8.*p2*p3/den2 03306 - 8.*p3*p2*zeta/den3; 03307 03308 case 12: 03309 return -2.*zeta/den 03310 + 4.*p3/den 03311 - 4.*p3*zeta/den2 03312 + 4.*p4/den 03313 - 4.*p4*zeta/den2 03314 + 8.*p3*p4/den2 03315 - 8.*p4*p3*zeta/den3; 03316 03317 case 13: 03318 return 8.*p3*p4/den2 03319 + 8.*p2*p4/den2 03320 + 8.*p2*p3/den2 03321 + 32.*p2*p3*p4/den3 03322 + 8.*p4*p1/den2 03323 + 8.*p1*p3/den2 03324 + 32.*p3*p4*p1/den3 03325 + 8.*p1*p2/den2 03326 + 32.*p2*p1*p4/den3 03327 + 32.*p1*p2*p3/den3 03328 + 96.*p1*p2*p3*p4/den4; 03329 03330 default: 03331 libmesh_error_msg("Invalid i = " << i); 03332 } 03333 } 03334 03335 default: 03336 libmesh_error_msg("Invalid j = " << j); 03337 } 03338 } 03339 03340 // G. Bedrosian, "Shape functions and integration formulas for 03341 // three-dimensional finite element analysis", Int. J. Numerical 03342 // Methods Engineering, vol 35, p. 95-108, 1992. 03343 case PYRAMID13: 03344 { 03345 libmesh_assert_less (i, 13); 03346 03347 const Real xi = p(0); 03348 const Real eta = p(1); 03349 const Real zeta = p(2); 03350 const Real eps = 1.e-35; 03351 03352 // Denominators are perturbed by epsilon to avoid 03353 // divide-by-zero issues. 03354 Real 03355 den = (-1. + zeta + eps), 03356 den2 = den*den, 03357 den3 = den2*den, 03358 xi2 = xi*xi, 03359 eta2 = eta*eta, 03360 zeta2 = zeta*zeta, 03361 zeta3 = zeta2*zeta; 03362 03363 switch (j) 03364 { 03365 case 0: // d^2()/dxi^2 03366 { 03367 switch(i) 03368 { 03369 case 0: 03370 case 1: 03371 return 0.5*(-1. + zeta + eta)/den; 03372 03373 case 2: 03374 case 3: 03375 return 0.5*(-1. + zeta - eta)/den; 03376 03377 case 4: 03378 case 6: 03379 case 8: 03380 case 9: 03381 case 10: 03382 case 11: 03383 case 12: 03384 return 0.; 03385 03386 case 5: 03387 return (1. - eta - zeta)/den; 03388 03389 case 7: 03390 return (1. + eta - zeta)/den; 03391 03392 default: 03393 libmesh_error_msg("Invalid i = " << i); 03394 } 03395 } 03396 03397 case 1: // d^2()/dxideta 03398 { 03399 switch(i) 03400 { 03401 case 0: 03402 return 0.25*(-1. + 2.*zeta + 2.*xi + 2.*eta)/den; 03403 03404 case 1: 03405 return -0.25*(-1. + 2.*zeta - 2.*xi + 2.*eta)/den; 03406 03407 case 2: 03408 return -0.25*(1. - 2.*zeta + 2.*xi + 2.*eta)/den; 03409 03410 case 3: 03411 return 0.25*(1. - 2.*zeta - 2.*xi + 2.*eta)/den; 03412 03413 case 4: 03414 return 0.; 03415 03416 case 5: 03417 return -xi/den; 03418 03419 case 6: 03420 return eta/den; 03421 03422 case 7: 03423 return xi/den; 03424 03425 case 8: 03426 return -eta/den; 03427 03428 case 9: 03429 return -zeta/den; 03430 03431 case 10: 03432 return zeta/den; 03433 03434 case 11: 03435 return -zeta/den; 03436 03437 case 12: 03438 return zeta/den; 03439 03440 default: 03441 libmesh_error_msg("Invalid i = " << i); 03442 } 03443 } 03444 03445 03446 case 2: // d^2()/deta^2 03447 { 03448 switch(i) 03449 { 03450 case 0: 03451 case 3: 03452 return 0.5*(-1. + zeta + xi)/den; 03453 03454 case 1: 03455 case 2: 03456 return 0.5*(-1. + zeta - xi)/den; 03457 03458 case 4: 03459 case 5: 03460 case 7: 03461 case 9: 03462 case 10: 03463 case 11: 03464 case 12: 03465 return 0.; 03466 03467 case 6: 03468 return (1. + xi - zeta)/den; 03469 03470 case 8: 03471 return (1. - xi - zeta)/den; 03472 03473 default: 03474 libmesh_error_msg("Invalid i = " << i); 03475 } 03476 } 03477 03478 03479 case 3: // d^2()/dxidzeta 03480 { 03481 switch(i) 03482 { 03483 case 0: 03484 return -0.25*(-1. + 2.*zeta - zeta2 + eta + 2.*eta*xi + eta2)/den2; 03485 03486 case 1: 03487 return 0.25*(-1. + 2.*zeta - zeta2 + eta - 2.*eta*xi + eta2)/den2; 03488 03489 case 2: 03490 return 0.25*(-1. + 2.*zeta - zeta2 - eta + 2.*eta*xi + eta2)/den2; 03491 03492 case 3: 03493 return -0.25*(-1. + 2.*zeta - zeta2 - eta - 2.*eta*xi + eta2)/den2; 03494 03495 case 4: 03496 return 0.; 03497 03498 case 5: 03499 return eta*xi/den2; 03500 03501 case 6: 03502 return -0.5*(1. + zeta2 + eta2 - 2.*zeta)/den2; 03503 03504 case 7: 03505 return -eta*xi/den2; 03506 03507 case 8: 03508 return 0.5*(1. + zeta2 + eta2 - 2.*zeta)/den2; 03509 03510 case 9: 03511 return (-1. - zeta2 + eta + 2.*zeta)/den2; 03512 03513 case 10: 03514 return -(-1. - zeta2 + eta + 2.*zeta)/den2; 03515 03516 case 11: 03517 return (1. + zeta2 + eta - 2.*zeta)/den2; 03518 03519 case 12: 03520 return -(1. + zeta2 + eta - 2.*zeta)/den2; 03521 03522 default: 03523 libmesh_error_msg("Invalid i = " << i); 03524 } 03525 } 03526 03527 case 4: // d^2()/detadzeta 03528 { 03529 switch(i) 03530 { 03531 case 0: 03532 return -0.25*(-1. + 2.*zeta - zeta2 + xi + 2.*eta*xi + xi2)/den2; 03533 03534 case 1: 03535 return 0.25*(1. - 2.*zeta + zeta2 + xi + 2.*eta*xi - xi2)/den2; 03536 03537 case 2: 03538 return 0.25*(-1. + 2.*zeta - zeta2 - xi + 2.*eta*xi + xi2)/den2; 03539 03540 case 3: 03541 return -0.25*(1. - 2.*zeta + zeta2 - xi + 2.*eta*xi - xi2)/den2; 03542 03543 case 4: 03544 return 0.; 03545 03546 case 5: 03547 return 0.5*(1. + xi2 + zeta2 - 2.*zeta)/den2; 03548 03549 case 6: 03550 return -eta*xi/den2; 03551 03552 case 7: 03553 return -0.5*(1. + xi2 + zeta2 - 2.*zeta)/den2; 03554 03555 case 8: 03556 return eta*xi/den2; 03557 03558 case 9: 03559 return (-1. - zeta2 + xi + 2.*zeta)/den2; 03560 03561 case 10: 03562 return -(1. + zeta2 + xi - 2.*zeta)/den2; 03563 03564 case 11: 03565 return (1. + zeta2 + xi - 2.*zeta)/den2; 03566 03567 case 12: 03568 return -(-1. - zeta2 + xi + 2.*zeta)/den2; 03569 03570 default: 03571 libmesh_error_msg("Invalid i = " << i); 03572 } 03573 } 03574 03575 case 5: // d^2()/dzeta^2 03576 { 03577 switch(i) 03578 { 03579 case 0: 03580 return 0.5*(xi + eta + 1.)*eta*xi/den3; 03581 03582 case 1: 03583 return -0.5*(eta - xi + 1.)*eta*xi/den3; 03584 03585 case 2: 03586 return -0.5*(xi + eta - 1.)*eta*xi/den3; 03587 03588 case 3: 03589 return 0.5*(eta - xi - 1.)*eta*xi/den3; 03590 03591 case 4: 03592 return 4.; 03593 03594 case 5: 03595 return -(1. - 3.*zeta + 3.*zeta2 - zeta3 + eta*xi2)/den3; 03596 03597 case 6: 03598 return (-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta2*xi)/den3; 03599 03600 case 7: 03601 return (-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta*xi2)/den3; 03602 03603 case 8: 03604 return -(1. - 3.*zeta + 3.*zeta2 - zeta3 + eta2*xi)/den3; 03605 03606 case 9: 03607 return -2.*(-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta*xi)/den3; 03608 03609 case 10: 03610 return 2.*(1. - 3.*zeta + 3.*zeta2 - zeta3 + eta*xi)/den3; 03611 03612 case 11: 03613 return -2.*(-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta*xi)/den3; 03614 03615 case 12: 03616 return 2.*(1. - 3.*zeta + 3.*zeta2 - zeta3 + eta*xi)/den3; 03617 03618 default: 03619 libmesh_error_msg("Invalid i = " << i); 03620 } 03621 } 03622 03623 default: 03624 libmesh_error_msg("Invalid j = " << j); 03625 } 03626 } 03627 03628 default: 03629 libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type); 03630 } 03631 } 03632 03633 03634 // unsupported order 03635 default: 03636 libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order); 03637 } 03638 03639 #endif 03640 03641 libmesh_error_msg("We'll never get here!"); 03642 return 0.; 03643 } 03644 03645 03646 03647 template <> 03648 Real FE<3,LAGRANGE>::shape_second_deriv(const Elem* elem, 03649 const Order order, 03650 const unsigned int i, 03651 const unsigned int j, 03652 const Point& p) 03653 { 03654 libmesh_assert(elem); 03655 03656 // call the orientation-independent shape function derivatives 03657 return FE<3,LAGRANGE>::shape_second_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p); 03658 } 03659 03660 } // namespace libMesh
generated by
