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