fe_nedelec_one.C
Go to the documentation of this file.00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2012 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 00003 00004 // This library is free software; you can redistribute it and/or 00005 // modify it under the terms of the GNU Lesser General Public 00006 // License as published by the Free Software Foundation; either 00007 // version 2.1 of the License, or (at your option) any later version. 00008 00009 // This library is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 // Lesser General Public License for more details. 00013 00014 // You should have received a copy of the GNU Lesser General Public 00015 // License along with this library; if not, write to the Free Software 00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00017 00018 00019 00020 // Local includes 00021 #include "libmesh/dof_map.h" 00022 #include "libmesh/fe.h" 00023 #include "libmesh/fe_interface.h" 00024 #include "libmesh/elem.h" 00025 #include "libmesh/threads.h" 00026 #include "libmesh/tensor_value.h" 00027 00028 namespace libMesh 00029 { 00030 00031 // ------------------------------------------------------------ 00032 // Nedelec first kind specific implementations 00033 00034 00035 // Anonymous namespace for local helper functions 00036 namespace { 00037 void nedelec_one_nodal_soln(const Elem* elem, 00038 const Order order, 00039 const std::vector<Number>& elem_soln, 00040 const int dim, 00041 std::vector<Number>& nodal_soln) 00042 { 00043 const unsigned int n_nodes = elem->n_nodes(); 00044 const ElemType elem_type = elem->type(); 00045 00046 const Order totalorder = static_cast<Order>(order+elem->p_level()); 00047 00048 nodal_soln.resize(n_nodes*dim); 00049 00050 FEType fe_type(totalorder, NEDELEC_ONE); 00051 00052 switch (totalorder) 00053 { 00054 case FIRST: 00055 { 00056 switch (elem_type) 00057 { 00058 case TRI6: 00059 { 00060 libmesh_assert_equal_to (elem_soln.size(), 3); 00061 libmesh_assert_equal_to (nodal_soln.size(), 6*2); 00062 break; 00063 } 00064 case QUAD8: 00065 case QUAD9: 00066 { 00067 libmesh_assert_equal_to (elem_soln.size(), 4); 00068 00069 if (elem_type == QUAD8) 00070 libmesh_assert_equal_to (nodal_soln.size(), 8*2); 00071 else 00072 libmesh_assert_equal_to (nodal_soln.size(), 9*2); 00073 break; 00074 } 00075 case TET10: 00076 { 00077 libmesh_assert_equal_to (elem_soln.size(), 6); 00078 libmesh_assert_equal_to (nodal_soln.size(), 10*3); 00079 00080 libmesh_not_implemented(); 00081 00082 break; 00083 } 00084 00085 00086 case HEX20: 00087 case HEX27: 00088 { 00089 libmesh_assert_equal_to (elem_soln.size(), 12); 00090 00091 if (elem_type == HEX20) 00092 libmesh_assert_equal_to (nodal_soln.size(), 20*3); 00093 else 00094 libmesh_assert_equal_to (nodal_soln.size(), 27*3); 00095 00096 libmesh_not_implemented(); 00097 00098 break; 00099 } 00100 00101 default: 00102 { 00103 libmesh_error(); 00104 00105 break; 00106 } 00107 00108 } // switch(elem_type) 00109 00110 const unsigned int n_sf = 00111 FEInterface::n_shape_functions(dim, fe_type, elem_type); 00112 00113 std::vector<Point> refspace_nodes; 00114 FEVectorBase::get_refspace_nodes(elem_type,refspace_nodes); 00115 libmesh_assert_equal_to (refspace_nodes.size(), n_nodes); 00116 00117 00118 // Need to create new fe object so the shape function as the FETransformation 00119 // applied to it. 00120 AutoPtr<FEVectorBase> vis_fe = FEVectorBase::build(dim,fe_type); 00121 00122 const std::vector<std::vector<RealGradient> >& vis_phi = vis_fe->get_phi(); 00123 00124 vis_fe->reinit(elem,&refspace_nodes); 00125 00126 for( unsigned int n = 0; n < n_nodes; n++ ) 00127 { 00128 libmesh_assert_equal_to (elem_soln.size(), n_sf); 00129 00130 // Zero before summation 00131 nodal_soln[2*n] = 0; 00132 nodal_soln[2*n+1] = 0; 00133 00134 // u = Sum (u_i phi_i) 00135 for (unsigned int i=0; i<n_sf; i++) 00136 { 00137 nodal_soln[2*n] += elem_soln[i]*(vis_phi[i][n](0)); 00138 nodal_soln[2*n+1] += elem_soln[i]*(vis_phi[i][n](1)); 00139 } 00140 } 00141 00142 return; 00143 } // case FIRST 00144 00145 default: 00146 { 00147 libmesh_error(); 00148 } 00149 00150 }//switch (totalorder) 00151 00152 return; 00153 } // nedelec_one_nodal_soln 00154 00155 00156 unsigned int nedelec_one_n_dofs(const ElemType t, const Order o) 00157 { 00158 switch (o) 00159 { 00160 case FIRST: 00161 { 00162 switch (t) 00163 { 00164 case TRI6: 00165 return 3; 00166 00167 case QUAD8: 00168 case QUAD9: 00169 return 4; 00170 00171 case TET10: 00172 return 6; 00173 00174 case HEX20: 00175 case HEX27: 00176 return 12; 00177 00178 default: 00179 { 00180 #ifdef DEBUG 00181 libMesh::err << "ERROR: Bad ElemType = " << t 00182 << " for FIRST order approximation!" 00183 << std::endl; 00184 #endif 00185 libmesh_error(); 00186 } 00187 } 00188 } 00189 00190 default: 00191 libmesh_error(); 00192 } 00193 00194 libmesh_error(); 00195 return 0; 00196 } 00197 00198 00199 00200 00201 unsigned int nedelec_one_n_dofs_at_node(const ElemType t, 00202 const Order o, 00203 const unsigned int n) 00204 { 00205 switch (o) 00206 { 00207 case FIRST: 00208 { 00209 switch (t) 00210 { 00211 case TRI6: 00212 { 00213 switch (n) 00214 { 00215 case 0: 00216 case 1: 00217 case 2: 00218 return 0; 00219 case 3: 00220 case 4: 00221 case 5: 00222 return 1; 00223 00224 default: 00225 libmesh_error(); 00226 } 00227 } 00228 case QUAD8: 00229 { 00230 switch (n) 00231 { 00232 case 0: 00233 case 1: 00234 case 2: 00235 case 3: 00236 return 0; 00237 case 4: 00238 case 5: 00239 case 6: 00240 case 7: 00241 return 1; 00242 00243 default: 00244 libmesh_error(); 00245 } 00246 } 00247 case QUAD9: 00248 { 00249 switch (n) 00250 { 00251 case 0: 00252 case 1: 00253 case 2: 00254 case 3: 00255 case 8: 00256 return 0; 00257 case 4: 00258 case 5: 00259 case 6: 00260 case 7: 00261 return 1; 00262 00263 default: 00264 libmesh_error(); 00265 } 00266 } 00267 case TET10: 00268 { 00269 switch (n) 00270 { 00271 case 0: 00272 case 1: 00273 case 2: 00274 case 3: 00275 return 0; 00276 case 4: 00277 case 5: 00278 case 6: 00279 case 7: 00280 case 8: 00281 case 9: 00282 return 1; 00283 00284 default: 00285 libmesh_error(); 00286 } 00287 } 00288 00289 case HEX20: 00290 { 00291 switch (n) 00292 { 00293 case 0: 00294 case 1: 00295 case 2: 00296 case 3: 00297 case 4: 00298 case 5: 00299 case 6: 00300 case 7: 00301 return 0; 00302 case 8: 00303 case 9: 00304 case 10: 00305 case 11: 00306 case 12: 00307 case 13: 00308 case 14: 00309 case 15: 00310 case 16: 00311 case 17: 00312 case 18: 00313 case 19: 00314 return 1; 00315 00316 default: 00317 libmesh_error(); 00318 } 00319 } 00320 case HEX27: 00321 { 00322 switch (n) 00323 { 00324 case 0: 00325 case 1: 00326 case 2: 00327 case 3: 00328 case 4: 00329 case 5: 00330 case 6: 00331 case 7: 00332 case 20: 00333 case 21: 00334 case 22: 00335 case 23: 00336 case 24: 00337 case 25: 00338 case 26: 00339 return 0; 00340 case 8: 00341 case 9: 00342 case 10: 00343 case 11: 00344 case 12: 00345 case 13: 00346 case 14: 00347 case 15: 00348 case 16: 00349 case 17: 00350 case 18: 00351 case 19: 00352 return 1; 00353 00354 default: 00355 libmesh_error(); 00356 } 00357 } 00358 default: 00359 { 00360 #ifdef DEBUG 00361 libMesh::err << "ERROR: Bad ElemType = " << t 00362 << " for FIRST order approximation!" 00363 << std::endl; 00364 #endif 00365 libmesh_error(); 00366 } 00367 } 00368 } 00369 00370 default: 00371 libmesh_error(); 00372 } 00373 00374 libmesh_error(); 00375 return 0; 00376 00377 } 00378 00379 00380 00381 #ifdef LIBMESH_ENABLE_AMR 00382 void nedelec_one_compute_constraints (DofConstraints &/*constraints*/, 00383 DofMap &/*dof_map*/, 00384 const unsigned int /*variable_number*/, 00385 const Elem* elem, 00386 const unsigned Dim) 00387 { 00388 // Only constrain elements in 2,3D. 00389 if (Dim == 1) 00390 return; 00391 00392 libmesh_assert(elem); 00393 00394 libmesh_not_implemented(); 00395 00396 /* 00397 // Only constrain active and ancestor elements 00398 if (elem->subactive()) 00399 return; 00400 00401 FEType fe_type = dof_map.variable_type(variable_number); 00402 fe_type.order = static_cast<Order>(fe_type.order + elem->p_level()); 00403 00404 std::vector<unsigned int> my_dof_indices, parent_dof_indices; 00405 00406 // Look at the element faces. Check to see if we need to 00407 // build constraints. 00408 for (unsigned int s=0; s<elem->n_sides(); s++) 00409 if (elem->neighbor(s) != NULL) 00410 if (elem->neighbor(s)->level() < elem->level()) // constrain dofs shared between 00411 { // this element and ones coarser 00412 // than this element. 00413 // Get pointers to the elements of interest and its parent. 00414 const Elem* parent = elem->parent(); 00415 00416 // This can't happen... Only level-0 elements have NULL 00417 // parents, and no level-0 elements can be at a higher 00418 // level than their neighbors! 00419 libmesh_assert(parent); 00420 00421 const AutoPtr<Elem> my_side (elem->build_side(s)); 00422 const AutoPtr<Elem> parent_side (parent->build_side(s)); 00423 00424 // This function gets called element-by-element, so there 00425 // will be a lot of memory allocation going on. We can 00426 // at least minimize this for the case of the dof indices 00427 // by efficiently preallocating the requisite storage. 00428 my_dof_indices.reserve (my_side->n_nodes()); 00429 parent_dof_indices.reserve (parent_side->n_nodes()); 00430 00431 dof_map.dof_indices (my_side.get(), my_dof_indices, 00432 variable_number); 00433 dof_map.dof_indices (parent_side.get(), parent_dof_indices, 00434 variable_number); 00435 00436 for (unsigned int my_dof=0; 00437 my_dof<FEInterface::n_dofs(Dim-1, fe_type, my_side->type()); 00438 my_dof++) 00439 { 00440 libmesh_assert_less (my_dof, my_side->n_nodes()); 00441 00442 // My global dof index. 00443 const unsigned int my_dof_g = my_dof_indices[my_dof]; 00444 00445 // The support point of the DOF 00446 const Point& support_point = my_side->point(my_dof); 00447 00448 // Figure out where my node lies on their reference element. 00449 const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type, 00450 parent_side.get(), 00451 support_point); 00452 00453 // Compute the parent's side shape function values. 00454 for (unsigned int their_dof=0; 00455 their_dof<FEInterface::n_dofs(Dim-1, fe_type, parent_side->type()); 00456 their_dof++) 00457 { 00458 libmesh_assert_less (their_dof, parent_side->n_nodes()); 00459 00460 // Their global dof index. 00461 const unsigned int their_dof_g = 00462 parent_dof_indices[their_dof]; 00463 00464 const Real their_dof_value = FEInterface::shape(Dim-1, 00465 fe_type, 00466 parent_side->type(), 00467 their_dof, 00468 mapped_point); 00469 00470 // Only add non-zero and non-identity values 00471 // for Lagrange basis functions. 00472 if ((std::abs(their_dof_value) > 1.e-5) && 00473 (std::abs(their_dof_value) < .999)) 00474 { 00475 // since we may be running this method concurretly 00476 // on multiple threads we need to acquire a lock 00477 // before modifying the shared constraint_row object. 00478 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx); 00479 00480 // A reference to the constraint row. 00481 DofConstraintRow& constraint_row = constraints[my_dof_g].first; 00482 00483 constraint_row.insert(std::make_pair (their_dof_g, 00484 their_dof_value)); 00485 } 00486 #ifdef DEBUG 00487 // Protect for the case u_i = 0.999 u_j, 00488 // in which case i better equal j. 00489 else if (their_dof_value >= .999) 00490 libmesh_assert_equal_to (my_dof_g, their_dof_g); 00491 #endif 00492 } 00493 } 00494 } 00495 */ 00496 } // nedelec_one_compute_constrants() 00497 #endif // #ifdef LIBMESH_ENABLE_AMR 00498 00499 } // anonymous namespace 00500 00501 #define NEDELEC_LOW_D_ERROR_MESSAGE \ 00502 libMesh::err << "ERROR: This method makes no sense for low-D elements!" \ 00503 << std::endl; \ 00504 libmesh_error(); 00505 00506 00507 // Do full-specialization for every dimension, instead 00508 // of explicit instantiation at the end of this file. 00509 template <> 00510 void FE<0,NEDELEC_ONE>::nodal_soln(const Elem*, 00511 const Order, 00512 const std::vector<Number>&, 00513 std::vector<Number>&) 00514 { NEDELEC_LOW_D_ERROR_MESSAGE } 00515 00516 template <> 00517 void FE<1,NEDELEC_ONE>::nodal_soln(const Elem*, 00518 const Order, 00519 const std::vector<Number>&, 00520 std::vector<Number>&) 00521 { NEDELEC_LOW_D_ERROR_MESSAGE } 00522 00523 template <> 00524 void FE<2,NEDELEC_ONE>::nodal_soln(const Elem* elem, 00525 const Order order, 00526 const std::vector<Number>& elem_soln, 00527 std::vector<Number>& nodal_soln) 00528 { nedelec_one_nodal_soln(elem, order, elem_soln, 2 /*dim*/, nodal_soln); } 00529 00530 template <> 00531 void FE<3,NEDELEC_ONE>::nodal_soln(const Elem* elem, 00532 const Order order, 00533 const std::vector<Number>& elem_soln, 00534 std::vector<Number>& nodal_soln) 00535 { nedelec_one_nodal_soln(elem, order, elem_soln, 3 /*dim*/, nodal_soln); } 00536 00537 00538 // Do full-specialization for every dimension, instead 00539 // of explicit instantiation at the end of this function. 00540 // This could be macro-ified. 00541 template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE } 00542 template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE } 00543 template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs(const ElemType t, const Order o) { return nedelec_one_n_dofs(t, o); } 00544 template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs(const ElemType t, const Order o) { return nedelec_one_n_dofs(t, o); } 00545 00546 00547 // Do full-specialization for every dimension, instead 00548 // of explicit instantiation at the end of this function. 00549 template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE } 00550 template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE } 00551 template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return nedelec_one_n_dofs_at_node(t, o, n); } 00552 template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return nedelec_one_n_dofs_at_node(t, o, n); } 00553 00554 00555 // Nedelec first type elements have no dofs per element 00556 // FIXME: Only for first order! 00557 template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE } 00558 template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE } 00559 template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { return 0; } 00560 template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { return 0; } 00561 00562 // Nedelec first type FEMs are always tangentially continuous 00563 template <> FEContinuity FE<0,NEDELEC_ONE>::get_continuity() const { NEDELEC_LOW_D_ERROR_MESSAGE } 00564 template <> FEContinuity FE<1,NEDELEC_ONE>::get_continuity() const { NEDELEC_LOW_D_ERROR_MESSAGE } 00565 template <> FEContinuity FE<2,NEDELEC_ONE>::get_continuity() const { return H_CURL; } 00566 template <> FEContinuity FE<3,NEDELEC_ONE>::get_continuity() const { return H_CURL; } 00567 00568 // Nedelec first type FEMs are not hierarchic 00569 template <> bool FE<0,NEDELEC_ONE>::is_hierarchic() const { NEDELEC_LOW_D_ERROR_MESSAGE } 00570 template <> bool FE<1,NEDELEC_ONE>::is_hierarchic() const { NEDELEC_LOW_D_ERROR_MESSAGE } 00571 template <> bool FE<2,NEDELEC_ONE>::is_hierarchic() const { return false; } 00572 template <> bool FE<3,NEDELEC_ONE>::is_hierarchic() const { return false; } 00573 00574 // Nedelec first type FEM shapes always need to be reinit'ed (because of orientation dependence) 00575 template <> bool FE<0,NEDELEC_ONE>::shapes_need_reinit() const { NEDELEC_LOW_D_ERROR_MESSAGE } 00576 template <> bool FE<1,NEDELEC_ONE>::shapes_need_reinit() const { NEDELEC_LOW_D_ERROR_MESSAGE } 00577 template <> bool FE<2,NEDELEC_ONE>::shapes_need_reinit() const { return true; } 00578 template <> bool FE<3,NEDELEC_ONE>::shapes_need_reinit() const { return true; } 00579 00580 #ifdef LIBMESH_ENABLE_AMR 00581 template <> 00582 void FE<0,NEDELEC_ONE>::compute_constraints (DofConstraints &, 00583 DofMap &, 00584 const unsigned int, 00585 const Elem*) 00586 { NEDELEC_LOW_D_ERROR_MESSAGE } 00587 00588 template <> 00589 void FE<1,NEDELEC_ONE>::compute_constraints (DofConstraints &, 00590 DofMap &, 00591 const unsigned int, 00592 const Elem*) 00593 { NEDELEC_LOW_D_ERROR_MESSAGE } 00594 00595 template <> 00596 void FE<2,NEDELEC_ONE>::compute_constraints (DofConstraints &constraints, 00597 DofMap &dof_map, 00598 const unsigned int variable_number, 00599 const Elem* elem) 00600 { nedelec_one_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/2); } 00601 00602 template <> 00603 void FE<3,NEDELEC_ONE>::compute_constraints (DofConstraints &constraints, 00604 DofMap &dof_map, 00605 const unsigned int variable_number, 00606 const Elem* elem) 00607 { nedelec_one_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/3); } 00608 #endif // LIBMESH_ENABLE_AMR 00609 00610 // Specialize useless shape function methods 00611 template <> 00612 RealGradient FE<0,NEDELEC_ONE>::shape(const ElemType, const Order,const unsigned int,const Point&) 00613 { NEDELEC_LOW_D_ERROR_MESSAGE } 00614 template <> 00615 RealGradient FE<0,NEDELEC_ONE>::shape(const Elem*,const Order,const unsigned int,const Point&) 00616 { NEDELEC_LOW_D_ERROR_MESSAGE } 00617 template <> 00618 RealGradient FE<0,NEDELEC_ONE>::shape_deriv(const ElemType, const Order,const unsigned int, 00619 const unsigned int,const Point&) 00620 { NEDELEC_LOW_D_ERROR_MESSAGE } 00621 template <> 00622 RealGradient FE<0,NEDELEC_ONE>::shape_deriv(const Elem*,const Order,const unsigned int, 00623 const unsigned int,const Point&) 00624 { NEDELEC_LOW_D_ERROR_MESSAGE } 00625 template <> 00626 RealGradient FE<0,NEDELEC_ONE>::shape_second_deriv(const ElemType, const Order,const unsigned int, 00627 const unsigned int,const Point&) 00628 { NEDELEC_LOW_D_ERROR_MESSAGE } 00629 template <> 00630 RealGradient FE<0,NEDELEC_ONE>::shape_second_deriv(const Elem*,const Order,const unsigned int, 00631 const unsigned int,const Point&) 00632 { NEDELEC_LOW_D_ERROR_MESSAGE } 00633 00634 template <> 00635 RealGradient FE<1,NEDELEC_ONE>::shape(const ElemType, const Order,const unsigned int,const Point&) 00636 { NEDELEC_LOW_D_ERROR_MESSAGE } 00637 template <> 00638 RealGradient FE<1,NEDELEC_ONE>::shape(const Elem*,const Order,const unsigned int,const Point&) 00639 { NEDELEC_LOW_D_ERROR_MESSAGE } 00640 template <> 00641 RealGradient FE<1,NEDELEC_ONE>::shape_deriv(const ElemType, const Order,const unsigned int, 00642 const unsigned int,const Point&) 00643 { NEDELEC_LOW_D_ERROR_MESSAGE } 00644 template <> 00645 RealGradient FE<1,NEDELEC_ONE>::shape_deriv(const Elem*,const Order,const unsigned int, 00646 const unsigned int,const Point&) 00647 { NEDELEC_LOW_D_ERROR_MESSAGE } 00648 template <> 00649 RealGradient FE<1,NEDELEC_ONE>::shape_second_deriv(const ElemType, const Order,const unsigned int, 00650 const unsigned int,const Point&) 00651 { NEDELEC_LOW_D_ERROR_MESSAGE } 00652 template <> 00653 RealGradient FE<1,NEDELEC_ONE>::shape_second_deriv(const Elem*,const Order,const unsigned int, 00654 const unsigned int,const Point&) 00655 { NEDELEC_LOW_D_ERROR_MESSAGE } 00656 00657 } // namespace libMesh
Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:46 UTC
Hosted By: