fe_bernstein_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++ includes 00020 00021 00022 // Local includes 00023 #include "libmesh/libmesh_config.h" 00024 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES 00025 00026 #include "libmesh/fe.h" 00027 #include "libmesh/elem.h" 00028 00029 00030 namespace libMesh 00031 { 00032 00033 00034 00035 template <> 00036 Real FE<3,BERNSTEIN>::shape(const ElemType, 00037 const Order, 00038 const unsigned int, 00039 const Point&) 00040 { 00041 libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge and face orientation is needed."); 00042 return 0.; 00043 } 00044 00045 00046 00047 template <> 00048 Real FE<3,BERNSTEIN>::shape(const Elem* elem, 00049 const Order order, 00050 const unsigned int i, 00051 const Point& p) 00052 { 00053 00054 #if LIBMESH_DIM == 3 00055 00056 libmesh_assert(elem); 00057 const ElemType type = elem->type(); 00058 00059 const Order totalorder = static_cast<Order>(order + elem->p_level()); 00060 00061 switch (totalorder) 00062 { 00063 // 1st order Bernstein. 00064 case FIRST: 00065 { 00066 switch (type) 00067 { 00068 00069 // Bernstein shape functions on the tetrahedron. 00070 case TET4: 00071 case TET10: 00072 { 00073 libmesh_assert_less (i, 4); 00074 00075 // Area coordinates 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: return zeta0; 00084 case 1: return zeta1; 00085 case 2: return zeta2; 00086 case 3: return zeta3; 00087 00088 default: 00089 libmesh_error_msg("Invalid shape function index i = " << i); 00090 } 00091 } 00092 00093 // Bernstein shape functions on the hexahedral. 00094 case HEX20: 00095 case HEX27: 00096 { 00097 libmesh_assert_less (i, 8); 00098 00099 // Compute hex shape functions as a tensor-product 00100 const Real xi = p(0); 00101 const Real eta = p(1); 00102 const Real zeta = p(2); 00103 00104 // The only way to make any sense of this 00105 // is to look at the mgflo/mg2/mgf documentation 00106 // and make the cut-out cube! 00107 // 0 1 2 3 4 5 6 7 00108 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0}; 00109 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1}; 00110 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1}; 00111 00112 return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)* 00113 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)* 00114 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta)); 00115 } 00116 00117 00118 default: 00119 libmesh_error_msg("Invalid element type = " << type); 00120 } 00121 } 00122 00123 00124 00125 00126 case SECOND: 00127 { 00128 switch (type) 00129 { 00130 00131 // Bernstein shape functions on the tetrahedron. 00132 case TET10: 00133 { 00134 libmesh_assert_less (i, 10); 00135 00136 // Area coordinates 00137 const Real zeta1 = p(0); 00138 const Real zeta2 = p(1); 00139 const Real zeta3 = p(2); 00140 const Real zeta0 = 1. - zeta1 - zeta2 - zeta3; 00141 00142 switch(i) 00143 { 00144 case 0: return zeta0*zeta0; 00145 case 1: return zeta1*zeta1; 00146 case 2: return zeta2*zeta2; 00147 case 3: return zeta3*zeta3; 00148 case 4: return 2.*zeta0*zeta1; 00149 case 5: return 2.*zeta1*zeta2; 00150 case 6: return 2.*zeta0*zeta2; 00151 case 7: return 2.*zeta3*zeta0; 00152 case 8: return 2.*zeta1*zeta3; 00153 case 9: return 2.*zeta2*zeta3; 00154 00155 default: 00156 libmesh_error_msg("Invalid shape function index i = " << i); 00157 } 00158 } 00159 00160 // Bernstein shape functions on the 20-noded hexahedral. 00161 case HEX20: 00162 { 00163 libmesh_assert_less (i, 20); 00164 00165 // Compute hex shape functions as a tensor-product 00166 const Real xi = p(0); 00167 const Real eta = p(1); 00168 const Real zeta = p(2); 00169 00170 // The only way to make any sense of this 00171 // is to look at the mgflo/mg2/mgf documentation 00172 // and make the cut-out cube! 00173 // 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 00174 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}; 00175 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}; 00176 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}; 00177 //To compute the hex20 shape functions the original shape functions for hex27 are used. 00178 //scalx[i] tells how often the original x-th shape function has to be added to the original i-th shape function 00179 //to compute the new i-th shape function for hex20 00180 //example: B_0^HEX20 = B_0^HEX27 - 0.25*B_20^HEX27 - 0.25*B_21^HEX27 + 0*B_22^HEX27 + 0*B_23^HEX27 - 0.25*B_24^HEX27 + 0*B_25^HEX27 - 0.25*B_26^HEX27 00181 // B_0^HEX20 = B_0^HEX27 + scal20[0]*B_20^HEX27 + scal21[0]*B_21^HEX27 + ... 00182 static const Real scal20[] = {-0.25, -0.25, -0.25, -0.25, 0, 0, 0, 0, 0.5, 0.5, 0.5, 0.5, 0, 0, 0, 0, 0, 0, 0, 0}; 00183 static const Real scal21[] = {-0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0, 0, 0}; 00184 static const Real scal22[] = {0, -0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0, 0}; 00185 static const Real scal23[] = {0, 0, -0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0}; 00186 static const Real scal24[] = {-0.25, 0, 0, -0.25, -0.25, 0, 0, -0.25, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0, 0, 0, 0.5}; 00187 static const Real scal25[] = {0, 0, 0, 0, -0.25, -0.25, -0.25, -0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0.5, 0.5}; 00188 static const Real scal26[] = {-0.25, -0.25, -0.25, -0.25, -0.25, -0.25, -0.25, -0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}; 00189 00190 return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)* 00191 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)* 00192 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta) 00193 +scal20[i]* 00194 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[20], xi)* 00195 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[20], eta)* 00196 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[20], zeta) 00197 +scal21[i]* 00198 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[21], xi)* 00199 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[21], eta)* 00200 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[21], zeta) 00201 +scal22[i]* 00202 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[22], xi)* 00203 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[22], eta)* 00204 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[22], zeta) 00205 +scal23[i]* 00206 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[23], xi)* 00207 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[23], eta)* 00208 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[23], zeta) 00209 +scal24[i]* 00210 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[24], xi)* 00211 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[24], eta)* 00212 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[24], zeta) 00213 +scal25[i]* 00214 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[25], xi)* 00215 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[25], eta)* 00216 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[25], zeta) 00217 +scal26[i]* 00218 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[26], xi)* 00219 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[26], eta)* 00220 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[26], zeta)); 00221 } 00222 00223 // Bernstein shape functions on the hexahedral. 00224 case HEX27: 00225 { 00226 libmesh_assert_less (i, 27); 00227 00228 // Compute hex shape functions as a tensor-product 00229 const Real xi = p(0); 00230 const Real eta = p(1); 00231 const Real zeta = p(2); 00232 00233 // The only way to make any sense of this 00234 // is to look at the mgflo/mg2/mgf documentation 00235 // and make the cut-out cube! 00236 // 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 00237 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}; 00238 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}; 00239 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}; 00240 00241 return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)* 00242 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)* 00243 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta)); 00244 } 00245 00246 00247 default: 00248 libmesh_error_msg("Invalid element type = " << type); 00249 } 00250 00251 } 00252 00253 00254 00255 // 3rd-order Bernstein. 00256 case THIRD: 00257 { 00258 switch (type) 00259 { 00260 00261 // // Bernstein shape functions on the tetrahedron. 00262 // case TET10: 00263 // { 00264 // libmesh_assert_less (i, 20); 00265 00266 // // Area coordinates 00267 // const Real zeta1 = p(0); 00268 // const Real zeta2 = p(1); 00269 // const Real zeta3 = p(2); 00270 // const Real zeta0 = 1. - zeta1 - zeta2 - zeta3; 00271 00272 00273 // unsigned int shape=i; 00274 00275 // // handle the edge orientation 00276 00277 // if ((i== 4||i== 5) && elem->node(0) > elem->node(1))shape= 9-i; //Edge 0 00278 // if ((i== 6||i== 7) && elem->node(1) > elem->node(2))shape=13-i; //Edge 1 00279 // if ((i== 8||i== 9) && elem->node(0) > elem->node(2))shape=17-i; //Edge 2 00280 // if ((i==10||i==11) && elem->node(0) > elem->node(3))shape=21-i; //Edge 3 00281 // if ((i==12||i==13) && elem->node(1) > elem->node(3))shape=25-i; //Edge 4 00282 // if ((i==14||i==15) && elem->node(2) > elem->node(3))shape=29-i; //Edge 5 00283 00284 // // No need to handle face orientation in 3rd order. 00285 00286 00287 // switch(shape) 00288 // { 00289 // //point function 00290 // case 0: return zeta0*zeta0*zeta0; 00291 // case 1: return zeta1*zeta1*zeta1; 00292 // case 2: return zeta2*zeta2*zeta2; 00293 // case 3: return zeta3*zeta3*zeta3; 00294 00295 // //edge functions 00296 // case 4: return 3.*zeta0*zeta0*zeta1; 00297 // case 5: return 3.*zeta1*zeta1*zeta0; 00298 00299 // case 6: return 3.*zeta1*zeta1*zeta2; 00300 // case 7: return 3.*zeta2*zeta2*zeta1; 00301 00302 // case 8: return 3.*zeta0*zeta0*zeta2; 00303 // case 9: return 3.*zeta2*zeta2*zeta0; 00304 00305 // case 10: return 3.*zeta0*zeta0*zeta3; 00306 // case 11: return 3.*zeta3*zeta3*zeta0; 00307 00308 // case 12: return 3.*zeta1*zeta1*zeta3; 00309 // case 13: return 3.*zeta3*zeta3*zeta1; 00310 00311 // case 14: return 3.*zeta2*zeta2*zeta3; 00312 // case 15: return 3.*zeta3*zeta3*zeta2; 00313 00314 // //face functions 00315 // case 16: return 6.*zeta0*zeta1*zeta2; 00316 // case 17: return 6.*zeta0*zeta1*zeta3; 00317 // case 18: return 6.*zeta1*zeta2*zeta3; 00318 // case 19: return 6.*zeta2*zeta0*zeta3; 00319 00320 // default: 00321 // libmesh_error_msg("Invalid shape function index i = " << i); 00322 // } 00323 // } 00324 00325 00326 // Bernstein shape functions on the hexahedral. 00327 case HEX27: 00328 { 00329 libmesh_assert_less (i, 64); 00330 00331 // Compute hex shape functions as a tensor-product 00332 const Real xi = p(0); 00333 const Real eta = p(1); 00334 const Real zeta = p(2); 00335 Real xi_mapped = p(0); 00336 Real eta_mapped = p(1); 00337 Real zeta_mapped = p(2); 00338 00339 // The only way to make any sense of this 00340 // is to look at the mgflo/mg2/mgf documentation 00341 // and make the cut-out cube! 00342 // Nodes 0 1 2 3 4 5 6 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 20 20 21 21 21 21 22 22 22 22 23 23 23 23 24 24 24 24 25 25 25 25 26 26 26 26 26 26 26 26 00343 // DOFS 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 27 18 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 60 62 63 00344 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3}; 00345 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3, 0, 0, 0, 0, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3}; 00346 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3}; 00347 00348 00349 00350 // handle the edge orientation 00351 { 00352 // Edge 0 00353 if ( (i0[i] >= 2) && (i1[i] == 0) && (i2[i] == 0)) 00354 { 00355 if (elem->point(0) != std::min(elem->point(0), elem->point(1))) 00356 xi_mapped = -xi; 00357 } 00358 // Edge 1 00359 else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] == 0)) 00360 { 00361 if (elem->point(1) != std::min(elem->point(1), elem->point(2))) 00362 eta_mapped = -eta; 00363 } 00364 // Edge 2 00365 else if ((i0[i] >= 2) && (i1[i] == 1) && (i2[i] == 0)) 00366 { 00367 if (elem->point(3) != std::min(elem->point(3), elem->point(2))) 00368 xi_mapped = -xi; 00369 } 00370 // Edge 3 00371 else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] == 0)) 00372 { 00373 if (elem->point(0) != std::min(elem->point(0), elem->point(3))) 00374 eta_mapped = -eta; 00375 } 00376 // Edge 4 00377 else if ((i0[i] == 0) && (i1[i] == 0) && (i2[i] >=2 )) 00378 { 00379 if (elem->point(0) != std::min(elem->point(0), elem->point(4))) 00380 zeta_mapped = -zeta; 00381 } 00382 // Edge 5 00383 else if ((i0[i] == 1) && (i1[i] == 0) && (i2[i] >=2 )) 00384 { 00385 if (elem->point(1) != std::min(elem->point(1), elem->point(5))) 00386 zeta_mapped = -zeta; 00387 } 00388 // Edge 6 00389 else if ((i0[i] == 1) && (i1[i] == 1) && (i2[i] >=2 )) 00390 { 00391 if (elem->point(2) != std::min(elem->point(2), elem->point(6))) 00392 zeta_mapped = -zeta; 00393 } 00394 // Edge 7 00395 else if ((i0[i] == 0) && (i1[i] == 1) && (i2[i] >=2 )) 00396 { 00397 if (elem->point(3) != std::min(elem->point(3), elem->point(7))) 00398 zeta_mapped = -zeta; 00399 } 00400 // Edge 8 00401 else if ((i0[i] >=2 ) && (i1[i] == 0) && (i2[i] == 1)) 00402 { 00403 if (elem->point(4) != std::min(elem->point(4), elem->point(5))) 00404 xi_mapped = -xi; 00405 } 00406 // Edge 9 00407 else if ((i0[i] == 1) && (i1[i] >=2 ) && (i2[i] == 1)) 00408 { 00409 if (elem->point(5) != std::min(elem->point(5), elem->point(6))) 00410 eta_mapped = -eta; 00411 } 00412 // Edge 10 00413 else if ((i0[i] >=2 ) && (i1[i] == 1) && (i2[i] == 1)) 00414 { 00415 if (elem->point(7) != std::min(elem->point(7), elem->point(6))) 00416 xi_mapped = -xi; 00417 } 00418 // Edge 11 00419 else if ((i0[i] == 0) && (i1[i] >=2 ) && (i2[i] == 1)) 00420 { 00421 if (elem->point(4) != std::min(elem->point(4), elem->point(7))) 00422 eta_mapped = -eta; 00423 } 00424 } 00425 00426 00427 // handle the face orientation 00428 { 00429 // Face 0 00430 if ( (i2[i] == 0) && (i0[i] >= 2) && (i1[i] >= 2)) 00431 { 00432 const Point min_point = std::min(elem->point(1), 00433 std::min(elem->point(2), 00434 std::min(elem->point(0), 00435 elem->point(3)))); 00436 if (elem->point(0) == min_point) 00437 if (elem->point(1) == std::min(elem->point(1), elem->point(3))) 00438 { 00439 // Case 1 00440 xi_mapped = xi; 00441 eta_mapped = eta; 00442 } 00443 else 00444 { 00445 // Case 2 00446 xi_mapped = eta; 00447 eta_mapped = xi; 00448 } 00449 00450 else if (elem->point(3) == min_point) 00451 if (elem->point(0) == std::min(elem->point(0), elem->point(2))) 00452 { 00453 // Case 3 00454 xi_mapped = -eta; 00455 eta_mapped = xi; 00456 } 00457 else 00458 { 00459 // Case 4 00460 xi_mapped = xi; 00461 eta_mapped = -eta; 00462 } 00463 00464 else if (elem->point(2) == min_point) 00465 if (elem->point(3) == std::min(elem->point(3), elem->point(1))) 00466 { 00467 // Case 5 00468 xi_mapped = -xi; 00469 eta_mapped = -eta; 00470 } 00471 else 00472 { 00473 // Case 6 00474 xi_mapped = -eta; 00475 eta_mapped = -xi; 00476 } 00477 00478 else if (elem->point(1) == min_point) 00479 { 00480 if (elem->point(2) == std::min(elem->point(2), elem->point(0))) 00481 { 00482 // Case 7 00483 xi_mapped = eta; 00484 eta_mapped = -xi; 00485 } 00486 else 00487 { 00488 // Case 8 00489 xi_mapped = -xi; 00490 eta_mapped = eta; 00491 } 00492 } 00493 } 00494 00495 00496 // Face 1 00497 else if ((i1[i] == 0) && (i0[i] >= 2) && (i2[i] >= 2)) 00498 { 00499 const Point min_point = std::min(elem->point(0), 00500 std::min(elem->point(1), 00501 std::min(elem->point(5), 00502 elem->point(4)))); 00503 if (elem->point(0) == min_point) 00504 if (elem->point(1) == std::min(elem->point(1), elem->point(4))) 00505 { 00506 // Case 1 00507 xi_mapped = xi; 00508 zeta_mapped = zeta; 00509 } 00510 else 00511 { 00512 // Case 2 00513 xi_mapped = zeta; 00514 zeta_mapped = xi; 00515 } 00516 00517 else if (elem->point(1) == min_point) 00518 if (elem->point(5) == std::min(elem->point(5), elem->point(0))) 00519 { 00520 // Case 3 00521 xi_mapped = zeta; 00522 zeta_mapped = -xi; 00523 } 00524 else 00525 { 00526 // Case 4 00527 xi_mapped = -xi; 00528 zeta_mapped = zeta; 00529 } 00530 00531 else if (elem->point(5) == min_point) 00532 if (elem->point(4) == std::min(elem->point(4), elem->point(1))) 00533 { 00534 // Case 5 00535 xi_mapped = -xi; 00536 zeta_mapped = -zeta; 00537 } 00538 else 00539 { 00540 // Case 6 00541 xi_mapped = -zeta; 00542 zeta_mapped = -xi; 00543 } 00544 00545 else if (elem->point(4) == min_point) 00546 { 00547 if (elem->point(0) == std::min(elem->point(0), elem->point(5))) 00548 { 00549 // Case 7 00550 xi_mapped = -xi; 00551 zeta_mapped = zeta; 00552 } 00553 else 00554 { 00555 // Case 8 00556 xi_mapped = xi; 00557 zeta_mapped = -zeta; 00558 } 00559 } 00560 } 00561 00562 00563 // Face 2 00564 else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] >= 2)) 00565 { 00566 const Point min_point = std::min(elem->point(1), 00567 std::min(elem->point(2), 00568 std::min(elem->point(6), 00569 elem->point(5)))); 00570 if (elem->point(1) == min_point) 00571 if (elem->point(2) == std::min(elem->point(2), elem->point(5))) 00572 { 00573 // Case 1 00574 eta_mapped = eta; 00575 zeta_mapped = zeta; 00576 } 00577 else 00578 { 00579 // Case 2 00580 eta_mapped = zeta; 00581 zeta_mapped = eta; 00582 } 00583 00584 else if (elem->point(2) == min_point) 00585 if (elem->point(6) == std::min(elem->point(6), elem->point(1))) 00586 { 00587 // Case 3 00588 eta_mapped = zeta; 00589 zeta_mapped = -eta; 00590 } 00591 else 00592 { 00593 // Case 4 00594 eta_mapped = -eta; 00595 zeta_mapped = zeta; 00596 } 00597 00598 else if (elem->point(6) == min_point) 00599 if (elem->point(5) == std::min(elem->point(5), elem->point(2))) 00600 { 00601 // Case 5 00602 eta_mapped = -eta; 00603 zeta_mapped = -zeta; 00604 } 00605 else 00606 { 00607 // Case 6 00608 eta_mapped = -zeta; 00609 zeta_mapped = -eta; 00610 } 00611 00612 else if (elem->point(5) == min_point) 00613 { 00614 if (elem->point(1) == std::min(elem->point(1), elem->point(6))) 00615 { 00616 // Case 7 00617 eta_mapped = -zeta; 00618 zeta_mapped = eta; 00619 } 00620 else 00621 { 00622 // Case 8 00623 eta_mapped = eta; 00624 zeta_mapped = -zeta; 00625 } 00626 } 00627 } 00628 00629 00630 // Face 3 00631 else if ((i1[i] == 1) && (i0[i] >= 2) && (i2[i] >= 2)) 00632 { 00633 const Point min_point = std::min(elem->point(2), 00634 std::min(elem->point(3), 00635 std::min(elem->point(7), 00636 elem->point(6)))); 00637 if (elem->point(3) == min_point) 00638 if (elem->point(2) == std::min(elem->point(2), elem->point(7))) 00639 { 00640 // Case 1 00641 xi_mapped = xi; 00642 zeta_mapped = zeta; 00643 } 00644 else 00645 { 00646 // Case 2 00647 xi_mapped = zeta; 00648 zeta_mapped = xi; 00649 } 00650 00651 else if (elem->point(7) == min_point) 00652 if (elem->point(3) == std::min(elem->point(3), elem->point(6))) 00653 { 00654 // Case 3 00655 xi_mapped = -zeta; 00656 zeta_mapped = xi; 00657 } 00658 else 00659 { 00660 // Case 4 00661 xi_mapped = xi; 00662 zeta_mapped = -zeta; 00663 } 00664 00665 else if (elem->point(6) == min_point) 00666 if (elem->point(7) == std::min(elem->point(7), elem->point(2))) 00667 { 00668 // Case 5 00669 xi_mapped = -xi; 00670 zeta_mapped = -zeta; 00671 } 00672 else 00673 { 00674 // Case 6 00675 xi_mapped = -zeta; 00676 zeta_mapped = -xi; 00677 } 00678 00679 else if (elem->point(2) == min_point) 00680 { 00681 if (elem->point(6) == std::min(elem->point(3), elem->point(6))) 00682 { 00683 // Case 7 00684 xi_mapped = zeta; 00685 zeta_mapped = -xi; 00686 } 00687 else 00688 { 00689 // Case 8 00690 xi_mapped = -xi; 00691 zeta_mapped = zeta; 00692 } 00693 } 00694 } 00695 00696 00697 // Face 4 00698 else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] >= 2)) 00699 { 00700 const Point min_point = std::min(elem->point(3), 00701 std::min(elem->point(0), 00702 std::min(elem->point(4), 00703 elem->point(7)))); 00704 if (elem->point(0) == min_point) 00705 if (elem->point(3) == std::min(elem->point(3), elem->point(4))) 00706 { 00707 // Case 1 00708 eta_mapped = eta; 00709 zeta_mapped = zeta; 00710 } 00711 else 00712 { 00713 // Case 2 00714 eta_mapped = zeta; 00715 zeta_mapped = eta; 00716 } 00717 00718 else if (elem->point(4) == min_point) 00719 if (elem->point(0) == std::min(elem->point(0), elem->point(7))) 00720 { 00721 // Case 3 00722 eta_mapped = -zeta; 00723 zeta_mapped = eta; 00724 } 00725 else 00726 { 00727 // Case 4 00728 eta_mapped = eta; 00729 zeta_mapped = -zeta; 00730 } 00731 00732 else if (elem->point(7) == min_point) 00733 if (elem->point(4) == std::min(elem->point(4), elem->point(3))) 00734 { 00735 // Case 5 00736 eta_mapped = -eta; 00737 zeta_mapped = -zeta; 00738 } 00739 else 00740 { 00741 // Case 6 00742 eta_mapped = -zeta; 00743 zeta_mapped = -eta; 00744 } 00745 00746 else if (elem->point(3) == min_point) 00747 { 00748 if (elem->point(7) == std::min(elem->point(7), elem->point(0))) 00749 { 00750 // Case 7 00751 eta_mapped = zeta; 00752 zeta_mapped = -eta; 00753 } 00754 else 00755 { 00756 // Case 8 00757 eta_mapped = -eta; 00758 zeta_mapped = zeta; 00759 } 00760 } 00761 } 00762 00763 00764 // Face 5 00765 else if ((i2[i] == 1) && (i0[i] >= 2) && (i1[i] >= 2)) 00766 { 00767 const Point min_point = std::min(elem->point(4), 00768 std::min(elem->point(5), 00769 std::min(elem->point(6), 00770 elem->point(7)))); 00771 if (elem->point(4) == min_point) 00772 if (elem->point(5) == std::min(elem->point(5), elem->point(7))) 00773 { 00774 // Case 1 00775 xi_mapped = xi; 00776 eta_mapped = eta; 00777 } 00778 else 00779 { 00780 // Case 2 00781 xi_mapped = eta; 00782 eta_mapped = xi; 00783 } 00784 00785 else if (elem->point(5) == min_point) 00786 if (elem->point(6) == std::min(elem->point(6), elem->point(4))) 00787 { 00788 // Case 3 00789 xi_mapped = eta; 00790 eta_mapped = -xi; 00791 } 00792 else 00793 { 00794 // Case 4 00795 xi_mapped = -xi; 00796 eta_mapped = eta; 00797 } 00798 00799 else if (elem->point(6) == min_point) 00800 if (elem->point(7) == std::min(elem->point(7), elem->point(5))) 00801 { 00802 // Case 5 00803 xi_mapped = -xi; 00804 eta_mapped = -eta; 00805 } 00806 else 00807 { 00808 // Case 6 00809 xi_mapped = -eta; 00810 eta_mapped = -xi; 00811 } 00812 00813 else if (elem->point(7) == min_point) 00814 { 00815 if (elem->point(4) == std::min(elem->point(4), elem->point(6))) 00816 { 00817 // Case 7 00818 xi_mapped = -eta; 00819 eta_mapped = xi; 00820 } 00821 else 00822 { 00823 // Case 8 00824 xi_mapped = xi; 00825 eta_mapped = eta; 00826 } 00827 } 00828 } 00829 } 00830 00831 return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi_mapped)* 00832 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta_mapped)* 00833 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta_mapped)); 00834 } 00835 00836 00837 default: 00838 libmesh_error_msg("Invalid element type = " << type); 00839 } //case HEX27 00840 00841 }//case THIRD 00842 00843 00844 // 4th-order Bernstein. 00845 case FOURTH: 00846 { 00847 switch (type) 00848 { 00849 00850 // Bernstein shape functions on the hexahedral. 00851 case HEX27: 00852 { 00853 libmesh_assert_less (i, 125); 00854 00855 // Compute hex shape functions as a tensor-product 00856 const Real xi = p(0); 00857 const Real eta = p(1); 00858 const Real zeta = p(2); 00859 Real xi_mapped = p(0); 00860 Real eta_mapped = p(1); 00861 Real zeta_mapped = p(2); 00862 00863 // The only way to make any sense of this 00864 // is to look at the mgflo/mg2/mgf documentation 00865 // and make the cut-out cube! 00866 // Nodes 0 1 2 3 4 5 6 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 20 20 21 21 21 21 22 22 22 22 23 23 23 23 24 24 24 24 25 25 25 25 26 26 26 26 26 26 26 26 00867 // DOFS 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 27 18 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 60 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 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 00868 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 2, 3, 4, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4}; 00869 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4}; 00870 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4}; 00871 00872 00873 00874 // handle the edge orientation 00875 { 00876 // Edge 0 00877 if ( (i0[i] >= 2) && (i1[i] == 0) && (i2[i] == 0)) 00878 { 00879 if (elem->point(0) != std::min(elem->point(0), elem->point(1))) 00880 xi_mapped = -xi; 00881 } 00882 // Edge 1 00883 else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] == 0)) 00884 { 00885 if (elem->point(1) != std::min(elem->point(1), elem->point(2))) 00886 eta_mapped = -eta; 00887 } 00888 // Edge 2 00889 else if ((i0[i] >= 2) && (i1[i] == 1) && (i2[i] == 0)) 00890 { 00891 if (elem->point(3) != std::min(elem->point(3), elem->point(2))) 00892 xi_mapped = -xi; 00893 } 00894 // Edge 3 00895 else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] == 0)) 00896 { 00897 if (elem->point(0) != std::min(elem->point(0), elem->point(3))) 00898 eta_mapped = -eta; 00899 } 00900 // Edge 4 00901 else if ((i0[i] == 0) && (i1[i] == 0) && (i2[i] >=2 )) 00902 { 00903 if (elem->point(0) != std::min(elem->point(0), elem->point(4))) 00904 zeta_mapped = -zeta; 00905 } 00906 // Edge 5 00907 else if ((i0[i] == 1) && (i1[i] == 0) && (i2[i] >=2 )) 00908 { 00909 if (elem->point(1) != std::min(elem->point(1), elem->point(5))) 00910 zeta_mapped = -zeta; 00911 } 00912 // Edge 6 00913 else if ((i0[i] == 1) && (i1[i] == 1) && (i2[i] >=2 )) 00914 { 00915 if (elem->point(2) != std::min(elem->point(2), elem->point(6))) 00916 zeta_mapped = -zeta; 00917 } 00918 // Edge 7 00919 else if ((i0[i] == 0) && (i1[i] == 1) && (i2[i] >=2 )) 00920 { 00921 if (elem->point(3) != std::min(elem->point(3), elem->point(7))) 00922 zeta_mapped = -zeta; 00923 } 00924 // Edge 8 00925 else if ((i0[i] >=2 ) && (i1[i] == 0) && (i2[i] == 1)) 00926 { 00927 if (elem->point(4) != std::min(elem->point(4), elem->point(5))) 00928 xi_mapped = -xi; 00929 } 00930 // Edge 9 00931 else if ((i0[i] == 1) && (i1[i] >=2 ) && (i2[i] == 1)) 00932 { 00933 if (elem->point(5) != std::min(elem->point(5), elem->point(6))) 00934 eta_mapped = -eta; 00935 } 00936 // Edge 10 00937 else if ((i0[i] >=2 ) && (i1[i] == 1) && (i2[i] == 1)) 00938 { 00939 if (elem->point(7) != std::min(elem->point(7), elem->point(6))) 00940 xi_mapped = -xi; 00941 } 00942 // Edge 11 00943 else if ((i0[i] == 0) && (i1[i] >=2 ) && (i2[i] == 1)) 00944 { 00945 if (elem->point(4) != std::min(elem->point(4), elem->point(7))) 00946 eta_mapped = -eta; 00947 } 00948 } 00949 00950 00951 // handle the face orientation 00952 { 00953 // Face 0 00954 if ( (i2[i] == 0) && (i0[i] >= 2) && (i1[i] >= 2)) 00955 { 00956 const Point min_point = std::min(elem->point(1), 00957 std::min(elem->point(2), 00958 std::min(elem->point(0), 00959 elem->point(3)))); 00960 if (elem->point(0) == min_point) 00961 if (elem->point(1) == std::min(elem->point(1), elem->point(3))) 00962 { 00963 // Case 1 00964 xi_mapped = xi; 00965 eta_mapped = eta; 00966 } 00967 else 00968 { 00969 // Case 2 00970 xi_mapped = eta; 00971 eta_mapped = xi; 00972 } 00973 00974 else if (elem->point(3) == min_point) 00975 if (elem->point(0) == std::min(elem->point(0), elem->point(2))) 00976 { 00977 // Case 3 00978 xi_mapped = -eta; 00979 eta_mapped = xi; 00980 } 00981 else 00982 { 00983 // Case 4 00984 xi_mapped = xi; 00985 eta_mapped = -eta; 00986 } 00987 00988 else if (elem->point(2) == min_point) 00989 if (elem->point(3) == std::min(elem->point(3), elem->point(1))) 00990 { 00991 // Case 5 00992 xi_mapped = -xi; 00993 eta_mapped = -eta; 00994 } 00995 else 00996 { 00997 // Case 6 00998 xi_mapped = -eta; 00999 eta_mapped = -xi; 01000 } 01001 01002 else if (elem->point(1) == min_point) 01003 { 01004 if (elem->point(2) == std::min(elem->point(2), elem->point(0))) 01005 { 01006 // Case 7 01007 xi_mapped = eta; 01008 eta_mapped = -xi; 01009 } 01010 else 01011 { 01012 // Case 8 01013 xi_mapped = -xi; 01014 eta_mapped = eta; 01015 } 01016 } 01017 } 01018 01019 01020 // Face 1 01021 else if ((i1[i] == 0) && (i0[i] >= 2) && (i2[i] >= 2)) 01022 { 01023 const Point min_point = std::min(elem->point(0), 01024 std::min(elem->point(1), 01025 std::min(elem->point(5), 01026 elem->point(4)))); 01027 if (elem->point(0) == min_point) 01028 if (elem->point(1) == std::min(elem->point(1), elem->point(4))) 01029 { 01030 // Case 1 01031 xi_mapped = xi; 01032 zeta_mapped = zeta; 01033 } 01034 else 01035 { 01036 // Case 2 01037 xi_mapped = zeta; 01038 zeta_mapped = xi; 01039 } 01040 01041 else if (elem->point(1) == min_point) 01042 if (elem->point(5) == std::min(elem->point(5), elem->point(0))) 01043 { 01044 // Case 3 01045 xi_mapped = zeta; 01046 zeta_mapped = -xi; 01047 } 01048 else 01049 { 01050 // Case 4 01051 xi_mapped = -xi; 01052 zeta_mapped = zeta; 01053 } 01054 01055 else if (elem->point(5) == min_point) 01056 if (elem->point(4) == std::min(elem->point(4), elem->point(1))) 01057 { 01058 // Case 5 01059 xi_mapped = -xi; 01060 zeta_mapped = -zeta; 01061 } 01062 else 01063 { 01064 // Case 6 01065 xi_mapped = -zeta; 01066 zeta_mapped = -xi; 01067 } 01068 01069 else if (elem->point(4) == min_point) 01070 { 01071 if (elem->point(0) == std::min(elem->point(0), elem->point(5))) 01072 { 01073 // Case 7 01074 xi_mapped = -xi; 01075 zeta_mapped = zeta; 01076 } 01077 else 01078 { 01079 // Case 8 01080 xi_mapped = xi; 01081 zeta_mapped = -zeta; 01082 } 01083 } 01084 } 01085 01086 01087 // Face 2 01088 else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] >= 2)) 01089 { 01090 const Point min_point = std::min(elem->point(1), 01091 std::min(elem->point(2), 01092 std::min(elem->point(6), 01093 elem->point(5)))); 01094 if (elem->point(1) == min_point) 01095 if (elem->point(2) == std::min(elem->point(2), elem->point(5))) 01096 { 01097 // Case 1 01098 eta_mapped = eta; 01099 zeta_mapped = zeta; 01100 } 01101 else 01102 { 01103 // Case 2 01104 eta_mapped = zeta; 01105 zeta_mapped = eta; 01106 } 01107 01108 else if (elem->point(2) == min_point) 01109 if (elem->point(6) == std::min(elem->point(6), elem->point(1))) 01110 { 01111 // Case 3 01112 eta_mapped = zeta; 01113 zeta_mapped = -eta; 01114 } 01115 else 01116 { 01117 // Case 4 01118 eta_mapped = -eta; 01119 zeta_mapped = zeta; 01120 } 01121 01122 else if (elem->point(6) == min_point) 01123 if (elem->point(5) == std::min(elem->point(5), elem->point(2))) 01124 { 01125 // Case 5 01126 eta_mapped = -eta; 01127 zeta_mapped = -zeta; 01128 } 01129 else 01130 { 01131 // Case 6 01132 eta_mapped = -zeta; 01133 zeta_mapped = -eta; 01134 } 01135 01136 else if (elem->point(5) == min_point) 01137 { 01138 if (elem->point(1) == std::min(elem->point(1), elem->point(6))) 01139 { 01140 // Case 7 01141 eta_mapped = -zeta; 01142 zeta_mapped = eta; 01143 } 01144 else 01145 { 01146 // Case 8 01147 eta_mapped = eta; 01148 zeta_mapped = -zeta; 01149 } 01150 } 01151 } 01152 01153 01154 // Face 3 01155 else if ((i1[i] == 1) && (i0[i] >= 2) && (i2[i] >= 2)) 01156 { 01157 const Point min_point = std::min(elem->point(2), 01158 std::min(elem->point(3), 01159 std::min(elem->point(7), 01160 elem->point(6)))); 01161 if (elem->point(3) == min_point) 01162 if (elem->point(2) == std::min(elem->point(2), elem->point(7))) 01163 { 01164 // Case 1 01165 xi_mapped = xi; 01166 zeta_mapped = zeta; 01167 } 01168 else 01169 { 01170 // Case 2 01171 xi_mapped = zeta; 01172 zeta_mapped = xi; 01173 } 01174 01175 else if (elem->point(7) == min_point) 01176 if (elem->point(3) == std::min(elem->point(3), elem->point(6))) 01177 { 01178 // Case 3 01179 xi_mapped = -zeta; 01180 zeta_mapped = xi; 01181 } 01182 else 01183 { 01184 // Case 4 01185 xi_mapped = xi; 01186 zeta_mapped = -zeta; 01187 } 01188 01189 else if (elem->point(6) == min_point) 01190 if (elem->point(7) == std::min(elem->point(7), elem->point(2))) 01191 { 01192 // Case 5 01193 xi_mapped = -xi; 01194 zeta_mapped = -zeta; 01195 } 01196 else 01197 { 01198 // Case 6 01199 xi_mapped = -zeta; 01200 zeta_mapped = -xi; 01201 } 01202 01203 else if (elem->point(2) == min_point) 01204 { 01205 if (elem->point(6) == std::min(elem->point(3), elem->point(6))) 01206 { 01207 // Case 7 01208 xi_mapped = zeta; 01209 zeta_mapped = -xi; 01210 } 01211 else 01212 { 01213 // Case 8 01214 xi_mapped = -xi; 01215 zeta_mapped = zeta; 01216 } 01217 } 01218 } 01219 01220 01221 // Face 4 01222 else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] >= 2)) 01223 { 01224 const Point min_point = std::min(elem->point(3), 01225 std::min(elem->point(0), 01226 std::min(elem->point(4), 01227 elem->point(7)))); 01228 if (elem->point(0) == min_point) 01229 if (elem->point(3) == std::min(elem->point(3), elem->point(4))) 01230 { 01231 // Case 1 01232 eta_mapped = eta; 01233 zeta_mapped = zeta; 01234 } 01235 else 01236 { 01237 // Case 2 01238 eta_mapped = zeta; 01239 zeta_mapped = eta; 01240 } 01241 01242 else if (elem->point(4) == min_point) 01243 if (elem->point(0) == std::min(elem->point(0), elem->point(7))) 01244 { 01245 // Case 3 01246 eta_mapped = -zeta; 01247 zeta_mapped = eta; 01248 } 01249 else 01250 { 01251 // Case 4 01252 eta_mapped = eta; 01253 zeta_mapped = -zeta; 01254 } 01255 01256 else if (elem->point(7) == min_point) 01257 if (elem->point(4) == std::min(elem->point(4), elem->point(3))) 01258 { 01259 // Case 5 01260 eta_mapped = -eta; 01261 zeta_mapped = -zeta; 01262 } 01263 else 01264 { 01265 // Case 6 01266 eta_mapped = -zeta; 01267 zeta_mapped = -eta; 01268 } 01269 01270 else if (elem->point(3) == min_point) 01271 { 01272 if (elem->point(7) == std::min(elem->point(7), elem->point(0))) 01273 { 01274 // Case 7 01275 eta_mapped = zeta; 01276 zeta_mapped = -eta; 01277 } 01278 else 01279 { 01280 // Case 8 01281 eta_mapped = -eta; 01282 zeta_mapped = zeta; 01283 } 01284 } 01285 } 01286 01287 01288 // Face 5 01289 else if ((i2[i] == 1) && (i0[i] >= 2) && (i1[i] >= 2)) 01290 { 01291 const Point min_point = std::min(elem->point(4), 01292 std::min(elem->point(5), 01293 std::min(elem->point(6), 01294 elem->point(7)))); 01295 if (elem->point(4) == min_point) 01296 if (elem->point(5) == std::min(elem->point(5), elem->point(7))) 01297 { 01298 // Case 1 01299 xi_mapped = xi; 01300 eta_mapped = eta; 01301 } 01302 else 01303 { 01304 // Case 2 01305 xi_mapped = eta; 01306 eta_mapped = xi; 01307 } 01308 01309 else if (elem->point(5) == min_point) 01310 if (elem->point(6) == std::min(elem->point(6), elem->point(4))) 01311 { 01312 // Case 3 01313 xi_mapped = eta; 01314 eta_mapped = -xi; 01315 } 01316 else 01317 { 01318 // Case 4 01319 xi_mapped = -xi; 01320 eta_mapped = eta; 01321 } 01322 01323 else if (elem->point(6) == min_point) 01324 if (elem->point(7) == std::min(elem->point(7), elem->point(5))) 01325 { 01326 // Case 5 01327 xi_mapped = -xi; 01328 eta_mapped = -eta; 01329 } 01330 else 01331 { 01332 // Case 6 01333 xi_mapped = -eta; 01334 eta_mapped = -xi; 01335 } 01336 01337 else if (elem->point(7) == min_point) 01338 { 01339 if (elem->point(4) == std::min(elem->point(4), elem->point(6))) 01340 { 01341 // Case 7 01342 xi_mapped = -eta; 01343 eta_mapped = xi; 01344 } 01345 else 01346 { 01347 // Case 8 01348 xi_mapped = xi; 01349 eta_mapped = eta; 01350 } 01351 } 01352 } 01353 01354 01355 } 01356 01357 01358 return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi_mapped)* 01359 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta_mapped)* 01360 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta_mapped)); 01361 } 01362 01363 01364 default: 01365 libmesh_error_msg("Invalid element type = " << type); 01366 } 01367 } 01368 01369 01370 default: 01371 libmesh_error_msg("Invalid totalorder = " << totalorder); 01372 } 01373 01374 #endif 01375 01376 libmesh_error_msg("We'll never get here!"); 01377 return 0.; 01378 } 01379 01380 01381 01382 01383 template <> 01384 Real FE<3,BERNSTEIN>::shape_deriv(const ElemType, 01385 const Order, 01386 const unsigned int, 01387 const unsigned int, 01388 const Point& ) 01389 { 01390 libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge and face orientation is needed."); 01391 return 0.; 01392 } 01393 01394 01395 01396 template <> 01397 Real FE<3,BERNSTEIN>::shape_deriv(const Elem* elem, 01398 const Order order, 01399 const unsigned int i, 01400 const unsigned int j, 01401 const Point& p) 01402 { 01403 01404 #if LIBMESH_DIM == 3 01405 libmesh_assert(elem); 01406 const ElemType type = elem->type(); 01407 01408 const Order totalorder = static_cast<Order>(order + elem->p_level()); 01409 01410 libmesh_assert_less (j, 3); 01411 01412 switch (totalorder) 01413 { 01414 // 1st order Bernstein. 01415 case FIRST: 01416 { 01417 switch (type) 01418 { 01419 // Bernstein shape functions on the tetrahedron. 01420 case TET4: 01421 case TET10: 01422 { 01423 // I have been lazy here and am using finite differences 01424 // to compute the derivatives! 01425 const Real eps = 1.e-6; 01426 01427 libmesh_assert_less (i, 4); 01428 libmesh_assert_less (j, 3); 01429 01430 01431 switch (j) 01432 { 01433 // d()/dxi 01434 case 0: 01435 { 01436 const Point pp(p(0)+eps, p(1), p(2)); 01437 const Point pm(p(0)-eps, p(1), p(2)); 01438 01439 return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) - 01440 FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01441 } 01442 01443 // d()/deta 01444 case 1: 01445 { 01446 const Point pp(p(0), p(1)+eps, p(2)); 01447 const Point pm(p(0), p(1)-eps, p(2)); 01448 01449 return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) - 01450 FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01451 } 01452 // d()/dzeta 01453 case 2: 01454 { 01455 const Point pp(p(0), p(1), p(2)+eps); 01456 const Point pm(p(0), p(1), p(2)-eps); 01457 01458 return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) - 01459 FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01460 } 01461 default: 01462 libmesh_error_msg("Invalid derivative index j = " << j); 01463 } 01464 } 01465 01466 01467 01468 01469 // Bernstein shape functions on the hexahedral. 01470 case HEX20: 01471 case HEX27: 01472 { 01473 libmesh_assert_less (i, 8); 01474 01475 // Compute hex shape functions as a tensor-product 01476 const Real xi = p(0); 01477 const Real eta = p(1); 01478 const Real zeta = p(2); 01479 01480 // The only way to make any sense of this 01481 // is to look at the mgflo/mg2/mgf documentation 01482 // and make the cut-out cube! 01483 // 0 1 2 3 4 5 6 7 01484 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0}; 01485 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1}; 01486 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1}; 01487 01488 switch (j) 01489 { 01490 // d()/dxi 01491 case 0: 01492 return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)* 01493 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)* 01494 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta)); 01495 01496 // d()/deta 01497 case 1: 01498 return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)* 01499 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)* 01500 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta)); 01501 01502 // d()/dzeta 01503 case 2: 01504 return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)* 01505 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)* 01506 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta)); 01507 01508 default: 01509 libmesh_error_msg("Invalid derivative index j = " << j); 01510 } 01511 } 01512 01513 default: 01514 libmesh_error_msg("Invalid element type = " << type); 01515 } 01516 } 01517 01518 01519 01520 01521 case SECOND: 01522 { 01523 switch (type) 01524 { 01525 // Bernstein shape functions on the tetrahedron. 01526 case TET10: 01527 { 01528 // I have been lazy here and am using finite differences 01529 // to compute the derivatives! 01530 const Real eps = 1.e-6; 01531 01532 libmesh_assert_less (i, 10); 01533 libmesh_assert_less (j, 3); 01534 01535 01536 switch (j) 01537 { 01538 // d()/dxi 01539 case 0: 01540 { 01541 const Point pp(p(0)+eps, p(1), p(2)); 01542 const Point pm(p(0)-eps, p(1), p(2)); 01543 01544 return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) - 01545 FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01546 } 01547 01548 // d()/deta 01549 case 1: 01550 { 01551 const Point pp(p(0), p(1)+eps, p(2)); 01552 const Point pm(p(0), p(1)-eps, p(2)); 01553 01554 return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) - 01555 FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01556 } 01557 // d()/dzeta 01558 case 2: 01559 { 01560 const Point pp(p(0), p(1), p(2)+eps); 01561 const Point pm(p(0), p(1), p(2)-eps); 01562 01563 return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) - 01564 FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01565 } 01566 default: 01567 libmesh_error_msg("Invalid derivative index j = " << j); 01568 } 01569 } 01570 01571 // Bernstein shape functions on the hexahedral. 01572 case HEX20: 01573 { 01574 libmesh_assert_less (i, 20); 01575 01576 // Compute hex shape functions as a tensor-product 01577 const Real xi = p(0); 01578 const Real eta = p(1); 01579 const Real zeta = p(2); 01580 01581 // The only way to make any sense of this 01582 // is to look at the mgflo/mg2/mgf documentation 01583 // and make the cut-out cube! 01584 // 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 01585 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}; 01586 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}; 01587 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}; 01588 static const Real scal20[] = {-0.25, -0.25, -0.25, -0.25, 0, 0, 0, 0, 0.5, 0.5, 0.5, 0.5, 0, 0, 0, 0, 0, 0, 0, 0}; 01589 static const Real scal21[] = {-0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0, 0, 0}; 01590 static const Real scal22[] = {0, -0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0, 0}; 01591 static const Real scal23[] = {0, 0, -0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0}; 01592 static const Real scal24[] = {-0.25, 0, 0, -0.25, -0.25, 0, 0, -0.25, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0, 0, 0, 0.5}; 01593 static const Real scal25[] = {0, 0, 0, 0, -0.25, -0.25, -0.25, -0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0.5, 0.5}; 01594 static const Real scal26[] = {-0.25, -0.25, -0.25, -0.25, -0.25, -0.25, -0.25, -0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}; 01595 01596 switch (j) 01597 { 01598 // d()/dxi 01599 case 0: 01600 return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)* 01601 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)* 01602 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta) 01603 +scal20[i]* 01604 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[20], 0, xi)* 01605 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[20], eta)* 01606 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[20], zeta) 01607 +scal21[i]* 01608 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[21], 0, xi)* 01609 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[21], eta)* 01610 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[21], zeta) 01611 +scal22[i]* 01612 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[22], 0, xi)* 01613 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[22], eta)* 01614 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[22], zeta) 01615 +scal23[i]* 01616 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[23], 0, xi)* 01617 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[23], eta)* 01618 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[23], zeta) 01619 +scal24[i]* 01620 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[24], 0, xi)* 01621 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[24], eta)* 01622 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[24], zeta) 01623 +scal25[i]* 01624 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[25], 0, xi)* 01625 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[25], eta)* 01626 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[25], zeta) 01627 +scal26[i]* 01628 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[26], 0, xi)* 01629 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[26], eta)* 01630 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[26], zeta)); 01631 01632 // d()/deta 01633 case 1: 01634 return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)* 01635 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)* 01636 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta) 01637 +scal20[i]* 01638 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[20], xi)* 01639 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[20], 0, eta)* 01640 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[20], zeta) 01641 +scal21[i]* 01642 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[21], xi)* 01643 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[21], 0, eta)* 01644 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[21], zeta) 01645 +scal22[i]* 01646 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[22], xi)* 01647 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[22], 0, eta)* 01648 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[22], zeta) 01649 +scal23[i]* 01650 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[23], xi)* 01651 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[23], 0, eta)* 01652 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[23], zeta) 01653 +scal24[i]* 01654 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[24], xi)* 01655 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[24], 0, eta)* 01656 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[24], zeta) 01657 +scal25[i]* 01658 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[25], xi)* 01659 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[25], 0, eta)* 01660 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[25], zeta) 01661 +scal26[i]* 01662 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[26], xi)* 01663 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[26], 0, eta)* 01664 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[26], zeta)); 01665 01666 // d()/dzeta 01667 case 2: 01668 return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)* 01669 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)* 01670 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta) 01671 +scal20[i]* 01672 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[20], xi)* 01673 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[20], eta)* 01674 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[20], 0, zeta) 01675 +scal21[i]* 01676 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[21], xi)* 01677 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[21], eta)* 01678 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[21], 0, zeta) 01679 +scal22[i]* 01680 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[22], xi)* 01681 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[22], eta)* 01682 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[22], 0, zeta) 01683 +scal23[i]* 01684 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[23], xi)* 01685 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[23], eta)* 01686 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[23], 0, zeta) 01687 +scal24[i]* 01688 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[24], xi)* 01689 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[24], eta)* 01690 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[24], 0, zeta) 01691 +scal25[i]* 01692 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[25], xi)* 01693 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[25], eta)* 01694 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[25], 0, zeta) 01695 +scal26[i]* 01696 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[26], xi)* 01697 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[26], eta)* 01698 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[26], 0, zeta)); 01699 01700 default: 01701 libmesh_error_msg("Invalid derivative index j = " << j); 01702 } 01703 } 01704 01705 // Bernstein shape functions on the hexahedral. 01706 case HEX27: 01707 { 01708 libmesh_assert_less (i, 27); 01709 01710 // Compute hex shape functions as a tensor-product 01711 const Real xi = p(0); 01712 const Real eta = p(1); 01713 const Real zeta = p(2); 01714 01715 // The only way to make any sense of this 01716 // is to look at the mgflo/mg2/mgf documentation 01717 // and make the cut-out cube! 01718 // 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 01719 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}; 01720 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}; 01721 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}; 01722 01723 switch (j) 01724 { 01725 // d()/dxi 01726 case 0: 01727 return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)* 01728 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)* 01729 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta)); 01730 01731 // d()/deta 01732 case 1: 01733 return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)* 01734 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)* 01735 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta)); 01736 01737 // d()/dzeta 01738 case 2: 01739 return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)* 01740 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)* 01741 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta)); 01742 01743 default: 01744 libmesh_error_msg("Invalid derivative index j = " << j); 01745 } 01746 } 01747 01748 01749 default: 01750 libmesh_error_msg("Invalid element type = " << type); 01751 } 01752 } 01753 01754 01755 01756 // 3rd-order Bernstein. 01757 case THIRD: 01758 { 01759 switch (type) 01760 { 01761 01762 // // Bernstein shape functions derivatives. 01763 // case TET10: 01764 // { 01765 // // I have been lazy here and am using finite differences 01766 // // to compute the derivatives! 01767 // const Real eps = 1.e-6; 01768 01769 // libmesh_assert_less (i, 20); 01770 // libmesh_assert_less (j, 3); 01771 01772 // switch (j) 01773 // { 01774 // // d()/dxi 01775 // case 0: 01776 // { 01777 // const Point pp(p(0)+eps, p(1), p(2)); 01778 // const Point pm(p(0)-eps, p(1), p(2)); 01779 01780 // return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) - 01781 // FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01782 // } 01783 01784 // // d()/deta 01785 // case 1: 01786 // { 01787 // const Point pp(p(0), p(1)+eps, p(2)); 01788 // const Point pm(p(0), p(1)-eps, p(2)); 01789 01790 // return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) - 01791 // FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01792 // } 01793 // // d()/dzeta 01794 // case 2: 01795 // { 01796 // const Point pp(p(0), p(1), p(2)+eps); 01797 // const Point pm(p(0), p(1), p(2)-eps); 01798 01799 // return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) - 01800 // FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01801 // } 01802 // default: 01803 // libmesh_error_msg("Invalid derivative index j = " << j); 01804 // } 01805 01806 01807 // } 01808 01809 01810 // Bernstein shape functions on the hexahedral. 01811 case HEX27: 01812 { 01813 // I have been lazy here and am using finite differences 01814 // to compute the derivatives! 01815 const Real eps = 1.e-6; 01816 01817 libmesh_assert_less (i, 64); 01818 libmesh_assert_less (j, 3); 01819 01820 switch (j) 01821 { 01822 // d()/dxi 01823 case 0: 01824 { 01825 const Point pp(p(0)+eps, p(1), p(2)); 01826 const Point pm(p(0)-eps, p(1), p(2)); 01827 01828 return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) - 01829 FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01830 } 01831 01832 // d()/deta 01833 case 1: 01834 { 01835 const Point pp(p(0), p(1)+eps, p(2)); 01836 const Point pm(p(0), p(1)-eps, p(2)); 01837 01838 return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) - 01839 FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01840 } 01841 // d()/dzeta 01842 case 2: 01843 { 01844 const Point pp(p(0), p(1), p(2)+eps); 01845 const Point pm(p(0), p(1), p(2)-eps); 01846 01847 return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) - 01848 FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01849 } 01850 default: 01851 libmesh_error_msg("Invalid derivative index j = " << j); 01852 } 01853 01854 } 01855 01856 // // Compute hex shape functions as a tensor-product 01857 // const Real xi = p(0); 01858 // const Real eta = p(1); 01859 // const Real zeta = p(2); 01860 // Real xi_mapped = p(0); 01861 // Real eta_mapped = p(1); 01862 // Real zeta_mapped = p(2); 01863 01864 // // The only way to make any sense of this 01865 // // is to look at the mgflo/mg2/mgf documentation 01866 // // and make the cut-out cube! 01867 // // Nodes 0 1 2 3 4 5 6 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 20 20 21 21 21 21 22 22 22 22 23 23 23 23 24 24 24 24 25 25 25 25 26 26 26 26 26 26 26 26 01868 // // DOFS 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 27 18 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 60 62 63 01869 // static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3}; 01870 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3, 0, 0, 0, 0, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3}; 01871 // static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3}; 01872 01873 01874 01875 // // handle the edge orientation 01876 // { 01877 // // Edge 0 01878 // if ((i1[i] == 0) && (i2[i] == 0)) 01879 // { 01880 // if (elem->node(0) != std::min(elem->node(0), elem->node(1))) 01881 // xi_mapped = -xi; 01882 // } 01883 // // Edge 1 01884 // else if ((i0[i] == 1) && (i2[i] == 0)) 01885 // { 01886 // if (elem->node(1) != std::min(elem->node(1), elem->node(2))) 01887 // eta_mapped = -eta; 01888 // } 01889 // // Edge 2 01890 // else if ((i1[i] == 1) && (i2[i] == 0)) 01891 // { 01892 // if (elem->node(3) != std::min(elem->node(3), elem->node(2))) 01893 // xi_mapped = -xi; 01894 // } 01895 // // Edge 3 01896 // else if ((i0[i] == 0) && (i2[i] == 0)) 01897 // { 01898 // if (elem->node(0) != std::min(elem->node(0), elem->node(3))) 01899 // eta_mapped = -eta; 01900 // } 01901 // // Edge 4 01902 // else if ((i0[i] == 0) && (i1[i] == 0)) 01903 // { 01904 // if (elem->node(0) != std::min(elem->node(0), elem->node(4))) 01905 // zeta_mapped = -zeta; 01906 // } 01907 // // Edge 5 01908 // else if ((i0[i] == 1) && (i1[i] == 0)) 01909 // { 01910 // if (elem->node(1) != std::min(elem->node(1), elem->node(5))) 01911 // zeta_mapped = -zeta; 01912 // } 01913 // // Edge 6 01914 // else if ((i0[i] == 1) && (i1[i] == 1)) 01915 // { 01916 // if (elem->node(2) != std::min(elem->node(2), elem->node(6))) 01917 // zeta_mapped = -zeta; 01918 // } 01919 // // Edge 7 01920 // else if ((i0[i] == 0) && (i1[i] == 1)) 01921 // { 01922 // if (elem->node(3) != std::min(elem->node(3), elem->node(7))) 01923 // zeta_mapped = -zeta; 01924 // } 01925 // // Edge 8 01926 // else if ((i1[i] == 0) && (i2[i] == 1)) 01927 // { 01928 // if (elem->node(4) != std::min(elem->node(4), elem->node(5))) 01929 // xi_mapped = -xi; 01930 // } 01931 // // Edge 9 01932 // else if ((i0[i] == 1) && (i2[i] == 1)) 01933 // { 01934 // if (elem->node(5) != std::min(elem->node(5), elem->node(6))) 01935 // eta_mapped = -eta; 01936 // } 01937 // // Edge 10 01938 // else if ((i1[i] == 1) && (i2[i] == 1)) 01939 // { 01940 // if (elem->node(7) != std::min(elem->node(7), elem->node(6))) 01941 // xi_mapped = -xi; 01942 // } 01943 // // Edge 11 01944 // else if ((i0[i] == 0) && (i2[i] == 1)) 01945 // { 01946 // if (elem->node(4) != std::min(elem->node(4), elem->node(7))) 01947 // eta_mapped = -eta; 01948 // } 01949 // } 01950 01951 01952 // // handle the face orientation 01953 // { 01954 // // Face 0 01955 // if ( (i2[i] == 0) && (i0[i] >= 2) && (i1[i] >= 2)) 01956 // { 01957 // const unsigned int min_node = std::min(elem->node(1), 01958 // std::min(elem->node(2), 01959 // std::min(elem->node(0), 01960 // elem->node(3)))); 01961 // if (elem->node(0) == min_node) 01962 // if (elem->node(1) == std::min(elem->node(1), elem->node(3))) 01963 // { 01964 // // Case 1 01965 // xi_mapped = xi; 01966 // eta_mapped = eta; 01967 // } 01968 // else 01969 // { 01970 // // Case 2 01971 // xi_mapped = eta; 01972 // eta_mapped = xi; 01973 // } 01974 01975 // else if (elem->node(3) == min_node) 01976 // if (elem->node(0) == std::min(elem->node(0), elem->node(2))) 01977 // { 01978 // // Case 3 01979 // xi_mapped = -eta; 01980 // eta_mapped = xi; 01981 // } 01982 // else 01983 // { 01984 // // Case 4 01985 // xi_mapped = xi; 01986 // eta_mapped = -eta; 01987 // } 01988 01989 // else if (elem->node(2) == min_node) 01990 // if (elem->node(3) == std::min(elem->node(3), elem->node(1))) 01991 // { 01992 // // Case 5 01993 // xi_mapped = -xi; 01994 // eta_mapped = -eta; 01995 // } 01996 // else 01997 // { 01998 // // Case 6 01999 // xi_mapped = -eta; 02000 // eta_mapped = -xi; 02001 // } 02002 02003 // else if (elem->node(1) == min_node) 02004 // if (elem->node(2) == std::min(elem->node(2), elem->node(0))) 02005 // { 02006 // // Case 7 02007 // xi_mapped = eta; 02008 // eta_mapped = -xi; 02009 // } 02010 // else 02011 // { 02012 // // Case 8 02013 // xi_mapped = -xi; 02014 // eta_mapped = eta; 02015 // } 02016 // } 02017 02018 02019 // // Face 1 02020 // else if ((i1[i] == 0) && (i0[i] >= 2) && (i2[i] >= 2)) 02021 // { 02022 // const unsigned int min_node = std::min(elem->node(0), 02023 // std::min(elem->node(1), 02024 // std::min(elem->node(5), 02025 // elem->node(4)))); 02026 // if (elem->node(0) == min_node) 02027 // if (elem->node(1) == std::min(elem->node(1), elem->node(4))) 02028 // { 02029 // // Case 1 02030 // xi_mapped = xi; 02031 // zeta_mapped = zeta; 02032 // } 02033 // else 02034 // { 02035 // // Case 2 02036 // xi_mapped = zeta; 02037 // zeta_mapped = xi; 02038 // } 02039 02040 // else if (elem->node(1) == min_node) 02041 // if (elem->node(5) == std::min(elem->node(5), elem->node(0))) 02042 // { 02043 // // Case 3 02044 // xi_mapped = zeta; 02045 // zeta_mapped = -xi; 02046 // } 02047 // else 02048 // { 02049 // // Case 4 02050 // xi_mapped = -xi; 02051 // zeta_mapped = zeta; 02052 // } 02053 02054 // else if (elem->node(5) == min_node) 02055 // if (elem->node(4) == std::min(elem->node(4), elem->node(1))) 02056 // { 02057 // // Case 5 02058 // xi_mapped = -xi; 02059 // zeta_mapped = -zeta; 02060 // } 02061 // else 02062 // { 02063 // // Case 6 02064 // xi_mapped = -zeta; 02065 // zeta_mapped = -xi; 02066 // } 02067 02068 // else if (elem->node(4) == min_node) 02069 // if (elem->node(0) == std::min(elem->node(0), elem->node(5))) 02070 // { 02071 // // Case 7 02072 // xi_mapped = -xi; 02073 // zeta_mapped = zeta; 02074 // } 02075 // else 02076 // { 02077 // // Case 8 02078 // xi_mapped = xi; 02079 // zeta_mapped = -zeta; 02080 // } 02081 // } 02082 02083 02084 // // Face 2 02085 // else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] >= 2)) 02086 // { 02087 // const unsigned int min_node = std::min(elem->node(1), 02088 // std::min(elem->node(2), 02089 // std::min(elem->node(6), 02090 // elem->node(5)))); 02091 // if (elem->node(1) == min_node) 02092 // if (elem->node(2) == std::min(elem->node(2), elem->node(5))) 02093 // { 02094 // // Case 1 02095 // eta_mapped = eta; 02096 // zeta_mapped = zeta; 02097 // } 02098 // else 02099 // { 02100 // // Case 2 02101 // eta_mapped = zeta; 02102 // zeta_mapped = eta; 02103 // } 02104 02105 // else if (elem->node(2) == min_node) 02106 // if (elem->node(6) == std::min(elem->node(6), elem->node(1))) 02107 // { 02108 // // Case 3 02109 // eta_mapped = zeta; 02110 // zeta_mapped = -eta; 02111 // } 02112 // else 02113 // { 02114 // // Case 4 02115 // eta_mapped = -eta; 02116 // zeta_mapped = zeta; 02117 // } 02118 02119 // else if (elem->node(6) == min_node) 02120 // if (elem->node(5) == std::min(elem->node(5), elem->node(2))) 02121 // { 02122 // // Case 5 02123 // eta_mapped = -eta; 02124 // zeta_mapped = -zeta; 02125 // } 02126 // else 02127 // { 02128 // // Case 6 02129 // eta_mapped = -zeta; 02130 // zeta_mapped = -eta; 02131 // } 02132 02133 // else if (elem->node(5) == min_node) 02134 // if (elem->node(1) == std::min(elem->node(1), elem->node(6))) 02135 // { 02136 // // Case 7 02137 // eta_mapped = -zeta; 02138 // zeta_mapped = eta; 02139 // } 02140 // else 02141 // { 02142 // // Case 8 02143 // eta_mapped = eta; 02144 // zeta_mapped = -zeta; 02145 // } 02146 // } 02147 02148 02149 // // Face 3 02150 // else if ((i1[i] == 1) && (i0[i] >= 2) && (i2[i] >= 2)) 02151 // { 02152 // const unsigned int min_node = std::min(elem->node(2), 02153 // std::min(elem->node(3), 02154 // std::min(elem->node(7), 02155 // elem->node(6)))); 02156 // if (elem->node(3) == min_node) 02157 // if (elem->node(2) == std::min(elem->node(2), elem->node(7))) 02158 // { 02159 // // Case 1 02160 // xi_mapped = xi; 02161 // zeta_mapped = zeta; 02162 // } 02163 // else 02164 // { 02165 // // Case 2 02166 // xi_mapped = zeta; 02167 // zeta_mapped = xi; 02168 // } 02169 02170 // else if (elem->node(7) == min_node) 02171 // if (elem->node(3) == std::min(elem->node(3), elem->node(6))) 02172 // { 02173 // // Case 3 02174 // xi_mapped = -zeta; 02175 // zeta_mapped = xi; 02176 // } 02177 // else 02178 // { 02179 // // Case 4 02180 // xi_mapped = xi; 02181 // zeta_mapped = -zeta; 02182 // } 02183 02184 // else if (elem->node(6) == min_node) 02185 // if (elem->node(7) == std::min(elem->node(7), elem->node(2))) 02186 // { 02187 // // Case 5 02188 // xi_mapped = -xi; 02189 // zeta_mapped = -zeta; 02190 // } 02191 // else 02192 // { 02193 // // Case 6 02194 // xi_mapped = -zeta; 02195 // zeta_mapped = -xi; 02196 // } 02197 02198 // else if (elem->node(2) == min_node) 02199 // if (elem->node(6) == std::min(elem->node(3), elem->node(6))) 02200 // { 02201 // // Case 7 02202 // xi_mapped = zeta; 02203 // zeta_mapped = -xi; 02204 // } 02205 // else 02206 // { 02207 // // Case 8 02208 // xi_mapped = -xi; 02209 // zeta_mapped = zeta; 02210 // } 02211 // } 02212 02213 02214 // // Face 4 02215 // else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] >= 2)) 02216 // { 02217 // const unsigned int min_node = std::min(elem->node(3), 02218 // std::min(elem->node(0), 02219 // std::min(elem->node(4), 02220 // elem->node(7)))); 02221 // if (elem->node(0) == min_node) 02222 // if (elem->node(3) == std::min(elem->node(3), elem->node(4))) 02223 // { 02224 // // Case 1 02225 // eta_mapped = eta; 02226 // zeta_mapped = zeta; 02227 // } 02228 // else 02229 // { 02230 // // Case 2 02231 // eta_mapped = zeta; 02232 // zeta_mapped = eta; 02233 // } 02234 02235 // else if (elem->node(4) == min_node) 02236 // if (elem->node(0) == std::min(elem->node(0), elem->node(7))) 02237 // { 02238 // // Case 3 02239 // eta_mapped = -zeta; 02240 // zeta_mapped = eta; 02241 // } 02242 // else 02243 // { 02244 // // Case 4 02245 // eta_mapped = eta; 02246 // zeta_mapped = -zeta; 02247 // } 02248 02249 // else if (elem->node(7) == min_node) 02250 // if (elem->node(4) == std::min(elem->node(4), elem->node(3))) 02251 // { 02252 // // Case 5 02253 // eta_mapped = -eta; 02254 // zeta_mapped = -zeta; 02255 // } 02256 // else 02257 // { 02258 // // Case 6 02259 // eta_mapped = -zeta; 02260 // zeta_mapped = -eta; 02261 // } 02262 02263 // else if (elem->node(3) == min_node) 02264 // if (elem->node(7) == std::min(elem->node(7), elem->node(0))) 02265 // { 02266 // // Case 7 02267 // eta_mapped = zeta; 02268 // zeta_mapped = -eta; 02269 // } 02270 // else 02271 // { 02272 // // Case 8 02273 // eta_mapped = -eta; 02274 // zeta_mapped = zeta; 02275 // } 02276 // } 02277 02278 02279 // // Face 5 02280 // else if ((i2[i] == 1) && (i0[i] >= 2) && (i1[i] >= 2)) 02281 // { 02282 // const unsigned int min_node = std::min(elem->node(4), 02283 // std::min(elem->node(5), 02284 // std::min(elem->node(6), 02285 // elem->node(7)))); 02286 // if (elem->node(4) == min_node) 02287 // if (elem->node(5) == std::min(elem->node(5), elem->node(7))) 02288 // { 02289 // // Case 1 02290 // xi_mapped = xi; 02291 // eta_mapped = eta; 02292 // } 02293 // else 02294 // { 02295 // // Case 2 02296 // xi_mapped = eta; 02297 // eta_mapped = xi; 02298 // } 02299 02300 // else if (elem->node(5) == min_node) 02301 // if (elem->node(6) == std::min(elem->node(6), elem->node(4))) 02302 // { 02303 // // Case 3 02304 // xi_mapped = eta; 02305 // eta_mapped = -xi; 02306 // } 02307 // else 02308 // { 02309 // // Case 4 02310 // xi_mapped = -xi; 02311 // eta_mapped = eta; 02312 // } 02313 02314 // else if (elem->node(6) == min_node) 02315 // if (elem->node(7) == std::min(elem->node(7), elem->node(5))) 02316 // { 02317 // // Case 5 02318 // xi_mapped = -xi; 02319 // eta_mapped = -eta; 02320 // } 02321 // else 02322 // { 02323 // // Case 6 02324 // xi_mapped = -eta; 02325 // eta_mapped = -xi; 02326 // } 02327 02328 // else if (elem->node(7) == min_node) 02329 // if (elem->node(4) == std::min(elem->node(4), elem->node(6))) 02330 // { 02331 // // Case 7 02332 // xi_mapped = -eta; 02333 // eta_mapped = xi; 02334 // } 02335 // else 02336 // { 02337 // // Case 8 02338 // xi_mapped = xi; 02339 // eta_mapped = eta; 02340 // } 02341 // } 02342 02343 02344 // } 02345 02346 02347 02348 // libmesh_assert_less (j, 3); 02349 02350 // switch (j) 02351 // { 02352 // // d()/dxi 02353 // case 0: 02354 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi_mapped)* 02355 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)* 02356 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped)); 02357 02358 // // d()/deta 02359 // case 1: 02360 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)* 02361 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta_mapped)* 02362 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped)); 02363 02364 // // d()/dzeta 02365 // case 2: 02366 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)* 02367 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)* 02368 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta_mapped)); 02369 02370 // default: 02371 // libmesh_error_msg("Invalid derivative index j = " << j); 02372 // } 02373 // } 02374 02375 02376 default: 02377 libmesh_error_msg("Invalid element type = " << type); 02378 } 02379 } 02380 02381 // 4th-order Bernstein. 02382 case FOURTH: 02383 { 02384 switch (type) 02385 { 02386 02387 // Bernstein shape functions derivatives on the hexahedral. 02388 case HEX27: 02389 { 02390 const Real eps = 1.e-6; 02391 02392 libmesh_assert_less (i, 125); 02393 libmesh_assert_less (j, 3); 02394 02395 switch (j) 02396 { 02397 // d()/dxi 02398 case 0: 02399 { 02400 const Point pp(p(0)+eps, p(1), p(2)); 02401 const Point pm(p(0)-eps, p(1), p(2)); 02402 02403 return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) - 02404 FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 02405 } 02406 02407 // d()/deta 02408 case 1: 02409 { 02410 const Point pp(p(0), p(1)+eps, p(2)); 02411 const Point pm(p(0), p(1)-eps, p(2)); 02412 02413 return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) - 02414 FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 02415 } 02416 // d()/dzeta 02417 case 2: 02418 { 02419 const Point pp(p(0), p(1), p(2)+eps); 02420 const Point pm(p(0), p(1), p(2)-eps); 02421 02422 return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) - 02423 FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 02424 } 02425 default: 02426 libmesh_error_msg("Invalid derivative index j = " << j); 02427 } 02428 } 02429 02430 // // Compute hex shape functions as a tensor-product 02431 // const Real xi = p(0); 02432 // const Real eta = p(1); 02433 // const Real zeta = p(2); 02434 // Real xi_mapped = p(0); 02435 // Real eta_mapped = p(1); 02436 // Real zeta_mapped = p(2); 02437 02438 // // The only way to make any sense of this 02439 // // is to look at the mgflo/mg2/mgf documentation 02440 // // and make the cut-out cube! 02441 // // Nodes 0 1 2 3 4 5 6 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 20 20 21 21 21 21 22 22 22 22 23 23 23 23 24 24 24 24 25 25 25 25 26 26 26 26 26 26 26 26 02442 // // DOFS 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 27 18 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 60 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 02443 // static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 2, 3, 4, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4}; 02444 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4}; 02445 // static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4}; 02446 02447 02448 02449 // // handle the edge orientation 02450 // { 02451 // // Edge 0 02452 // if ((i1[i] == 0) && (i2[i] == 0)) 02453 // { 02454 // if (elem->node(0) != std::min(elem->node(0), elem->node(1))) 02455 // xi_mapped = -xi; 02456 // } 02457 // // Edge 1 02458 // else if ((i0[i] == 1) && (i2[i] == 0)) 02459 // { 02460 // if (elem->node(1) != std::min(elem->node(1), elem->node(2))) 02461 // eta_mapped = -eta; 02462 // } 02463 // // Edge 2 02464 // else if ((i1[i] == 1) && (i2[i] == 0)) 02465 // { 02466 // if (elem->node(3) != std::min(elem->node(3), elem->node(2))) 02467 // xi_mapped = -xi; 02468 // } 02469 // // Edge 3 02470 // else if ((i0[i] == 0) && (i2[i] == 0)) 02471 // { 02472 // if (elem->node(0) != std::min(elem->node(0), elem->node(3))) 02473 // eta_mapped = -eta; 02474 // } 02475 // // Edge 4 02476 // else if ((i0[i] == 0) && (i1[i] == 0)) 02477 // { 02478 // if (elem->node(0) != std::min(elem->node(0), elem->node(4))) 02479 // zeta_mapped = -zeta; 02480 // } 02481 // // Edge 5 02482 // else if ((i0[i] == 1) && (i1[i] == 0)) 02483 // { 02484 // if (elem->node(1) != std::min(elem->node(1), elem->node(5))) 02485 // zeta_mapped = -zeta; 02486 // } 02487 // // Edge 6 02488 // else if ((i0[i] == 1) && (i1[i] == 1)) 02489 // { 02490 // if (elem->node(2) != std::min(elem->node(2), elem->node(6))) 02491 // zeta_mapped = -zeta; 02492 // } 02493 // // Edge 7 02494 // else if ((i0[i] == 0) && (i1[i] == 1)) 02495 // { 02496 // if (elem->node(3) != std::min(elem->node(3), elem->node(7))) 02497 // zeta_mapped = -zeta; 02498 // } 02499 // // Edge 8 02500 // else if ((i1[i] == 0) && (i2[i] == 1)) 02501 // { 02502 // if (elem->node(4) != std::min(elem->node(4), elem->node(5))) 02503 // xi_mapped = -xi; 02504 // } 02505 // // Edge 9 02506 // else if ((i0[i] == 1) && (i2[i] == 1)) 02507 // { 02508 // if (elem->node(5) != std::min(elem->node(5), elem->node(6))) 02509 // eta_mapped = -eta; 02510 // } 02511 // // Edge 10 02512 // else if ((i1[i] == 1) && (i2[i] == 1)) 02513 // { 02514 // if (elem->node(7) != std::min(elem->node(7), elem->node(6))) 02515 // xi_mapped = -xi; 02516 // } 02517 // // Edge 11 02518 // else if ((i0[i] == 0) && (i2[i] == 1)) 02519 // { 02520 // if (elem->node(4) != std::min(elem->node(4), elem->node(7))) 02521 // eta_mapped = -eta; 02522 // } 02523 // } 02524 02525 02526 // // handle the face orientation 02527 // { 02528 // // Face 0 02529 // if ( (i2[i] == 0) && (i0[i] >= 2) && (i1[i] >= 2)) 02530 // { 02531 // const unsigned int min_node = std::min(elem->node(1), 02532 // std::min(elem->node(2), 02533 // std::min(elem->node(0), 02534 // elem->node(3)))); 02535 // if (elem->node(0) == min_node) 02536 // if (elem->node(1) == std::min(elem->node(1), elem->node(3))) 02537 // { 02538 // // Case 1 02539 // xi_mapped = xi; 02540 // eta_mapped = eta; 02541 // } 02542 // else 02543 // { 02544 // // Case 2 02545 // xi_mapped = eta; 02546 // eta_mapped = xi; 02547 // } 02548 02549 // else if (elem->node(3) == min_node) 02550 // if (elem->node(0) == std::min(elem->node(0), elem->node(2))) 02551 // { 02552 // // Case 3 02553 // xi_mapped = -eta; 02554 // eta_mapped = xi; 02555 // } 02556 // else 02557 // { 02558 // // Case 4 02559 // xi_mapped = xi; 02560 // eta_mapped = -eta; 02561 // } 02562 02563 // else if (elem->node(2) == min_node) 02564 // if (elem->node(3) == std::min(elem->node(3), elem->node(1))) 02565 // { 02566 // // Case 5 02567 // xi_mapped = -xi; 02568 // eta_mapped = -eta; 02569 // } 02570 // else 02571 // { 02572 // // Case 6 02573 // xi_mapped = -eta; 02574 // eta_mapped = -xi; 02575 // } 02576 02577 // else if (elem->node(1) == min_node) 02578 // if (elem->node(2) == std::min(elem->node(2), elem->node(0))) 02579 // { 02580 // // Case 7 02581 // xi_mapped = eta; 02582 // eta_mapped = -xi; 02583 // } 02584 // else 02585 // { 02586 // // Case 8 02587 // xi_mapped = -xi; 02588 // eta_mapped = eta; 02589 // } 02590 // } 02591 02592 02593 // // Face 1 02594 // else if ((i1[i] == 0) && (i0[i] >= 2) && (i2[i] >= 2)) 02595 // { 02596 // const unsigned int min_node = std::min(elem->node(0), 02597 // std::min(elem->node(1), 02598 // std::min(elem->node(5), 02599 // elem->node(4)))); 02600 // if (elem->node(0) == min_node) 02601 // if (elem->node(1) == std::min(elem->node(1), elem->node(4))) 02602 // { 02603 // // Case 1 02604 // xi_mapped = xi; 02605 // zeta_mapped = zeta; 02606 // } 02607 // else 02608 // { 02609 // // Case 2 02610 // xi_mapped = zeta; 02611 // zeta_mapped = xi; 02612 // } 02613 02614 // else if (elem->node(1) == min_node) 02615 // if (elem->node(5) == std::min(elem->node(5), elem->node(0))) 02616 // { 02617 // // Case 3 02618 // xi_mapped = zeta; 02619 // zeta_mapped = -xi; 02620 // } 02621 // else 02622 // { 02623 // // Case 4 02624 // xi_mapped = -xi; 02625 // zeta_mapped = zeta; 02626 // } 02627 02628 // else if (elem->node(5) == min_node) 02629 // if (elem->node(4) == std::min(elem->node(4), elem->node(1))) 02630 // { 02631 // // Case 5 02632 // xi_mapped = -xi; 02633 // zeta_mapped = -zeta; 02634 // } 02635 // else 02636 // { 02637 // // Case 6 02638 // xi_mapped = -zeta; 02639 // zeta_mapped = -xi; 02640 // } 02641 02642 // else if (elem->node(4) == min_node) 02643 // if (elem->node(0) == std::min(elem->node(0), elem->node(5))) 02644 // { 02645 // // Case 7 02646 // xi_mapped = -xi; 02647 // zeta_mapped = zeta; 02648 // } 02649 // else 02650 // { 02651 // // Case 8 02652 // xi_mapped = xi; 02653 // zeta_mapped = -zeta; 02654 // } 02655 // } 02656 02657 02658 // // Face 2 02659 // else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] >= 2)) 02660 // { 02661 // const unsigned int min_node = std::min(elem->node(1), 02662 // std::min(elem->node(2), 02663 // std::min(elem->node(6), 02664 // elem->node(5)))); 02665 // if (elem->node(1) == min_node) 02666 // if (elem->node(2) == std::min(elem->node(2), elem->node(5))) 02667 // { 02668 // // Case 1 02669 // eta_mapped = eta; 02670 // zeta_mapped = zeta; 02671 // } 02672 // else 02673 // { 02674 // // Case 2 02675 // eta_mapped = zeta; 02676 // zeta_mapped = eta; 02677 // } 02678 02679 // else if (elem->node(2) == min_node) 02680 // if (elem->node(6) == std::min(elem->node(6), elem->node(1))) 02681 // { 02682 // // Case 3 02683 // eta_mapped = zeta; 02684 // zeta_mapped = -eta; 02685 // } 02686 // else 02687 // { 02688 // // Case 4 02689 // eta_mapped = -eta; 02690 // zeta_mapped = zeta; 02691 // } 02692 02693 // else if (elem->node(6) == min_node) 02694 // if (elem->node(5) == std::min(elem->node(5), elem->node(2))) 02695 // { 02696 // // Case 5 02697 // eta_mapped = -eta; 02698 // zeta_mapped = -zeta; 02699 // } 02700 // else 02701 // { 02702 // // Case 6 02703 // eta_mapped = -zeta; 02704 // zeta_mapped = -eta; 02705 // } 02706 02707 // else if (elem->node(5) == min_node) 02708 // if (elem->node(1) == std::min(elem->node(1), elem->node(6))) 02709 // { 02710 // // Case 7 02711 // eta_mapped = -zeta; 02712 // zeta_mapped = eta; 02713 // } 02714 // else 02715 // { 02716 // // Case 8 02717 // eta_mapped = eta; 02718 // zeta_mapped = -zeta; 02719 // } 02720 // } 02721 02722 02723 // // Face 3 02724 // else if ((i1[i] == 1) && (i0[i] >= 2) && (i2[i] >= 2)) 02725 // { 02726 // const unsigned int min_node = std::min(elem->node(2), 02727 // std::min(elem->node(3), 02728 // std::min(elem->node(7), 02729 // elem->node(6)))); 02730 // if (elem->node(3) == min_node) 02731 // if (elem->node(2) == std::min(elem->node(2), elem->node(7))) 02732 // { 02733 // // Case 1 02734 // xi_mapped = xi; 02735 // zeta_mapped = zeta; 02736 // } 02737 // else 02738 // { 02739 // // Case 2 02740 // xi_mapped = zeta; 02741 // zeta_mapped = xi; 02742 // } 02743 02744 // else if (elem->node(7) == min_node) 02745 // if (elem->node(3) == std::min(elem->node(3), elem->node(6))) 02746 // { 02747 // // Case 3 02748 // xi_mapped = -zeta; 02749 // zeta_mapped = xi; 02750 // } 02751 // else 02752 // { 02753 // // Case 4 02754 // xi_mapped = xi; 02755 // zeta_mapped = -zeta; 02756 // } 02757 02758 // else if (elem->node(6) == min_node) 02759 // if (elem->node(7) == std::min(elem->node(7), elem->node(2))) 02760 // { 02761 // // Case 5 02762 // xi_mapped = -xi; 02763 // zeta_mapped = -zeta; 02764 // } 02765 // else 02766 // { 02767 // // Case 6 02768 // xi_mapped = -zeta; 02769 // zeta_mapped = -xi; 02770 // } 02771 02772 // else if (elem->node(2) == min_node) 02773 // if (elem->node(6) == std::min(elem->node(3), elem->node(6))) 02774 // { 02775 // // Case 7 02776 // xi_mapped = zeta; 02777 // zeta_mapped = -xi; 02778 // } 02779 // else 02780 // { 02781 // // Case 8 02782 // xi_mapped = -xi; 02783 // zeta_mapped = zeta; 02784 // } 02785 // } 02786 02787 02788 // // Face 4 02789 // else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] >= 2)) 02790 // { 02791 // const unsigned int min_node = std::min(elem->node(3), 02792 // std::min(elem->node(0), 02793 // std::min(elem->node(4), 02794 // elem->node(7)))); 02795 // if (elem->node(0) == min_node) 02796 // if (elem->node(3) == std::min(elem->node(3), elem->node(4))) 02797 // { 02798 // // Case 1 02799 // eta_mapped = eta; 02800 // zeta_mapped = zeta; 02801 // } 02802 // else 02803 // { 02804 // // Case 2 02805 // eta_mapped = zeta; 02806 // zeta_mapped = eta; 02807 // } 02808 02809 // else if (elem->node(4) == min_node) 02810 // if (elem->node(0) == std::min(elem->node(0), elem->node(7))) 02811 // { 02812 // // Case 3 02813 // eta_mapped = -zeta; 02814 // zeta_mapped = eta; 02815 // } 02816 // else 02817 // { 02818 // // Case 4 02819 // eta_mapped = eta; 02820 // zeta_mapped = -zeta; 02821 // } 02822 02823 // else if (elem->node(7) == min_node) 02824 // if (elem->node(4) == std::min(elem->node(4), elem->node(3))) 02825 // { 02826 // // Case 5 02827 // eta_mapped = -eta; 02828 // zeta_mapped = -zeta; 02829 // } 02830 // else 02831 // { 02832 // // Case 6 02833 // eta_mapped = -zeta; 02834 // zeta_mapped = -eta; 02835 // } 02836 02837 // else if (elem->node(3) == min_node) 02838 // if (elem->node(7) == std::min(elem->node(7), elem->node(0))) 02839 // { 02840 // // Case 7 02841 // eta_mapped = zeta; 02842 // zeta_mapped = -eta; 02843 // } 02844 // else 02845 // { 02846 // // Case 8 02847 // eta_mapped = -eta; 02848 // zeta_mapped = zeta; 02849 // } 02850 // } 02851 02852 02853 // // Face 5 02854 // else if ((i2[i] == 1) && (i0[i] >= 2) && (i1[i] >= 2)) 02855 // { 02856 // const unsigned int min_node = std::min(elem->node(4), 02857 // std::min(elem->node(5), 02858 // std::min(elem->node(6), 02859 // elem->node(7)))); 02860 // if (elem->node(4) == min_node) 02861 // if (elem->node(5) == std::min(elem->node(5), elem->node(7))) 02862 // { 02863 // // Case 1 02864 // xi_mapped = xi; 02865 // eta_mapped = eta; 02866 // } 02867 // else 02868 // { 02869 // // Case 2 02870 // xi_mapped = eta; 02871 // eta_mapped = xi; 02872 // } 02873 02874 // else if (elem->node(5) == min_node) 02875 // if (elem->node(6) == std::min(elem->node(6), elem->node(4))) 02876 // { 02877 // // Case 3 02878 // xi_mapped = eta; 02879 // eta_mapped = -xi; 02880 // } 02881 // else 02882 // { 02883 // // Case 4 02884 // xi_mapped = -xi; 02885 // eta_mapped = eta; 02886 // } 02887 02888 // else if (elem->node(6) == min_node) 02889 // if (elem->node(7) == std::min(elem->node(7), elem->node(5))) 02890 // { 02891 // // Case 5 02892 // xi_mapped = -xi; 02893 // eta_mapped = -eta; 02894 // } 02895 // else 02896 // { 02897 // // Case 6 02898 // xi_mapped = -eta; 02899 // eta_mapped = -xi; 02900 // } 02901 02902 // else if (elem->node(7) == min_node) 02903 // if (elem->node(4) == std::min(elem->node(4), elem->node(6))) 02904 // { 02905 // // Case 7 02906 // xi_mapped = -eta; 02907 // eta_mapped = xi; 02908 // } 02909 // else 02910 // { 02911 // // Case 8 02912 // xi_mapped = xi; 02913 // eta_mapped = eta; 02914 // } 02915 // } 02916 02917 02918 // } 02919 02920 02921 02922 // libmesh_assert_less (j, 3); 02923 02924 // switch (j) 02925 // { 02926 // // d()/dxi 02927 // case 0: 02928 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi_mapped)* 02929 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)* 02930 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped)); 02931 02932 // // d()/deta 02933 // case 1: 02934 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)* 02935 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta_mapped)* 02936 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped)); 02937 02938 // // d()/dzeta 02939 // case 2: 02940 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)* 02941 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)* 02942 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta_mapped)); 02943 02944 // default: 02945 // libmesh_error_msg("Invalid derivative index j = " << j); 02946 // } 02947 // } 02948 02949 02950 default: 02951 libmesh_error_msg("Invalid element type = " << type); 02952 } 02953 } 02954 02955 02956 default: 02957 libmesh_error_msg("Invalid totalorder = " << totalorder); 02958 } 02959 02960 #endif 02961 02962 libmesh_error_msg("We'll never get here!"); 02963 return 0.; 02964 } 02965 02966 02967 02968 template <> 02969 Real FE<3,BERNSTEIN>::shape_second_deriv(const ElemType, 02970 const Order, 02971 const unsigned int, 02972 const unsigned int, 02973 const Point&) 02974 { 02975 static bool warning_given = false; 02976 02977 if (!warning_given) 02978 libMesh::err << "Second derivatives for Bernstein elements " 02979 << "are not yet implemented!" 02980 << std::endl; 02981 02982 warning_given = true; 02983 return 0.; 02984 } 02985 02986 02987 02988 template <> 02989 Real FE<3,BERNSTEIN>::shape_second_deriv(const Elem*, 02990 const Order, 02991 const unsigned int, 02992 const unsigned int, 02993 const Point&) 02994 { 02995 static bool warning_given = false; 02996 02997 if (!warning_given) 02998 libMesh::err << "Second derivatives for Bernstein elements " 02999 << "are not yet implemented!" 03000 << std::endl; 03001 03002 warning_given = true; 03003 return 0.; 03004 } 03005 03006 } // namespace libMesh 03007 03008 03009 03010 #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
generated by
