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