fe_lagrange_vec.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/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 // Lagrange-specific implementations 00033 00034 00035 // Anonymous namespace for local helper functions 00036 namespace { 00037 void lagrange_vec_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 type = elem->type(); 00045 00046 const Order totalorder = static_cast<Order>(order+elem->p_level()); 00047 00048 nodal_soln.resize(dim*n_nodes); 00049 00050 switch (totalorder) 00051 { 00052 // linear Lagrange shape functions 00053 case FIRST: 00054 { 00055 switch (type) 00056 { 00057 case TRI6: 00058 { 00059 libmesh_assert_equal_to (elem_soln.size(), 2*3); 00060 libmesh_assert_equal_to (nodal_soln.size(), 2*6); 00061 00062 // node 0 components 00063 nodal_soln[0] = elem_soln[0]; 00064 nodal_soln[1] = elem_soln[1]; 00065 00066 // node 1 components 00067 nodal_soln[2] = elem_soln[2]; 00068 nodal_soln[3] = elem_soln[3]; 00069 00070 // node 2 components 00071 nodal_soln[4] = elem_soln[4]; 00072 nodal_soln[5] = elem_soln[5]; 00073 00074 // node 3 components 00075 nodal_soln[6] = .5*(elem_soln[0] + elem_soln[2]); 00076 nodal_soln[7] = .5*(elem_soln[1] + elem_soln[3]); 00077 00078 // node 4 components 00079 nodal_soln[8] = .5*(elem_soln[2] + elem_soln[4]); 00080 nodal_soln[9] = .5*(elem_soln[3] + elem_soln[5]); 00081 00082 // node 5 components 00083 nodal_soln[10] = .5*(elem_soln[0] + elem_soln[4]); 00084 nodal_soln[11] = .5*(elem_soln[1] + elem_soln[5]); 00085 00086 return; 00087 } 00088 00089 00090 case QUAD8: 00091 case QUAD9: 00092 { 00093 libmesh_assert_equal_to (elem_soln.size(), 2*4); 00094 00095 if (type == QUAD8) 00096 libmesh_assert_equal_to (nodal_soln.size(), 2*8); 00097 else 00098 libmesh_assert_equal_to (nodal_soln.size(), 2*9); 00099 00100 // node 0 components 00101 nodal_soln[0] = elem_soln[0]; 00102 nodal_soln[1] = elem_soln[1]; 00103 00104 // node 1 components 00105 nodal_soln[2] = elem_soln[2]; 00106 nodal_soln[3] = elem_soln[3]; 00107 00108 // node 2 components 00109 nodal_soln[4] = elem_soln[4]; 00110 nodal_soln[5] = elem_soln[5]; 00111 00112 // node 3 components 00113 nodal_soln[6] = elem_soln[6]; 00114 nodal_soln[7] = elem_soln[7]; 00115 00116 // node 4 components 00117 nodal_soln[8] = .5*(elem_soln[0] + elem_soln[2]); 00118 nodal_soln[9] = .5*(elem_soln[1] + elem_soln[3]); 00119 00120 // node 5 components 00121 nodal_soln[10] = .5*(elem_soln[2] + elem_soln[4]); 00122 nodal_soln[11] = .5*(elem_soln[3] + elem_soln[5]); 00123 00124 // node 6 components 00125 nodal_soln[12] = .5*(elem_soln[4] + elem_soln[6]); 00126 nodal_soln[13] = .5*(elem_soln[5] + elem_soln[7]); 00127 00128 // node 7 components 00129 nodal_soln[14] = .5*(elem_soln[6] + elem_soln[0]); 00130 nodal_soln[15] = .5*(elem_soln[7] + elem_soln[1]); 00131 00132 if (type == QUAD9) 00133 { 00134 // node 8 components 00135 nodal_soln[16] = .25*(elem_soln[0] + elem_soln[2] + elem_soln[4] + elem_soln[6]); 00136 nodal_soln[17] = .25*(elem_soln[1] + elem_soln[3] + elem_soln[5] + elem_soln[7]); 00137 } 00138 00139 return; 00140 } 00141 00142 00143 case TET10: 00144 { 00145 libmesh_assert_equal_to (elem_soln.size(), 3*4); 00146 libmesh_assert_equal_to (nodal_soln.size(), 3*10); 00147 00148 // node 0 components 00149 nodal_soln[0] = elem_soln[0]; 00150 nodal_soln[1] = elem_soln[1]; 00151 nodal_soln[2] = elem_soln[2]; 00152 00153 // node 1 components 00154 nodal_soln[3] = elem_soln[3]; 00155 nodal_soln[4] = elem_soln[4]; 00156 nodal_soln[5] = elem_soln[5]; 00157 00158 // node 2 components 00159 nodal_soln[6] = elem_soln[6]; 00160 nodal_soln[7] = elem_soln[7]; 00161 nodal_soln[8] = elem_soln[8]; 00162 00163 // node 3 components 00164 nodal_soln[9] = elem_soln[9]; 00165 nodal_soln[10] = elem_soln[10]; 00166 nodal_soln[11] = elem_soln[11]; 00167 00168 // node 4 components 00169 nodal_soln[12] = .5*(elem_soln[0] + elem_soln[3]); 00170 nodal_soln[13] = .5*(elem_soln[1] + elem_soln[4]); 00171 nodal_soln[14] = .5*(elem_soln[2] + elem_soln[5]); 00172 00173 // node 5 components 00174 nodal_soln[15] = .5*(elem_soln[3] + elem_soln[6]); 00175 nodal_soln[16] = .5*(elem_soln[4] + elem_soln[7]); 00176 nodal_soln[17] = .5*(elem_soln[5] + elem_soln[8]); 00177 00178 // node 6 components 00179 nodal_soln[18] = .5*(elem_soln[6] + elem_soln[0]); 00180 nodal_soln[19] = .5*(elem_soln[7] + elem_soln[1]); 00181 nodal_soln[20] = .5*(elem_soln[8] + elem_soln[2]); 00182 00183 // node 7 components 00184 nodal_soln[21] = .5*(elem_soln[9] + elem_soln[0]); 00185 nodal_soln[22] = .5*(elem_soln[10] + elem_soln[1]); 00186 nodal_soln[23] = .5*(elem_soln[11] + elem_soln[2]); 00187 00188 // node 8 components 00189 nodal_soln[24] = .5*(elem_soln[9] + elem_soln[3]); 00190 nodal_soln[25] = .5*(elem_soln[10] + elem_soln[4]); 00191 nodal_soln[26] = .5*(elem_soln[11] + elem_soln[5]); 00192 00193 // node 9 components 00194 nodal_soln[27] = .5*(elem_soln[9] + elem_soln[6]); 00195 nodal_soln[28] = .5*(elem_soln[10] + elem_soln[7]); 00196 nodal_soln[29] = .5*(elem_soln[11] + elem_soln[8]); 00197 00198 return; 00199 } 00200 00201 00202 case HEX20: 00203 case HEX27: 00204 { 00205 libmesh_assert_equal_to (elem_soln.size(), 3*8); 00206 00207 if (type == HEX20) 00208 libmesh_assert_equal_to (nodal_soln.size(), 3*20); 00209 else 00210 libmesh_assert_equal_to (nodal_soln.size(), 3*27); 00211 00212 // node 0 components 00213 nodal_soln[0] = elem_soln[0]; 00214 nodal_soln[1] = elem_soln[1]; 00215 nodal_soln[2] = elem_soln[2]; 00216 00217 // node 1 components 00218 nodal_soln[3] = elem_soln[3]; 00219 nodal_soln[4] = elem_soln[4]; 00220 nodal_soln[5] = elem_soln[5]; 00221 00222 // node 2 components 00223 nodal_soln[6] = elem_soln[6]; 00224 nodal_soln[7] = elem_soln[7]; 00225 nodal_soln[8] = elem_soln[8]; 00226 00227 // node 3 components 00228 nodal_soln[9] = elem_soln[9]; 00229 nodal_soln[10] = elem_soln[10]; 00230 nodal_soln[11] = elem_soln[11]; 00231 00232 // node 4 components 00233 nodal_soln[12] = elem_soln[12]; 00234 nodal_soln[13] = elem_soln[13]; 00235 nodal_soln[14] = elem_soln[14]; 00236 00237 // node 5 components 00238 nodal_soln[15] = elem_soln[15]; 00239 nodal_soln[16] = elem_soln[16]; 00240 nodal_soln[17] = elem_soln[17]; 00241 00242 // node 6 components 00243 nodal_soln[18] = elem_soln[18]; 00244 nodal_soln[19] = elem_soln[19]; 00245 nodal_soln[20] = elem_soln[20]; 00246 00247 // node 7 components 00248 nodal_soln[21] = elem_soln[21]; 00249 nodal_soln[22] = elem_soln[22]; 00250 nodal_soln[23] = elem_soln[23]; 00251 00252 // node 8 components 00253 nodal_soln[24] = .5*(elem_soln[0] + elem_soln[3]); 00254 nodal_soln[25] = .5*(elem_soln[1] + elem_soln[4]); 00255 nodal_soln[26] = .5*(elem_soln[2] + elem_soln[5]); 00256 00257 // node 9 components 00258 nodal_soln[27] = .5*(elem_soln[3] + elem_soln[6]); 00259 nodal_soln[28] = .5*(elem_soln[4] + elem_soln[7]); 00260 nodal_soln[29] = .5*(elem_soln[5] + elem_soln[8]); 00261 00262 // node 10 components 00263 nodal_soln[30] = .5*(elem_soln[6] + elem_soln[9]); 00264 nodal_soln[31] = .5*(elem_soln[7] + elem_soln[10]); 00265 nodal_soln[32] = .5*(elem_soln[8] + elem_soln[11]); 00266 00267 // node 11 components 00268 nodal_soln[33] = .5*(elem_soln[9] + elem_soln[0]); 00269 nodal_soln[34] = .5*(elem_soln[10] + elem_soln[1]); 00270 nodal_soln[35] = .5*(elem_soln[11] + elem_soln[2]); 00271 00272 // node 12 components 00273 nodal_soln[36] = .5*(elem_soln[0] + elem_soln[12]); 00274 nodal_soln[37] = .5*(elem_soln[1] + elem_soln[13]); 00275 nodal_soln[38] = .5*(elem_soln[2] + elem_soln[14]); 00276 00277 // node 13 components 00278 nodal_soln[39] = .5*(elem_soln[3] + elem_soln[15]); 00279 nodal_soln[40] = .5*(elem_soln[4] + elem_soln[16]); 00280 nodal_soln[41] = .5*(elem_soln[5] + elem_soln[17]); 00281 00282 // node 14 components 00283 nodal_soln[42] = .5*(elem_soln[6] + elem_soln[18]); 00284 nodal_soln[43] = .5*(elem_soln[7] + elem_soln[19]); 00285 nodal_soln[44] = .5*(elem_soln[8] + elem_soln[20]); 00286 00287 // node 15 components 00288 nodal_soln[45] = .5*(elem_soln[9] + elem_soln[21]); 00289 nodal_soln[46] = .5*(elem_soln[10] + elem_soln[22]); 00290 nodal_soln[47] = .5*(elem_soln[11] + elem_soln[23]); 00291 00292 // node 16 components 00293 nodal_soln[48] = .5*(elem_soln[12] + elem_soln[15]); 00294 nodal_soln[49] = .5*(elem_soln[13] + elem_soln[16]); 00295 nodal_soln[50] = .5*(elem_soln[14] + elem_soln[17]); 00296 00297 // node 17 components 00298 nodal_soln[51] = .5*(elem_soln[15] + elem_soln[18]); 00299 nodal_soln[52] = .5*(elem_soln[16] + elem_soln[19]); 00300 nodal_soln[53] = .5*(elem_soln[17] + elem_soln[20]); 00301 00302 // node 18 components 00303 nodal_soln[54] = .5*(elem_soln[18] + elem_soln[21]); 00304 nodal_soln[55] = .5*(elem_soln[19] + elem_soln[22]); 00305 nodal_soln[56] = .5*(elem_soln[20] + elem_soln[23]); 00306 00307 // node 19 components 00308 nodal_soln[57] = .5*(elem_soln[12] + elem_soln[21]); 00309 nodal_soln[58] = .5*(elem_soln[13] + elem_soln[22]); 00310 nodal_soln[59] = .5*(elem_soln[14] + elem_soln[23]); 00311 00312 if (type == HEX27) 00313 { 00314 // node 20 components 00315 nodal_soln[60] = .25*(elem_soln[0] + elem_soln[3] + elem_soln[6] + elem_soln[9]); 00316 nodal_soln[61] = .25*(elem_soln[1] + elem_soln[4] + elem_soln[7] + elem_soln[10]); 00317 nodal_soln[62] = .25*(elem_soln[2] + elem_soln[5] + elem_soln[8] + elem_soln[11]); 00318 00319 // node 21 components 00320 nodal_soln[63] = .25*(elem_soln[0] + elem_soln[3] + elem_soln[12] + elem_soln[15]); 00321 nodal_soln[64] = .25*(elem_soln[1] + elem_soln[4] + elem_soln[13] + elem_soln[16]); 00322 nodal_soln[65] = .25*(elem_soln[2] + elem_soln[5] + elem_soln[14] + elem_soln[17]); 00323 00324 // node 22 components 00325 nodal_soln[66] = .25*(elem_soln[3] + elem_soln[6] + elem_soln[15] + elem_soln[18]); 00326 nodal_soln[67] = .25*(elem_soln[4] + elem_soln[7] + elem_soln[16] + elem_soln[19]); 00327 nodal_soln[68] = .25*(elem_soln[5] + elem_soln[8] + elem_soln[17] + elem_soln[20]); 00328 00329 // node 23 components 00330 nodal_soln[69] = .25*(elem_soln[6] + elem_soln[9] + elem_soln[18] + elem_soln[21]); 00331 nodal_soln[70] = .25*(elem_soln[7] + elem_soln[10] + elem_soln[19] + elem_soln[22]); 00332 nodal_soln[71] = .25*(elem_soln[8] + elem_soln[11] + elem_soln[20] + elem_soln[23]); 00333 00334 // node 24 components 00335 nodal_soln[72] = .25*(elem_soln[9] + elem_soln[0] + elem_soln[21] + elem_soln[12]); 00336 nodal_soln[73] = .25*(elem_soln[10] + elem_soln[1] + elem_soln[22] + elem_soln[13]); 00337 nodal_soln[74] = .25*(elem_soln[11] + elem_soln[2] + elem_soln[23] + elem_soln[14]); 00338 00339 // node 25 components 00340 nodal_soln[75] = .25*(elem_soln[12] + elem_soln[15] + elem_soln[18] + elem_soln[21]); 00341 nodal_soln[76] = .25*(elem_soln[13] + elem_soln[16] + elem_soln[19] + elem_soln[22]); 00342 nodal_soln[77] = .25*(elem_soln[14] + elem_soln[17] + elem_soln[20] + elem_soln[23]); 00343 00344 // node 26 components 00345 nodal_soln[78] = .125*(elem_soln[0] + elem_soln[3] + elem_soln[6] + elem_soln[9] + 00346 elem_soln[12] + elem_soln[15] + elem_soln[18] + elem_soln[21]); 00347 00348 nodal_soln[79] = .125*(elem_soln[1] + elem_soln[4] + elem_soln[7] + elem_soln[10] + 00349 elem_soln[13] + elem_soln[16] + elem_soln[19] + elem_soln[22]); 00350 00351 nodal_soln[80] = .125*(elem_soln[2] + elem_soln[5] + elem_soln[8] + elem_soln[11] + 00352 elem_soln[14] + elem_soln[17] + elem_soln[20] + elem_soln[23]); 00353 } 00354 00355 return; 00356 } 00357 00358 00359 case PRISM15: 00360 case PRISM18: 00361 { 00362 libmesh_assert_equal_to (elem_soln.size(), 3*6); 00363 00364 if (type == PRISM15) 00365 libmesh_assert_equal_to (nodal_soln.size(), 3*15); 00366 else 00367 libmesh_assert_equal_to (nodal_soln.size(), 3*18); 00368 00369 // node 0 components 00370 nodal_soln[0] = elem_soln[0]; 00371 nodal_soln[1] = elem_soln[1]; 00372 nodal_soln[2] = elem_soln[2]; 00373 00374 // node 1 components 00375 nodal_soln[3] = elem_soln[3]; 00376 nodal_soln[4] = elem_soln[4]; 00377 nodal_soln[5] = elem_soln[5]; 00378 00379 // node 2 components 00380 nodal_soln[6] = elem_soln[6]; 00381 nodal_soln[7] = elem_soln[7]; 00382 nodal_soln[8] = elem_soln[8]; 00383 00384 // node 3 components 00385 nodal_soln[9] = elem_soln[9]; 00386 nodal_soln[10] = elem_soln[10]; 00387 nodal_soln[11] = elem_soln[11]; 00388 00389 // node 4 components 00390 nodal_soln[12] = elem_soln[12]; 00391 nodal_soln[13] = elem_soln[13]; 00392 nodal_soln[14] = elem_soln[14]; 00393 00394 // node 5 components 00395 nodal_soln[15] = elem_soln[15]; 00396 nodal_soln[16] = elem_soln[16]; 00397 nodal_soln[17] = elem_soln[17]; 00398 00399 // node 6 components 00400 nodal_soln[18] = .5*(elem_soln[0] + elem_soln[3]); 00401 nodal_soln[19] = .5*(elem_soln[1] + elem_soln[4]); 00402 nodal_soln[20] = .5*(elem_soln[2] + elem_soln[5]); 00403 00404 // node 7 components 00405 nodal_soln[21] = .5*(elem_soln[3] + elem_soln[6]); 00406 nodal_soln[22] = .5*(elem_soln[4] + elem_soln[7]); 00407 nodal_soln[23] = .5*(elem_soln[5] + elem_soln[8]); 00408 00409 // node 8 components 00410 nodal_soln[24] = .5*(elem_soln[0] + elem_soln[6]); 00411 nodal_soln[25] = .5*(elem_soln[1] + elem_soln[7]); 00412 nodal_soln[26] = .5*(elem_soln[2] + elem_soln[8]); 00413 00414 // node 9 components 00415 nodal_soln[27] = .5*(elem_soln[0] + elem_soln[9]); 00416 nodal_soln[28] = .5*(elem_soln[1] + elem_soln[10]); 00417 nodal_soln[29] = .5*(elem_soln[2] + elem_soln[11]); 00418 00419 // node 10 components 00420 nodal_soln[30] = .5*(elem_soln[3] + elem_soln[12]); 00421 nodal_soln[31] = .5*(elem_soln[4] + elem_soln[13]); 00422 nodal_soln[32] = .5*(elem_soln[5] + elem_soln[14]); 00423 00424 // node 11 components 00425 nodal_soln[33] = .5*(elem_soln[6] + elem_soln[15]); 00426 nodal_soln[34] = .5*(elem_soln[7] + elem_soln[16]); 00427 nodal_soln[35] = .5*(elem_soln[8] + elem_soln[17]); 00428 00429 // node 12 components 00430 nodal_soln[36] = .5*(elem_soln[9] + elem_soln[12]); 00431 nodal_soln[37] = .5*(elem_soln[10] + elem_soln[13]); 00432 nodal_soln[38] = .5*(elem_soln[11] + elem_soln[14]); 00433 00434 // node 13 components 00435 nodal_soln[39] = .5*(elem_soln[12] + elem_soln[15]); 00436 nodal_soln[40] = .5*(elem_soln[13] + elem_soln[16]); 00437 nodal_soln[41] = .5*(elem_soln[14] + elem_soln[17]); 00438 00439 // node 14 components 00440 nodal_soln[42] = .5*(elem_soln[12] + elem_soln[15]); 00441 nodal_soln[43] = .5*(elem_soln[13] + elem_soln[16]); 00442 nodal_soln[44] = .5*(elem_soln[14] + elem_soln[17]); 00443 00444 if (type == PRISM18) 00445 { 00446 // node 15 components 00447 nodal_soln[45] = .25*(elem_soln[0] + elem_soln[3] + elem_soln[12] + elem_soln[9]); 00448 nodal_soln[46] = .25*(elem_soln[1] + elem_soln[4] + elem_soln[13] + elem_soln[10]); 00449 nodal_soln[47] = .25*(elem_soln[2] + elem_soln[5] + elem_soln[14] + elem_soln[11]); 00450 00451 // node 16 components 00452 nodal_soln[48] = .25*(elem_soln[3] + elem_soln[6] + elem_soln[15] + elem_soln[12]); 00453 nodal_soln[49] = .25*(elem_soln[4] + elem_soln[7] + elem_soln[16] + elem_soln[13]); 00454 nodal_soln[50] = .25*(elem_soln[5] + elem_soln[8] + elem_soln[17] + elem_soln[14]); 00455 00456 // node 17 components 00457 nodal_soln[51] = .25*(elem_soln[6] + elem_soln[0] + elem_soln[9] + elem_soln[15]); 00458 nodal_soln[52] = .25*(elem_soln[7] + elem_soln[1] + elem_soln[10] + elem_soln[16]); 00459 nodal_soln[53] = .25*(elem_soln[8] + elem_soln[2] + elem_soln[11] + elem_soln[17]); 00460 } 00461 00462 return; 00463 } 00464 00465 default: 00466 { 00467 // By default the element solution _is_ nodal, 00468 // so just copy it. 00469 nodal_soln = elem_soln; 00470 00471 return; 00472 } 00473 } 00474 } 00475 00476 case SECOND: 00477 { 00478 switch (type) 00479 { 00480 default: 00481 { 00482 // By default the element solution _is_ nodal, 00483 // so just copy it. 00484 nodal_soln = elem_soln; 00485 00486 return; 00487 } 00488 } 00489 } 00490 00491 default: 00492 { 00493 00494 } 00495 00496 } // switch(totalorder) 00497 00498 }// void lagrange_vec_nodal_soln 00499 00500 } // anonymous namespace 00501 00502 00503 // Do full-specialization for every dimension, instead 00504 // of explicit instantiation at the end of this file. 00505 // This could be macro-ified so that it fits on one line... 00506 template <> 00507 void FE<0,LAGRANGE_VEC>::nodal_soln(const Elem* elem, 00508 const Order order, 00509 const std::vector<Number>& elem_soln, 00510 std::vector<Number>& nodal_soln) 00511 { FE<0,LAGRANGE>::nodal_soln(elem, order, elem_soln, nodal_soln); } 00512 00513 template <> 00514 void FE<1,LAGRANGE_VEC>::nodal_soln(const Elem* elem, 00515 const Order order, 00516 const std::vector<Number>& elem_soln, 00517 std::vector<Number>& nodal_soln) 00518 { FE<1,LAGRANGE>::nodal_soln(elem, order, elem_soln, nodal_soln); } 00519 00520 template <> 00521 void FE<2,LAGRANGE_VEC>::nodal_soln(const Elem* elem, 00522 const Order order, 00523 const std::vector<Number>& elem_soln, 00524 std::vector<Number>& nodal_soln) 00525 { lagrange_vec_nodal_soln(elem, order, elem_soln, 2 /*dimension*/, nodal_soln); } 00526 00527 template <> 00528 void FE<3,LAGRANGE_VEC>::nodal_soln(const Elem* elem, 00529 const Order order, 00530 const std::vector<Number>& elem_soln, 00531 std::vector<Number>& nodal_soln) 00532 { lagrange_vec_nodal_soln(elem, order, elem_soln, 3 /*dimension*/, nodal_soln); } 00533 00534 00535 // Specialize for shape function routines by leveraging scalar LAGRANGE elements 00536 00537 // 0-D 00538 template <> RealGradient FE<0,LAGRANGE_VEC>::shape(const ElemType type, const Order order, 00539 const unsigned int i, const Point& p) 00540 { 00541 Real value = FE<0,LAGRANGE>::shape( type, order, i, p ); 00542 return libMesh::RealGradient( value ); 00543 } 00544 template <> RealGradient FE<0,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order, 00545 const unsigned int i, const unsigned int j, 00546 const Point& p) 00547 { 00548 Real value = FE<0,LAGRANGE>::shape_deriv( type, order, i, j, p ); 00549 return libMesh::RealGradient( value ); 00550 } 00551 template <> RealGradient FE<0,LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order, 00552 const unsigned int i, const unsigned int j, 00553 const Point& p) 00554 { 00555 Real value = FE<0,LAGRANGE>::shape_second_deriv( type, order, i, j, p ); 00556 return libMesh::RealGradient( value ); 00557 } 00558 00559 // 1-D 00560 template <> RealGradient FE<1,LAGRANGE_VEC>::shape(const ElemType type, const Order order, 00561 const unsigned int i, const Point& p) 00562 { 00563 Real value = FE<1,LAGRANGE>::shape( type, order, i, p ); 00564 return libMesh::RealGradient( value ); 00565 } 00566 template <> RealGradient FE<1,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order, 00567 const unsigned int i, const unsigned int j, 00568 const Point& p) 00569 { 00570 Real value = FE<1,LAGRANGE>::shape_deriv( type, order, i, j, p ); 00571 return libMesh::RealGradient( value ); 00572 } 00573 template <> RealGradient FE<1,LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order, 00574 const unsigned int i, const unsigned int j, 00575 const Point& p) 00576 { 00577 Real value = FE<1,LAGRANGE>::shape_second_deriv( type, order, i, j, p ); 00578 return libMesh::RealGradient( value ); 00579 } 00580 00581 // 2-D 00582 template <> RealGradient FE<2,LAGRANGE_VEC>::shape(const ElemType type, const Order order, 00583 const unsigned int i, const Point& p) 00584 { 00585 Real value = FE<2,LAGRANGE>::shape( type, order, i/2, p ); 00586 00587 switch( i%2 ) 00588 { 00589 case 0: 00590 return libMesh::RealGradient( value ); 00591 00592 case 1: 00593 return libMesh::RealGradient( Real(0), value ); 00594 00595 default: 00596 libmesh_error_msg("i%2 must be either 0 or 1!"); 00597 } 00598 00599 //dummy 00600 return libMesh::RealGradient(); 00601 } 00602 template <> RealGradient FE<2,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order, 00603 const unsigned int i, const unsigned int j, 00604 const Point& p) 00605 { 00606 Real value = FE<2,LAGRANGE>::shape_deriv( type, order, i/2, j, p ); 00607 00608 switch( i%2 ) 00609 { 00610 case 0: 00611 return libMesh::RealGradient( value ); 00612 00613 case 1: 00614 return libMesh::RealGradient( Real(0), value ); 00615 00616 default: 00617 libmesh_error_msg("i%2 must be either 0 or 1!"); 00618 } 00619 00620 //dummy 00621 return libMesh::RealGradient(); 00622 } 00623 template <> RealGradient FE<2,LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order, 00624 const unsigned int i, const unsigned int j, 00625 const Point& p) 00626 { 00627 Real value = FE<2,LAGRANGE>::shape_second_deriv( type, order, i/2, j, p ); 00628 00629 switch( i%2 ) 00630 { 00631 case 0: 00632 return libMesh::RealGradient( value ); 00633 00634 case 1: 00635 return libMesh::RealGradient( Real(0), value ); 00636 00637 default: 00638 libmesh_error_msg("i%2 must be either 0 or 1!"); 00639 } 00640 00641 //dummy 00642 return libMesh::RealGradient(); 00643 } 00644 00645 00646 // 3-D 00647 template <> RealGradient FE<3,LAGRANGE_VEC>::shape(const ElemType type, const Order order, 00648 const unsigned int i, const Point& p) 00649 { 00650 Real value = FE<3,LAGRANGE>::shape( type, order, i/3, p ); 00651 00652 switch( i%3 ) 00653 { 00654 case 0: 00655 return libMesh::RealGradient( value ); 00656 00657 case 1: 00658 return libMesh::RealGradient( Real(0), value ); 00659 00660 case 2: 00661 return libMesh::RealGradient( Real(0), Real(0), value ); 00662 00663 default: 00664 libmesh_error_msg("i%3 must be 0, 1, or 2!"); 00665 } 00666 00667 //dummy 00668 return libMesh::RealGradient(); 00669 } 00670 template <> RealGradient FE<3,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order, 00671 const unsigned int i, const unsigned int j, 00672 const Point& p) 00673 { 00674 Real value = FE<3,LAGRANGE>::shape_deriv( type, order, i/3, j, p ); 00675 00676 switch( i%3 ) 00677 { 00678 case 0: 00679 return libMesh::RealGradient( value ); 00680 00681 case 1: 00682 return libMesh::RealGradient( Real(0), value ); 00683 00684 case 2: 00685 return libMesh::RealGradient( Real(0), Real(0), value ); 00686 00687 default: 00688 libmesh_error_msg("i%3 must be 0, 1, or 2!"); 00689 } 00690 00691 //dummy 00692 return libMesh::RealGradient(); 00693 } 00694 template <> RealGradient FE<3,LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order, 00695 const unsigned int i, const unsigned int j, 00696 const Point& p) 00697 { 00698 Real value = FE<3,LAGRANGE>::shape_second_deriv( type, order, i/3, j, p ); 00699 00700 switch( i%3 ) 00701 { 00702 case 0: 00703 return libMesh::RealGradient( value ); 00704 00705 case 1: 00706 return libMesh::RealGradient( Real(0), value ); 00707 00708 case 2: 00709 return libMesh::RealGradient( Real(0), Real(0), value ); 00710 00711 default: 00712 libmesh_error_msg("i%3 must be 0, 1, or 2!"); 00713 } 00714 00715 //dummy 00716 return libMesh::RealGradient(); 00717 } 00718 00719 00720 00721 // 0-D 00722 template <> RealGradient FE<0,LAGRANGE_VEC>::shape(const Elem* elem, const Order order, 00723 const unsigned int i, const Point& p) 00724 { 00725 Real value = FE<0,LAGRANGE>::shape( elem->type(), static_cast<Order>(order + elem->p_level()), i, p); 00726 return libMesh::RealGradient( value ); 00727 } 00728 template <> RealGradient FE<0,LAGRANGE_VEC>::shape_deriv(const Elem* elem, const Order order, 00729 const unsigned int i, const unsigned int j, 00730 const Point& p) 00731 { 00732 Real value = FE<0,LAGRANGE>::shape_deriv( elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p); 00733 return libMesh::RealGradient( value ); 00734 } 00735 template <> RealGradient FE<0,LAGRANGE_VEC>::shape_second_deriv(const Elem* elem, const Order order, 00736 const unsigned int i, const unsigned int j, 00737 const Point& p) 00738 { 00739 Real value = FE<0,LAGRANGE>::shape_second_deriv( elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p); 00740 return libMesh::RealGradient( value ); 00741 } 00742 00743 // 1-D 00744 template <> RealGradient FE<1,LAGRANGE_VEC>::shape(const Elem* elem, const Order order, 00745 const unsigned int i, const Point& p) 00746 { 00747 Real value = FE<1,LAGRANGE>::shape( elem->type(), static_cast<Order>(order + elem->p_level()), i, p); 00748 return libMesh::RealGradient( value ); 00749 } 00750 template <> RealGradient FE<1,LAGRANGE_VEC>::shape_deriv(const Elem* elem, const Order order, 00751 const unsigned int i, const unsigned int j, 00752 const Point& p) 00753 { 00754 Real value = FE<1,LAGRANGE>::shape_deriv( elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p); 00755 return libMesh::RealGradient( value ); 00756 } 00757 template <> RealGradient FE<1,LAGRANGE_VEC>::shape_second_deriv(const Elem* elem, const Order order, 00758 const unsigned int i, const unsigned int j, 00759 const Point& p) 00760 { 00761 Real value = FE<1,LAGRANGE>::shape_second_deriv( elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p); 00762 return libMesh::RealGradient( value ); 00763 } 00764 00765 // 2-D 00766 template <> RealGradient FE<2,LAGRANGE_VEC>::shape(const Elem* elem, const Order order, 00767 const unsigned int i, const Point& p) 00768 { 00769 Real value = FE<2,LAGRANGE>::shape( elem->type(), static_cast<Order>(order + elem->p_level()), i/2, p ); 00770 00771 switch( i%2 ) 00772 { 00773 case 0: 00774 return libMesh::RealGradient( value ); 00775 00776 case 1: 00777 return libMesh::RealGradient( Real(0), value ); 00778 00779 default: 00780 libmesh_error_msg("i%2 must be either 0 or 1!"); 00781 } 00782 00783 //dummy 00784 return libMesh::RealGradient(); 00785 } 00786 template <> RealGradient FE<2,LAGRANGE_VEC>::shape_deriv(const Elem* elem, const Order order, 00787 const unsigned int i, const unsigned int j, 00788 const Point& p) 00789 { 00790 Real value = FE<2,LAGRANGE>::shape_deriv( elem->type(), static_cast<Order>(order + elem->p_level()), i/2, j, p ); 00791 00792 switch( i%2 ) 00793 { 00794 case 0: 00795 return libMesh::RealGradient( value ); 00796 00797 case 1: 00798 return libMesh::RealGradient( Real(0), value ); 00799 00800 default: 00801 libmesh_error_msg("i%2 must be either 0 or 1!"); 00802 } 00803 00804 //dummy 00805 return libMesh::RealGradient(); 00806 } 00807 template <> RealGradient FE<2,LAGRANGE_VEC>::shape_second_deriv(const Elem* elem, const Order order, 00808 const unsigned int i, const unsigned int j, 00809 const Point& p) 00810 { 00811 Real value = FE<2,LAGRANGE>::shape_second_deriv( elem->type(), static_cast<Order>(order + elem->p_level()), i/2, j, p ); 00812 00813 switch( i%2 ) 00814 { 00815 case 0: 00816 return libMesh::RealGradient( value ); 00817 00818 case 1: 00819 return libMesh::RealGradient( Real(0), value ); 00820 00821 default: 00822 libmesh_error_msg("i%2 must be either 0 or 1!"); 00823 } 00824 00825 //dummy 00826 return libMesh::RealGradient(); 00827 } 00828 00829 // 3-D 00830 template <> RealGradient FE<3,LAGRANGE_VEC>::shape(const Elem* elem, const Order order, 00831 const unsigned int i, const Point& p) 00832 { 00833 Real value = FE<3,LAGRANGE>::shape( elem->type(), static_cast<Order>(order + elem->p_level()), i/3, p ); 00834 00835 switch( i%3 ) 00836 { 00837 case 0: 00838 return libMesh::RealGradient( value ); 00839 00840 case 1: 00841 return libMesh::RealGradient( Real(0), value ); 00842 00843 case 2: 00844 return libMesh::RealGradient( Real(0), Real(0), value ); 00845 00846 default: 00847 libmesh_error_msg("i%3 must be 0, 1, or 2!"); 00848 } 00849 00850 //dummy 00851 return libMesh::RealGradient(); 00852 } 00853 template <> RealGradient FE<3,LAGRANGE_VEC>::shape_deriv(const Elem* elem, const Order order, 00854 const unsigned int i, const unsigned int j, 00855 const Point& p) 00856 { 00857 Real value = FE<3,LAGRANGE>::shape_deriv( elem->type(), static_cast<Order>(order + elem->p_level()), i/3, j, p ); 00858 00859 switch( i%3 ) 00860 { 00861 case 0: 00862 return libMesh::RealGradient( value ); 00863 00864 case 1: 00865 return libMesh::RealGradient( Real(0), value ); 00866 00867 case 2: 00868 return libMesh::RealGradient( Real(0), Real(0), value ); 00869 00870 default: 00871 libmesh_error_msg("i%3 must be 0, 1, or 2!"); 00872 } 00873 00874 //dummy 00875 return libMesh::RealGradient(); 00876 } 00877 template <> RealGradient FE<3,LAGRANGE_VEC>::shape_second_deriv(const Elem* elem, const Order order, 00878 const unsigned int i, const unsigned int j, 00879 const Point& p) 00880 { 00881 Real value = FE<3,LAGRANGE>::shape_second_deriv( elem->type(), static_cast<Order>(order + elem->p_level()), i/3, j, p ); 00882 00883 switch( i%3 ) 00884 { 00885 case 0: 00886 return libMesh::RealGradient( value ); 00887 00888 case 1: 00889 return libMesh::RealGradient( Real(0), value ); 00890 00891 case 2: 00892 return libMesh::RealGradient( Real(0), Real(0), value ); 00893 00894 default: 00895 libmesh_error_msg("i%3 must be 0, 1, or 2!"); 00896 } 00897 00898 //dummy 00899 return libMesh::RealGradient(); 00900 } 00901 00902 // Do full-specialization for every dimension, instead 00903 // of explicit instantiation at the end of this function. 00904 // This could be macro-ified. 00905 template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<0,LAGRANGE>::n_dofs(t,o); } 00906 template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<1,LAGRANGE>::n_dofs(t,o); } 00907 template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 2*FE<2,LAGRANGE>::n_dofs(t,o); } 00908 template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 3*FE<3,LAGRANGE>::n_dofs(t,o); } 00909 00910 00911 // Do full-specialization for every dimension, instead 00912 // of explicit instantiation at the end of this function. 00913 template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<0,LAGRANGE>::n_dofs_at_node(t,o,n); } 00914 template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<1,LAGRANGE>::n_dofs_at_node(t,o,n); } 00915 template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return 2*FE<2,LAGRANGE>::n_dofs_at_node(t,o,n); } 00916 template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return 3*FE<2,LAGRANGE>::n_dofs_at_node(t,o,n); } 00917 00918 00919 // Lagrange elements have no dofs per element 00920 // (just at the nodes) 00921 template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; } 00922 template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; } 00923 template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; } 00924 template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; } 00925 00926 // Lagrange FEMs are always C^0 continuous 00927 template <> FEContinuity FE<0,LAGRANGE_VEC>::get_continuity() const { return C_ZERO; } 00928 template <> FEContinuity FE<1,LAGRANGE_VEC>::get_continuity() const { return C_ZERO; } 00929 template <> FEContinuity FE<2,LAGRANGE_VEC>::get_continuity() const { return C_ZERO; } 00930 template <> FEContinuity FE<3,LAGRANGE_VEC>::get_continuity() const { return C_ZERO; } 00931 00932 // Lagrange FEMs are not hierarchic 00933 template <> bool FE<0,LAGRANGE_VEC>::is_hierarchic() const { return false; } 00934 template <> bool FE<1,LAGRANGE_VEC>::is_hierarchic() const { return false; } 00935 template <> bool FE<2,LAGRANGE_VEC>::is_hierarchic() const { return false; } 00936 template <> bool FE<3,LAGRANGE_VEC>::is_hierarchic() const { return false; } 00937 00938 // Lagrange FEM shapes do not need reinit (is this always true?) 00939 template <> bool FE<0,LAGRANGE_VEC>::shapes_need_reinit() const { return false; } 00940 template <> bool FE<1,LAGRANGE_VEC>::shapes_need_reinit() const { return false; } 00941 template <> bool FE<2,LAGRANGE_VEC>::shapes_need_reinit() const { return false; } 00942 template <> bool FE<3,LAGRANGE_VEC>::shapes_need_reinit() const { return false; } 00943 00944 // Methods for computing Lagrange constraints. Note: we pass the 00945 // dimension as the last argument to the anonymous helper function. 00946 // Also note: we only need instantiations of this function for 00947 // Dim==2 and 3. 00948 #ifdef LIBMESH_ENABLE_AMR 00949 template <> 00950 void FE<2,LAGRANGE_VEC>::compute_constraints (DofConstraints &constraints, 00951 DofMap &dof_map, 00952 const unsigned int variable_number, 00953 const Elem* elem) 00954 { //libmesh_not_implemented(); 00955 FEVectorBase::compute_proj_constraints(constraints, dof_map, variable_number, elem); 00956 } 00957 00958 template <> 00959 void FE<3,LAGRANGE_VEC>::compute_constraints (DofConstraints &constraints, 00960 DofMap &dof_map, 00961 const unsigned int variable_number, 00962 const Elem* elem) 00963 { //libmesh_not_implemented(); 00964 FEVectorBase::compute_proj_constraints(constraints, dof_map, variable_number, elem); 00965 } 00966 #endif // LIBMESH_ENABLE_AMR 00967 00968 } // namespace libMesh
generated by
