fe_hermite_shape_3D.C
Go to the documentation of this file.
00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 00003 00004 // This library is free software; you can redistribute it and/or 00005 // modify it under the terms of the GNU Lesser General Public 00006 // License as published by the Free Software Foundation; either 00007 // version 2.1 of the License, or (at your option) any later version. 00008 00009 // This library is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 // Lesser General Public License for more details. 00013 00014 // You should have received a copy of the GNU Lesser General Public 00015 // License along with this library; if not, write to the Free Software 00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00017 00018 00019 // C++ inlcludes 00020 00021 // Local includes 00022 #include "libmesh/fe.h" 00023 #include "libmesh/elem.h" 00024 #include "libmesh/number_lookups.h" 00025 00026 namespace 00027 { 00028 using namespace libMesh; 00029 00030 // Compute the static coefficients for an element 00031 void hermite_compute_coefs(const Elem* elem, std::vector<std::vector<Real> > & dxdxi 00032 00033 #ifdef DEBUG 00034 , std::vector<Real> & dydxi, std::vector<Real> & dzdeta, std::vector<Real> & dxdzeta, 00035 std::vector<Real> & dzdxi, std::vector<Real> & dxdeta, std::vector<Real> & dydzeta 00036 #endif 00037 ) 00038 { 00039 00040 const Order mapping_order (elem->default_order()); 00041 const ElemType mapping_elem_type (elem->type()); 00042 const int n_mapping_shape_functions = 00043 FE<3,LAGRANGE>::n_shape_functions(mapping_elem_type, 00044 mapping_order); 00045 00046 static const Point dofpt[2] = {Point(-1,-1,-1), Point(1,1,1)}; 00047 00048 for (int p = 0; p != 2; ++p) 00049 { 00050 dxdxi[0][p] = 0; 00051 dxdxi[1][p] = 0; 00052 dxdxi[2][p] = 0; 00053 #ifdef DEBUG 00054 dydxi[p] = 0; 00055 dzdeta[p] = 0; 00056 dxdzeta[p] = 0; 00057 dzdxi[p] = 0; 00058 dxdeta[p] = 0; 00059 dydzeta[p] = 0; 00060 #endif 00061 for (int i = 0; i != n_mapping_shape_functions; ++i) 00062 { 00063 const Real ddxi = FE<3,LAGRANGE>::shape_deriv 00064 (mapping_elem_type, mapping_order, i, 0, dofpt[p]); 00065 const Real ddeta = FE<3,LAGRANGE>::shape_deriv 00066 (mapping_elem_type, mapping_order, i, 1, dofpt[p]); 00067 const Real ddzeta = FE<3,LAGRANGE>::shape_deriv 00068 (mapping_elem_type, mapping_order, i, 2, dofpt[p]); 00069 00070 // dxdeta, dxdzeta, dydxi, dydzeta, dzdxi, dzdeta should all 00071 // be 0! 00072 const Point &point_i = elem->point(i); 00073 dxdxi[0][p] += point_i(0) * ddxi; 00074 dxdxi[1][p] += point_i(1) * ddeta; 00075 dxdxi[2][p] += point_i(2) * ddzeta; 00076 #ifdef DEBUG 00077 dydxi[p] += point_i(1) * ddxi; 00078 dzdeta[p] += point_i(2) * ddeta; 00079 dxdzeta[p] += point_i(0) * ddzeta; 00080 dzdxi[p] += point_i(2) * ddxi; 00081 dxdeta[p] += point_i(0) * ddeta; 00082 dydzeta[p] += point_i(1) * ddzeta; 00083 #endif 00084 } 00085 00086 // No singular elements! 00087 libmesh_assert(dxdxi[0][p]); 00088 libmesh_assert(dxdxi[1][p]); 00089 libmesh_assert(dxdxi[2][p]); 00090 // No non-rectilinear or non-axis-aligned elements! 00091 #ifdef DEBUG 00092 libmesh_assert_less (std::abs(dydxi[p]), TOLERANCE); 00093 libmesh_assert_less (std::abs(dzdeta[p]), TOLERANCE); 00094 libmesh_assert_less (std::abs(dxdzeta[p]), TOLERANCE); 00095 libmesh_assert_less (std::abs(dzdxi[p]), TOLERANCE); 00096 libmesh_assert_less (std::abs(dxdeta[p]), TOLERANCE); 00097 libmesh_assert_less (std::abs(dydzeta[p]), TOLERANCE); 00098 #endif 00099 } 00100 } 00101 00102 00103 00104 Real hermite_bases_3D 00105 (std::vector<unsigned int> &bases1D, 00106 const std::vector<std::vector<Real> > &dxdxi, 00107 const Order &o, 00108 unsigned int i) 00109 { 00110 bases1D.clear(); 00111 bases1D.resize(3,0); 00112 Real coef = 1.0; 00113 00114 unsigned int e = o-2; 00115 00116 // Nodes 00117 if (i < 64) 00118 { 00119 switch (i / 8) 00120 { 00121 case 0: 00122 break; 00123 case 1: 00124 bases1D[0] = 1; 00125 break; 00126 case 2: 00127 bases1D[0] = 1; 00128 bases1D[1] = 1; 00129 break; 00130 case 3: 00131 bases1D[1] = 1; 00132 break; 00133 case 4: 00134 bases1D[2] = 1; 00135 break; 00136 case 5: 00137 bases1D[0] = 1; 00138 bases1D[2] = 1; 00139 break; 00140 case 6: 00141 bases1D[0] = 1; 00142 bases1D[1] = 1; 00143 bases1D[2] = 1; 00144 break; 00145 case 7: 00146 bases1D[1] = 1; 00147 bases1D[2] = 1; 00148 break; 00149 default: 00150 libmesh_error_msg("Invalid basis node = " << i/8); 00151 } 00152 00153 unsigned int basisnum = i%8; 00154 switch (basisnum) // DoF type 00155 { 00156 case 0: // DoF = value at node 00157 coef = 1.0; 00158 break; 00159 case 1: // DoF = x derivative at node 00160 coef = dxdxi[0][bases1D[0]]; 00161 bases1D[0] += 2; break; 00162 case 2: // DoF = y derivative at node 00163 coef = dxdxi[1][bases1D[1]]; 00164 bases1D[1] += 2; break; 00165 case 3: // DoF = xy derivative at node 00166 coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]]; 00167 bases1D[0] += 2; bases1D[1] += 2; break; 00168 case 4: // DoF = z derivative at node 00169 coef = dxdxi[2][bases1D[2]]; 00170 bases1D[2] += 2; break; 00171 case 5: // DoF = xz derivative at node 00172 coef = dxdxi[0][bases1D[0]] * dxdxi[2][bases1D[2]]; 00173 bases1D[0] += 2; bases1D[2] += 2; break; 00174 case 6: // DoF = yz derivative at node 00175 coef = dxdxi[1][bases1D[1]] * dxdxi[2][bases1D[2]]; 00176 bases1D[1] += 2; bases1D[2] += 2; break; 00177 case 7: // DoF = xyz derivative at node 00178 coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]] * dxdxi[2][bases1D[2]]; 00179 bases1D[0] += 2; bases1D[1] += 2; bases1D[2] += 2; break; 00180 default: 00181 libmesh_error_msg("Invalid basisnum = " << basisnum); 00182 } 00183 } 00184 // Edges 00185 else if (i < 64 + 12*4*e) 00186 { 00187 unsigned int basisnum = (i - 64) % (4*e); 00188 switch ((i - 64) / (4*e)) 00189 { 00190 case 0: 00191 bases1D[0] = basisnum / 4 + 4; 00192 bases1D[1] = basisnum % 4 / 2 * 2; 00193 bases1D[2] = basisnum % 2 * 2; 00194 if (basisnum % 4 / 2) 00195 coef *= dxdxi[1][0]; 00196 if (basisnum % 2) 00197 coef *= dxdxi[2][0]; 00198 break; 00199 case 1: 00200 bases1D[0] = basisnum % 4 / 2 * 2 + 1; 00201 bases1D[1] = basisnum / 4 + 4; 00202 bases1D[2] = basisnum % 2 * 2; 00203 if (basisnum % 4 / 2) 00204 coef *= dxdxi[0][1]; 00205 if (basisnum % 2) 00206 coef *= dxdxi[2][0]; 00207 break; 00208 case 2: 00209 bases1D[0] = basisnum / 4 + 4; 00210 bases1D[1] = basisnum % 4 / 2 * 2 + 1; 00211 bases1D[2] = basisnum % 2 * 2; 00212 if (basisnum % 4 / 2) 00213 coef *= dxdxi[1][1]; 00214 if (basisnum % 2) 00215 coef *= dxdxi[2][0]; 00216 break; 00217 case 3: 00218 bases1D[0] = basisnum % 4 / 2 * 2; 00219 bases1D[1] = basisnum / 4 + 4; 00220 bases1D[2] = basisnum % 2 * 2; 00221 if (basisnum % 4 / 2) 00222 coef *= dxdxi[0][0]; 00223 if (basisnum % 2) 00224 coef *= dxdxi[2][0]; 00225 break; 00226 case 4: 00227 bases1D[0] = basisnum % 4 / 2 * 2; 00228 bases1D[1] = basisnum % 2 * 2; 00229 bases1D[2] = basisnum / 4 + 4; 00230 if (basisnum % 4 / 2) 00231 coef *= dxdxi[0][0]; 00232 if (basisnum % 2) 00233 coef *= dxdxi[1][0]; 00234 break; 00235 case 5: 00236 bases1D[0] = basisnum % 4 / 2 * 2 + 1; 00237 bases1D[1] = basisnum % 2 * 2; 00238 bases1D[2] = basisnum / 4 + 4; 00239 if (basisnum % 4 / 2) 00240 coef *= dxdxi[0][1]; 00241 if (basisnum % 2) 00242 coef *= dxdxi[1][0]; 00243 break; 00244 case 6: 00245 bases1D[0] = basisnum % 4 / 2 * 2 + 1; 00246 bases1D[1] = basisnum % 2 * 2 + 1; 00247 bases1D[2] = basisnum / 4 + 4; 00248 if (basisnum % 4 / 2) 00249 coef *= dxdxi[0][1]; 00250 if (basisnum % 2) 00251 coef *= dxdxi[1][1]; 00252 break; 00253 case 7: 00254 bases1D[0] = basisnum % 4 / 2 * 2; 00255 bases1D[1] = basisnum % 2 * 2 + 1; 00256 bases1D[2] = basisnum / 4 + 4; 00257 if (basisnum % 4 / 2) 00258 coef *= dxdxi[0][0]; 00259 if (basisnum % 2) 00260 coef *= dxdxi[1][1]; 00261 break; 00262 case 8: 00263 bases1D[0] = basisnum / 4 + 4; 00264 bases1D[1] = basisnum % 4 / 2 * 2; 00265 bases1D[2] = basisnum % 2 * 2 + 1; 00266 if (basisnum % 4 / 2) 00267 coef *= dxdxi[1][0]; 00268 if (basisnum % 2) 00269 coef *= dxdxi[2][1]; 00270 break; 00271 case 9: 00272 bases1D[0] = basisnum % 4 / 2 * 2 + 1; 00273 bases1D[1] = basisnum / 4 + 4; 00274 bases1D[2] = basisnum % 2 * 2; 00275 if (basisnum % 4 / 2) 00276 coef *= dxdxi[0][1]; 00277 if (basisnum % 2) 00278 coef *= dxdxi[2][1]; 00279 break; 00280 case 10: 00281 bases1D[0] = basisnum / 4 + 4; 00282 bases1D[1] = basisnum % 4 / 2 * 2 + 1; 00283 bases1D[2] = basisnum % 2 * 2 + 1; 00284 if (basisnum % 4 / 2) 00285 coef *= dxdxi[1][1]; 00286 if (basisnum % 2) 00287 coef *= dxdxi[2][1]; 00288 break; 00289 case 11: 00290 bases1D[0] = basisnum % 4 / 2 * 2; 00291 bases1D[1] = basisnum / 4 + 4; 00292 bases1D[2] = basisnum % 2 * 2 + 1; 00293 if (basisnum % 4 / 2) 00294 coef *= dxdxi[0][0]; 00295 if (basisnum % 2) 00296 coef *= dxdxi[2][1]; 00297 break; 00298 default: 00299 libmesh_error_msg("Invalid basis node = " << (i - 64) / (4*e)); 00300 } 00301 } 00302 // Faces 00303 else if (i < 64 + 12*4*e + 6*2*e*e) 00304 { 00305 unsigned int basisnum = (i - 64 - 12*4*e) % (2*e*e); 00306 switch ((i - 64 - 12*4*e) / (2*e*e)) 00307 { 00308 case 0: 00309 bases1D[0] = square_number_column[basisnum / 2]; 00310 bases1D[1] = square_number_row[basisnum / 2]; 00311 bases1D[2] = basisnum % 2 * 2; 00312 if (basisnum % 2) 00313 coef *= dxdxi[2][0]; 00314 break; 00315 case 1: 00316 bases1D[0] = square_number_column[basisnum / 2]; 00317 bases1D[1] = basisnum % 2 * 2; 00318 bases1D[2] = square_number_row[basisnum / 2]; 00319 if (basisnum % 2) 00320 coef *= dxdxi[1][0]; 00321 break; 00322 case 2: 00323 bases1D[0] = basisnum % 2 * 2 + 1; 00324 bases1D[1] = square_number_column[basisnum / 2]; 00325 bases1D[2] = square_number_row[basisnum / 2]; 00326 if (basisnum % 2) 00327 coef *= dxdxi[0][1]; 00328 break; 00329 case 3: 00330 bases1D[0] = square_number_column[basisnum / 2]; 00331 bases1D[1] = basisnum % 2 * 2 + 1; 00332 bases1D[2] = square_number_row[basisnum / 2]; 00333 if (basisnum % 2) 00334 coef *= dxdxi[1][1]; 00335 break; 00336 case 4: 00337 bases1D[0] = basisnum % 2 * 2; 00338 bases1D[1] = square_number_column[basisnum / 2]; 00339 bases1D[2] = square_number_row[basisnum / 2]; 00340 if (basisnum % 2) 00341 coef *= dxdxi[0][0]; 00342 break; 00343 case 5: 00344 bases1D[0] = square_number_column[basisnum / 2]; 00345 bases1D[1] = square_number_row[basisnum / 2]; 00346 bases1D[2] = basisnum % 2 * 2 + 1; 00347 if (basisnum % 2) 00348 coef *= dxdxi[2][1]; 00349 break; 00350 default: 00351 libmesh_error_msg("Invalid basis node = " << (i - 64 - 12*4*e) / (2*e*e)); 00352 } 00353 } 00354 // Interior 00355 else 00356 { 00357 unsigned int basisnum = i - 64 - 12*4*e; 00358 bases1D[0] = cube_number_column[basisnum] + 4; 00359 bases1D[1] = cube_number_row[basisnum] + 4; 00360 bases1D[2] = cube_number_page[basisnum] + 4; 00361 } 00362 00363 // No singular elements 00364 libmesh_assert(coef); 00365 return coef; 00366 } 00367 00368 00369 } // end anonymous namespace 00370 00371 00372 namespace libMesh 00373 { 00374 00375 00376 template <> 00377 Real FE<3,HERMITE>::shape(const ElemType, 00378 const Order, 00379 const unsigned int, 00380 const Point&) 00381 { 00382 libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom."); 00383 return 0.; 00384 } 00385 00386 00387 00388 template <> 00389 Real FE<3,HERMITE>::shape(const Elem* elem, 00390 const Order order, 00391 const unsigned int i, 00392 const Point& p) 00393 { 00394 libmesh_assert(elem); 00395 00396 std::vector<std::vector<Real> > dxdxi(3, std::vector<Real>(2, 0)); 00397 00398 #ifdef DEBUG 00399 std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2); 00400 std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2); 00401 #endif //DEBUG 00402 00403 hermite_compute_coefs(elem, dxdxi 00404 #ifdef DEBUG 00405 , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta 00406 #endif 00407 ); 00408 00409 const ElemType type = elem->type(); 00410 00411 const Order totalorder = static_cast<Order>(order + elem->p_level()); 00412 00413 switch (totalorder) 00414 { 00415 // 3rd-order tricubic Hermite functions 00416 case THIRD: 00417 { 00418 switch (type) 00419 { 00420 case HEX8: 00421 case HEX20: 00422 case HEX27: 00423 { 00424 libmesh_assert_less (i, 64); 00425 00426 std::vector<unsigned int> bases1D; 00427 00428 Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i); 00429 00430 return coef * 00431 FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) * 00432 FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) * 00433 FEHermite<1>::hermite_raw_shape(bases1D[2],p(2)); 00434 } 00435 default: 00436 libmesh_error_msg("ERROR: Unsupported element type " << type); 00437 } 00438 } 00439 // by default throw an error 00440 default: 00441 libmesh_error_msg("ERROR: Unsupported polynomial order " << totalorder); 00442 } 00443 00444 libmesh_error_msg("We'll never get here!"); 00445 return 0.; 00446 } 00447 00448 00449 00450 template <> 00451 Real FE<3,HERMITE>::shape_deriv(const ElemType, 00452 const Order, 00453 const unsigned int, 00454 const unsigned int, 00455 const Point&) 00456 { 00457 libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom."); 00458 return 0.; 00459 } 00460 00461 00462 00463 template <> 00464 Real FE<3,HERMITE>::shape_deriv(const Elem* elem, 00465 const Order order, 00466 const unsigned int i, 00467 const unsigned int j, 00468 const Point& p) 00469 { 00470 libmesh_assert(elem); 00471 libmesh_assert (j == 0 || j == 1 || j == 2); 00472 00473 std::vector<std::vector<Real> > dxdxi(3, std::vector<Real>(2, 0)); 00474 00475 #ifdef DEBUG 00476 std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2); 00477 std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2); 00478 #endif //DEBUG 00479 00480 hermite_compute_coefs(elem, dxdxi 00481 #ifdef DEBUG 00482 , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta 00483 #endif 00484 ); 00485 00486 const ElemType type = elem->type(); 00487 00488 const Order totalorder = static_cast<Order>(order + elem->p_level()); 00489 00490 switch (totalorder) 00491 { 00492 // 3rd-order tricubic Hermite functions 00493 case THIRD: 00494 { 00495 switch (type) 00496 { 00497 case HEX8: 00498 case HEX20: 00499 case HEX27: 00500 { 00501 libmesh_assert_less (i, 64); 00502 00503 std::vector<unsigned int> bases1D; 00504 00505 Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i); 00506 00507 switch (j) // Derivative type 00508 { 00509 case 0: 00510 return coef * 00511 FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) * 00512 FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) * 00513 FEHermite<1>::hermite_raw_shape(bases1D[2],p(2)); 00514 break; 00515 case 1: 00516 return coef * 00517 FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) * 00518 FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) * 00519 FEHermite<1>::hermite_raw_shape(bases1D[2],p(2)); 00520 break; 00521 case 2: 00522 return coef * 00523 FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) * 00524 FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) * 00525 FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2)); 00526 break; 00527 default: 00528 libmesh_error_msg("Invalid shape function derivative j = " << j); 00529 } 00530 00531 } 00532 default: 00533 libmesh_error_msg("ERROR: Unsupported element type " << type); 00534 } 00535 } 00536 // by default throw an error 00537 default: 00538 libmesh_error_msg("ERROR: Unsupported polynomial order " << totalorder); 00539 } 00540 00541 libmesh_error_msg("We'll never get here!"); 00542 return 0.; 00543 } 00544 00545 00546 00547 template <> 00548 Real FE<3,HERMITE>::shape_second_deriv(const Elem* elem, 00549 const Order order, 00550 const unsigned int i, 00551 const unsigned int j, 00552 const Point& p) 00553 { 00554 libmesh_assert(elem); 00555 00556 std::vector<std::vector<Real> > dxdxi(3, std::vector<Real>(2, 0)); 00557 00558 #ifdef DEBUG 00559 std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2); 00560 std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2); 00561 #endif //DEBUG 00562 00563 hermite_compute_coefs(elem, dxdxi 00564 #ifdef DEBUG 00565 , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta 00566 #endif 00567 ); 00568 00569 const ElemType type = elem->type(); 00570 00571 const Order totalorder = static_cast<Order>(order + elem->p_level()); 00572 00573 switch (totalorder) 00574 { 00575 // 3rd-order tricubic Hermite functions 00576 case THIRD: 00577 { 00578 switch (type) 00579 { 00580 case HEX8: 00581 case HEX20: 00582 case HEX27: 00583 { 00584 libmesh_assert_less (i, 64); 00585 00586 std::vector<unsigned int> bases1D; 00587 00588 Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i); 00589 00590 switch (j) // Derivative type 00591 { 00592 case 0: 00593 return coef * 00594 FEHermite<1>::hermite_raw_shape_second_deriv(bases1D[0],p(0)) * 00595 FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) * 00596 FEHermite<1>::hermite_raw_shape(bases1D[2],p(2)); 00597 break; 00598 case 1: 00599 return coef * 00600 FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) * 00601 FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) * 00602 FEHermite<1>::hermite_raw_shape(bases1D[2],p(2)); 00603 break; 00604 case 2: 00605 return coef * 00606 FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) * 00607 FEHermite<1>::hermite_raw_shape_second_deriv(bases1D[1],p(1)) * 00608 FEHermite<1>::hermite_raw_shape(bases1D[2],p(2)); 00609 break; 00610 case 3: 00611 return coef * 00612 FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) * 00613 FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) * 00614 FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2)); 00615 break; 00616 case 4: 00617 return coef * 00618 FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) * 00619 FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) * 00620 FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2)); 00621 break; 00622 case 5: 00623 return coef * 00624 FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) * 00625 FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) * 00626 FEHermite<1>::hermite_raw_shape_second_deriv(bases1D[2],p(2)); 00627 break; 00628 default: 00629 libmesh_error_msg("Invalid shape function derivative j = " << j); 00630 } 00631 00632 } 00633 default: 00634 libmesh_error_msg("ERROR: Unsupported element type " << type); 00635 } 00636 } 00637 // by default throw an error 00638 default: 00639 libmesh_error_msg("ERROR: Unsupported polynomial order " << totalorder); 00640 } 00641 00642 libmesh_error_msg("We'll never get here!"); 00643 return 0.; 00644 } 00645 00646 } // namespace libMesh
generated by
