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