fe_bernstein_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 // Local includes 00020 #include "libmesh/libmesh_config.h" 00021 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES 00022 00023 #include "libmesh/fe.h" 00024 #include "libmesh/elem.h" 00025 #include "libmesh/number_lookups.h" 00026 #include "libmesh/utility.h" 00027 00028 00029 namespace libMesh 00030 { 00031 00032 00033 template <> 00034 Real FE<2,BERNSTEIN>::shape(const ElemType, 00035 const Order, 00036 const unsigned int, 00037 const Point&) 00038 { 00039 libMesh::err << "Bernstein polynomials require the element type\n" 00040 << "because edge orientation is needed." 00041 << std::endl; 00042 00043 libmesh_error(); 00044 return 0.; 00045 } 00046 00047 00048 00049 template <> 00050 Real FE<2,BERNSTEIN>::shape(const Elem* elem, 00051 const Order order, 00052 const unsigned int i, 00053 const Point& p) 00054 { 00055 libmesh_assert(elem); 00056 00057 const ElemType type = elem->type(); 00058 00059 const Order totalorder = static_cast<Order>(order + elem->p_level()); 00060 00061 // Declare that we are using our own special power function 00062 // from the Utility namespace. This saves typing later. 00063 using Utility::pow; 00064 00065 switch (type) 00066 { 00067 // Hierarchic shape functions on the quadrilateral. 00068 case QUAD4: 00069 case QUAD9: 00070 { 00071 // Compute quad shape functions as a tensor-product 00072 const Real xi = p(0); 00073 const Real eta = p(1); 00074 00075 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)); 00076 00077 // Example i, i0, i1 values for totalorder = 5: 00078 // 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 00079 // 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, 3, 2, 4, 4, 4, 3, 2, 5, 5, 5, 5, 4, 3, 2}; 00080 // 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, 3, 3, 2, 3, 4, 4, 4, 2, 3, 4, 5, 5, 5, 5}; 00081 00082 unsigned int i0, i1; 00083 00084 // Vertex DoFs 00085 if (i == 0) 00086 { i0 = 0; i1 = 0; } 00087 else if (i == 1) 00088 { i0 = 1; i1 = 0; } 00089 else if (i == 2) 00090 { i0 = 1; i1 = 1; } 00091 else if (i == 3) 00092 { i0 = 0; i1 = 1; } 00093 00094 00095 // Edge DoFs 00096 else if (i < totalorder + 3u) 00097 { i0 = i - 2; i1 = 0; } 00098 else if (i < 2u*totalorder + 2) 00099 { i0 = 1; i1 = i - totalorder - 1; } 00100 else if (i < 3u*totalorder + 1) 00101 { i0 = i - 2u*totalorder; i1 = 1; } 00102 else if (i < 4u*totalorder) 00103 { i0 = 0; i1 = i - 3u*totalorder + 1; } 00104 // Interior DoFs. Use Roy's number look up 00105 else 00106 { 00107 unsigned int basisnum = i - 4*totalorder; 00108 i0 = square_number_column[basisnum] + 2; 00109 i1 = square_number_row[basisnum] + 2; 00110 } 00111 00112 00113 // Flip odd degree of freedom values if necessary 00114 // to keep continuity on sides. 00115 if ((i>= 4 && i<= 4+ totalorder-2u) && elem->point(0) > elem->point(1)) i0=totalorder+2-i0; // 00116 else if((i>= 4+ totalorder-1u && i<= 4+2*totalorder-3u) && elem->point(1) > elem->point(2)) i1=totalorder+2-i1; 00117 else if((i>= 4+2*totalorder-2u && i<= 4+3*totalorder-4u) && elem->point(3) > elem->point(2)) i0=totalorder+2-i0; 00118 else if((i>= 4+3*totalorder-3u && i<= 4+4*totalorder-5u) && elem->point(0) > elem->point(3)) i1=totalorder+2-i1; 00119 00120 00121 return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0, xi)* 00122 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1, eta)); 00123 } 00124 // handle serendipity QUAD8 element separately 00125 case QUAD8: 00126 { 00127 libmesh_assert_less (totalorder, 3); 00128 00129 const Real xi = p(0); 00130 const Real eta = p(1); 00131 00132 libmesh_assert_less (i, 8); 00133 00134 // 0 1 2 3 4 5 6 7 8 00135 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2}; 00136 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2}; 00137 static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5}; 00138 00139 //B_t,i0(i)|xi * B_s,i1(i)|eta + scal(i) * B_t,2|xi * B_t,2|eta 00140 return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)* 00141 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta) 00142 +scal[i]* 00143 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[8], xi)* 00144 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[8], eta)); 00145 00146 } 00147 00148 case TRI3: 00149 libmesh_assert_less (totalorder, 2); 00150 case TRI6: 00151 switch (totalorder) 00152 { 00153 case FIRST: 00154 { 00155 const Real x=p(0); 00156 const Real y=p(1); 00157 const Real r=1.-x-y; 00158 00159 libmesh_assert_less (i, 3); 00160 00161 switch(i) 00162 { 00163 case 0: return r; //f0,0,1 00164 case 1: return x; //f0,1,1 00165 case 2: return y; //f1,0,1 00166 00167 default: libmesh_error(); return 0; 00168 } 00169 } 00170 case SECOND: 00171 { 00172 const Real x=p(0); 00173 const Real y=p(1); 00174 const Real r=1.-x-y; 00175 00176 libmesh_assert_less (i, 6); 00177 00178 switch(i) 00179 { 00180 case 0: return r*r; 00181 case 1: return x*x; 00182 case 2: return y*y; 00183 00184 case 3: return 2.*x*r; 00185 case 4: return 2.*x*y; 00186 case 5: return 2.*r*y; 00187 00188 default: libmesh_error(); return 0; 00189 } 00190 } 00191 case THIRD: 00192 { 00193 const Real x=p(0); 00194 const Real y=p(1); 00195 const Real r=1.-x-y; 00196 libmesh_assert_less (i, 10); 00197 00198 unsigned int shape=i; 00199 00200 00201 if((i==3||i==4) && elem->point(0) > elem->point(1)) shape=7-i; 00202 if((i==5||i==6) && elem->point(1) > elem->point(2)) shape=11-i; 00203 if((i==7||i==8) && elem->point(0) > elem->point(2)) shape=15-i; 00204 00205 switch(shape) 00206 { 00207 case 0: return r*r*r; 00208 case 1: return x*x*x; 00209 case 2: return y*y*y; 00210 00211 case 3: return 3.*x*r*r; 00212 case 4: return 3.*x*x*r; 00213 00214 case 5: return 3.*y*x*x; 00215 case 6: return 3.*y*y*x; 00216 00217 case 7: return 3.*y*r*r; 00218 case 8: return 3.*y*y*r; 00219 00220 case 9: return 6.*x*y*r; 00221 00222 default: libmesh_error(); return 0; 00223 } 00224 } 00225 case FOURTH: 00226 { 00227 const Real x=p(0); 00228 const Real y=p(1); 00229 const Real r=1-x-y; 00230 unsigned int shape=i; 00231 00232 libmesh_assert_less (i, 15); 00233 00234 if((i==3||i== 5) && elem->point(0) > elem->point(1))shape=8-i; 00235 if((i==6||i== 8) && elem->point(1) > elem->point(2))shape=14-i; 00236 if((i==9||i==11) && elem->point(0) > elem->point(2))shape=20-i; 00237 00238 00239 switch(shape) 00240 { 00241 // point functions 00242 case 0: return r*r*r*r; 00243 case 1: return x*x*x*x; 00244 case 2: return y*y*y*y; 00245 00246 // edge functions 00247 case 3: return 4.*x*r*r*r; 00248 case 4: return 6.*x*x*r*r; 00249 case 5: return 4.*x*x*x*r; 00250 00251 case 6: return 4.*y*x*x*x; 00252 case 7: return 6.*y*y*x*x; 00253 case 8: return 4.*y*y*y*x; 00254 00255 case 9: return 4.*y*r*r*r; 00256 case 10: return 6.*y*y*r*r; 00257 case 11: return 4.*y*y*y*r; 00258 00259 // inner functions 00260 case 12: return 12.*x*y*r*r; 00261 case 13: return 12.*x*x*y*r; 00262 case 14: return 12.*x*y*y*r; 00263 00264 default: libmesh_error(); return 0; 00265 } 00266 } 00267 case FIFTH: 00268 { 00269 const Real x=p(0); 00270 const Real y=p(1); 00271 const Real r=1-x-y; 00272 unsigned int shape=i; 00273 00274 libmesh_assert_less (i, 21); 00275 00276 if((i>= 3&&i<= 6) && elem->point(0) > elem->point(1))shape=9-i; 00277 if((i>= 7&&i<=10) && elem->point(1) > elem->point(2))shape=17-i; 00278 if((i>=11&&i<=14) && elem->point(0) > elem->point(2))shape=25-i; 00279 00280 switch(shape) 00281 { 00282 //point functions 00283 case 0: return pow<5>(r); 00284 case 1: return pow<5>(x); 00285 case 2: return pow<5>(y); 00286 00287 //edge functions 00288 case 3: return 5.*x *pow<4>(r); 00289 case 4: return 10.*pow<2>(x)*pow<3>(r); 00290 case 5: return 10.*pow<3>(x)*pow<2>(r); 00291 case 6: return 5.*pow<4>(x)*r; 00292 00293 case 7: return 5.*y *pow<4>(x); 00294 case 8: return 10.*pow<2>(y)*pow<3>(x); 00295 case 9: return 10.*pow<3>(y)*pow<2>(x); 00296 case 10: return 5.*pow<4>(y)*x; 00297 00298 case 11: return 5.*y *pow<4>(r); 00299 case 12: return 10.*pow<2>(y)*pow<3>(r); 00300 case 13: return 10.*pow<3>(y)*pow<2>(r); 00301 case 14: return 5.*pow<4>(y)*r; 00302 00303 //inner functions 00304 case 15: return 20.*x*y*pow<3>(r); 00305 case 16: return 30.*x*pow<2>(y)*pow<2>(r); 00306 case 17: return 30.*pow<2>(x)*y*pow<2>(r); 00307 case 18: return 20.*x*pow<3>(y)*r; 00308 case 19: return 20.*pow<3>(x)*y*r; 00309 case 20: return 30.*pow<2>(x)*pow<2>(y)*r; 00310 00311 default: libmesh_error(); return 0; 00312 } 00313 } 00314 case SIXTH: 00315 { 00316 const Real x=p(0); 00317 const Real y=p(1); 00318 const Real r=1-x-y; 00319 unsigned int shape=i; 00320 00321 libmesh_assert_less (i, 28); 00322 00323 if((i>= 3&&i<= 7) && elem->point(0) > elem->point(1))shape=10-i; 00324 if((i>= 8&&i<=12) && elem->point(1) > elem->point(2))shape=20-i; 00325 if((i>=13&&i<=17) && elem->point(0) > elem->point(2))shape=30-i; 00326 00327 switch(shape) 00328 { 00329 //point functions 00330 case 0: return pow<6>(r); 00331 case 1: return pow<6>(x); 00332 case 2: return pow<6>(y); 00333 00334 //edge functions 00335 case 3: return 6.*x *pow<5>(r); 00336 case 4: return 15.*pow<2>(x)*pow<4>(r); 00337 case 5: return 20.*pow<3>(x)*pow<3>(r); 00338 case 6: return 15.*pow<4>(x)*pow<2>(r); 00339 case 7: return 6.*pow<5>(x)*r; 00340 00341 case 8: return 6.*y *pow<5>(x); 00342 case 9: return 15.*pow<2>(y)*pow<4>(x); 00343 case 10: return 20.*pow<3>(y)*pow<3>(x); 00344 case 11: return 15.*pow<4>(y)*pow<2>(x); 00345 case 12: return 6.*pow<5>(y)*x; 00346 00347 case 13: return 6.*y *pow<5>(r); 00348 case 14: return 15.*pow<2>(y)*pow<4>(r); 00349 case 15: return 20.*pow<3>(y)*pow<3>(r); 00350 case 16: return 15.*pow<4>(y)*pow<2>(r); 00351 case 17: return 6.*pow<5>(y)*r; 00352 00353 //inner functions 00354 case 18: return 30.*x*y*pow<4>(r); 00355 case 19: return 60.*x*pow<2>(y)*pow<3>(r); 00356 case 20: return 60.* pow<2>(x)*y*pow<3>(r); 00357 case 21: return 60.*x*pow<3>(y)*pow<2>(r); 00358 case 22: return 60.*pow<3>(x)*y*pow<2>(r); 00359 case 23: return 90.*pow<2>(x)*pow<2>(y)*pow<2>(r); 00360 case 24: return 30.*x*pow<4>(y)*r; 00361 case 25: return 60.*pow<2>(x)*pow<3>(y)*r; 00362 case 26: return 60.*pow<3>(x)*pow<2>(y)*r; 00363 case 27: return 30.*pow<4>(x)*y*r; 00364 00365 default: libmesh_error(); return 0; 00366 } // switch shape 00367 } // case TRI6 00368 default: 00369 { 00370 libMesh::err << "ERROR: element order!" << std::endl; 00371 libmesh_error(); 00372 } 00373 00374 00375 } // switch order 00376 00377 00378 default: 00379 { 00380 libMesh::err << "ERROR: Unsupported element type!" << std::endl; 00381 libmesh_error(); 00382 } 00383 00384 } // switch type 00385 00386 00387 // old code 00388 // switch (totalorder) 00389 // { 00390 00391 // case FIRST: 00392 00393 // switch(type) 00394 // { 00395 00396 // case TRI6: 00397 // { 00398 // const Real x=p(0); 00399 // const Real y=p(1); 00400 // const Real r=1.-x-y; 00401 00402 // libmesh_assert_less (i, 3); 00403 00404 // switch(i) 00405 // { 00406 // case 0: return r; //f0,0,1 00407 // case 1: return x; //f0,1,1 00408 // case 2: return y; //f1,0,1 00409 // } 00410 // } 00411 00412 // case QUAD8: 00413 // case QUAD9: 00414 // { 00415 // // Compute quad shape functions as a tensor-product 00416 // const Real xi = p(0); 00417 // const Real eta = p(1); 00418 00419 // libmesh_assert_less (i, 4); 00420 00421 // // 0 1 2 3 00422 // static const unsigned int i0[] = {0, 1, 1, 0}; 00423 // static const unsigned int i1[] = {0, 0, 1, 1}; 00424 00425 // return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)* 00426 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)); 00427 00428 // } 00429 // default: 00430 // libmesh_error(); 00431 // } 00432 00433 // case SECOND: 00434 // switch(type) 00435 // { 00436 // case TRI6: 00437 // { 00438 // const Real x=p(0); 00439 // const Real y=p(1); 00440 // const Real r=1.-x-y; 00441 00442 // libmesh_assert_less (i, 6); 00443 00444 // switch(i) 00445 // { 00446 // case 0: return r*r; 00447 // case 1: return x*x; 00448 // case 2: return y*y; 00449 00450 // case 3: return 2.*x*r; 00451 // case 4: return 2.*x*y; 00452 // case 5: return 2.*r*y; 00453 // } 00454 // } 00455 00456 // // Bernstein shape functions on the 8-noded quadrilateral. 00457 // case QUAD8: 00458 // { 00459 // const Real xi = p(0); 00460 // const Real eta = p(1); 00461 00462 // libmesh_assert_less (i, 8); 00463 00464 // // 0 1 2 3 4 5 6 7 8 00465 // static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2}; 00466 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2}; 00467 // static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5}; 00468 00469 // return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)* 00470 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta) 00471 // +scal[i]* 00472 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[8], xi)* 00473 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[8], eta)); 00474 00475 // //B_t,i0(i)|xi * B_s,i1(i)|eta + scal(i) * B_t,2|xi * B_t,2|eta 00476 // } 00477 00478 // // Bernstein shape functions on the 9-noded quadrilateral. 00479 // case QUAD9: 00480 // { 00481 // // Compute quad shape functions as a tensor-product 00482 // const Real xi = p(0); 00483 // const Real eta = p(1); 00484 00485 // libmesh_assert_less (i, 9); 00486 00487 // // 0 1 2 3 4 5 6 7 8 00488 // static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2}; 00489 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2}; 00490 00491 // return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)* 00492 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)); 00493 00494 // } 00495 // default: 00496 // libmesh_error(); 00497 // } 00498 00499 // case THIRD: 00500 // switch(type) 00501 // { 00502 // case TRI6: 00503 // { 00504 // const Real x=p(0); 00505 // const Real y=p(1); 00506 // const Real r=1.-x-y; 00507 // libmesh_assert_less (i, 10); 00508 00509 // unsigned int shape=i; 00510 00511 00512 // if((i==3||i==4) && elem->node(0) > elem->node(1))shape=7-i; 00513 // if((i==5||i==6) && elem->node(1) > elem->node(2))shape=11-i; 00514 // if((i==7||i==8) && elem->node(0) > elem->node(2))shape=15-i; 00515 00516 00517 // switch(shape) 00518 // { 00519 // case 0: return r*r*r; 00520 // case 1: return x*x*x; 00521 // case 2: return y*y*y; 00522 00523 // case 3: return 3.*x*r*r; 00524 // case 4: return 3.*x*x*r; 00525 00526 // case 5: return 3.*y*x*x; 00527 // case 6: return 3.*y*y*x; 00528 00529 // case 7: return 3.*y*r*r; 00530 // case 8: return 3.*y*y*r; 00531 00532 // case 9: return 6.*x*y*r; 00533 // } 00534 // } 00535 00536 // // Bernstein shape functions on the quadrilateral. 00537 // case QUAD8: 00538 // case QUAD9: 00539 // { 00540 // // Compute quad shape functions as a tensor-product 00541 // Real xi = p(0); 00542 // Real eta = p(1); 00543 00544 // libmesh_assert_less (i, 16); 00545 00546 // // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 00547 // static const unsigned int i0_reg[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3}; 00548 // static const unsigned int i1_reg[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3}; 00549 00550 // unsigned int i0=i0_reg[i]; 00551 // unsigned int i1=i1_reg[i]; 00552 00553 // if((i== 4||i== 5) && elem->node(0) > elem->node(1)) i0=5-i0; // 2->3 3->2 00554 // if((i== 6||i== 7) && elem->node(1) > elem->node(2)) i1=5-i1; 00555 // if((i== 8||i== 9) && elem->node(3) > elem->node(2)) i0=5-i0; 00556 // if((i==10||i==11) && elem->node(0) > elem->node(3)) i1=5-i1; 00557 00558 // // element dof orientation is needed when used with ifems 00559 // // if(i > 11) 00560 // // { 00561 // // const unsigned int min_node = std::min(elem->node(1), 00562 // // std::min(elem->node(2), 00563 // // std::min(elem->node(0), 00564 // // elem->node(3)))); 00565 // // if (elem->node(0) == min_node) 00566 // // if (elem->node(1) == std::min(elem->node(1), elem->node(3))) 00567 // // { 00568 // // // Case 1 00569 // // xi = xi; 00570 // // eta = eta; 00571 // // } 00572 // // else 00573 // // { 00574 // // // Case 2 00575 // // xi = eta; 00576 // // eta = xi; 00577 // // } 00578 00579 // // else if (elem->node(3) == min_node) 00580 // // if (elem->node(0) == std::min(elem->node(0), elem->node(2))) 00581 // // { 00582 // // // Case 3 00583 // // xi = -eta; 00584 // // eta = xi; 00585 // // } 00586 // // else 00587 // // { 00588 // // // Case 4 00589 // // xi = xi; 00590 // // eta = -eta; 00591 // // } 00592 00593 // // else if (elem->node(2) == min_node) 00594 // // if (elem->node(3) == std::min(elem->node(3), elem->node(1))) 00595 // // { 00596 // // // Case 5 00597 // // xi = -xi; 00598 // // eta = -eta; 00599 // // } 00600 // // else 00601 // // { 00602 // // // Case 6 00603 // // xi = -eta; 00604 // // eta = -xi; 00605 // // } 00606 00607 // // else if (elem->node(1) == min_node) 00608 // // if (elem->node(2) == std::min(elem->node(2), elem->node(0))) 00609 // // { 00610 // // // Case 7 00611 // // xi = eta; 00612 // // eta = -xi; 00613 // // } 00614 // // else 00615 // // { 00616 // // // Case 8 00617 // // xi = -xi; 00618 // // eta = eta; 00619 // // } 00620 // // } 00621 00622 00623 // return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0, xi)* 00624 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1, eta)); 00625 00626 // } 00627 00628 // default: 00629 // libmesh_error(); 00630 00631 // } 00632 00633 // case FOURTH: 00634 // switch(type) 00635 // { 00636 // case TRI6: 00637 // { 00638 // const Real x=p(0); 00639 // const Real y=p(1); 00640 // const Real r=1-x-y; 00641 // unsigned int shape=i; 00642 00643 // libmesh_assert_less (i, 15); 00644 00645 // if((i==3||i== 5) && elem->node(0) > elem->node(1))shape=8-i; 00646 // if((i==6||i== 8) && elem->node(1) > elem->node(2))shape=14-i; 00647 // if((i==9||i==11) && elem->node(0) > elem->node(2))shape=20-i; 00648 00649 00650 // switch(shape) 00651 // { 00652 // // point functions 00653 // case 0: return r*r*r*r; 00654 // case 1: return x*x*x*x; 00655 // case 2: return y*y*y*y; 00656 00657 // // edge functions 00658 // case 3: return 4.*x*r*r*r; 00659 // case 4: return 6.*x*x*r*r; 00660 // case 5: return 4.*x*x*x*r; 00661 00662 // case 6: return 4.*y*x*x*x; 00663 // case 7: return 6.*y*y*x*x; 00664 // case 8: return 4.*y*y*y*x; 00665 00666 // case 9: return 4.*y*r*r*r; 00667 // case 10: return 6.*y*y*r*r; 00668 // case 11: return 4.*y*y*y*r; 00669 00670 // // inner functions 00671 // case 12: return 12.*x*y*r*r; 00672 // case 13: return 12.*x*x*y*r; 00673 // case 14: return 12.*x*y*y*r; 00674 00675 // } 00676 // } 00677 00678 00679 // // Bernstein shape functions on the quadrilateral. 00680 // case QUAD8: 00681 // case QUAD9: 00682 // { 00683 // // Compute quad shape functions as a tensor-product 00684 // const Real xi = p(0); 00685 // const Real eta = p(1); 00686 00687 // libmesh_assert_less (i, 25); 00688 00689 // // 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 00690 // static const unsigned int i0_reg[] = {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}; 00691 // static const unsigned int i1_reg[] = {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}; 00692 00693 // unsigned int i0=i0_reg[i]; 00694 // unsigned int i1=i1_reg[i]; 00695 00696 // if((i>= 4&&i<= 6) && elem->node(0) > elem->node(1)) i0=6-i0; // 2->4, 4->2 00697 // if((i>= 7&&i<= 9) && elem->node(1) > elem->node(2)) i1=6-i1; 00698 // if((i>=10&&i<=12) && elem->node(3) > elem->node(2)) i0=6-i0; 00699 // if((i>=13&&i<=15) && elem->node(0) > elem->node(3)) i1=6-i1; 00700 00701 // return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0, xi)* 00702 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1, eta)); 00703 00704 // } 00705 00706 // default: 00707 // libmesh_error(); 00708 00709 // } 00710 00711 // case FIFTH: 00712 // switch(type) 00713 // { 00714 // case TRI6: 00715 // { 00716 // const Real x=p(0); 00717 // const Real y=p(1); 00718 // const Real r=1-x-y; 00719 // unsigned int shape=i; 00720 00721 // libmesh_assert_less (i, 21); 00722 00723 // if((i>= 3&&i<= 6) && elem->node(0) > elem->node(1))shape=9-i; 00724 // if((i>= 7&&i<=10) && elem->node(1) > elem->node(2))shape=17-i; 00725 // if((i>=11&&i<=14) && elem->node(0) > elem->node(2))shape=25-i; 00726 00727 // switch(shape) 00728 // { 00729 // //point functions 00730 // case 0: return pow<5>(r); 00731 // case 1: return pow<5>(x); 00732 // case 2: return pow<5>(y); 00733 00734 // //edge functions 00735 // case 3: return 5.*x *pow<4>(r); 00736 // case 4: return 10.*pow<2>(x)*pow<3>(r); 00737 // case 5: return 10.*pow<3>(x)*pow<2>(r); 00738 // case 6: return 5.*pow<4>(x)*r; 00739 00740 // case 7: return 5.*y *pow<4>(x); 00741 // case 8: return 10.*pow<2>(y)*pow<3>(x); 00742 // case 9: return 10.*pow<3>(y)*pow<2>(x); 00743 // case 10: return 5.*pow<4>(y)*x; 00744 00745 // case 11: return 5.*y *pow<4>(r); 00746 // case 12: return 10.*pow<2>(y)*pow<3>(r); 00747 // case 13: return 10.*pow<3>(y)*pow<2>(r); 00748 // case 14: return 5.*pow<4>(y)*r; 00749 00750 // //inner functions 00751 // case 15: return 20.*x*y*pow<3>(r); 00752 // case 16: return 30.*x*pow<2>(y)*pow<2>(r); 00753 // case 17: return 30.*pow<2>(x)*y*pow<2>(r); 00754 // case 18: return 20.*x*pow<3>(y)*r; 00755 // case 19: return 20.*pow<3>(x)*y*r; 00756 // case 20: return 30.*pow<2>(x)*pow<2>(y)*r; 00757 00758 // } 00759 // } 00760 00761 00762 // // Bernstein shape functions on the quadrilateral. 00763 // case QUAD8: 00764 // case QUAD9: 00765 // { 00766 // // Compute quad shape functions as a tensor-product 00767 // const Real xi = p(0); 00768 // const Real eta = p(1); 00769 00770 // libmesh_assert_less (i, 36); 00771 00772 // // 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 00773 // static const unsigned int i0_reg[] = {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}; 00774 // static const unsigned int i1_reg[] = {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}; 00775 00776 // unsigned int i0=i0_reg[i]; 00777 // unsigned int i1=i1_reg[i]; 00778 00779 // if((i>= 4&&i<= 7) && elem->node(0) > elem->node(1)) i0=7-i0; 00780 // if((i>= 8&&i<=11) && elem->node(1) > elem->node(2)) i1=7-i1; 00781 // if((i>=12&&i<=15) && elem->node(3) > elem->node(2)) i0=7-i0; 00782 // if((i>=16&&i<=19) && elem->node(0) > elem->node(3)) i1=7-i1; 00783 00784 // return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0, xi)* 00785 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1, eta)); 00786 00787 // } 00788 // default: 00789 // libmesh_error(); 00790 // } 00791 00792 // case SIXTH: 00793 // switch(type) 00794 // { 00795 // case TRI6: 00796 // { 00797 // const Real x=p(0); 00798 // const Real y=p(1); 00799 // const Real r=1-x-y; 00800 // unsigned int shape=i; 00801 00802 // libmesh_assert_less (i, 28); 00803 00804 // if((i>= 3&&i<= 7) && elem->node(0) > elem->node(1))shape=10-i; 00805 // if((i>= 8&&i<=12) && elem->node(1) > elem->node(2))shape=20-i; 00806 // if((i>=13&&i<=17) && elem->node(0) > elem->node(2))shape=30-i; 00807 00808 // switch(shape) 00809 // { 00810 // //point functions 00811 // case 0: return pow<6>(r); 00812 // case 1: return pow<6>(x); 00813 // case 2: return pow<6>(y); 00814 00815 // //edge functions 00816 // case 3: return 6.*x *pow<5>(r); 00817 // case 4: return 15.*pow<2>(x)*pow<4>(r); 00818 // case 5: return 20.*pow<3>(x)*pow<3>(r); 00819 // case 6: return 15.*pow<4>(x)*pow<2>(r); 00820 // case 7: return 6.*pow<5>(x)*r; 00821 00822 // case 8: return 6.*y *pow<5>(x); 00823 // case 9: return 15.*pow<2>(y)*pow<4>(x); 00824 // case 10: return 20.*pow<3>(y)*pow<3>(x); 00825 // case 11: return 15.*pow<4>(y)*pow<2>(x); 00826 // case 12: return 6.*pow<5>(y)*x; 00827 00828 // case 13: return 6.*y *pow<5>(r); 00829 // case 14: return 15.*pow<2>(y)*pow<4>(r); 00830 // case 15: return 20.*pow<3>(y)*pow<3>(r); 00831 // case 16: return 15.*pow<4>(y)*pow<2>(r); 00832 // case 17: return 6.*pow<5>(y)*r; 00833 00834 // //inner functions 00835 // case 18: return 30.*x*y*pow<4>(r); 00836 // case 19: return 60.*x*pow<2>(y)*pow<3>(r); 00837 // case 20: return 60.* pow<2>(x)*y*pow<3>(r); 00838 // case 21: return 60.*x*pow<3>(y)*pow<2>(r); 00839 // case 22: return 60.*pow<3>(x)*y*pow<2>(r); 00840 // case 23: return 90.*pow<2>(x)*pow<2>(y)*pow<2>(r); 00841 // case 24: return 30.*x*pow<4>(y)*r; 00842 // case 25: return 60.*pow<2>(x)*pow<3>(y)*r; 00843 // case 26: return 60.*pow<3>(x)*pow<2>(y)*r; 00844 // case 27: return 30.*pow<4>(x)*y*r; 00845 00846 // } // switch shape 00847 // } // case TRI6 00848 00849 00850 00851 // // Bernstein shape functions on the quadrilateral. 00852 // case QUAD8: 00853 // case QUAD9: 00854 // { 00855 // // Compute quad shape functions as a tensor-product 00856 // const Real xi = p(0); 00857 // const Real eta = p(1); 00858 00859 // libmesh_assert_less (i, 49); 00860 00861 // // 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 00862 // static const unsigned int i0_reg[] = {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}; 00863 // static const unsigned int i1_reg[] = {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}; 00864 00865 // unsigned int i0=i0_reg[i]; 00866 // unsigned int i1=i1_reg[i]; 00867 00868 // if((i>= 4&&i<= 8) && elem->node(0) > elem->node(1)) i0=8-i0; 00869 // if((i>= 9&&i<=13) && elem->node(1) > elem->node(2)) i1=8-i1; 00870 // if((i>=14&&i<=18) && elem->node(3) > elem->node(2)) i0=8-i0; 00871 // if((i>=19&&i<=23) && elem->node(0) > elem->node(3)) i1=8-i1; 00872 00873 // return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0, xi)* 00874 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1, eta)); 00875 00876 // } // case QUAD8/9 00877 00878 // default: 00879 // libmesh_error(); 00880 00881 // } 00882 // // 7th-order Bernstein. 00883 // case SEVENTH: 00884 // { 00885 // switch (type) 00886 // { 00887 00888 // // Szabo-Babuska shape functions on the quadrilateral. 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, 64); 00896 00897 // // 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 00898 // 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}; 00899 // 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}; 00900 00901 // Real f=1.; 00902 00903 // switch(i) 00904 // { 00905 // case 5: // edge 0 nodes 00906 // case 7: 00907 // case 9: 00908 // if (elem->node(0) > elem->node(1))f = -1.; 00909 // break; 00910 // case 11: // edge 1 nodes 00911 // case 13: 00912 // case 15: 00913 // if (elem->node(1) > elem->node(2))f = -1.; 00914 // break; 00915 // case 17: // edge 2 nodes 00916 // case 19: 00917 // case 21: 00918 // if (elem->node(3) > elem->node(2))f = -1.; 00919 // break; 00920 // case 23: // edge 3 nodes 00921 // case 25: 00922 // case 27: 00923 // if (elem->node(0) > elem->node(3))f = -1.; 00924 // break; 00925 // } 00926 00927 // return f*(FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)* 00928 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)); 00929 00930 // } // case QUAD8/QUAD9 00931 00932 // default: 00933 // libmesh_error(); 00934 00935 // } // switch type 00936 00937 // } // case SEVENTH 00938 00939 00940 // // by default throw an error 00941 // default: 00942 // libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl; 00943 // libmesh_error(); 00944 // } 00945 00946 libmesh_error(); 00947 return 0.; 00948 } 00949 00950 00951 00952 template <> 00953 Real FE<2,BERNSTEIN>::shape_deriv(const ElemType, 00954 const Order, 00955 const unsigned int, 00956 const unsigned int, 00957 const Point&) 00958 { 00959 libMesh::err << "Bernstein polynomials require the element type\n" 00960 << "because edge orientation is needed." 00961 << std::endl; 00962 00963 libmesh_error(); 00964 return 0.; 00965 } 00966 00967 00968 00969 template <> 00970 Real FE<2,BERNSTEIN>::shape_deriv(const Elem* elem, 00971 const Order order, 00972 const unsigned int i, 00973 const unsigned int j, 00974 const Point& p) 00975 { 00976 libmesh_assert(elem); 00977 00978 const ElemType type = elem->type(); 00979 00980 const Order totalorder = static_cast<Order>(order + elem->p_level()); 00981 00982 switch (type) 00983 { 00984 // Hierarchic shape functions on the quadrilateral. 00985 case QUAD4: 00986 case QUAD9: 00987 { 00988 // Compute quad shape functions as a tensor-product 00989 const Real xi = p(0); 00990 const Real eta = p(1); 00991 00992 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)); 00993 00994 unsigned int i0, i1; 00995 00996 // Vertex DoFs 00997 if (i == 0) 00998 { i0 = 0; i1 = 0; } 00999 else if (i == 1) 01000 { i0 = 1; i1 = 0; } 01001 else if (i == 2) 01002 { i0 = 1; i1 = 1; } 01003 else if (i == 3) 01004 { i0 = 0; i1 = 1; } 01005 01006 01007 // Edge DoFs 01008 else if (i < totalorder + 3u) 01009 { i0 = i - 2; i1 = 0; } 01010 else if (i < 2u*totalorder + 2) 01011 { i0 = 1; i1 = i - totalorder - 1; } 01012 else if (i < 3u*totalorder + 1) 01013 { i0 = i - 2u*totalorder; i1 = 1; } 01014 else if (i < 4u*totalorder) 01015 { i0 = 0; i1 = i - 3u*totalorder + 1; } 01016 // Interior DoFs 01017 else 01018 { 01019 unsigned int basisnum = i - 4*totalorder; 01020 i0 = square_number_column[basisnum] + 2; 01021 i1 = square_number_row[basisnum] + 2; 01022 } 01023 01024 01025 // Flip odd degree of freedom values if necessary 01026 // to keep continuity on sides 01027 if ((i>= 4 && i<= 4+ totalorder-2u) && elem->point(0) > elem->point(1)) i0=totalorder+2-i0; // 01028 else if((i>= 4+ totalorder-1u && i<= 4+2*totalorder-3u) && elem->point(1) > elem->point(2)) i1=totalorder+2-i1; 01029 else if((i>= 4+2*totalorder-2u && i<= 4+3*totalorder-4u) && elem->point(3) > elem->point(2)) i0=totalorder+2-i0; 01030 else if((i>= 4+3*totalorder-3u && i<= 4+4*totalorder-5u) && elem->point(0) > elem->point(3)) i1=totalorder+2-i1; 01031 01032 switch (j) 01033 { 01034 // d()/dxi 01035 case 0: 01036 return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)* 01037 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, eta)); 01038 01039 // d()/deta 01040 case 1: 01041 return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, xi)* 01042 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta)); 01043 01044 } 01045 } 01046 01047 // Bernstein shape functions on the 8-noded quadrilateral 01048 // is handled separately. 01049 case QUAD8: 01050 { 01051 libmesh_assert_less (totalorder, 3); 01052 01053 const Real xi = p(0); 01054 const Real eta = p(1); 01055 01056 libmesh_assert_less (i, 8); 01057 01058 // 0 1 2 3 4 5 6 7 8 01059 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2}; 01060 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2}; 01061 static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5}; 01062 switch (j) 01063 { 01064 // d()/dxi 01065 case 0: 01066 return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)* 01067 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta) 01068 +scal[i]* 01069 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[8], 0, xi)* 01070 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[8], eta)); 01071 01072 // d()/deta 01073 case 1: 01074 return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)* 01075 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta) 01076 +scal[i]* 01077 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[8], xi)* 01078 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[8], 0, eta)); 01079 01080 default: 01081 libmesh_error(); 01082 } 01083 } 01084 01085 case TRI3: 01086 libmesh_assert_less (totalorder, 2); 01087 case TRI6: 01088 { 01089 // I have been lazy here and am using finite differences 01090 // to compute the derivatives! 01091 const Real eps = 1.e-6; 01092 01093 switch (j) 01094 { 01095 // d()/dxi 01096 case 0: 01097 { 01098 const Point pp(p(0)+eps, p(1)); 01099 const Point pm(p(0)-eps, p(1)); 01100 01101 return (FE<2,BERNSTEIN>::shape(elem, totalorder, i, pp) - 01102 FE<2,BERNSTEIN>::shape(elem, totalorder, i, pm))/2./eps; 01103 } 01104 01105 // d()/deta 01106 case 1: 01107 { 01108 const Point pp(p(0), p(1)+eps); 01109 const Point pm(p(0), p(1)-eps); 01110 01111 return (FE<2,BERNSTEIN>::shape(elem, totalorder, i, pp) - 01112 FE<2,BERNSTEIN>::shape(elem, totalorder, i, pm))/2./eps; 01113 } 01114 01115 01116 default: 01117 libmesh_error(); 01118 } 01119 } 01120 01121 default: 01122 { 01123 libMesh::err << "ERROR: Unsupported element type!" << std::endl; 01124 libmesh_error(); 01125 } 01126 01127 } 01128 01129 // old code 01130 // switch (totalorder) 01131 // { 01132 01133 // // 1st order Bernsteins are aquivalent to Lagrange elements. 01134 // case FIRST: 01135 // { 01136 // switch (type) 01137 // { 01138 // // Bernstein shape functions on the triangle. 01139 // case TRI6: 01140 // { 01141 // // I have been lazy here and am using finite differences 01142 // // to compute the derivatives! 01143 // const Real eps = 1.e-6; 01144 01145 // libmesh_assert_less (i, 3); 01146 // libmesh_assert_less (j, 2); 01147 01148 // switch (j) 01149 // { 01150 // // d()/dxi 01151 // case 0: 01152 // { 01153 // const Point pp(p(0)+eps, p(1)); 01154 // const Point pm(p(0)-eps, p(1)); 01155 01156 // return (FE<2,BERNSTEIN>::shape(elem, order, i, pp) - 01157 // FE<2,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01158 // } 01159 01160 // // d()/deta 01161 // case 1: 01162 // { 01163 // const Point pp(p(0), p(1)+eps); 01164 // const Point pm(p(0), p(1)-eps); 01165 01166 // return (FE<2,BERNSTEIN>::shape(elem, order, i, pp) - 01167 // FE<2,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01168 // } 01169 01170 01171 // default: 01172 // libmesh_error(); 01173 // } 01174 // } 01175 01176 // // Bernstein shape functions on the quadrilateral. 01177 // case QUAD8: 01178 // case QUAD9: 01179 // { 01180 // // Compute quad shape functions as a tensor-product 01181 // const Real xi = p(0); 01182 // const Real eta = p(1); 01183 01184 // libmesh_assert_less (i, 4); 01185 01186 // // 0 1 2 3 01187 // static const unsigned int i0[] = {0, 1, 1, 0}; 01188 // static const unsigned int i1[] = {0, 0, 1, 1}; 01189 01190 // switch (j) 01191 // { 01192 // // d()/dxi 01193 // case 0: 01194 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)* 01195 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)); 01196 01197 // // d()/deta 01198 // case 1: 01199 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)* 01200 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)); 01201 01202 // default: 01203 // libmesh_error(); 01204 // } 01205 // } 01206 01207 // default: 01208 // libmesh_error(); 01209 // } 01210 // } 01211 01212 // // second order Bernsteins. 01213 // case SECOND: 01214 // { 01215 // switch (type) 01216 // { 01217 // // Bernstein shape functions on the triangle. 01218 // case TRI6: 01219 // { 01220 // // I have been lazy here and am using finite differences 01221 // // to compute the derivatives! 01222 // const Real eps = 1.e-6; 01223 01224 // libmesh_assert_less (i, 6); 01225 // libmesh_assert_less (j, 2); 01226 01227 // switch (j) 01228 // { 01229 // // d()/dxi 01230 // case 0: 01231 // { 01232 // const Point pp(p(0)+eps, p(1)); 01233 // const Point pm(p(0)-eps, p(1)); 01234 01235 // return (FE<2,BERNSTEIN>::shape(elem, order, i, pp) - 01236 // FE<2,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01237 // } 01238 01239 // // d()/deta 01240 // case 1: 01241 // { 01242 // const Point pp(p(0), p(1)+eps); 01243 // const Point pm(p(0), p(1)-eps); 01244 01245 // return (FE<2,BERNSTEIN>::shape(elem, order, i, pp) - 01246 // FE<2,BERNSTEIN>::shape(elem, order, i, pm))/2./eps; 01247 // } 01248 01249 01250 // default: 01251 // libmesh_error(); 01252 // } 01253 // } 01254 01255 // // Bernstein shape functions on the 8-noded quadrilateral. 01256 // case QUAD8: 01257 // { 01258 // const Real xi = p(0); 01259 // const Real eta = p(1); 01260 01261 // libmesh_assert_less (i, 8); 01262 01263 // // 0 1 2 3 4 5 6 7 8 01264 // static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2}; 01265 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2}; 01266 // static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5}; 01267 // switch (j) 01268 // { 01269 // // d()/dxi 01270 // case 0: 01271 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)* 01272 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta) 01273 // +scal[i]* 01274 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[8], 0, xi)* 01275 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[8], eta)); 01276 01277 // // d()/deta 01278 // case 1: 01279 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)* 01280 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta) 01281 // +scal[i]* 01282 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[8], xi)* 01283 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[8], 0, eta)); 01284 01285 // default: 01286 // libmesh_error(); 01287 // } 01288 // } 01289 01290 01291 // // Bernstein shape functions on the 9-noded quadrilateral. 01292 // case QUAD9: 01293 // { 01294 // // Compute quad shape functions as a tensor-product 01295 // const Real xi = p(0); 01296 // const Real eta = p(1); 01297 01298 // libmesh_assert_less (i, 9); 01299 01300 // // 0 1 2 3 4 5 6 7 8 01301 // static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2}; 01302 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2}; 01303 01304 // switch (j) 01305 // { 01306 // // d()/dxi 01307 // case 0: 01308 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)* 01309 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)); 01310 01311 // // d()/deta 01312 // case 1: 01313 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)* 01314 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)); 01315 01316 // default: 01317 // libmesh_error(); 01318 // } 01319 // } 01320 01321 // default: 01322 // libmesh_error(); 01323 // } 01324 // } 01325 01326 01327 // // 3rd-order Bernsteins. 01328 // case THIRD: 01329 // { 01330 // switch (type) 01331 // { 01332 // // Bernstein shape functions on the triangle. 01333 // case TRI6: 01334 // { 01335 // // I have been lazy here and am using finite differences 01336 // // to compute the derivatives! 01337 // const Real eps = 1.e-6; 01338 01339 // libmesh_assert_less (i, 10); 01340 // libmesh_assert_less (j, 2); 01341 01342 01343 // unsigned int shape=i; 01344 01345 // /** 01346 // * Note, since the derivatives are computed using FE::shape 01347 // * the continuity along neighboring elements (shape_flip) 01348 // * is already considered. 01349 // */ 01350 01351 01352 // switch (j) 01353 // { 01354 // // d()/dxi 01355 // case 0: 01356 // { 01357 // const Point pp(p(0)+eps, p(1)); 01358 // const Point pm(p(0)-eps, p(1)); 01359 01360 // return (FE<2,BERNSTEIN>::shape(elem, totalorder, shape, pp) - 01361 // FE<2,BERNSTEIN>::shape(elem, totalorder, shape, pm))/2./eps; 01362 // } 01363 01364 // // d()/deta 01365 // case 1: 01366 // { 01367 // const Point pp(p(0), p(1)+eps); 01368 // const Point pm(p(0), p(1)-eps); 01369 01370 // return (FE<2,BERNSTEIN>::shape(elem, totalorder, shape, pp) - 01371 // FE<2,BERNSTEIN>::shape(elem, totalorder, shape, pm))/2./eps; 01372 // } 01373 01374 01375 // default: 01376 // libmesh_error(); 01377 // } 01378 // } 01379 01380 01381 01382 // // Bernstein shape functions on the quadrilateral. 01383 // case QUAD8: 01384 // case QUAD9: 01385 // { 01386 // // Compute quad shape functions as a tensor-product 01387 // const Real xi = p(0); 01388 // const Real eta = p(1); 01389 01390 // libmesh_assert_less (i, 16); 01391 01392 // // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 01393 // static const unsigned int i0_reg[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3}; 01394 // static const unsigned int i1_reg[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3}; 01395 01396 // unsigned int i0=i0_reg[i]; 01397 // unsigned int i1=i1_reg[i]; 01398 01399 // if((i== 4||i== 5) && elem->node(0) > elem->node(1)) i0=5-i0; 01400 // if((i== 6||i== 7) && elem->node(1) > elem->node(2)) i1=5-i1; 01401 // if((i== 8||i== 9) && elem->node(3) > elem->node(2)) i0=5-i0; 01402 // if((i==10||i==11) && elem->node(0) > elem->node(3)) i1=5-i1; 01403 01404 // switch (j) 01405 // { 01406 // // d()/dxi 01407 // case 0: 01408 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)* 01409 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, eta)); 01410 01411 // // d()/deta 01412 // case 1: 01413 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, xi)* 01414 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta)); 01415 01416 // default: 01417 // libmesh_error(); 01418 // } 01419 // } 01420 01421 // default: 01422 // libmesh_error(); 01423 // } 01424 // } 01425 01426 01427 01428 01429 // // 4th-order Bernsteins. 01430 // case FOURTH: 01431 // { 01432 // switch (type) 01433 // { 01434 // // Bernstein shape functions on the triangle. 01435 // case TRI6: 01436 // { 01437 // // I have been lazy here and am using finite differences 01438 // // to compute the derivatives! 01439 // const Real eps = 1.e-6; 01440 01441 // libmesh_assert_less (i, 15); 01442 // libmesh_assert_less (j, 2); 01443 01444 // unsigned int shape=i; 01445 01446 // switch (j) 01447 // { 01448 // // d()/dxi 01449 // case 0: 01450 // { 01451 // const Point pp(p(0)+eps, p(1)); 01452 // const Point pm(p(0)-eps, p(1)); 01453 01454 // return (FE<2,BERNSTEIN>::shape(elem, order, shape, pp) - 01455 // FE<2,BERNSTEIN>::shape(elem, order, shape, pm))/2./eps; 01456 // } 01457 01458 // // d()/deta 01459 // case 1: 01460 // { 01461 // const Point pp(p(0), p(1)+eps); 01462 // const Point pm(p(0), p(1)-eps); 01463 01464 // return (FE<2,BERNSTEIN>::shape(elem, order, shape, pp) - 01465 // FE<2,BERNSTEIN>::shape(elem, order, shape, pm))/2./eps; 01466 // } 01467 01468 01469 // default: 01470 // libmesh_error(); 01471 // } 01472 // } 01473 01474 01475 01476 // // Bernstein shape functions on the quadrilateral. 01477 // case QUAD8: 01478 // case QUAD9: 01479 // { 01480 // // Compute quad shape functions as a tensor-product 01481 // const Real xi = p(0); 01482 // const Real eta = p(1); 01483 01484 // libmesh_assert_less (i, 25); 01485 01486 // // 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 01487 // static const unsigned int i0_reg[] = {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}; 01488 // static const unsigned int i1_reg[] = {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}; 01489 01490 // unsigned int i0=i0_reg[i]; 01491 // unsigned int i1=i1_reg[i]; 01492 01493 // if((i== 4||i== 6) && elem->node(0) > elem->node(1)) i0=6-i0; 01494 // if((i== 7||i== 9) && elem->node(1) > elem->node(2)) i1=6-i1; 01495 // if((i==10||i==12) && elem->node(3) > elem->node(2)) i0=6-i0; 01496 // if((i==13||i==15) && elem->node(0) > elem->node(3)) i1=6-i1; 01497 01498 01499 // switch (j) 01500 // { 01501 // // d()/dxi 01502 // case 0: 01503 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)* 01504 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, eta)); 01505 01506 // // d()/deta 01507 // case 1: 01508 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, xi)* 01509 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta)); 01510 01511 // default: 01512 // libmesh_error(); 01513 // } 01514 // } 01515 01516 // default: 01517 // libmesh_error(); 01518 // } 01519 // } 01520 01521 01522 01523 01524 // // 5th-order Bernsteins. 01525 // case FIFTH: 01526 // { 01527 // // Bernstein shape functions on the quadrilateral. 01528 // switch (type) 01529 // { 01530 // // Bernstein shape functions on the triangle. 01531 // case TRI6: 01532 // { 01533 // // I have been lazy here and am using finite differences 01534 // // to compute the derivatives! 01535 // const Real eps = 1.e-6; 01536 01537 // libmesh_assert_less (i, 21); 01538 // libmesh_assert_less (j, 2); 01539 01540 // unsigned int shape=i; 01541 01542 // switch (j) 01543 // { 01544 // // d()/dxi 01545 // case 0: 01546 // { 01547 // const Point pp(p(0)+eps, p(1)); 01548 // const Point pm(p(0)-eps, p(1)); 01549 01550 // return (FE<2,BERNSTEIN>::shape(elem, order, shape, pp) - 01551 // FE<2,BERNSTEIN>::shape(elem, order, shape, pm))/2./eps; 01552 // } 01553 01554 // // d()/deta 01555 // case 1: 01556 // { 01557 // const Point pp(p(0), p(1)+eps); 01558 // const Point pm(p(0), p(1)-eps); 01559 01560 // return (FE<2,BERNSTEIN>::shape(elem, order, shape, pp) - 01561 // FE<2,BERNSTEIN>::shape(elem, order, shape, pm))/2./eps; 01562 // } 01563 01564 01565 // default: 01566 // libmesh_error(); 01567 // } 01568 // } // case TRI6 01569 01570 01571 01572 // case QUAD8: 01573 // case QUAD9: 01574 // { 01575 // // Compute quad shape functions as a tensor-product 01576 // const Real xi = p(0); 01577 // const Real eta = p(1); 01578 01579 // libmesh_assert_less (i, 36); 01580 01581 // // 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 01582 // static const unsigned int i0_reg[] = {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}; 01583 // static const unsigned int i1_reg[] = {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}; 01584 01585 // unsigned int i0=i0_reg[i]; 01586 // unsigned int i1=i1_reg[i]; 01587 01588 // if((i>= 4&&i<= 7) && elem->node(0) > elem->node(1)) i0=7-i0; 01589 // if((i>= 8&&i<=11) && elem->node(1) > elem->node(2)) i1=7-i1; 01590 // if((i>=12&&i<=15) && elem->node(3) > elem->node(2)) i0=7-i0; 01591 // if((i>=16&&i<=19) && elem->node(0) > elem->node(3)) i1=7-i1; 01592 01593 01594 // switch (j) 01595 // { 01596 // // d()/dxi 01597 // case 0: 01598 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)* 01599 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, eta)); 01600 01601 // // d()/deta 01602 // case 1: 01603 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, xi)* 01604 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta)); 01605 01606 // default: 01607 // libmesh_error(); 01608 // } 01609 // } // case QUAD8/9 01610 01611 // default: 01612 // libmesh_error(); 01613 // } 01614 01615 // } 01616 01617 01618 // // 6th-order Bernsteins. 01619 // case SIXTH: 01620 // { 01621 // // Bernstein shape functions on the quadrilateral. 01622 // switch (type) 01623 // { 01624 // // Bernstein shape functions on the triangle. 01625 // case TRI6: 01626 // { 01627 // // I have been lazy here and am using finite differences 01628 // // to compute the derivatives! 01629 // const Real eps = 1.e-6; 01630 01631 // libmesh_assert_less (i, 28); 01632 // libmesh_assert_less (j, 2); 01633 01634 // unsigned int shape=i; 01635 01636 // switch (j) 01637 // { 01638 // // d()/dxi 01639 // case 0: 01640 // { 01641 // const Point pp(p(0)+eps, p(1)); 01642 // const Point pm(p(0)-eps, p(1)); 01643 01644 // return (FE<2,BERNSTEIN>::shape(elem, order, shape, pp) - 01645 // FE<2,BERNSTEIN>::shape(elem, order, shape, pm))/2./eps; 01646 // } 01647 01648 // // d()/deta 01649 // case 1: 01650 // { 01651 // const Point pp(p(0), p(1)+eps); 01652 // const Point pm(p(0), p(1)-eps); 01653 01654 // return (FE<2,BERNSTEIN>::shape(elem, order, shape, pp) - 01655 // FE<2,BERNSTEIN>::shape(elem, order, shape, pm))/2./eps; 01656 // } 01657 01658 01659 // default: 01660 // libmesh_error(); 01661 // } 01662 // } // case TRI6 01663 01664 01665 01666 // case QUAD8: 01667 // case QUAD9: 01668 // { 01669 // // Compute quad shape functions as a tensor-product 01670 // const Real xi = p(0); 01671 // const Real eta = p(1); 01672 01673 // libmesh_assert_less (i, 49); 01674 01675 // // 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 01676 // static const unsigned int i0_reg[] = {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}; 01677 // static const unsigned int i1_reg[] = {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}; 01678 01679 // unsigned int i0=i0_reg[i]; 01680 // unsigned int i1=i1_reg[i]; 01681 01682 // if((i>= 4&&i<= 8) && elem->node(0) > elem->node(1)) i0=8-i0; 01683 // if((i>= 9&&i<=13) && elem->node(1) > elem->node(2)) i1=8-i1; 01684 // if((i>=14&&i<=18) && elem->node(3) > elem->node(2)) i0=8-i0; 01685 // if((i>=19&&i<=23) && elem->node(0) > elem->node(3)) i1=8-i1; 01686 01687 01688 // switch (j) 01689 // { 01690 // // d()/dxi 01691 // case 0: 01692 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)* 01693 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, eta)); 01694 01695 // // d()/deta 01696 // case 1: 01697 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, xi)* 01698 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta)); 01699 01700 // default: 01701 // libmesh_error(); 01702 // } 01703 // } // case QUAD8/9 01704 01705 // default: 01706 // libmesh_error(); 01707 // } 01708 01709 // // 7th-order Bernstein. 01710 // case SEVENTH: 01711 // { 01712 // // Szabo-Babuska shape functions on the quadrilateral. 01713 // switch (type) 01714 // { 01715 01716 01717 // case QUAD9: 01718 // { 01719 // // Compute quad shape functions as a tensor-product 01720 // const Real xi = p(0); 01721 // const Real eta = p(1); 01722 01723 // libmesh_assert_less (i, 64); 01724 01725 // // 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 01726 // 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}; 01727 // 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}; 01728 01729 // Real f=1.; 01730 01731 // switch(i) 01732 // { 01733 // case 5: // edge 0 nodes 01734 // case 7: 01735 // case 9: 01736 // if (elem->node(0) > elem->node(1))f = -1.; 01737 // break; 01738 // case 11: // edge 1 nodes 01739 // case 13: 01740 // case 15: 01741 // if (elem->node(1) > elem->node(2))f = -1.; 01742 // break; 01743 // case 17: // edge 2 nodes 01744 // case 19: 01745 // case 21: 01746 // if (elem->node(3) > elem->node(2))f = -1.; 01747 // break; 01748 // case 23: // edge 3 nodes 01749 // case 25: 01750 // case 27: 01751 // if (elem->node(0) > elem->node(3))f = -1.; 01752 // break; 01753 // } 01754 01755 01756 // switch (j) 01757 // { 01758 // // d()/dxi 01759 // case 0: 01760 // return f*(FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)* 01761 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)); 01762 01763 // // d()/deta 01764 // case 1: 01765 // return f*(FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)* 01766 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)); 01767 01768 // default: 01769 // libmesh_error(); 01770 // } 01771 // } 01772 01773 // default: 01774 // libmesh_error(); 01775 // } 01776 // } 01777 01778 // } 01779 // // by default throw an error;call the orientation-independent shape functions 01780 // default: 01781 // libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl; 01782 // libmesh_error(); 01783 // } 01784 01785 01786 libmesh_error(); 01787 return 0.; 01788 } 01789 01790 01791 01792 template <> 01793 Real FE<2,BERNSTEIN>::shape_second_deriv(const ElemType, 01794 const Order, 01795 const unsigned int, 01796 const unsigned int, 01797 const Point&) 01798 { 01799 static bool warning_given = false; 01800 01801 if (!warning_given) 01802 libMesh::err << "Second derivatives for Bernstein elements " 01803 << "are not yet implemented!" 01804 << std::endl; 01805 01806 warning_given = true; 01807 return 0.; 01808 } 01809 01810 01811 01812 template <> 01813 Real FE<2,BERNSTEIN>::shape_second_deriv(const Elem*, 01814 const Order, 01815 const unsigned int, 01816 const unsigned int, 01817 const Point&) 01818 { 01819 static bool warning_given = false; 01820 01821 if (!warning_given) 01822 libMesh::err << "Second derivatives for Bernstein elements " 01823 << "are not yet implemented!" 01824 << std::endl; 01825 01826 warning_given = true; 01827 return 0.; 01828 } 01829 01830 } // namespace libMesh 01831 01832 01833 #endif // LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:46 UTC
Hosted By: