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