fe_nedelec_one_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++ inlcludes 00020 00021 // Local includes 00022 #include "libmesh/fe.h" 00023 #include "libmesh/elem.h" 00024 00025 namespace libMesh 00026 { 00027 00028 template <> 00029 RealGradient FE<2,NEDELEC_ONE>::shape(const ElemType, 00030 const Order, 00031 const unsigned int, 00032 const Point&) 00033 { 00034 #if LIBMESH_DIM > 1 00035 libMesh::err << "Nedelec elements require the element type\n" 00036 << "because edge orientation is needed." 00037 << std::endl; 00038 libmesh_error(); 00039 #endif // LIBMESH_DIM > 1 00040 00041 libmesh_error(); 00042 return RealGradient(); 00043 } 00044 00045 00046 // An excellent discussion of Nedelec shape functions is given in 00047 // http://www.dealii.org/developer/reports/nedelec/nedelec.pdf 00048 template <> 00049 RealGradient FE<2,NEDELEC_ONE>::shape(const Elem* elem, 00050 const Order order, 00051 const unsigned int i, 00052 const Point& p) 00053 { 00054 #if LIBMESH_DIM > 1 00055 libmesh_assert(elem); 00056 00057 const Order total_order = static_cast<Order>(order + elem->p_level()); 00058 00059 switch (total_order) 00060 { 00061 case FIRST: 00062 { 00063 switch (elem->type()) 00064 { 00065 case QUAD8: 00066 case QUAD9: 00067 { 00068 libmesh_assert_less (i, 4); 00069 00070 const Real xi = p(0); 00071 const Real eta = p(1); 00072 00073 libmesh_assert_less_equal ( std::fabs(xi), 1.0 ); 00074 libmesh_assert_less_equal ( std::fabs(eta), 1.0 ); 00075 00076 switch(i) 00077 { 00078 case 0: 00079 { 00080 if( elem->point(0) > elem->point(1) ) 00081 return RealGradient( -0.25*(1.0-eta), 0.0 ); 00082 else 00083 return RealGradient( 0.25*(1.0-eta), 0.0 ); 00084 } 00085 case 1: 00086 { 00087 if( elem->point(1) > elem->point(2) ) 00088 return RealGradient( 0.0, -0.25*(1.0+xi) ); 00089 else 00090 return RealGradient( 0.0, 0.25*(1.0+xi) ); 00091 } 00092 00093 case 2: 00094 { 00095 if( elem->point(2) > elem->point(3) ) 00096 return RealGradient( 0.25*(1.0+eta), 0.0 ); 00097 else 00098 return RealGradient( -0.25*(1.0+eta), 0.0 ); 00099 } 00100 case 3: 00101 { 00102 if( elem->point(3) > elem->point(0) ) 00103 return RealGradient( 0.0, -0.25*(xi-1.0) ); 00104 else 00105 return RealGradient( 0.0, 0.25*(xi-1.0) ); 00106 } 00107 00108 default: 00109 libmesh_error(); 00110 00111 } 00112 00113 return RealGradient(); 00114 } 00115 00116 case TRI6: 00117 { 00118 const Real xi = p(0); 00119 const Real eta = p(1); 00120 00121 libmesh_assert_less (i, 3); 00122 00123 switch(i) 00124 { 00125 case 0: 00126 { 00127 if( elem->point(0) > elem->point(1) ) 00128 return RealGradient( -1.0+eta, -xi ); 00129 else 00130 return RealGradient( 1.0-eta, xi ); 00131 } 00132 case 1: 00133 { 00134 if( elem->point(1) > elem->point(2) ) 00135 return RealGradient( eta, -xi ); 00136 else 00137 return RealGradient( -eta, xi ); 00138 } 00139 00140 case 2: 00141 { 00142 if( elem->point(2) > elem->point(0) ) 00143 return RealGradient( eta, -xi+1.0 ); 00144 else 00145 return RealGradient( -eta, xi-1.0 ); 00146 } 00147 00148 default: 00149 libmesh_error(); 00150 00151 } 00152 } 00153 00154 default: 00155 { 00156 libMesh::err << "ERROR: Unsupported 2D element type!: " << elem->type() 00157 << std::endl; 00158 libmesh_error(); 00159 } 00160 } 00161 } 00162 00163 // unsupported order 00164 default: 00165 { 00166 libMesh::err << "ERROR: Unsupported 2D FE order!: " << total_order 00167 << std::endl; 00168 libmesh_error(); 00169 } 00170 } 00171 #endif // LIBMESH_DIM > 1 00172 00173 libmesh_error(); 00174 return RealGradient(); 00175 } 00176 00177 00178 00179 template <> 00180 RealGradient FE<2,NEDELEC_ONE>::shape_deriv(const ElemType, 00181 const Order, 00182 const unsigned int, 00183 const unsigned int, 00184 const Point&) 00185 { 00186 #if LIBMESH_DIM > 1 00187 libMesh::err << "Nedelec elements require the element type\n" 00188 << "because edge orientation is needed." 00189 << std::endl; 00190 libmesh_error(); 00191 #endif // LIBMESH_DIM > 1 00192 00193 libmesh_error(); 00194 return RealGradient(); 00195 } 00196 00197 00198 00199 template <> 00200 RealGradient FE<2,NEDELEC_ONE>::shape_deriv(const Elem* elem, 00201 const Order order, 00202 const unsigned int i, 00203 const unsigned int j, 00204 const Point&) 00205 { 00206 #if LIBMESH_DIM > 1 00207 libmesh_assert(elem); 00208 libmesh_assert_less (j, 2); 00209 00210 const Order total_order = static_cast<Order>(order + elem->p_level()); 00211 00212 switch (total_order) 00213 { 00214 // linear Lagrange shape functions 00215 case FIRST: 00216 { 00217 switch (elem->type()) 00218 { 00219 case QUAD8: 00220 case QUAD9: 00221 { 00222 libmesh_assert_less (i, 4); 00223 00224 switch (j) 00225 { 00226 // d()/dxi 00227 case 0: 00228 { 00229 switch(i) 00230 { 00231 case 0: 00232 case 2: 00233 return RealGradient(); 00234 case 1: 00235 { 00236 if( elem->point(1) > elem->point(2) ) 00237 return RealGradient( 0.0, -0.25 ); 00238 else 00239 return RealGradient( 0.0, 0.25 ); 00240 } 00241 case 3: 00242 { 00243 if( elem->point(3) > elem->point(0) ) 00244 return RealGradient( 0.0, -0.25 ); 00245 else 00246 return RealGradient( 0.0, 0.25 ); 00247 } 00248 default: 00249 libmesh_error(); 00250 } 00251 } // j=0 00252 00253 // d()/deta 00254 case 1: 00255 { 00256 switch(i) 00257 { 00258 case 1: 00259 case 3: 00260 return RealGradient(); 00261 case 0: 00262 { 00263 if( elem->point(0) > elem->point(1) ) 00264 return RealGradient( 0.25 ); 00265 else 00266 return RealGradient( -0.25 ); 00267 } 00268 case 2: 00269 { 00270 if( elem->point(2) > elem->point(3) ) 00271 return RealGradient( 0.25 ); 00272 else 00273 return RealGradient( -0.25 ); 00274 } 00275 default: 00276 libmesh_error(); 00277 } 00278 } // j=1 00279 00280 default: 00281 libmesh_error(); 00282 } 00283 00284 return RealGradient(); 00285 } 00286 00287 case TRI6: 00288 { 00289 libmesh_assert_less (i, 3); 00290 00291 // Account for edge flipping 00292 Real f = 1.0; 00293 00294 switch(i) 00295 { 00296 case 0: 00297 { 00298 if( elem->point(0) > elem->point(1) ) 00299 f = -1.0; 00300 break; 00301 } 00302 case 1: 00303 { 00304 if( elem->point(1) > elem->point(2) ) 00305 f = -1.0; 00306 break; 00307 } 00308 case 2: 00309 { 00310 if( elem->point(2) > elem->point(0) ) 00311 f = -1.0; 00312 break; 00313 } 00314 default: 00315 libmesh_error(); 00316 } 00317 00318 switch (j) 00319 { 00320 // d()/dxi 00321 case 0: 00322 { 00323 return RealGradient( 0.0, f*1.0); 00324 } 00325 // d()/deta 00326 case 1: 00327 { 00328 return RealGradient( f*(-1.0) ); 00329 } 00330 default: 00331 libmesh_error(); 00332 } 00333 } 00334 00335 default: 00336 { 00337 libMesh::err << "ERROR: Unsupported 2D element type!: " << elem->type() 00338 << std::endl; 00339 libmesh_error(); 00340 } 00341 } 00342 } 00343 // unsupported order 00344 default: 00345 { 00346 libMesh::err << "ERROR: Unsupported 2D FE order!: " << total_order 00347 << std::endl; 00348 libmesh_error(); 00349 } 00350 } 00351 #endif // LIBMESH_DIM > 1 00352 00353 libmesh_error(); 00354 return RealGradient(); 00355 } 00356 00357 00358 00359 00360 template <> 00361 RealGradient FE<2,NEDELEC_ONE>::shape_second_deriv(const ElemType, 00362 const Order, 00363 const unsigned int, 00364 const unsigned int, 00365 const Point&) 00366 { 00367 #if LIBMESH_DIM > 1 00368 libMesh::err << "Nedelec elements require the element type\n" 00369 << "because edge orientation is needed." 00370 << std::endl; 00371 libmesh_error(); 00372 #endif // LIBMESH_DIM > 1 00373 00374 libmesh_error(); 00375 return RealGradient(); 00376 } 00377 00378 00379 00380 template <> 00381 RealGradient FE<2,NEDELEC_ONE>::shape_second_deriv(const Elem* elem, 00382 const Order order, 00383 const unsigned int i, 00384 const unsigned int j, 00385 const Point&) 00386 { 00387 #if LIBMESH_DIM > 1 00388 libmesh_assert(elem); 00389 00390 // j = 0 ==> d^2 phi / dxi^2 00391 // j = 1 ==> d^2 phi / dxi deta 00392 // j = 2 ==> d^2 phi / deta^2 00393 libmesh_assert_less (j, 3); 00394 00395 const Order total_order = static_cast<Order>(order + elem->p_level()); 00396 00397 switch (total_order) 00398 { 00399 // linear Lagrange shape functions 00400 case FIRST: 00401 { 00402 switch (elem->type()) 00403 { 00404 case QUAD8: 00405 case QUAD9: 00406 { 00407 libmesh_assert_less (i, 4); 00408 // All second derivatives for linear quads are zero. 00409 return RealGradient(); 00410 } 00411 00412 case TRI6: 00413 { 00414 libmesh_assert_less (i, 3); 00415 // All second derivatives for linear triangles are zero. 00416 return RealGradient(); 00417 } 00418 00419 default: 00420 { 00421 libMesh::err << "ERROR: Unsupported 2D element type!: " << elem->type() 00422 << std::endl; 00423 libmesh_error(); 00424 } 00425 00426 } // end switch (type) 00427 } // end case FIRST 00428 00429 // unsupported order 00430 default: 00431 { 00432 libMesh::err << "ERROR: Unsupported 2D FE order!: " << total_order 00433 << std::endl; 00434 libmesh_error(); 00435 } 00436 00437 } // end switch (order) 00438 00439 #endif // LIBMESH_DIM > 1 00440 00441 libmesh_error(); 00442 return RealGradient(); 00443 } 00444 00445 } // namespace libMesh
Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:46 UTC
Hosted By: