fe_l2_hierarchic_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 // C++ includes 00020 00021 // Local includes 00022 #include "libmesh/fe.h" 00023 #include "libmesh/elem.h" 00024 #include "libmesh/number_lookups.h" 00025 #include "libmesh/utility.h" 00026 00027 00028 namespace libMesh 00029 { 00030 00031 template <> 00032 Real FE<2,L2_HIERARCHIC>::shape(const ElemType, 00033 const Order, 00034 const unsigned int, 00035 const Point&) 00036 { 00037 libMesh::err << "Hierarchic polynomials require the element type\n" 00038 << "because edge orientation is needed." 00039 << std::endl; 00040 00041 libmesh_error(); 00042 return 0.; 00043 } 00044 00045 00046 00047 template <> 00048 Real FE<2,L2_HIERARCHIC>::shape(const Elem* elem, 00049 const Order order, 00050 const unsigned int i, 00051 const Point& p) 00052 { 00053 libmesh_assert(elem); 00054 00055 const Order totalorder = static_cast<Order>(order+elem->p_level()); 00056 libmesh_assert_greater (totalorder, 0); 00057 00058 switch (elem->type()) 00059 { 00060 case TRI3: 00061 case TRI6: 00062 { 00063 const Real zeta1 = p(0); 00064 const Real zeta2 = p(1); 00065 const Real zeta0 = 1. - zeta1 - zeta2; 00066 00067 libmesh_assert_less (i, (totalorder+1u)*(totalorder+2u)/2); 00068 libmesh_assert (elem->type() == TRI6 || totalorder < 2); 00069 00070 // Vertex DoFs 00071 if (i == 0) 00072 return zeta0; 00073 else if (i == 1) 00074 return zeta1; 00075 else if (i == 2) 00076 return zeta2; 00077 // Edge DoFs 00078 else if (i < totalorder + 2u) 00079 { 00080 // Avoid returning NaN on vertices! 00081 if (zeta0 + zeta1 == 0.) 00082 return 0.; 00083 00084 const unsigned int basisorder = i - 1; 00085 // Get factors to account for edge-flipping 00086 Real f0 = 1; 00087 if (basisorder%2 && (elem->point(0) > elem->point(1))) 00088 f0 = -1.; 00089 00090 Real edgeval = (zeta1 - zeta0) / (zeta1 + zeta0); 00091 Real crossfunc = zeta0 + zeta1; 00092 for (unsigned int n=1; n != basisorder; ++n) 00093 crossfunc *= (zeta0 + zeta1); 00094 00095 return f0 * crossfunc * 00096 FE<1,L2_HIERARCHIC>::shape(EDGE3, totalorder, 00097 basisorder, edgeval); 00098 } 00099 else if (i < 2u*totalorder + 1) 00100 { 00101 // Avoid returning NaN on vertices! 00102 if (zeta1 + zeta2 == 0.) 00103 return 0.; 00104 00105 const unsigned int basisorder = i - totalorder; 00106 // Get factors to account for edge-flipping 00107 Real f1 = 1; 00108 if (basisorder%2 && (elem->point(1) > elem->point(2))) 00109 f1 = -1.; 00110 00111 Real edgeval = (zeta2 - zeta1) / (zeta2 + zeta1); 00112 Real crossfunc = zeta2 + zeta1; 00113 for (unsigned int n=1; n != basisorder; ++n) 00114 crossfunc *= (zeta2 + zeta1); 00115 00116 return f1 * crossfunc * 00117 FE<1,L2_HIERARCHIC>::shape(EDGE3, totalorder, 00118 basisorder, edgeval); 00119 } 00120 else if (i < 3u*totalorder) 00121 { 00122 // Avoid returning NaN on vertices! 00123 if (zeta0 + zeta2 == 0.) 00124 return 0.; 00125 00126 const unsigned int basisorder = i - (2u*totalorder) + 1; 00127 // Get factors to account for edge-flipping 00128 Real f2 = 1; 00129 if (basisorder%2 && (elem->point(2) > elem->point(0))) 00130 f2 = -1.; 00131 00132 Real edgeval = (zeta0 - zeta2) / (zeta0 + zeta2); 00133 Real crossfunc = zeta0 + zeta2; 00134 for (unsigned int n=1; n != basisorder; ++n) 00135 crossfunc *= (zeta0 + zeta2); 00136 00137 return f2 * crossfunc * 00138 FE<1,L2_HIERARCHIC>::shape(EDGE3, totalorder, 00139 basisorder, edgeval); 00140 } 00141 // Interior DoFs 00142 else 00143 { 00144 const unsigned int basisnum = i - (3u*totalorder); 00145 unsigned int exp0 = triangular_number_column[basisnum] + 1; 00146 unsigned int exp1 = triangular_number_row[basisnum] + 1 - 00147 triangular_number_column[basisnum]; 00148 00149 Real returnval = 1; 00150 for (unsigned int n = 0; n != exp0; ++n) 00151 returnval *= zeta0; 00152 for (unsigned int n = 0; n != exp1; ++n) 00153 returnval *= zeta1; 00154 returnval *= zeta2; 00155 return returnval; 00156 } 00157 } 00158 00159 // Hierarchic shape functions on the quadrilateral. 00160 case QUAD4: 00161 libmesh_assert_less (totalorder, 2); 00162 case QUAD8: 00163 case QUAD9: 00164 { 00165 // Compute quad shape functions as a tensor-product 00166 const Real xi = p(0); 00167 const Real eta = p(1); 00168 00169 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)); 00170 00171 // Example i, i0, i1 values for totalorder = 5: 00172 // 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 00173 // 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}; 00174 // 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}; 00175 00176 unsigned int i0, i1; 00177 00178 // Vertex DoFs 00179 if (i == 0) 00180 { i0 = 0; i1 = 0; } 00181 else if (i == 1) 00182 { i0 = 1; i1 = 0; } 00183 else if (i == 2) 00184 { i0 = 1; i1 = 1; } 00185 else if (i == 3) 00186 { i0 = 0; i1 = 1; } 00187 // Edge DoFs 00188 else if (i < totalorder + 3u) 00189 { i0 = i - 2; i1 = 0; } 00190 else if (i < 2u*totalorder + 2) 00191 { i0 = 1; i1 = i - totalorder - 1; } 00192 else if (i < 3u*totalorder + 1) 00193 { i0 = i - 2u*totalorder; i1 = 1; } 00194 else if (i < 4u*totalorder) 00195 { i0 = 0; i1 = i - 3u*totalorder + 1; } 00196 // Interior DoFs 00197 else 00198 { 00199 unsigned int basisnum = i - 4*totalorder; 00200 i0 = square_number_column[basisnum] + 2; 00201 i1 = square_number_row[basisnum] + 2; 00202 } 00203 00204 // Flip odd degree of freedom values if necessary 00205 // to keep continuity on sides 00206 Real f = 1.; 00207 00208 if ((i0%2) && (i0 > 2) && (i1 == 0)) 00209 f = (elem->point(0) > elem->point(1))?-1.:1.; 00210 else if ((i0%2) && (i0>2) && (i1 == 1)) 00211 f = (elem->point(3) > elem->point(2))?-1.:1.; 00212 else if ((i0 == 0) && (i1%2) && (i1>2)) 00213 f = (elem->point(0) > elem->point(3))?-1.:1.; 00214 else if ((i0 == 1) && (i1%2) && (i1>2)) 00215 f = (elem->point(1) > elem->point(2))?-1.:1.; 00216 00217 return f*(FE<1,L2_HIERARCHIC>::shape(EDGE3, totalorder, i0, xi)* 00218 FE<1,L2_HIERARCHIC>::shape(EDGE3, totalorder, i1, eta)); 00219 } 00220 00221 default: 00222 libMesh::err << "ERROR: Unsupported element type!" << std::endl; 00223 libmesh_error(); 00224 } 00225 00226 return 0.; 00227 } 00228 00229 00230 00231 template <> 00232 Real FE<2,L2_HIERARCHIC>::shape_deriv(const ElemType, 00233 const Order, 00234 const unsigned int, 00235 const unsigned int, 00236 const Point&) 00237 { 00238 libMesh::err << "Hierarchic polynomials require the element type\n" 00239 << "because edge orientation is needed." 00240 << std::endl; 00241 00242 libmesh_error(); 00243 return 0.; 00244 } 00245 00246 00247 00248 template <> 00249 Real FE<2,L2_HIERARCHIC>::shape_deriv(const Elem* elem, 00250 const Order order, 00251 const unsigned int i, 00252 const unsigned int j, 00253 const Point& p) 00254 { 00255 libmesh_assert(elem); 00256 00257 const ElemType type = elem->type(); 00258 00259 const Order totalorder = static_cast<Order>(order+elem->p_level()); 00260 00261 libmesh_assert_greater (totalorder, 0); 00262 00263 switch (type) 00264 { 00265 // 1st & 2nd-order Hierarchics. 00266 case TRI3: 00267 case TRI6: 00268 { 00269 const Real eps = 1.e-6; 00270 00271 libmesh_assert_less (j, 2); 00272 00273 switch (j) 00274 { 00275 // d()/dxi 00276 case 0: 00277 { 00278 const Point pp(p(0)+eps, p(1)); 00279 const Point pm(p(0)-eps, p(1)); 00280 00281 return (FE<2,L2_HIERARCHIC>::shape(elem, order, i, pp) - 00282 FE<2,L2_HIERARCHIC>::shape(elem, order, i, pm))/2./eps; 00283 } 00284 00285 // d()/deta 00286 case 1: 00287 { 00288 const Point pp(p(0), p(1)+eps); 00289 const Point pm(p(0), p(1)-eps); 00290 00291 return (FE<2,L2_HIERARCHIC>::shape(elem, order, i, pp) - 00292 FE<2,L2_HIERARCHIC>::shape(elem, order, i, pm))/2./eps; 00293 } 00294 00295 00296 default: 00297 libmesh_error(); 00298 } 00299 } 00300 00301 case QUAD4: 00302 libmesh_assert_less (totalorder, 2); 00303 case QUAD8: 00304 case QUAD9: 00305 { 00306 // Compute quad shape functions as a tensor-product 00307 const Real xi = p(0); 00308 const Real eta = p(1); 00309 00310 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)); 00311 00312 // Example i, i0, i1 values for totalorder = 5: 00313 // 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 00314 // 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}; 00315 // 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}; 00316 00317 unsigned int i0, i1; 00318 00319 // Vertex DoFs 00320 if (i == 0) 00321 { i0 = 0; i1 = 0; } 00322 else if (i == 1) 00323 { i0 = 1; i1 = 0; } 00324 else if (i == 2) 00325 { i0 = 1; i1 = 1; } 00326 else if (i == 3) 00327 { i0 = 0; i1 = 1; } 00328 // Edge DoFs 00329 else if (i < totalorder + 3u) 00330 { i0 = i - 2; i1 = 0; } 00331 else if (i < 2u*totalorder + 2) 00332 { i0 = 1; i1 = i - totalorder - 1; } 00333 else if (i < 3u*totalorder + 1u) 00334 { i0 = i - 2u*totalorder; i1 = 1; } 00335 else if (i < 4u*totalorder) 00336 { i0 = 0; i1 = i - 3u*totalorder + 1; } 00337 // Interior DoFs 00338 else 00339 { 00340 unsigned int basisnum = i - 4*totalorder; 00341 i0 = square_number_column[basisnum] + 2; 00342 i1 = square_number_row[basisnum] + 2; 00343 } 00344 00345 // Flip odd degree of freedom values if necessary 00346 // to keep continuity on sides 00347 Real f = 1.; 00348 00349 if ((i0%2) && (i0 > 2) && (i1 == 0)) 00350 f = (elem->point(0) > elem->point(1))?-1.:1.; 00351 else if ((i0%2) && (i0>2) && (i1 == 1)) 00352 f = (elem->point(3) > elem->point(2))?-1.:1.; 00353 else if ((i0 == 0) && (i1%2) && (i1>2)) 00354 f = (elem->point(0) > elem->point(3))?-1.:1.; 00355 else if ((i0 == 1) && (i1%2) && (i1>2)) 00356 f = (elem->point(1) > elem->point(2))?-1.:1.; 00357 00358 switch (j) 00359 { 00360 // d()/dxi 00361 case 0: 00362 return f*(FE<1,L2_HIERARCHIC>::shape_deriv(EDGE3, totalorder, i0, 0, xi)* 00363 FE<1,L2_HIERARCHIC>::shape (EDGE3, totalorder, i1, eta)); 00364 00365 // d()/deta 00366 case 1: 00367 return f*(FE<1,L2_HIERARCHIC>::shape (EDGE3, totalorder, i0, xi)* 00368 FE<1,L2_HIERARCHIC>::shape_deriv(EDGE3, totalorder, i1, 0, eta)); 00369 00370 default: 00371 libmesh_error(); 00372 } 00373 00374 } 00375 00376 default: 00377 libMesh::err << "ERROR: Unsupported element type!" << std::endl; 00378 libmesh_error(); 00379 } 00380 00381 return 0.; 00382 } 00383 00384 00385 00386 template <> 00387 Real FE<2,L2_HIERARCHIC>::shape_second_deriv(const ElemType, 00388 const Order, 00389 const unsigned int, 00390 const unsigned int, 00391 const Point&) 00392 { 00393 libMesh::err << "Hierarchic polynomials require the element type\n" 00394 << "because edge orientation is needed." 00395 << std::endl; 00396 00397 libmesh_error(); 00398 return 0.; 00399 } 00400 00401 00402 00403 template <> 00404 Real FE<2,L2_HIERARCHIC>::shape_second_deriv(const Elem* elem, 00405 const Order order, 00406 const unsigned int i, 00407 const unsigned int j, 00408 const Point& p) 00409 { 00410 libmesh_assert(elem); 00411 00412 // I have been lazy here and am using finite differences 00413 // to compute the derivatives! 00414 const Real eps = 1.e-6; 00415 Point pp, pm; 00416 unsigned int prevj = libMesh::invalid_uint; 00417 00418 switch (j) 00419 { 00420 // d^2()/dxi^2 00421 case 0: 00422 { 00423 pp = Point(p(0)+eps, p(1)); 00424 pm = Point(p(0)-eps, p(1)); 00425 prevj = 0; 00426 break; 00427 } 00428 00429 // d^2()/dxideta 00430 case 1: 00431 { 00432 pp = Point(p(0), p(1)+eps); 00433 pm = Point(p(0), p(1)-eps); 00434 prevj = 0; 00435 break; 00436 } 00437 00438 // d^2()/deta^2 00439 case 2: 00440 { 00441 pp = Point(p(0), p(1)+eps); 00442 pm = Point(p(0), p(1)-eps); 00443 prevj = 1; 00444 break; 00445 } 00446 default: 00447 libmesh_error(); 00448 } 00449 return (FE<2,L2_HIERARCHIC>::shape_deriv(elem, order, i, prevj, pp) - 00450 FE<2,L2_HIERARCHIC>::shape_deriv(elem, order, i, prevj, pm) 00451 )/2./eps; 00452 } 00453 00454 } // namespace libMesh
Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:46 UTC
Hosted By: