fe_interface.C
Go to the documentation of this file.
00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 00003 00004 // This library is free software; you can redistribute it and/or 00005 // modify it under the terms of the GNU Lesser General Public 00006 // License as published by the Free Software Foundation; either 00007 // version 2.1 of the License, or (at your option) any later version. 00008 00009 // This library is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 // Lesser General Public License for more details. 00013 00014 // You should have received a copy of the GNU Lesser General Public 00015 // License along with this library; if not, write to the Free Software 00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00017 00018 00019 00020 // Local includes 00021 #include "libmesh/fe_interface.h" 00022 #include "libmesh/elem.h" 00023 #include "libmesh/fe.h" 00024 #include "libmesh/fe_compute_data.h" 00025 #include "libmesh/dof_map.h" 00026 00027 namespace libMesh 00028 { 00029 00030 //------------------------------------------------------------ 00031 //FEInterface class members 00032 FEInterface::FEInterface() 00033 { 00034 libmesh_error_msg("ERROR: Do not define an object of this type."); 00035 } 00036 00037 00038 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES 00039 #define fe_family_switch(dim, func_and_args, prefix, suffix) \ 00040 do { \ 00041 switch (fe_t.family) \ 00042 { \ 00043 case CLOUGH: \ 00044 prefix FE<dim,CLOUGH>::func_and_args suffix \ 00045 case HERMITE: \ 00046 prefix FE<dim,HERMITE>::func_and_args suffix \ 00047 case HIERARCHIC: \ 00048 prefix FE<dim,HIERARCHIC>::func_and_args suffix \ 00049 case L2_HIERARCHIC: \ 00050 prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \ 00051 case LAGRANGE: \ 00052 prefix FE<dim,LAGRANGE>::func_and_args suffix \ 00053 case L2_LAGRANGE: \ 00054 prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \ 00055 case MONOMIAL: \ 00056 prefix FE<dim,MONOMIAL>::func_and_args suffix \ 00057 case SCALAR: \ 00058 prefix FE<dim,SCALAR>::func_and_args suffix \ 00059 case BERNSTEIN: \ 00060 prefix FE<dim,BERNSTEIN>::func_and_args suffix \ 00061 case SZABAB: \ 00062 prefix FE<dim,SZABAB>::func_and_args suffix \ 00063 case XYZ: \ 00064 prefix FEXYZ<dim>::func_and_args suffix \ 00065 case SUBDIVISION: \ 00066 libmesh_assert_equal_to (dim, 2); \ 00067 prefix FE<2,SUBDIVISION>::func_and_args suffix \ 00068 default: \ 00069 libmesh_error_msg("Unsupported family = " << fe_t.family); \ 00070 } \ 00071 } while (0) 00072 00073 #define fe_family_with_vec_switch(dim, func_and_args, prefix, suffix) \ 00074 do { \ 00075 switch (fe_t.family) \ 00076 { \ 00077 case CLOUGH: \ 00078 prefix FE<dim,CLOUGH>::func_and_args suffix \ 00079 case HERMITE: \ 00080 prefix FE<dim,HERMITE>::func_and_args suffix \ 00081 case HIERARCHIC: \ 00082 prefix FE<dim,HIERARCHIC>::func_and_args suffix \ 00083 case L2_HIERARCHIC: \ 00084 prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \ 00085 case LAGRANGE: \ 00086 prefix FE<dim,LAGRANGE>::func_and_args suffix \ 00087 case LAGRANGE_VEC: \ 00088 prefix FELagrangeVec<dim>::func_and_args suffix \ 00089 case L2_LAGRANGE: \ 00090 prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \ 00091 case MONOMIAL: \ 00092 prefix FE<dim,MONOMIAL>::func_and_args suffix \ 00093 case SCALAR: \ 00094 prefix FE<dim,SCALAR>::func_and_args suffix \ 00095 case BERNSTEIN: \ 00096 prefix FE<dim,BERNSTEIN>::func_and_args suffix \ 00097 case SZABAB: \ 00098 prefix FE<dim,SZABAB>::func_and_args suffix \ 00099 case XYZ: \ 00100 prefix FEXYZ<dim>::func_and_args suffix \ 00101 case SUBDIVISION: \ 00102 libmesh_assert_equal_to (dim, 2); \ 00103 prefix FE<2,SUBDIVISION>::func_and_args suffix \ 00104 case NEDELEC_ONE: \ 00105 prefix FENedelecOne<dim>::func_and_args suffix \ 00106 default: \ 00107 libmesh_error_msg("Unsupported family = " << fe_t.family); \ 00108 } \ 00109 } while (0) 00110 00111 #define fe_scalar_vec_error_switch(dim, func_and_args, prefix, suffix) \ 00112 do { \ 00113 switch (fe_t.family) \ 00114 { \ 00115 case CLOUGH: \ 00116 prefix FE<dim,CLOUGH>::func_and_args suffix \ 00117 case HERMITE: \ 00118 prefix FE<dim,HERMITE>::func_and_args suffix \ 00119 case HIERARCHIC: \ 00120 prefix FE<dim,HIERARCHIC>::func_and_args suffix \ 00121 case L2_HIERARCHIC: \ 00122 prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \ 00123 case LAGRANGE: \ 00124 prefix FE<dim,LAGRANGE>::func_and_args suffix \ 00125 case L2_LAGRANGE: \ 00126 prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \ 00127 case MONOMIAL: \ 00128 prefix FE<dim,MONOMIAL>::func_and_args suffix \ 00129 case SCALAR: \ 00130 prefix FE<dim,SCALAR>::func_and_args suffix \ 00131 case BERNSTEIN: \ 00132 prefix FE<dim,BERNSTEIN>::func_and_args suffix \ 00133 case SZABAB: \ 00134 prefix FE<dim,SZABAB>::func_and_args suffix \ 00135 case XYZ: \ 00136 prefix FEXYZ<dim>::func_and_args suffix \ 00137 case SUBDIVISION: \ 00138 libmesh_assert_equal_to (dim, 2); \ 00139 prefix FE<2,SUBDIVISION>::func_and_args suffix \ 00140 case LAGRANGE_VEC: \ 00141 case NEDELEC_ONE: \ 00142 libmesh_error_msg("Error: Can only request scalar valued elements for Real FEInterface::func_and_args"); \ 00143 default: \ 00144 libmesh_error_msg("Unsupported family = " << fe_t.family); \ 00145 } \ 00146 } while(0) 00147 00148 00149 #define fe_vector_scalar_error_switch(dim, func_and_args, prefix, suffix) \ 00150 do { \ 00151 switch (fe_t.family) \ 00152 { \ 00153 case LAGRANGE_VEC: \ 00154 prefix FELagrangeVec<dim>::func_and_args suffix \ 00155 case NEDELEC_ONE: \ 00156 prefix FENedelecOne<dim>::func_and_args suffix \ 00157 case HERMITE: \ 00158 case HIERARCHIC: \ 00159 case L2_HIERARCHIC: \ 00160 case LAGRANGE: \ 00161 case L2_LAGRANGE: \ 00162 case MONOMIAL: \ 00163 case SCALAR: \ 00164 case BERNSTEIN: \ 00165 case SZABAB: \ 00166 case XYZ: \ 00167 case SUBDIVISION: \ 00168 libmesh_error_msg("Error: Can only request vector valued elements for RealGradient FEInterface::shape"); \ 00169 default: \ 00170 libmesh_error_msg("Unsupported family = " << fe_t.family); \ 00171 } \ 00172 } while(0) 00173 00174 #else 00175 #define fe_family_switch(dim, func_and_args, prefix, suffix) \ 00176 do { \ 00177 switch (fe_t.family) \ 00178 { \ 00179 case CLOUGH: \ 00180 prefix FE<dim,CLOUGH>::func_and_args suffix \ 00181 case HERMITE: \ 00182 prefix FE<dim,HERMITE>::func_and_args suffix \ 00183 case HIERARCHIC: \ 00184 prefix FE<dim,HIERARCHIC>::func_and_args suffix \ 00185 case L2_HIERARCHIC: \ 00186 prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \ 00187 case LAGRANGE: \ 00188 prefix FE<dim,LAGRANGE>::func_and_args suffix \ 00189 case L2_LAGRANGE: \ 00190 prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \ 00191 case MONOMIAL: \ 00192 prefix FE<dim,MONOMIAL>::func_and_args suffix \ 00193 case SCALAR: \ 00194 prefix FE<dim,SCALAR>::func_and_args suffix \ 00195 case XYZ: \ 00196 prefix FEXYZ<dim>::func_and_args suffix \ 00197 case SUBDIVISION: \ 00198 libmesh_assert_equal_to (dim, 2); \ 00199 prefix FE<2,SUBDIVISION>::func_and_args suffix \ 00200 default: \ 00201 libmesh_error_msg("Unsupported family = " << fe_t.family); \ 00202 } \ 00203 } while (0) 00204 00205 #define fe_family_with_vec_switch(dim, func_and_args, prefix, suffix) \ 00206 do { \ 00207 switch (fe_t.family) \ 00208 { \ 00209 case CLOUGH: \ 00210 prefix FE<dim,CLOUGH>::func_and_args suffix \ 00211 case HERMITE: \ 00212 prefix FE<dim,HERMITE>::func_and_args suffix \ 00213 case HIERARCHIC: \ 00214 prefix FE<dim,HIERARCHIC>::func_and_args suffix \ 00215 case L2_HIERARCHIC: \ 00216 prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \ 00217 case LAGRANGE: \ 00218 prefix FE<dim,LAGRANGE>::func_and_args suffix \ 00219 case LAGRANGE_VEC: \ 00220 prefix FELagrangeVec<dim>::func_and_args suffix \ 00221 case L2_LAGRANGE: \ 00222 prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \ 00223 case MONOMIAL: \ 00224 prefix FE<dim,MONOMIAL>::func_and_args suffix \ 00225 case SCALAR: \ 00226 prefix FE<dim,SCALAR>::func_and_args suffix \ 00227 case XYZ: \ 00228 prefix FEXYZ<dim>::func_and_args suffix \ 00229 case SUBDIVISION: \ 00230 libmesh_assert_equal_to (dim, 2); \ 00231 prefix FE<2,SUBDIVISION>::func_and_args suffix \ 00232 case NEDELEC_ONE: \ 00233 prefix FENedelecOne<dim>::func_and_args suffix \ 00234 default: \ 00235 libmesh_error_msg("Unsupported family = " << fe_t.family); \ 00236 } \ 00237 } while (0) 00238 00239 #define fe_scalar_vec_error_switch(dim, func_and_args, prefix, suffix) \ 00240 do { \ 00241 switch (fe_t.family) \ 00242 { \ 00243 case CLOUGH: \ 00244 prefix FE<dim,CLOUGH>::func_and_args suffix \ 00245 case HERMITE: \ 00246 prefix FE<dim,HERMITE>::func_and_args suffix \ 00247 case HIERARCHIC: \ 00248 prefix FE<dim,HIERARCHIC>::func_and_args suffix \ 00249 case L2_HIERARCHIC: \ 00250 prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \ 00251 case LAGRANGE: \ 00252 prefix FE<dim,LAGRANGE>::func_and_args suffix \ 00253 case L2_LAGRANGE: \ 00254 prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \ 00255 case MONOMIAL: \ 00256 prefix FE<dim,MONOMIAL>::func_and_args suffix \ 00257 case SCALAR: \ 00258 prefix FE<dim,SCALAR>::func_and_args suffix \ 00259 case XYZ: \ 00260 prefix FEXYZ<dim>::func_and_args suffix \ 00261 case SUBDIVISION: \ 00262 libmesh_assert_equal_to (dim, 2); \ 00263 prefix FE<2,SUBDIVISION>::func_and_args suffix \ 00264 case LAGRANGE_VEC: \ 00265 case NEDELEC_ONE: \ 00266 libmesh_error_msg("Error: Can only request scalar valued elements for Real FEInterface::func_and_args"); \ 00267 default: \ 00268 libmesh_error_msg("Unsupported family = " << fe_t.family); \ 00269 } \ 00270 } while(0) 00271 00272 00273 #define fe_vector_scalar_error_switch(dim, func_and_args, prefix, suffix) \ 00274 do { \ 00275 switch (fe_t.family) \ 00276 { \ 00277 case LAGRANGE_VEC: \ 00278 prefix FELagrangeVec<dim>::func_and_args suffix \ 00279 case NEDELEC_ONE: \ 00280 prefix FENedelecOne<dim>::func_and_args suffix \ 00281 case HERMITE: \ 00282 case HIERARCHIC: \ 00283 case L2_HIERARCHIC: \ 00284 case LAGRANGE: \ 00285 case L2_LAGRANGE: \ 00286 case MONOMIAL: \ 00287 case SCALAR: \ 00288 case XYZ: \ 00289 case SUBDIVISION: \ 00290 libmesh_error_msg("Error: Can only request vector valued elements for RealGradient FEInterface::func_and_args"); \ 00291 default: \ 00292 libmesh_error_msg("Unsupported family = " << fe_t.family); \ 00293 } \ 00294 } while(0) 00295 #endif 00296 00297 00298 #define fe_switch(func_and_args) \ 00299 do { \ 00300 switch (dim) \ 00301 { \ 00302 /* 0D */ \ 00303 case 0: \ 00304 fe_family_switch (0, func_and_args, return, ;); \ 00305 /* 1D */ \ 00306 case 1: \ 00307 fe_family_switch (1, func_and_args, return, ;); \ 00308 /* 2D */ \ 00309 case 2: \ 00310 fe_family_switch (2, func_and_args, return, ;); \ 00311 /* 3D */ \ 00312 case 3: \ 00313 fe_family_switch (3, func_and_args, return, ;); \ 00314 default: \ 00315 libmesh_error_msg("Invalid dim = " << dim); \ 00316 } \ 00317 } while (0) 00318 00319 #define fe_with_vec_switch(func_and_args) \ 00320 do { \ 00321 switch (dim) \ 00322 { \ 00323 /* 0D */ \ 00324 case 0: \ 00325 fe_family_with_vec_switch (0, func_and_args, return, ;); \ 00326 /* 1D */ \ 00327 case 1: \ 00328 fe_family_with_vec_switch (1, func_and_args, return, ;); \ 00329 /* 2D */ \ 00330 case 2: \ 00331 fe_family_with_vec_switch (2, func_and_args, return, ;); \ 00332 /* 3D */ \ 00333 case 3: \ 00334 fe_family_with_vec_switch (3, func_and_args, return, ;); \ 00335 default: \ 00336 libmesh_error_msg("Invalid dim = " << dim); \ 00337 } \ 00338 } while (0) 00339 00340 00341 #define void_fe_switch(func_and_args) \ 00342 do { \ 00343 switch (dim) \ 00344 { \ 00345 /* 0D */ \ 00346 case 0: \ 00347 fe_family_switch (0, func_and_args, ;, ; return;); \ 00348 /* 1D */ \ 00349 case 1: \ 00350 fe_family_switch (1, func_and_args, ;, ; return;); \ 00351 /* 2D */ \ 00352 case 2: \ 00353 fe_family_switch (2, func_and_args, ;, ; return;); \ 00354 /* 3D */ \ 00355 case 3: \ 00356 fe_family_switch (3, func_and_args, ;, ; return;); \ 00357 default: \ 00358 libmesh_error_msg("Invalid dim = " << dim); \ 00359 } \ 00360 } while (0) 00361 00362 #define void_fe_with_vec_switch(func_and_args) \ 00363 do { \ 00364 switch (dim) \ 00365 { \ 00366 /* 0D */ \ 00367 case 0: \ 00368 fe_family_with_vec_switch (0, func_and_args, ;, ; return;); \ 00369 /* 1D */ \ 00370 case 1: \ 00371 fe_family_with_vec_switch (1, func_and_args, ;, ; return;); \ 00372 /* 2D */ \ 00373 case 2: \ 00374 fe_family_with_vec_switch (2, func_and_args, ;, ; return;); \ 00375 /* 3D */ \ 00376 case 3: \ 00377 fe_family_with_vec_switch (3, func_and_args, ;, ; return;); \ 00378 default: \ 00379 libmesh_error_msg("Invalid dim = " << dim); \ 00380 } \ 00381 } while (0) 00382 00383 00384 00385 unsigned int FEInterface::n_shape_functions(const unsigned int dim, 00386 const FEType& fe_t, 00387 const ElemType t) 00388 { 00389 00390 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 00391 /* 00392 * Since the FEType, stored in DofMap/(some System child), has to 00393 * be the _same_ for InfFE and FE, we have to catch calls 00394 * to infinite elements through the element type. 00395 */ 00396 00397 if ( is_InfFE_elem(t) ) 00398 return ifem_n_shape_functions(dim, fe_t, t); 00399 00400 #endif 00401 00402 const Order o = fe_t.order; 00403 00404 fe_with_vec_switch(n_shape_functions(t, o)); 00405 00406 libmesh_error_msg("We'll never get here!"); 00407 return 0; 00408 } 00409 00410 00411 00412 00413 00414 unsigned int FEInterface::n_dofs(const unsigned int dim, 00415 const FEType& fe_t, 00416 const ElemType t) 00417 { 00418 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 00419 00420 if ( is_InfFE_elem(t) ) 00421 return ifem_n_dofs(dim, fe_t, t); 00422 00423 #endif 00424 00425 const Order o = fe_t.order; 00426 00427 fe_with_vec_switch(n_dofs(t, o)); 00428 00429 libmesh_error_msg("We'll never get here!"); 00430 return 0; 00431 } 00432 00433 00434 00435 00436 unsigned int FEInterface::n_dofs_at_node(const unsigned int dim, 00437 const FEType& fe_t, 00438 const ElemType t, 00439 const unsigned int n) 00440 { 00441 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 00442 00443 if ( is_InfFE_elem(t) ) 00444 return ifem_n_dofs_at_node(dim, fe_t, t, n); 00445 00446 #endif 00447 00448 const Order o = fe_t.order; 00449 00450 fe_with_vec_switch(n_dofs_at_node(t, o, n)); 00451 00452 libmesh_error_msg("We'll never get here!"); 00453 return 0; 00454 } 00455 00456 00457 00458 00459 00460 unsigned int FEInterface::n_dofs_per_elem(const unsigned int dim, 00461 const FEType& fe_t, 00462 const ElemType t) 00463 { 00464 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 00465 00466 if ( is_InfFE_elem(t) ) 00467 return ifem_n_dofs_per_elem(dim, fe_t, t); 00468 00469 #endif 00470 00471 const Order o = fe_t.order; 00472 00473 fe_with_vec_switch(n_dofs_per_elem(t, o)); 00474 00475 libmesh_error_msg("We'll never get here!"); 00476 return 0; 00477 } 00478 00479 00480 00481 00482 void FEInterface::dofs_on_side(const Elem* const elem, 00483 const unsigned int dim, 00484 const FEType& fe_t, 00485 unsigned int s, 00486 std::vector<unsigned int>& di) 00487 { 00488 const Order o = fe_t.order; 00489 00490 void_fe_with_vec_switch(dofs_on_side(elem, o, s, di)); 00491 00492 libmesh_error_msg("We'll never get here!"); 00493 } 00494 00495 00496 00497 void FEInterface::dofs_on_edge(const Elem* const elem, 00498 const unsigned int dim, 00499 const FEType& fe_t, 00500 unsigned int e, 00501 std::vector<unsigned int>& di) 00502 { 00503 const Order o = fe_t.order; 00504 00505 void_fe_with_vec_switch(dofs_on_edge(elem, o, e, di)); 00506 00507 libmesh_error_msg("We'll never get here!"); 00508 } 00509 00510 00511 00512 00513 void FEInterface::nodal_soln(const unsigned int dim, 00514 const FEType& fe_t, 00515 const Elem* elem, 00516 const std::vector<Number>& elem_soln, 00517 std::vector<Number>& nodal_soln) 00518 { 00519 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 00520 00521 if ( is_InfFE_elem(elem->type()) ) 00522 { 00523 ifem_nodal_soln(dim, fe_t, elem, elem_soln, nodal_soln); 00524 return; 00525 } 00526 00527 #endif 00528 00529 const Order order = fe_t.order; 00530 00531 void_fe_with_vec_switch(nodal_soln(elem, order, elem_soln, nodal_soln)); 00532 } 00533 00534 00535 00536 00537 Point FEInterface::map(unsigned int dim, 00538 const FEType& fe_t, 00539 const Elem* elem, 00540 const Point& p) 00541 { 00542 fe_with_vec_switch(map(elem, p)); 00543 00544 libmesh_error_msg("We'll never get here!"); 00545 return Point(); 00546 } 00547 00548 00549 00550 00551 00552 Point FEInterface::inverse_map (const unsigned int dim, 00553 const FEType& fe_t, 00554 const Elem* elem, 00555 const Point& p, 00556 const Real tolerance, 00557 const bool secure) 00558 { 00559 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 00560 00561 if ( is_InfFE_elem(elem->type()) ) 00562 return ifem_inverse_map(dim, fe_t, elem, p,tolerance, secure); 00563 00564 #endif 00565 00566 fe_with_vec_switch(inverse_map(elem, p, tolerance, secure)); 00567 00568 libmesh_error_msg("We'll never get here!"); 00569 return Point(); 00570 } 00571 00572 00573 00574 00575 void FEInterface::inverse_map (const unsigned int dim, 00576 const FEType& fe_t, 00577 const Elem* elem, 00578 const std::vector<Point>& physical_points, 00579 std::vector<Point>& reference_points, 00580 const Real tolerance, 00581 const bool secure) 00582 { 00583 const std::size_t n_pts = physical_points.size(); 00584 00585 // Resize the vector 00586 reference_points.resize(n_pts); 00587 00588 if (n_pts == 0) 00589 { 00590 libMesh::err << "WARNING: empty vector physical_points!" 00591 << std::endl; 00592 libmesh_here(); 00593 return; 00594 } 00595 00596 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 00597 00598 if ( is_InfFE_elem(elem->type()) ) 00599 { 00600 ifem_inverse_map(dim, fe_t, elem, physical_points, reference_points, tolerance, secure); 00601 return; 00602 // libmesh_not_implemented(); 00603 } 00604 00605 #endif 00606 00607 void_fe_with_vec_switch(inverse_map(elem, physical_points, reference_points, tolerance, secure)); 00608 00609 libmesh_error_msg("We'll never get here!"); 00610 } 00611 00612 00613 00614 bool FEInterface::on_reference_element(const Point& p, 00615 const ElemType t, 00616 const Real eps) 00617 { 00618 return FEBase::on_reference_element(p,t,eps); 00619 } 00620 00621 00622 00623 00624 Real FEInterface::shape(const unsigned int dim, 00625 const FEType& fe_t, 00626 const ElemType t, 00627 const unsigned int i, 00628 const Point& p) 00629 { 00630 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 00631 00632 if ( is_InfFE_elem(t) ) 00633 return ifem_shape(dim, fe_t, t, i, p); 00634 00635 #endif 00636 00637 const Order o = fe_t.order; 00638 00639 fe_switch(shape(t,o,i,p)); 00640 00641 libmesh_error_msg("We'll never get here!"); 00642 return 0.; 00643 } 00644 00645 Real FEInterface::shape(const unsigned int dim, 00646 const FEType& fe_t, 00647 const Elem* elem, 00648 const unsigned int i, 00649 const Point& p) 00650 { 00651 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 00652 00653 if ( elem && is_InfFE_elem(elem->type()) ) 00654 return ifem_shape(dim, fe_t, elem, i, p); 00655 00656 #endif 00657 00658 const Order o = fe_t.order; 00659 00660 fe_switch(shape(elem,o,i,p)); 00661 00662 libmesh_error_msg("We'll never get here!"); 00663 return 0.; 00664 } 00665 00666 template<> 00667 void FEInterface::shape<Real>(const unsigned int dim, 00668 const FEType& fe_t, 00669 const ElemType t, 00670 const unsigned int i, 00671 const Point& p, 00672 Real& phi) 00673 { 00674 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 00675 00676 if ( is_InfFE_elem(t) ) 00677 phi = ifem_shape(dim, fe_t, t, i, p); 00678 00679 #endif 00680 00681 const Order o = fe_t.order; 00682 00683 switch(dim) 00684 { 00685 case 0: 00686 fe_scalar_vec_error_switch(0, shape(t,o,i,p), phi = , ; break;); 00687 break; 00688 case 1: 00689 fe_scalar_vec_error_switch(1, shape(t,o,i,p), phi = , ; break;); 00690 break; 00691 case 2: 00692 fe_scalar_vec_error_switch(2, shape(t,o,i,p), phi = , ; break;); 00693 break; 00694 case 3: 00695 fe_scalar_vec_error_switch(3, shape(t,o,i,p), phi = , ; break;); 00696 break; 00697 default: 00698 libmesh_error_msg("Invalid dimension = " << dim); 00699 } 00700 00701 return; 00702 } 00703 00704 template<> 00705 void FEInterface::shape<Real>(const unsigned int dim, 00706 const FEType& fe_t, 00707 const Elem* elem, 00708 const unsigned int i, 00709 const Point& p, 00710 Real& phi) 00711 { 00712 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 00713 00714 if ( elem && is_InfFE_elem(elem->type()) ) 00715 phi = ifem_shape(dim, fe_t, elem, i, p); 00716 00717 #endif 00718 00719 const Order o = fe_t.order; 00720 00721 switch(dim) 00722 { 00723 case 0: 00724 fe_scalar_vec_error_switch(0, shape(elem,o,i,p), phi = , ; break;); 00725 break; 00726 case 1: 00727 fe_scalar_vec_error_switch(1, shape(elem,o,i,p), phi = , ; break;); 00728 break; 00729 case 2: 00730 fe_scalar_vec_error_switch(2, shape(elem,o,i,p), phi = , ; break;); 00731 break; 00732 case 3: 00733 fe_scalar_vec_error_switch(3, shape(elem,o,i,p), phi = , ; break;); 00734 break; 00735 default: 00736 libmesh_error_msg("Invalid dimension = " << dim); 00737 } 00738 00739 return; 00740 } 00741 00742 template<> 00743 void FEInterface::shape<RealGradient>(const unsigned int dim, 00744 const FEType& fe_t, 00745 const ElemType t, 00746 const unsigned int i, 00747 const Point& p, 00748 RealGradient& phi) 00749 { 00750 const Order o = fe_t.order; 00751 00752 switch(dim) 00753 { 00754 case 0: 00755 fe_vector_scalar_error_switch(0, shape(t,o,i,p), phi = , ; break;); 00756 break; 00757 case 1: 00758 fe_vector_scalar_error_switch(1, shape(t,o,i,p), phi = , ; break;); 00759 break; 00760 case 2: 00761 fe_vector_scalar_error_switch(2, shape(t,o,i,p), phi = , ; break;); 00762 break; 00763 case 3: 00764 fe_vector_scalar_error_switch(3, shape(t,o,i,p), phi = , ; break;); 00765 break; 00766 default: 00767 libmesh_error_msg("Invalid dimension = " << dim); 00768 } 00769 00770 return; 00771 } 00772 00773 template<> 00774 void FEInterface::shape<RealGradient>(const unsigned int dim, 00775 const FEType& fe_t, 00776 const Elem* elem, 00777 const unsigned int i, 00778 const Point& p, 00779 RealGradient& phi) 00780 { 00781 const Order o = fe_t.order; 00782 00783 switch(dim) 00784 { 00785 case 0: 00786 fe_vector_scalar_error_switch(0, shape(elem,o,i,p), phi = , ; break;); 00787 break; 00788 case 1: 00789 fe_vector_scalar_error_switch(1, shape(elem,o,i,p), phi = , ; break;); 00790 break; 00791 case 2: 00792 fe_vector_scalar_error_switch(2, shape(elem,o,i,p), phi = , ; break;); 00793 break; 00794 case 3: 00795 fe_vector_scalar_error_switch(3, shape(elem,o,i,p), phi = , ; break;); 00796 break; 00797 default: 00798 libmesh_error_msg("Invalid dimension = " << dim); 00799 } 00800 00801 return; 00802 } 00803 00804 void FEInterface::compute_data(const unsigned int dim, 00805 const FEType& fe_t, 00806 const Elem* elem, 00807 FEComputeData& data) 00808 { 00809 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 00810 00811 if ( elem && is_InfFE_elem(elem->type()) ) 00812 { 00813 data.init(); 00814 ifem_compute_data(dim, fe_t, elem, data); 00815 return; 00816 } 00817 00818 #endif 00819 00820 FEType p_refined = fe_t; 00821 p_refined.order = static_cast<Order>(p_refined.order + elem->p_level()); 00822 00823 const unsigned int n_dof = n_dofs (dim, p_refined, elem->type()); 00824 const Point& p = data.p; 00825 data.shape.resize(n_dof); 00826 00827 // set default values for all the output fields 00828 data.init(); 00829 00830 for (unsigned int n=0; n<n_dof; n++) 00831 data.shape[n] = shape(dim, p_refined, elem, n, p); 00832 00833 return; 00834 } 00835 00836 00837 00838 #ifdef LIBMESH_ENABLE_AMR 00839 00840 void FEInterface::compute_constraints (DofConstraints &constraints, 00841 DofMap &dof_map, 00842 const unsigned int variable_number, 00843 const Elem* elem) 00844 { 00845 libmesh_assert(elem); 00846 00847 const FEType& fe_t = dof_map.variable_type(variable_number); 00848 00849 switch (elem->dim()) 00850 { 00851 case 0: 00852 case 1: 00853 { 00854 // No constraints in 0D/1D. 00855 return; 00856 } 00857 00858 00859 case 2: 00860 { 00861 switch (fe_t.family) 00862 { 00863 case CLOUGH: 00864 FE<2,CLOUGH>::compute_constraints (constraints, 00865 dof_map, 00866 variable_number, 00867 elem); return; 00868 00869 case HERMITE: 00870 FE<2,HERMITE>::compute_constraints (constraints, 00871 dof_map, 00872 variable_number, 00873 elem); return; 00874 00875 case LAGRANGE: 00876 FE<2,LAGRANGE>::compute_constraints (constraints, 00877 dof_map, 00878 variable_number, 00879 elem); return; 00880 00881 case HIERARCHIC: 00882 FE<2,HIERARCHIC>::compute_constraints (constraints, 00883 dof_map, 00884 variable_number, 00885 elem); return; 00886 00887 case L2_HIERARCHIC: 00888 FE<2,L2_HIERARCHIC>::compute_constraints (constraints, 00889 dof_map, 00890 variable_number, 00891 elem); return; 00892 00893 case LAGRANGE_VEC: 00894 FE<2,LAGRANGE_VEC>::compute_constraints (constraints, 00895 dof_map, 00896 variable_number, 00897 elem); return; 00898 00899 00900 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES 00901 case SZABAB: 00902 FE<2,SZABAB>::compute_constraints (constraints, 00903 dof_map, 00904 variable_number, 00905 elem); return; 00906 00907 case BERNSTEIN: 00908 FE<2,BERNSTEIN>::compute_constraints (constraints, 00909 dof_map, 00910 variable_number, 00911 elem); return; 00912 00913 #endif 00914 default: 00915 return; 00916 } 00917 } 00918 00919 00920 case 3: 00921 { 00922 switch (fe_t.family) 00923 { 00924 case HERMITE: 00925 FE<3,HERMITE>::compute_constraints (constraints, 00926 dof_map, 00927 variable_number, 00928 elem); return; 00929 00930 case LAGRANGE: 00931 FE<3,LAGRANGE>::compute_constraints (constraints, 00932 dof_map, 00933 variable_number, 00934 elem); return; 00935 00936 case HIERARCHIC: 00937 FE<3,HIERARCHIC>::compute_constraints (constraints, 00938 dof_map, 00939 variable_number, 00940 elem); return; 00941 00942 case L2_HIERARCHIC: 00943 FE<3,L2_HIERARCHIC>::compute_constraints (constraints, 00944 dof_map, 00945 variable_number, 00946 elem); return; 00947 00948 case LAGRANGE_VEC: 00949 FE<3,LAGRANGE_VEC>::compute_constraints (constraints, 00950 dof_map, 00951 variable_number, 00952 elem); return; 00953 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES 00954 case SZABAB: 00955 FE<3,SZABAB>::compute_constraints (constraints, 00956 dof_map, 00957 variable_number, 00958 elem); return; 00959 00960 case BERNSTEIN: 00961 FE<3,BERNSTEIN>::compute_constraints (constraints, 00962 dof_map, 00963 variable_number, 00964 elem); return; 00965 00966 #endif 00967 default: 00968 return; 00969 } 00970 } 00971 00972 00973 default: 00974 libmesh_error_msg("Invalid dimension = " << elem->dim()); 00975 } 00976 } 00977 00978 #endif // #ifdef LIBMESH_ENABLE_AMR 00979 00980 00981 00982 #ifdef LIBMESH_ENABLE_PERIODIC 00983 00984 void FEInterface::compute_periodic_constraints (DofConstraints &constraints, 00985 DofMap &dof_map, 00986 const PeriodicBoundaries &boundaries, 00987 const MeshBase &mesh, 00988 const PointLocatorBase* point_locator, 00989 const unsigned int variable_number, 00990 const Elem* elem) 00991 { 00992 // No element-specific optimizations currently exist 00993 FEBase::compute_periodic_constraints (constraints, 00994 dof_map, 00995 boundaries, 00996 mesh, 00997 point_locator, 00998 variable_number, 00999 elem); 01000 } 01001 01002 #endif // #ifdef LIBMESH_ENABLE_PERIODIC 01003 01004 01005 01006 unsigned int FEInterface::max_order(const FEType& fe_t, 01007 const ElemType& el_t) 01008 { 01009 // Yeah, I know, infinity is much larger than 11, but our 01010 // solvers don't seem to like high degree polynomials, and our 01011 // quadrature rules and number_lookups tables 01012 // need to go up higher. 01013 const unsigned int unlimited = 11; 01014 01015 // If we used 0 as a default, then elements missing from this 01016 // table (e.g. infinite elements) would be considered broken. 01017 const unsigned int unknown = unlimited; 01018 01019 switch (fe_t.family) 01020 { 01021 case LAGRANGE: 01022 case L2_LAGRANGE: // TODO: L2_LAGRANGE can have higher "max_order" than LAGRANGE 01023 case LAGRANGE_VEC: 01024 switch (el_t) 01025 { 01026 case EDGE2: 01027 case EDGE3: 01028 case EDGE4: 01029 return 3; 01030 case TRI3: 01031 return 1; 01032 case TRI6: 01033 return 2; 01034 case QUAD4: 01035 return 1; 01036 case QUAD8: 01037 case QUAD9: 01038 return 2; 01039 case TET4: 01040 return 1; 01041 case TET10: 01042 return 2; 01043 case HEX8: 01044 return 1; 01045 case HEX20: 01046 case HEX27: 01047 return 2; 01048 case PRISM6: 01049 return 1; 01050 case PRISM15: 01051 case PRISM18: 01052 return 2; 01053 case PYRAMID5: 01054 return 1; 01055 case PYRAMID13: 01056 case PYRAMID14: 01057 return 2; 01058 default: 01059 return unknown; 01060 } 01061 break; 01062 case MONOMIAL: 01063 switch (el_t) 01064 { 01065 case EDGE2: 01066 case EDGE3: 01067 case EDGE4: 01068 case TRI3: 01069 case TRI6: 01070 case QUAD4: 01071 case QUAD8: 01072 case QUAD9: 01073 case TET4: 01074 case TET10: 01075 case HEX8: 01076 case HEX20: 01077 case HEX27: 01078 case PRISM6: 01079 case PRISM15: 01080 case PRISM18: 01081 case PYRAMID5: 01082 case PYRAMID13: 01083 case PYRAMID14: 01084 return unlimited; 01085 default: 01086 return unknown; 01087 } 01088 break; 01089 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES 01090 case BERNSTEIN: 01091 switch (el_t) 01092 { 01093 case EDGE2: 01094 case EDGE3: 01095 case EDGE4: 01096 return unlimited; 01097 case TRI3: 01098 return 0; 01099 case TRI6: 01100 return 6; 01101 case QUAD4: 01102 return 0; 01103 case QUAD8: 01104 case QUAD9: 01105 return unlimited; 01106 case TET4: 01107 return 1; 01108 case TET10: 01109 return 2; 01110 case HEX8: 01111 return 0; 01112 case HEX20: 01113 return 2; 01114 case HEX27: 01115 return 4; 01116 case PRISM6: 01117 case PRISM15: 01118 case PRISM18: 01119 case PYRAMID5: 01120 case PYRAMID13: 01121 case PYRAMID14: 01122 return 0; 01123 default: 01124 return unknown; 01125 } 01126 break; 01127 case SZABAB: 01128 switch (el_t) 01129 { 01130 case EDGE2: 01131 case EDGE3: 01132 case EDGE4: 01133 return 7; 01134 case TRI3: 01135 return 0; 01136 case TRI6: 01137 return 7; 01138 case QUAD4: 01139 return 0; 01140 case QUAD8: 01141 case QUAD9: 01142 return 7; 01143 case TET4: 01144 case TET10: 01145 case HEX8: 01146 case HEX20: 01147 case HEX27: 01148 case PRISM6: 01149 case PRISM15: 01150 case PRISM18: 01151 case PYRAMID5: 01152 case PYRAMID13: 01153 case PYRAMID14: 01154 return 0; 01155 default: 01156 return unknown; 01157 } 01158 break; 01159 #endif 01160 case XYZ: 01161 switch (el_t) 01162 { 01163 case EDGE2: 01164 case EDGE3: 01165 case EDGE4: 01166 case TRI3: 01167 case TRI6: 01168 case QUAD4: 01169 case QUAD8: 01170 case QUAD9: 01171 case TET4: 01172 case TET10: 01173 case HEX8: 01174 case HEX20: 01175 case HEX27: 01176 case PRISM6: 01177 case PRISM15: 01178 case PRISM18: 01179 case PYRAMID5: 01180 case PYRAMID13: 01181 case PYRAMID14: 01182 return unlimited; 01183 default: 01184 return unknown; 01185 } 01186 break; 01187 case CLOUGH: 01188 switch (el_t) 01189 { 01190 case EDGE2: 01191 case EDGE3: 01192 return 3; 01193 case EDGE4: 01194 case TRI3: 01195 return 0; 01196 case TRI6: 01197 return 3; 01198 case QUAD4: 01199 case QUAD8: 01200 case QUAD9: 01201 case TET4: 01202 case TET10: 01203 case HEX8: 01204 case HEX20: 01205 case HEX27: 01206 case PRISM6: 01207 case PRISM15: 01208 case PRISM18: 01209 case PYRAMID5: 01210 case PYRAMID13: 01211 case PYRAMID14: 01212 return 0; 01213 default: 01214 return unknown; 01215 } 01216 break; 01217 case HERMITE: 01218 switch (el_t) 01219 { 01220 case EDGE2: 01221 case EDGE3: 01222 return unlimited; 01223 case EDGE4: 01224 case TRI3: 01225 case TRI6: 01226 return 0; 01227 case QUAD4: 01228 return 3; 01229 case QUAD8: 01230 case QUAD9: 01231 return unlimited; 01232 case TET4: 01233 case TET10: 01234 return 0; 01235 case HEX8: 01236 return 3; 01237 case HEX20: 01238 case HEX27: 01239 return unlimited; 01240 case PRISM6: 01241 case PRISM15: 01242 case PRISM18: 01243 case PYRAMID5: 01244 case PYRAMID13: 01245 case PYRAMID14: 01246 return 0; 01247 default: 01248 return unknown; 01249 } 01250 break; 01251 case HIERARCHIC: 01252 switch (el_t) 01253 { 01254 case EDGE2: 01255 case EDGE3: 01256 case EDGE4: 01257 return unlimited; 01258 case TRI3: 01259 return 1; 01260 case TRI6: 01261 return unlimited; 01262 case QUAD4: 01263 return 1; 01264 case QUAD8: 01265 case QUAD9: 01266 return unlimited; 01267 case TET4: 01268 case TET10: 01269 return 0; 01270 case HEX8: 01271 case HEX20: 01272 return 1; 01273 case HEX27: 01274 return unlimited; 01275 case PRISM6: 01276 case PRISM15: 01277 case PRISM18: 01278 case PYRAMID5: 01279 case PYRAMID13: 01280 case PYRAMID14: 01281 return 0; 01282 default: 01283 return unknown; 01284 } 01285 break; 01286 case L2_HIERARCHIC: 01287 switch (el_t) 01288 { 01289 case EDGE2: 01290 case EDGE3: 01291 case EDGE4: 01292 return unlimited; 01293 case TRI3: 01294 return 1; 01295 case TRI6: 01296 return unlimited; 01297 case QUAD4: 01298 return 1; 01299 case QUAD8: 01300 case QUAD9: 01301 return unlimited; 01302 case TET4: 01303 case TET10: 01304 return 0; 01305 case HEX8: 01306 case HEX20: 01307 return 1; 01308 case HEX27: 01309 return unlimited; 01310 case PRISM6: 01311 case PRISM15: 01312 case PRISM18: 01313 case PYRAMID5: 01314 case PYRAMID13: 01315 case PYRAMID14: 01316 return 0; 01317 default: 01318 return unknown; 01319 } 01320 break; 01321 case SUBDIVISION: 01322 switch (el_t) 01323 { 01324 case TRI3SUBDIVISION: 01325 return unlimited; 01326 default: 01327 return unknown; 01328 } 01329 break; 01330 case NEDELEC_ONE: 01331 switch (el_t) 01332 { 01333 case TRI6: 01334 case QUAD8: 01335 case QUAD9: 01336 case HEX20: 01337 case HEX27: 01338 return 1; 01339 default: 01340 return 0; 01341 } 01342 break; 01343 default: 01344 return 0; 01345 break; 01346 } 01347 } 01348 01349 01350 01351 bool FEInterface::extra_hanging_dofs(const FEType& fe_t) 01352 { 01353 switch (fe_t.family) 01354 { 01355 case LAGRANGE: 01356 case L2_LAGRANGE: 01357 case MONOMIAL: 01358 case L2_HIERARCHIC: 01359 case XYZ: 01360 case SUBDIVISION: 01361 case LAGRANGE_VEC: 01362 case NEDELEC_ONE: 01363 return false; 01364 case CLOUGH: 01365 case HERMITE: 01366 case HIERARCHIC: 01367 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES 01368 case BERNSTEIN: 01369 case SZABAB: 01370 #endif 01371 default: 01372 return true; 01373 } 01374 } 01375 01376 FEFieldType FEInterface::field_type( const FEType& fe_type ) 01377 { 01378 return FEInterface::field_type( fe_type.family ); 01379 } 01380 01381 FEFieldType FEInterface::field_type( const FEFamily& fe_family ) 01382 { 01383 switch (fe_family) 01384 { 01385 case LAGRANGE_VEC: 01386 case NEDELEC_ONE: 01387 return TYPE_VECTOR; 01388 default: 01389 return TYPE_SCALAR; 01390 } 01391 } 01392 01393 unsigned int FEInterface::n_vec_dim( const MeshBase& mesh, const FEType& fe_type ) 01394 { 01395 switch (fe_type.family) 01396 { 01397 //FIXME: We currently assume that the number of vector components is tied 01398 // to the mesh dimension. This will break for mixed-dimension meshes. 01399 case LAGRANGE_VEC: 01400 case NEDELEC_ONE: 01401 return mesh.mesh_dimension(); 01402 default: 01403 return 1; 01404 } 01405 } 01406 01407 } // namespace libMesh
generated by
