fe_l2_lagrange.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 00027 namespace libMesh 00028 { 00029 00030 // ------------------------------------------------------------ 00031 // Lagrange-specific implementations 00032 00033 00034 // Anonymous namespace for local helper functions 00035 namespace { 00036 void l2_lagrange_nodal_soln(const Elem* elem, 00037 const Order order, 00038 const std::vector<Number>& elem_soln, 00039 std::vector<Number>& nodal_soln) 00040 { 00041 const unsigned int n_nodes = elem->n_nodes(); 00042 const ElemType type = elem->type(); 00043 00044 const Order totalorder = static_cast<Order>(order+elem->p_level()); 00045 00046 nodal_soln.resize(n_nodes); 00047 00048 00049 00050 switch (totalorder) 00051 { 00052 // linear Lagrange shape functions 00053 case FIRST: 00054 { 00055 switch (type) 00056 { 00057 case EDGE3: 00058 { 00059 libmesh_assert_equal_to (elem_soln.size(), 2); 00060 libmesh_assert_equal_to (nodal_soln.size(), 3); 00061 00062 nodal_soln[0] = elem_soln[0]; 00063 nodal_soln[1] = elem_soln[1]; 00064 nodal_soln[2] = .5*(elem_soln[0] + elem_soln[1]); 00065 00066 return; 00067 } 00068 00069 case EDGE4: 00070 { 00071 libmesh_assert_equal_to (elem_soln.size(), 2); 00072 libmesh_assert_equal_to (nodal_soln.size(), 4); 00073 00074 nodal_soln[0] = elem_soln[0]; 00075 nodal_soln[1] = elem_soln[1]; 00076 nodal_soln[2] = (2.*elem_soln[0] + elem_soln[1])/3.; 00077 nodal_soln[3] = (elem_soln[0] + 2.*elem_soln[1])/3.; 00078 00079 return; 00080 } 00081 00082 00083 case TRI6: 00084 { 00085 libmesh_assert_equal_to (elem_soln.size(), 3); 00086 libmesh_assert_equal_to (nodal_soln.size(), 6); 00087 00088 nodal_soln[0] = elem_soln[0]; 00089 nodal_soln[1] = elem_soln[1]; 00090 nodal_soln[2] = elem_soln[2]; 00091 nodal_soln[3] = .5*(elem_soln[0] + elem_soln[1]); 00092 nodal_soln[4] = .5*(elem_soln[1] + elem_soln[2]); 00093 nodal_soln[5] = .5*(elem_soln[2] + elem_soln[0]); 00094 00095 return; 00096 } 00097 00098 00099 case QUAD8: 00100 case QUAD9: 00101 { 00102 libmesh_assert_equal_to (elem_soln.size(), 4); 00103 00104 if (type == QUAD8) 00105 libmesh_assert_equal_to (nodal_soln.size(), 8); 00106 else 00107 libmesh_assert_equal_to (nodal_soln.size(), 9); 00108 00109 00110 nodal_soln[0] = elem_soln[0]; 00111 nodal_soln[1] = elem_soln[1]; 00112 nodal_soln[2] = elem_soln[2]; 00113 nodal_soln[3] = elem_soln[3]; 00114 nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]); 00115 nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]); 00116 nodal_soln[6] = .5*(elem_soln[2] + elem_soln[3]); 00117 nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]); 00118 00119 if (type == QUAD9) 00120 nodal_soln[8] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]); 00121 00122 return; 00123 } 00124 00125 00126 case TET10: 00127 { 00128 libmesh_assert_equal_to (elem_soln.size(), 4); 00129 libmesh_assert_equal_to (nodal_soln.size(), 10); 00130 00131 nodal_soln[0] = elem_soln[0]; 00132 nodal_soln[1] = elem_soln[1]; 00133 nodal_soln[2] = elem_soln[2]; 00134 nodal_soln[3] = elem_soln[3]; 00135 nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]); 00136 nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]); 00137 nodal_soln[6] = .5*(elem_soln[2] + elem_soln[0]); 00138 nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]); 00139 nodal_soln[8] = .5*(elem_soln[3] + elem_soln[1]); 00140 nodal_soln[9] = .5*(elem_soln[3] + elem_soln[2]); 00141 00142 return; 00143 } 00144 00145 00146 case HEX20: 00147 case HEX27: 00148 { 00149 libmesh_assert_equal_to (elem_soln.size(), 8); 00150 00151 if (type == HEX20) 00152 libmesh_assert_equal_to (nodal_soln.size(), 20); 00153 else 00154 libmesh_assert_equal_to (nodal_soln.size(), 27); 00155 00156 nodal_soln[0] = elem_soln[0]; 00157 nodal_soln[1] = elem_soln[1]; 00158 nodal_soln[2] = elem_soln[2]; 00159 nodal_soln[3] = elem_soln[3]; 00160 nodal_soln[4] = elem_soln[4]; 00161 nodal_soln[5] = elem_soln[5]; 00162 nodal_soln[6] = elem_soln[6]; 00163 nodal_soln[7] = elem_soln[7]; 00164 nodal_soln[8] = .5*(elem_soln[0] + elem_soln[1]); 00165 nodal_soln[9] = .5*(elem_soln[1] + elem_soln[2]); 00166 nodal_soln[10] = .5*(elem_soln[2] + elem_soln[3]); 00167 nodal_soln[11] = .5*(elem_soln[3] + elem_soln[0]); 00168 nodal_soln[12] = .5*(elem_soln[0] + elem_soln[4]); 00169 nodal_soln[13] = .5*(elem_soln[1] + elem_soln[5]); 00170 nodal_soln[14] = .5*(elem_soln[2] + elem_soln[6]); 00171 nodal_soln[15] = .5*(elem_soln[3] + elem_soln[7]); 00172 nodal_soln[16] = .5*(elem_soln[4] + elem_soln[5]); 00173 nodal_soln[17] = .5*(elem_soln[5] + elem_soln[6]); 00174 nodal_soln[18] = .5*(elem_soln[6] + elem_soln[7]); 00175 nodal_soln[19] = .5*(elem_soln[4] + elem_soln[7]); 00176 00177 if (type == HEX27) 00178 { 00179 nodal_soln[20] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]); 00180 nodal_soln[21] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[5]); 00181 nodal_soln[22] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[6]); 00182 nodal_soln[23] = .25*(elem_soln[2] + elem_soln[3] + elem_soln[6] + elem_soln[7]); 00183 nodal_soln[24] = .25*(elem_soln[3] + elem_soln[0] + elem_soln[7] + elem_soln[4]); 00184 nodal_soln[25] = .25*(elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]); 00185 00186 nodal_soln[26] = .125*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3] + 00187 elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]); 00188 } 00189 00190 return; 00191 } 00192 00193 00194 case PRISM15: 00195 case PRISM18: 00196 { 00197 libmesh_assert_equal_to (elem_soln.size(), 6); 00198 00199 if (type == PRISM15) 00200 libmesh_assert_equal_to (nodal_soln.size(), 15); 00201 else 00202 libmesh_assert_equal_to (nodal_soln.size(), 18); 00203 00204 nodal_soln[0] = elem_soln[0]; 00205 nodal_soln[1] = elem_soln[1]; 00206 nodal_soln[2] = elem_soln[2]; 00207 nodal_soln[3] = elem_soln[3]; 00208 nodal_soln[4] = elem_soln[4]; 00209 nodal_soln[5] = elem_soln[5]; 00210 nodal_soln[6] = .5*(elem_soln[0] + elem_soln[1]); 00211 nodal_soln[7] = .5*(elem_soln[1] + elem_soln[2]); 00212 nodal_soln[8] = .5*(elem_soln[0] + elem_soln[2]); 00213 nodal_soln[9] = .5*(elem_soln[0] + elem_soln[3]); 00214 nodal_soln[10] = .5*(elem_soln[1] + elem_soln[4]); 00215 nodal_soln[11] = .5*(elem_soln[2] + elem_soln[5]); 00216 nodal_soln[12] = .5*(elem_soln[3] + elem_soln[4]); 00217 nodal_soln[13] = .5*(elem_soln[4] + elem_soln[5]); 00218 nodal_soln[14] = .5*(elem_soln[3] + elem_soln[5]); 00219 00220 if (type == PRISM18) 00221 { 00222 nodal_soln[15] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[3]); 00223 nodal_soln[16] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[4]); 00224 nodal_soln[17] = .25*(elem_soln[2] + elem_soln[0] + elem_soln[3] + elem_soln[5]); 00225 } 00226 00227 return; 00228 } 00229 00230 00231 00232 default: 00233 { 00234 // By default the element solution _is_ nodal, 00235 // so just copy it. 00236 nodal_soln = elem_soln; 00237 00238 return; 00239 } 00240 } 00241 } 00242 00243 case SECOND: 00244 { 00245 switch (type) 00246 { 00247 case EDGE4: 00248 { 00249 libmesh_assert_equal_to (elem_soln.size(), 3); 00250 libmesh_assert_equal_to (nodal_soln.size(), 4); 00251 00252 // Project quadratic solution onto cubic element nodes 00253 nodal_soln[0] = elem_soln[0]; 00254 nodal_soln[1] = elem_soln[1]; 00255 nodal_soln[2] = (2.*elem_soln[0] - elem_soln[1] + 00256 8.*elem_soln[2])/9.; 00257 nodal_soln[3] = (-elem_soln[0] + 2.*elem_soln[1] + 00258 8.*elem_soln[2])/9.; 00259 return; 00260 } 00261 00262 default: 00263 { 00264 // By default the element solution _is_ nodal, 00265 // so just copy it. 00266 nodal_soln = elem_soln; 00267 00268 return; 00269 } 00270 } 00271 } 00272 00273 00274 00275 00276 default: 00277 { 00278 // By default the element solution _is_ nodal, 00279 // so just copy it. 00280 nodal_soln = elem_soln; 00281 00282 return; 00283 } 00284 } 00285 } 00286 00287 00288 // TODO: We should make this work, for example, for SECOND on a TRI3 00289 // (this is valid with L2_LAGRANGE, but not with LAGRANGE) 00290 unsigned int l2_lagrange_n_dofs(const ElemType t, const Order o) 00291 { 00292 switch (o) 00293 { 00294 00295 // linear Lagrange shape functions 00296 case FIRST: 00297 { 00298 switch (t) 00299 { 00300 case NODEELEM: 00301 return 1; 00302 00303 case EDGE2: 00304 case EDGE3: 00305 case EDGE4: 00306 return 2; 00307 00308 case TRI3: 00309 case TRI6: 00310 return 3; 00311 00312 case QUAD4: 00313 case QUAD8: 00314 case QUAD9: 00315 return 4; 00316 00317 case TET4: 00318 case TET10: 00319 return 4; 00320 00321 case HEX8: 00322 case HEX20: 00323 case HEX27: 00324 return 8; 00325 00326 case PRISM6: 00327 case PRISM15: 00328 case PRISM18: 00329 return 6; 00330 00331 case PYRAMID5: 00332 return 5; 00333 00334 default: 00335 { 00336 #ifdef DEBUG 00337 libMesh::err << "ERROR: Bad ElemType = " << t 00338 << " for FIRST order approximation!" 00339 << std::endl; 00340 #endif 00341 libmesh_error(); 00342 } 00343 } 00344 } 00345 00346 00347 // quadratic Lagrange shape functions 00348 case SECOND: 00349 { 00350 switch (t) 00351 { 00352 case NODEELEM: 00353 return 1; 00354 00355 case EDGE3: 00356 return 3; 00357 00358 case TRI6: 00359 return 6; 00360 00361 case QUAD8: 00362 return 8; 00363 00364 case QUAD9: 00365 return 9; 00366 00367 case TET10: 00368 return 10; 00369 00370 case HEX20: 00371 return 20; 00372 00373 case HEX27: 00374 return 27; 00375 00376 case PRISM15: 00377 return 15; 00378 00379 case PRISM18: 00380 return 18; 00381 00382 default: 00383 { 00384 #ifdef DEBUG 00385 libMesh::err << "ERROR: Bad ElemType = " << t 00386 << " for SECOND order approximation!" 00387 << std::endl; 00388 #endif 00389 libmesh_error(); 00390 } 00391 } 00392 } 00393 00394 case THIRD: 00395 { 00396 switch (t) 00397 { 00398 case NODEELEM: 00399 return 1; 00400 00401 case EDGE4: 00402 return 4; 00403 00404 default: 00405 { 00406 #ifdef DEBUG 00407 libMesh::err << "ERROR: Bad ElemType = " << t 00408 << " for THIRD order approximation!" 00409 << std::endl; 00410 #endif 00411 libmesh_error(); 00412 } 00413 } 00414 } 00415 00416 default: 00417 libmesh_error(); 00418 } 00419 00420 libmesh_error(); 00421 return 0; 00422 } 00423 00424 } // anonymous namespace 00425 00426 00427 // Do full-specialization for every dimension, instead 00428 // of explicit instantiation at the end of this file. 00429 // This could be macro-ified so that it fits on one line... 00430 template <> 00431 void FE<0,L2_LAGRANGE>::nodal_soln(const Elem* elem, 00432 const Order order, 00433 const std::vector<Number>& elem_soln, 00434 std::vector<Number>& nodal_soln) 00435 { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); } 00436 00437 template <> 00438 void FE<1,L2_LAGRANGE>::nodal_soln(const Elem* elem, 00439 const Order order, 00440 const std::vector<Number>& elem_soln, 00441 std::vector<Number>& nodal_soln) 00442 { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); } 00443 00444 template <> 00445 void FE<2,L2_LAGRANGE>::nodal_soln(const Elem* elem, 00446 const Order order, 00447 const std::vector<Number>& elem_soln, 00448 std::vector<Number>& nodal_soln) 00449 { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); } 00450 00451 template <> 00452 void FE<3,L2_LAGRANGE>::nodal_soln(const Elem* elem, 00453 const Order order, 00454 const std::vector<Number>& elem_soln, 00455 std::vector<Number>& nodal_soln) 00456 { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); } 00457 00458 00459 // Do full-specialization for every dimension, instead 00460 // of explicit instantiation at the end of this function. 00461 // This could be macro-ified. 00462 template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); } 00463 template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); } 00464 template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); } 00465 template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); } 00466 00467 00468 // L2 Lagrange elements have all dofs on elements (hence none on nodes) 00469 template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; } 00470 template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; } 00471 template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; } 00472 template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; } 00473 00474 00475 // L2 Lagrange elements have all dofs on elements 00476 template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); } 00477 template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); } 00478 template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); } 00479 template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); } 00480 00481 // L2 Lagrange FEMs are DISCONTINUOUS 00482 template <> FEContinuity FE<0,L2_LAGRANGE>::get_continuity() const { return DISCONTINUOUS; } 00483 template <> FEContinuity FE<1,L2_LAGRANGE>::get_continuity() const { return DISCONTINUOUS; } 00484 template <> FEContinuity FE<2,L2_LAGRANGE>::get_continuity() const { return DISCONTINUOUS; } 00485 template <> FEContinuity FE<3,L2_LAGRANGE>::get_continuity() const { return DISCONTINUOUS; } 00486 00487 // Lagrange FEMs are not hierarchic 00488 template <> bool FE<0,L2_LAGRANGE>::is_hierarchic() const { return false; } 00489 template <> bool FE<1,L2_LAGRANGE>::is_hierarchic() const { return false; } 00490 template <> bool FE<2,L2_LAGRANGE>::is_hierarchic() const { return false; } 00491 template <> bool FE<3,L2_LAGRANGE>::is_hierarchic() const { return false; } 00492 00493 // Lagrange FEM shapes do not need reinit (is this always true?) 00494 template <> bool FE<0,L2_LAGRANGE>::shapes_need_reinit() const { return false; } 00495 template <> bool FE<1,L2_LAGRANGE>::shapes_need_reinit() const { return false; } 00496 template <> bool FE<2,L2_LAGRANGE>::shapes_need_reinit() const { return false; } 00497 template <> bool FE<3,L2_LAGRANGE>::shapes_need_reinit() const { return false; } 00498 00499 // We don't need any constraints for this DISCONTINUOUS basis, hence these 00500 // are NOOPs 00501 #ifdef LIBMESH_ENABLE_AMR 00502 template <> 00503 void FE<2,L2_LAGRANGE>::compute_constraints (DofConstraints &, 00504 DofMap &, 00505 const unsigned int , 00506 const Elem* ) 00507 { } 00508 00509 template <> 00510 void FE<3,L2_LAGRANGE>::compute_constraints (DofConstraints &, 00511 DofMap &, 00512 const unsigned int , 00513 const Elem* ) 00514 { } 00515 #endif // LIBMESH_ENABLE_AMR 00516 00517 } // namespace libMesh
Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:46 UTC
Hosted By: