libMesh
Loading...
Searching...
No Matches
fe_lagrange_vec.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3
4// This library is free software; you can redistribute it and/or
5// modify it under the terms of the GNU Lesser General Public
6// License as published by the Free Software Foundation; either
7// version 2.1 of the License, or (at your option) any later version.
8
9// This library is distributed in the hope that it will be useful,
10// but WITHOUT ANY WARRANTY; without even the implied warranty of
11// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12// Lesser General Public License for more details.
13
14// You should have received a copy of the GNU Lesser General Public
15// License along with this library; if not, write to the Free Software
16// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17
18
19
20// Local includes
21#include "libmesh/dof_map.h"
22#include "libmesh/elem.h"
23#include "libmesh/enum_to_string.h"
24#include "libmesh/fe.h"
25#include "libmesh/fe_interface.h"
26#include "libmesh/fe_macro.h"
27#include "libmesh/tensor_value.h"
28
29
30namespace libMesh
31{
32
33
34LIBMESH_DEFAULT_VECTORIZED_FE(0,LAGRANGE_VEC)
35LIBMESH_DEFAULT_VECTORIZED_FE(1,LAGRANGE_VEC)
36LIBMESH_DEFAULT_VECTORIZED_FE(2,LAGRANGE_VEC)
37LIBMESH_DEFAULT_VECTORIZED_FE(3,LAGRANGE_VEC)
38LIBMESH_DEFAULT_VECTORIZED_FE(0,L2_LAGRANGE_VEC)
39LIBMESH_DEFAULT_VECTORIZED_FE(1,L2_LAGRANGE_VEC)
40LIBMESH_DEFAULT_VECTORIZED_FE(2,L2_LAGRANGE_VEC)
41LIBMESH_DEFAULT_VECTORIZED_FE(3,L2_LAGRANGE_VEC)
42
43
44// ------------------------------------------------------------
45// Lagrange-specific implementations
46
47
48// Anonymous namespace for local helper functions
49namespace {
50void lagrange_vec_nodal_soln(const Elem * elem,
51 const Order order,
52 const std::vector<Number> & elem_soln,
53 const int dim,
54 std::vector<Number> & nodal_soln,
55 const bool add_p_level)
56{
57 const unsigned int n_nodes = elem->n_nodes();
58 const ElemType type = elem->type();
59
60 const Order totalorder = order + add_p_level*elem->p_level();
61
62 nodal_soln.resize(dim*n_nodes);
63
64 switch (totalorder)
65 {
66 // linear Lagrange shape functions
67 case FIRST:
68 {
69 switch (type)
70 {
71 case TRI7:
72 libmesh_assert_equal_to (nodal_soln.size(), 14);
73 nodal_soln[12] = (elem_soln[0] + elem_soln[2] + elem_soln[4])/3.;
74 nodal_soln[13] = (elem_soln[1] + elem_soln[3] + elem_soln[5])/3.;
75 libmesh_fallthrough();
76 case TRI6:
77 {
78 libmesh_assert (type == TRI7 || nodal_soln.size() == 12);
79 libmesh_assert_equal_to (elem_soln.size(), 2*3);
80
81 // node 0 components
82 nodal_soln[0] = elem_soln[0];
83 nodal_soln[1] = elem_soln[1];
84
85 // node 1 components
86 nodal_soln[2] = elem_soln[2];
87 nodal_soln[3] = elem_soln[3];
88
89 // node 2 components
90 nodal_soln[4] = elem_soln[4];
91 nodal_soln[5] = elem_soln[5];
92
93 // node 3 components
94 nodal_soln[6] = .5*(elem_soln[0] + elem_soln[2]);
95 nodal_soln[7] = .5*(elem_soln[1] + elem_soln[3]);
96
97 // node 4 components
98 nodal_soln[8] = .5*(elem_soln[2] + elem_soln[4]);
99 nodal_soln[9] = .5*(elem_soln[3] + elem_soln[5]);
100
101 // node 5 components
102 nodal_soln[10] = .5*(elem_soln[0] + elem_soln[4]);
103 nodal_soln[11] = .5*(elem_soln[1] + elem_soln[5]);
104
105 return;
106 }
107
108
109 case QUAD8:
110 case QUAD9:
111 {
112 libmesh_assert_equal_to (elem_soln.size(), 2*4);
113
114 if (type == QUAD8)
115 libmesh_assert_equal_to (nodal_soln.size(), 2*8);
116 else
117 libmesh_assert_equal_to (nodal_soln.size(), 2*9);
118
119 // node 0 components
120 nodal_soln[0] = elem_soln[0];
121 nodal_soln[1] = elem_soln[1];
122
123 // node 1 components
124 nodal_soln[2] = elem_soln[2];
125 nodal_soln[3] = elem_soln[3];
126
127 // node 2 components
128 nodal_soln[4] = elem_soln[4];
129 nodal_soln[5] = elem_soln[5];
130
131 // node 3 components
132 nodal_soln[6] = elem_soln[6];
133 nodal_soln[7] = elem_soln[7];
134
135 // node 4 components
136 nodal_soln[8] = .5*(elem_soln[0] + elem_soln[2]);
137 nodal_soln[9] = .5*(elem_soln[1] + elem_soln[3]);
138
139 // node 5 components
140 nodal_soln[10] = .5*(elem_soln[2] + elem_soln[4]);
141 nodal_soln[11] = .5*(elem_soln[3] + elem_soln[5]);
142
143 // node 6 components
144 nodal_soln[12] = .5*(elem_soln[4] + elem_soln[6]);
145 nodal_soln[13] = .5*(elem_soln[5] + elem_soln[7]);
146
147 // node 7 components
148 nodal_soln[14] = .5*(elem_soln[6] + elem_soln[0]);
149 nodal_soln[15] = .5*(elem_soln[7] + elem_soln[1]);
150
151 if (type == QUAD9)
152 {
153 // node 8 components
154 nodal_soln[16] = .25*(elem_soln[0] + elem_soln[2] + elem_soln[4] + elem_soln[6]);
155 nodal_soln[17] = .25*(elem_soln[1] + elem_soln[3] + elem_soln[5] + elem_soln[7]);
156 }
157
158 return;
159 }
160
161
162 case TET14:
163 libmesh_assert_equal_to (nodal_soln.size(), 3*14);
164
165 // node 10 components
166 nodal_soln[30] = 1./3. * (elem_soln[0] + elem_soln[3] + elem_soln[6]);
167 nodal_soln[31] = 1./3. * (elem_soln[1] + elem_soln[4] + elem_soln[7]);
168 nodal_soln[32] = 1./3. * (elem_soln[2] + elem_soln[5] + elem_soln[8]);
169
170 // node 11 components
171 nodal_soln[33] = 1./3. * (elem_soln[0] + elem_soln[3] + elem_soln[9]);
172 nodal_soln[34] = 1./3. * (elem_soln[1] + elem_soln[4] + elem_soln[10]);
173 nodal_soln[35] = 1./3. * (elem_soln[2] + elem_soln[5] + elem_soln[11]);
174
175 // node 12 components
176 nodal_soln[36] = 1./3. * (elem_soln[3] + elem_soln[6] + elem_soln[9]);
177 nodal_soln[37] = 1./3. * (elem_soln[4] + elem_soln[7] + elem_soln[10]);
178 nodal_soln[38] = 1./3. * (elem_soln[5] + elem_soln[8] + elem_soln[11]);
179
180 // node 13 components
181 nodal_soln[39] = 1./3. * (elem_soln[0] + elem_soln[6] + elem_soln[9]);
182 nodal_soln[40] = 1./3. * (elem_soln[1] + elem_soln[7] + elem_soln[10]);
183 nodal_soln[41] = 1./3. * (elem_soln[2] + elem_soln[8] + elem_soln[11]);
184
185 libmesh_fallthrough();
186 case TET10:
187 {
188 libmesh_assert (type == TET14 || nodal_soln.size() == 3*10);
189 libmesh_assert_equal_to (elem_soln.size(), 3*4);
190
191 // node 0 components
192 nodal_soln[0] = elem_soln[0];
193 nodal_soln[1] = elem_soln[1];
194 nodal_soln[2] = elem_soln[2];
195
196 // node 1 components
197 nodal_soln[3] = elem_soln[3];
198 nodal_soln[4] = elem_soln[4];
199 nodal_soln[5] = elem_soln[5];
200
201 // node 2 components
202 nodal_soln[6] = elem_soln[6];
203 nodal_soln[7] = elem_soln[7];
204 nodal_soln[8] = elem_soln[8];
205
206 // node 3 components
207 nodal_soln[9] = elem_soln[9];
208 nodal_soln[10] = elem_soln[10];
209 nodal_soln[11] = elem_soln[11];
210
211 // node 4 components
212 nodal_soln[12] = .5*(elem_soln[0] + elem_soln[3]);
213 nodal_soln[13] = .5*(elem_soln[1] + elem_soln[4]);
214 nodal_soln[14] = .5*(elem_soln[2] + elem_soln[5]);
215
216 // node 5 components
217 nodal_soln[15] = .5*(elem_soln[3] + elem_soln[6]);
218 nodal_soln[16] = .5*(elem_soln[4] + elem_soln[7]);
219 nodal_soln[17] = .5*(elem_soln[5] + elem_soln[8]);
220
221 // node 6 components
222 nodal_soln[18] = .5*(elem_soln[6] + elem_soln[0]);
223 nodal_soln[19] = .5*(elem_soln[7] + elem_soln[1]);
224 nodal_soln[20] = .5*(elem_soln[8] + elem_soln[2]);
225
226 // node 7 components
227 nodal_soln[21] = .5*(elem_soln[9] + elem_soln[0]);
228 nodal_soln[22] = .5*(elem_soln[10] + elem_soln[1]);
229 nodal_soln[23] = .5*(elem_soln[11] + elem_soln[2]);
230
231 // node 8 components
232 nodal_soln[24] = .5*(elem_soln[9] + elem_soln[3]);
233 nodal_soln[25] = .5*(elem_soln[10] + elem_soln[4]);
234 nodal_soln[26] = .5*(elem_soln[11] + elem_soln[5]);
235
236 // node 9 components
237 nodal_soln[27] = .5*(elem_soln[9] + elem_soln[6]);
238 nodal_soln[28] = .5*(elem_soln[10] + elem_soln[7]);
239 nodal_soln[29] = .5*(elem_soln[11] + elem_soln[8]);
240
241 return;
242 }
243
244
245 case HEX20:
246 case HEX27:
247 {
248 libmesh_assert_equal_to (elem_soln.size(), 3*8);
249
250 if (type == HEX20)
251 libmesh_assert_equal_to (nodal_soln.size(), 3*20);
252 else
253 libmesh_assert_equal_to (nodal_soln.size(), 3*27);
254
255 // node 0 components
256 nodal_soln[0] = elem_soln[0];
257 nodal_soln[1] = elem_soln[1];
258 nodal_soln[2] = elem_soln[2];
259
260 // node 1 components
261 nodal_soln[3] = elem_soln[3];
262 nodal_soln[4] = elem_soln[4];
263 nodal_soln[5] = elem_soln[5];
264
265 // node 2 components
266 nodal_soln[6] = elem_soln[6];
267 nodal_soln[7] = elem_soln[7];
268 nodal_soln[8] = elem_soln[8];
269
270 // node 3 components
271 nodal_soln[9] = elem_soln[9];
272 nodal_soln[10] = elem_soln[10];
273 nodal_soln[11] = elem_soln[11];
274
275 // node 4 components
276 nodal_soln[12] = elem_soln[12];
277 nodal_soln[13] = elem_soln[13];
278 nodal_soln[14] = elem_soln[14];
279
280 // node 5 components
281 nodal_soln[15] = elem_soln[15];
282 nodal_soln[16] = elem_soln[16];
283 nodal_soln[17] = elem_soln[17];
284
285 // node 6 components
286 nodal_soln[18] = elem_soln[18];
287 nodal_soln[19] = elem_soln[19];
288 nodal_soln[20] = elem_soln[20];
289
290 // node 7 components
291 nodal_soln[21] = elem_soln[21];
292 nodal_soln[22] = elem_soln[22];
293 nodal_soln[23] = elem_soln[23];
294
295 // node 8 components
296 nodal_soln[24] = .5*(elem_soln[0] + elem_soln[3]);
297 nodal_soln[25] = .5*(elem_soln[1] + elem_soln[4]);
298 nodal_soln[26] = .5*(elem_soln[2] + elem_soln[5]);
299
300 // node 9 components
301 nodal_soln[27] = .5*(elem_soln[3] + elem_soln[6]);
302 nodal_soln[28] = .5*(elem_soln[4] + elem_soln[7]);
303 nodal_soln[29] = .5*(elem_soln[5] + elem_soln[8]);
304
305 // node 10 components
306 nodal_soln[30] = .5*(elem_soln[6] + elem_soln[9]);
307 nodal_soln[31] = .5*(elem_soln[7] + elem_soln[10]);
308 nodal_soln[32] = .5*(elem_soln[8] + elem_soln[11]);
309
310 // node 11 components
311 nodal_soln[33] = .5*(elem_soln[9] + elem_soln[0]);
312 nodal_soln[34] = .5*(elem_soln[10] + elem_soln[1]);
313 nodal_soln[35] = .5*(elem_soln[11] + elem_soln[2]);
314
315 // node 12 components
316 nodal_soln[36] = .5*(elem_soln[0] + elem_soln[12]);
317 nodal_soln[37] = .5*(elem_soln[1] + elem_soln[13]);
318 nodal_soln[38] = .5*(elem_soln[2] + elem_soln[14]);
319
320 // node 13 components
321 nodal_soln[39] = .5*(elem_soln[3] + elem_soln[15]);
322 nodal_soln[40] = .5*(elem_soln[4] + elem_soln[16]);
323 nodal_soln[41] = .5*(elem_soln[5] + elem_soln[17]);
324
325 // node 14 components
326 nodal_soln[42] = .5*(elem_soln[6] + elem_soln[18]);
327 nodal_soln[43] = .5*(elem_soln[7] + elem_soln[19]);
328 nodal_soln[44] = .5*(elem_soln[8] + elem_soln[20]);
329
330 // node 15 components
331 nodal_soln[45] = .5*(elem_soln[9] + elem_soln[21]);
332 nodal_soln[46] = .5*(elem_soln[10] + elem_soln[22]);
333 nodal_soln[47] = .5*(elem_soln[11] + elem_soln[23]);
334
335 // node 16 components
336 nodal_soln[48] = .5*(elem_soln[12] + elem_soln[15]);
337 nodal_soln[49] = .5*(elem_soln[13] + elem_soln[16]);
338 nodal_soln[50] = .5*(elem_soln[14] + elem_soln[17]);
339
340 // node 17 components
341 nodal_soln[51] = .5*(elem_soln[15] + elem_soln[18]);
342 nodal_soln[52] = .5*(elem_soln[16] + elem_soln[19]);
343 nodal_soln[53] = .5*(elem_soln[17] + elem_soln[20]);
344
345 // node 18 components
346 nodal_soln[54] = .5*(elem_soln[18] + elem_soln[21]);
347 nodal_soln[55] = .5*(elem_soln[19] + elem_soln[22]);
348 nodal_soln[56] = .5*(elem_soln[20] + elem_soln[23]);
349
350 // node 19 components
351 nodal_soln[57] = .5*(elem_soln[12] + elem_soln[21]);
352 nodal_soln[58] = .5*(elem_soln[13] + elem_soln[22]);
353 nodal_soln[59] = .5*(elem_soln[14] + elem_soln[23]);
354
355 if (type == HEX27)
356 {
357 // node 20 components
358 nodal_soln[60] = .25*(elem_soln[0] + elem_soln[3] + elem_soln[6] + elem_soln[9]);
359 nodal_soln[61] = .25*(elem_soln[1] + elem_soln[4] + elem_soln[7] + elem_soln[10]);
360 nodal_soln[62] = .25*(elem_soln[2] + elem_soln[5] + elem_soln[8] + elem_soln[11]);
361
362 // node 21 components
363 nodal_soln[63] = .25*(elem_soln[0] + elem_soln[3] + elem_soln[12] + elem_soln[15]);
364 nodal_soln[64] = .25*(elem_soln[1] + elem_soln[4] + elem_soln[13] + elem_soln[16]);
365 nodal_soln[65] = .25*(elem_soln[2] + elem_soln[5] + elem_soln[14] + elem_soln[17]);
366
367 // node 22 components
368 nodal_soln[66] = .25*(elem_soln[3] + elem_soln[6] + elem_soln[15] + elem_soln[18]);
369 nodal_soln[67] = .25*(elem_soln[4] + elem_soln[7] + elem_soln[16] + elem_soln[19]);
370 nodal_soln[68] = .25*(elem_soln[5] + elem_soln[8] + elem_soln[17] + elem_soln[20]);
371
372 // node 23 components
373 nodal_soln[69] = .25*(elem_soln[6] + elem_soln[9] + elem_soln[18] + elem_soln[21]);
374 nodal_soln[70] = .25*(elem_soln[7] + elem_soln[10] + elem_soln[19] + elem_soln[22]);
375 nodal_soln[71] = .25*(elem_soln[8] + elem_soln[11] + elem_soln[20] + elem_soln[23]);
376
377 // node 24 components
378 nodal_soln[72] = .25*(elem_soln[9] + elem_soln[0] + elem_soln[21] + elem_soln[12]);
379 nodal_soln[73] = .25*(elem_soln[10] + elem_soln[1] + elem_soln[22] + elem_soln[13]);
380 nodal_soln[74] = .25*(elem_soln[11] + elem_soln[2] + elem_soln[23] + elem_soln[14]);
381
382 // node 25 components
383 nodal_soln[75] = .25*(elem_soln[12] + elem_soln[15] + elem_soln[18] + elem_soln[21]);
384 nodal_soln[76] = .25*(elem_soln[13] + elem_soln[16] + elem_soln[19] + elem_soln[22]);
385 nodal_soln[77] = .25*(elem_soln[14] + elem_soln[17] + elem_soln[20] + elem_soln[23]);
386
387 // node 26 components
388 nodal_soln[78] = .125*(elem_soln[0] + elem_soln[3] + elem_soln[6] + elem_soln[9] +
389 elem_soln[12] + elem_soln[15] + elem_soln[18] + elem_soln[21]);
390
391 nodal_soln[79] = .125*(elem_soln[1] + elem_soln[4] + elem_soln[7] + elem_soln[10] +
392 elem_soln[13] + elem_soln[16] + elem_soln[19] + elem_soln[22]);
393
394 nodal_soln[80] = .125*(elem_soln[2] + elem_soln[5] + elem_soln[8] + elem_soln[11] +
395 elem_soln[14] + elem_soln[17] + elem_soln[20] + elem_soln[23]);
396 }
397
398 return;
399 }
400
401
402 case PRISM21:
403 libmesh_assert_equal_to (nodal_soln.size(), 3*21);
404
405 // node 20 components
406 nodal_soln[60] = (elem_soln[27] + elem_soln[30] + elem_soln[33])/Real(3);
407 nodal_soln[61] = (elem_soln[28] + elem_soln[31] + elem_soln[34])/Real(3);
408 nodal_soln[62] = (elem_soln[29] + elem_soln[32] + elem_soln[35])/Real(3);
409 libmesh_fallthrough();
410 case PRISM20:
411 if (type == PRISM20)
412 libmesh_assert_equal_to (nodal_soln.size(), 3*20);
413
414 // node 18 components
415 nodal_soln[54] = (elem_soln[0] + elem_soln[3] + elem_soln[6])/Real(3);
416 nodal_soln[55] = (elem_soln[1] + elem_soln[4] + elem_soln[7])/Real(3);
417 nodal_soln[56] = (elem_soln[2] + elem_soln[5] + elem_soln[8])/Real(3);
418
419 // node 19 components
420 nodal_soln[57] = (elem_soln[9] + elem_soln[12] + elem_soln[15])/Real(3);
421 nodal_soln[58] = (elem_soln[10] + elem_soln[13] + elem_soln[16])/Real(3);
422 nodal_soln[59] = (elem_soln[11] + elem_soln[14] + elem_soln[17])/Real(3);
423
424 libmesh_fallthrough();
425 case PRISM18:
426 if (type == PRISM18)
427 libmesh_assert_equal_to (nodal_soln.size(), 3*18);
428
429 // node 15 components
430 nodal_soln[45] = .25*(elem_soln[0] + elem_soln[3] + elem_soln[12] + elem_soln[9]);
431 nodal_soln[46] = .25*(elem_soln[1] + elem_soln[4] + elem_soln[13] + elem_soln[10]);
432 nodal_soln[47] = .25*(elem_soln[2] + elem_soln[5] + elem_soln[14] + elem_soln[11]);
433
434 // node 16 components
435 nodal_soln[48] = .25*(elem_soln[3] + elem_soln[6] + elem_soln[15] + elem_soln[12]);
436 nodal_soln[49] = .25*(elem_soln[4] + elem_soln[7] + elem_soln[16] + elem_soln[13]);
437 nodal_soln[50] = .25*(elem_soln[5] + elem_soln[8] + elem_soln[17] + elem_soln[14]);
438
439 // node 17 components
440 nodal_soln[51] = .25*(elem_soln[6] + elem_soln[0] + elem_soln[9] + elem_soln[15]);
441 nodal_soln[52] = .25*(elem_soln[7] + elem_soln[1] + elem_soln[10] + elem_soln[16]);
442 nodal_soln[53] = .25*(elem_soln[8] + elem_soln[2] + elem_soln[11] + elem_soln[17]);
443
444 libmesh_fallthrough();
445 case PRISM15:
446 {
447 libmesh_assert_equal_to (elem_soln.size(), 3*6);
448
449 if (type == PRISM15)
450 libmesh_assert_equal_to (nodal_soln.size(), 3*15);
451
452 // node 0 components
453 nodal_soln[0] = elem_soln[0];
454 nodal_soln[1] = elem_soln[1];
455 nodal_soln[2] = elem_soln[2];
456
457 // node 1 components
458 nodal_soln[3] = elem_soln[3];
459 nodal_soln[4] = elem_soln[4];
460 nodal_soln[5] = elem_soln[5];
461
462 // node 2 components
463 nodal_soln[6] = elem_soln[6];
464 nodal_soln[7] = elem_soln[7];
465 nodal_soln[8] = elem_soln[8];
466
467 // node 3 components
468 nodal_soln[9] = elem_soln[9];
469 nodal_soln[10] = elem_soln[10];
470 nodal_soln[11] = elem_soln[11];
471
472 // node 4 components
473 nodal_soln[12] = elem_soln[12];
474 nodal_soln[13] = elem_soln[13];
475 nodal_soln[14] = elem_soln[14];
476
477 // node 5 components
478 nodal_soln[15] = elem_soln[15];
479 nodal_soln[16] = elem_soln[16];
480 nodal_soln[17] = elem_soln[17];
481
482 // node 6 components
483 nodal_soln[18] = .5*(elem_soln[0] + elem_soln[3]);
484 nodal_soln[19] = .5*(elem_soln[1] + elem_soln[4]);
485 nodal_soln[20] = .5*(elem_soln[2] + elem_soln[5]);
486
487 // node 7 components
488 nodal_soln[21] = .5*(elem_soln[3] + elem_soln[6]);
489 nodal_soln[22] = .5*(elem_soln[4] + elem_soln[7]);
490 nodal_soln[23] = .5*(elem_soln[5] + elem_soln[8]);
491
492 // node 8 components
493 nodal_soln[24] = .5*(elem_soln[0] + elem_soln[6]);
494 nodal_soln[25] = .5*(elem_soln[1] + elem_soln[7]);
495 nodal_soln[26] = .5*(elem_soln[2] + elem_soln[8]);
496
497 // node 9 components
498 nodal_soln[27] = .5*(elem_soln[0] + elem_soln[9]);
499 nodal_soln[28] = .5*(elem_soln[1] + elem_soln[10]);
500 nodal_soln[29] = .5*(elem_soln[2] + elem_soln[11]);
501
502 // node 10 components
503 nodal_soln[30] = .5*(elem_soln[3] + elem_soln[12]);
504 nodal_soln[31] = .5*(elem_soln[4] + elem_soln[13]);
505 nodal_soln[32] = .5*(elem_soln[5] + elem_soln[14]);
506
507 // node 11 components
508 nodal_soln[33] = .5*(elem_soln[6] + elem_soln[15]);
509 nodal_soln[34] = .5*(elem_soln[7] + elem_soln[16]);
510 nodal_soln[35] = .5*(elem_soln[8] + elem_soln[17]);
511
512 // node 12 components
513 nodal_soln[36] = .5*(elem_soln[9] + elem_soln[12]);
514 nodal_soln[37] = .5*(elem_soln[10] + elem_soln[13]);
515 nodal_soln[38] = .5*(elem_soln[11] + elem_soln[14]);
516
517 // node 13 components
518 nodal_soln[39] = .5*(elem_soln[12] + elem_soln[15]);
519 nodal_soln[40] = .5*(elem_soln[13] + elem_soln[16]);
520 nodal_soln[41] = .5*(elem_soln[14] + elem_soln[17]);
521
522 // node 14 components
523 nodal_soln[42] = .5*(elem_soln[12] + elem_soln[15]);
524 nodal_soln[43] = .5*(elem_soln[13] + elem_soln[16]);
525 nodal_soln[44] = .5*(elem_soln[14] + elem_soln[17]);
526
527 return;
528 }
529
530 default:
531 {
532 // By default the element solution _is_ nodal,
533 // so just copy it.
534 nodal_soln = elem_soln;
535
536 return;
537 }
538 }
539 }
540
541 case SECOND:
542 {
543 switch (type)
544 {
545 case TRI7:
546 {
547 libmesh_assert_equal_to (elem_soln.size(), 12);
548 libmesh_assert_equal_to (nodal_soln.size(), 14);
549
550 for (int i=0; i != 12; ++i)
551 nodal_soln[i] = elem_soln[i];
552
553 nodal_soln[12] = -1./9. * (elem_soln[0] + elem_soln[2] + elem_soln[4])
554 +4./9. * (elem_soln[6] + elem_soln[8] + elem_soln[10]);
555 nodal_soln[13] = -1./9. * (elem_soln[1] + elem_soln[3] + elem_soln[5])
556 +4./9. * (elem_soln[7] + elem_soln[9] + elem_soln[11]);
557
558 return;
559 }
560
561 case TET14:
562 {
563 libmesh_assert_equal_to (elem_soln.size(), 10*3);
564 libmesh_assert_equal_to (nodal_soln.size(), 14*3);
565
566 for (int i=0; i != 10*3; ++i)
567 nodal_soln[i] = elem_soln[i];
568
569 // node 10 components
570 nodal_soln[30] = -1./9. * (elem_soln[0] + elem_soln[3] + elem_soln[6])
571 +4./9. * (elem_soln[12] + elem_soln[15] + elem_soln[18]);
572 nodal_soln[31] = -1./9. * (elem_soln[1] + elem_soln[4] + elem_soln[7])
573 +4./9. * (elem_soln[13] + elem_soln[16] + elem_soln[19]);
574 nodal_soln[32] = -1./9. * (elem_soln[2] + elem_soln[5] + elem_soln[8])
575 +4./9. * (elem_soln[14] + elem_soln[17] + elem_soln[20]);
576
577 // node 11 components
578 nodal_soln[33] = -1./9. * (elem_soln[0] + elem_soln[3] + elem_soln[9])
579 +4./9. * (elem_soln[12] + elem_soln[21] + elem_soln[24]);
580 nodal_soln[34] = -1./9. * (elem_soln[1] + elem_soln[4] + elem_soln[10])
581 +4./9. * (elem_soln[13] + elem_soln[22] + elem_soln[25]);
582 nodal_soln[35] = -1./9. * (elem_soln[2] + elem_soln[5] + elem_soln[11])
583 +4./9. * (elem_soln[14] + elem_soln[23] + elem_soln[26]);
584
585 // node 12 components
586 nodal_soln[36] = -1./9. * (elem_soln[3] + elem_soln[6] + elem_soln[9])
587 +4./9. * (elem_soln[15] + elem_soln[24] + elem_soln[27]);
588 nodal_soln[37] = -1./9. * (elem_soln[4] + elem_soln[7] + elem_soln[10])
589 +4./9. * (elem_soln[16] + elem_soln[25] + elem_soln[28]);
590 nodal_soln[38] = -1./9. * (elem_soln[5] + elem_soln[8] + elem_soln[11])
591 +4./9. * (elem_soln[17] + elem_soln[26] + elem_soln[29]);
592
593 // node 13 components
594 nodal_soln[39] = -1./9. * (elem_soln[0] + elem_soln[6] + elem_soln[9])
595 +4./9. * (elem_soln[12] + elem_soln[21] + elem_soln[27]);
596 nodal_soln[40] = -1./9. * (elem_soln[1] + elem_soln[7] + elem_soln[10])
597 +4./9. * (elem_soln[13] + elem_soln[22] + elem_soln[28]);
598 nodal_soln[41] = -1./9. * (elem_soln[2] + elem_soln[8] + elem_soln[11])
599 +4./9. * (elem_soln[14] + elem_soln[23] + elem_soln[29]);
600
601 return;
602 }
603
604 default:
605 {
606 // By default the element solution _is_ nodal,
607 // so just copy it.
608 nodal_soln = elem_soln;
609
610 return;
611 }
612 }
613 }
614
615 case THIRD:
616 {
617 // By default the element solution _is_ nodal,
618 // so just copy it.
619 nodal_soln = elem_soln;
620
621 return;
622 }
623
624 default:
625 libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(totalorder) << " selected for LAGRANGE FE family!");
626 } // switch(totalorder)
627
628}// void lagrange_vec_nodal_soln
629
630} // anonymous namespace
631
632
633 // Do full-specialization for every dimension, instead
634 // of explicit instantiation at the end of this file.
635 // This could be macro-ified so that it fits on one line...
636LIBMESH_FE_NODAL_SOLN_DIM(LAGRANGE_VEC, (FE<0, LAGRANGE>::nodal_soln), 0)
637LIBMESH_FE_NODAL_SOLN_DIM(LAGRANGE_VEC, (FE<1, LAGRANGE>::nodal_soln), 1)
638
639template <>
640void FE<2,LAGRANGE_VEC>::nodal_soln(const Elem * elem,
641 const Order order,
642 const std::vector<Number> & elem_soln,
643 std::vector<Number> & nodal_soln,
644 const bool add_p_level,
645 const unsigned)
646{ lagrange_vec_nodal_soln(elem, order, elem_soln, 2 /*dim*/, nodal_soln, add_p_level); }
647
648template <>
649void FE<3,LAGRANGE_VEC>::nodal_soln(const Elem * elem,
650 const Order order,
651 const std::vector<Number> & elem_soln,
652 std::vector<Number> & nodal_soln,
653 const bool add_p_level,
654 const unsigned)
655{ lagrange_vec_nodal_soln(elem, order, elem_soln, 3 /*dim*/, nodal_soln, add_p_level); }
656
657LIBMESH_FE_SIDE_NODAL_SOLN(LAGRANGE_VEC)
658
659LIBMESH_FE_NODAL_SOLN_DIM(L2_LAGRANGE_VEC, (FE<0, LAGRANGE_VEC>::nodal_soln), 0)
660LIBMESH_FE_NODAL_SOLN_DIM(L2_LAGRANGE_VEC, (FE<1, LAGRANGE_VEC>::nodal_soln), 1)
661LIBMESH_FE_NODAL_SOLN_DIM(L2_LAGRANGE_VEC, (FE<2, LAGRANGE_VEC>::nodal_soln), 2)
662LIBMESH_FE_NODAL_SOLN_DIM(L2_LAGRANGE_VEC, (FE<3, LAGRANGE_VEC>::nodal_soln), 3)
663
664LIBMESH_FE_SIDE_NODAL_SOLN(L2_LAGRANGE_VEC)
665
666
667// Specialize for shape function routines by leveraging scalar LAGRANGE elements
668
669// 0-D
670template <> RealGradient FE<0,LAGRANGE_VEC>::shape(const ElemType type, const Order order,
671 const unsigned int i, const Point & p)
672{
673 Real value = FE<0,LAGRANGE>::shape( type, order, i, p );
674 return libMesh::RealGradient( value );
675}
676template <> RealGradient FE<0,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
677 const unsigned int i, const unsigned int j,
678 const Point & p)
679{
680 Real value = FE<0,LAGRANGE>::shape_deriv( type, order, i, j, p );
681 return libMesh::RealGradient( value );
682}
683
684#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
685
686template <> RealGradient FE<0,LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order,
687 const unsigned int i, const unsigned int j,
688 const Point & p)
689{
690 Real value = FE<0,LAGRANGE>::shape_second_deriv( type, order, i, j, p );
691 return libMesh::RealGradient( value );
692}
693#endif
694
695template <> RealGradient FE<0,L2_LAGRANGE_VEC>::shape(const ElemType type, const Order order,
696 const unsigned int i, const Point & p)
697{
698 return FE<0,LAGRANGE_VEC>::shape(type, order, i, p);
699}
700template <> RealGradient FE<0,L2_LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
701 const unsigned int i, const unsigned int j,
702 const Point & p)
703{
704 return FE<0,LAGRANGE_VEC>::shape_deriv(type, order, i, j, p);
705}
706
707#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
708
709template <> RealGradient FE<0,L2_LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order,
710 const unsigned int i, const unsigned int j,
711 const Point & p)
712{
713 return FE<0,LAGRANGE_VEC>::shape_second_deriv(type, order, i, j, p);
714}
715#endif
716
717// 1-D
718template <> RealGradient FE<1,LAGRANGE_VEC>::shape(const ElemType type, const Order order,
719 const unsigned int i, const Point & p)
720{
721 Real value = FE<1,LAGRANGE>::shape( type, order, i, p );
722 return libMesh::RealGradient( value );
723}
724template <> RealGradient FE<1,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
725 const unsigned int i, const unsigned int j,
726 const Point & p)
727{
728 Real value = FE<1,LAGRANGE>::shape_deriv( type, order, i, j, p );
729 return libMesh::RealGradient( value );
730}
731#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
732template <> RealGradient FE<1,LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order,
733 const unsigned int i, const unsigned int j,
734 const Point & p)
735{
736 Real value = FE<1,LAGRANGE>::shape_second_deriv( type, order, i, j, p );
737 return libMesh::RealGradient( value );
738}
739
740#endif
741
742template <> RealGradient FE<1,L2_LAGRANGE_VEC>::shape(const ElemType type, const Order order,
743 const unsigned int i, const Point & p)
744{
745 return FE<1,LAGRANGE_VEC>::shape(type, order, i, p);
746}
747template <> RealGradient FE<1,L2_LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
748 const unsigned int i, const unsigned int j,
749 const Point & p)
750{
751 return FE<1,LAGRANGE_VEC>::shape_deriv(type, order, i, j, p);
752}
753
754#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
755
756template <> RealGradient FE<1,L2_LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order,
757 const unsigned int i, const unsigned int j,
758 const Point & p)
759{
760 return FE<1,LAGRANGE_VEC>::shape_second_deriv(type, order, i, j, p);
761}
762#endif
763
764// 2-D
765template <> RealGradient FE<2,LAGRANGE_VEC>::shape(const ElemType type, const Order order,
766 const unsigned int i, const Point & p)
767{
768 Real value = FE<2,LAGRANGE>::shape( type, order, i/2, p );
769
770 switch( i%2 )
771 {
772 case 0:
773 return libMesh::RealGradient( value );
774
775 case 1:
776 return libMesh::RealGradient( Real(0), value );
777
778 default:
779 libmesh_error_msg("i%2 must be either 0 or 1!");
780 }
781
782 //dummy
783 return libMesh::RealGradient();
784}
785template <> RealGradient FE<2,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
786 const unsigned int i, const unsigned int j,
787 const Point & p)
788{
789 Real value = FE<2,LAGRANGE>::shape_deriv( type, order, i/2, j, p );
790
791 switch( i%2 )
792 {
793 case 0:
794 return libMesh::RealGradient( value );
795
796 case 1:
797 return libMesh::RealGradient( Real(0), value );
798
799 default:
800 libmesh_error_msg("i%2 must be either 0 or 1!");
801 }
802
803 //dummy
804 return libMesh::RealGradient();
805}
806
807#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
808template <> RealGradient FE<2,LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order,
809 const unsigned int i, const unsigned int j,
810 const Point & p)
811{
812 Real value = FE<2,LAGRANGE>::shape_second_deriv( type, order, i/2, j, p );
813
814 switch( i%2 )
815 {
816 case 0:
817 return libMesh::RealGradient( value );
818
819 case 1:
820 return libMesh::RealGradient( Real(0), value );
821
822 default:
823 libmesh_error_msg("i%2 must be either 0 or 1!");
824 }
825
826 //dummy
827 return libMesh::RealGradient();
828}
829
830#endif
831
832template <> RealGradient FE<2,L2_LAGRANGE_VEC>::shape(const ElemType type, const Order order,
833 const unsigned int i, const Point & p)
834{
835 return FE<2,LAGRANGE_VEC>::shape(type, order, i, p);
836}
837
838template <> RealGradient FE<2,L2_LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
839 const unsigned int i, const unsigned int j,
840 const Point & p)
841{
842 return FE<2,LAGRANGE_VEC>::shape_deriv(type, order, i, j, p);
843}
844
845#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
846template <> RealGradient FE<2,L2_LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order,
847 const unsigned int i, const unsigned int j,
848 const Point & p)
849{
850 return FE<2,LAGRANGE_VEC>::shape_second_deriv(type, order, i, j, p);
851}
852
853#endif
854
855// 3-D
856template <> RealGradient FE<3,LAGRANGE_VEC>::shape(const ElemType type, const Order order,
857 const unsigned int i, const Point & p)
858{
859 Real value = FE<3,LAGRANGE>::shape( type, order, i/3, p );
860
861 switch( i%3 )
862 {
863 case 0:
864 return libMesh::RealGradient( value );
865
866 case 1:
867 return libMesh::RealGradient( Real(0), value );
868
869 case 2:
870 return libMesh::RealGradient( Real(0), Real(0), value );
871
872 default:
873 libmesh_error_msg("i%3 must be 0, 1, or 2!");
874 }
875
876 //dummy
877 return libMesh::RealGradient();
878}
879template <> RealGradient FE<3,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
880 const unsigned int i, const unsigned int j,
881 const Point & p)
882{
883 Real value = FE<3,LAGRANGE>::shape_deriv( type, order, i/3, j, p );
884
885 switch( i%3 )
886 {
887 case 0:
888 return libMesh::RealGradient( value );
889
890 case 1:
891 return libMesh::RealGradient( Real(0), value );
892
893 case 2:
894 return libMesh::RealGradient( Real(0), Real(0), value );
895
896 default:
897 libmesh_error_msg("i%3 must be 0, 1, or 2!");
898 }
899
900 //dummy
901 return libMesh::RealGradient();
902}
903
904#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
905
906template <> RealGradient FE<3,LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order,
907 const unsigned int i, const unsigned int j,
908 const Point & p)
909{
910 Real value = FE<3,LAGRANGE>::shape_second_deriv( type, order, i/3, j, p );
911
912 switch( i%3 )
913 {
914 case 0:
915 return libMesh::RealGradient( value );
916
917 case 1:
918 return libMesh::RealGradient( Real(0), value );
919
920 case 2:
921 return libMesh::RealGradient( Real(0), Real(0), value );
922
923 default:
924 libmesh_error_msg("i%3 must be 0, 1, or 2!");
925 }
926
927 //dummy
928 return libMesh::RealGradient();
929}
930
931#endif
932
933template <> RealGradient FE<3,L2_LAGRANGE_VEC>::shape(const ElemType type, const Order order,
934 const unsigned int i, const Point & p)
935{
936 return FE<3,LAGRANGE_VEC>::shape(type, order, i, p);
937}
938
939template <> RealGradient FE<3,L2_LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
940 const unsigned int i, const unsigned int j,
941 const Point & p)
942{
943 return FE<3,LAGRANGE_VEC>::shape_deriv(type, order, i, j, p);
944}
945
946#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
947template <> RealGradient FE<3,L2_LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order,
948 const unsigned int i, const unsigned int j,
949 const Point & p)
950{
951 return FE<3,LAGRANGE_VEC>::shape_second_deriv(type, order, i, j, p);
952}
953
954#endif
955
956
957// 0-D
958template <> RealGradient FE<0,LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
959 const unsigned int i, const Point & p,
960 const bool add_p_level)
961{
962 Real value = FE<0,LAGRANGE>::shape( elem->type(), order + add_p_level*elem->p_level(), i, p);
963 return libMesh::RealGradient( value );
964}
965template <> RealGradient FE<0,LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
966 const unsigned int i, const unsigned int j,
967 const Point & p,
968 const bool add_p_level)
969{
970 Real value = FE<0,LAGRANGE>::shape_deriv( elem->type(), order + add_p_level*elem->p_level(), i, j, p);
971 return libMesh::RealGradient( value );
972}
973
974#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
975
976template <> RealGradient FE<0,LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
977 const unsigned int i, const unsigned int j,
978 const Point & p,
979 const bool add_p_level)
980{
981 Real value = FE<0,LAGRANGE>::shape_second_deriv( elem->type(), order + add_p_level*elem->p_level(), i, j, p);
982 return libMesh::RealGradient( value );
983}
984
985#endif
986
987template <> RealGradient FE<0,L2_LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
988 const unsigned int i, const Point & p,
989 const bool add_p_level)
990{
991 return FE<0,LAGRANGE_VEC>::shape(elem, order, i, p, add_p_level);
992}
993template <> RealGradient FE<0,L2_LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
994 const unsigned int i, const unsigned int j,
995 const Point & p,
996 const bool add_p_level)
997{
998 return FE<0,LAGRANGE_VEC>::shape_deriv(elem, order, i, j, p, add_p_level);
999}
1000
1001#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1002
1003template <> RealGradient FE<0,L2_LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
1004 const unsigned int i, const unsigned int j,
1005 const Point & p,
1006 const bool add_p_level)
1007{
1008 return FE<0,LAGRANGE_VEC>::shape_second_deriv(elem, order, i, j, p, add_p_level);
1009}
1010
1011#endif
1012
1013// 1-D
1014template <> RealGradient FE<1,LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
1015 const unsigned int i, const Point & p,
1016 const bool add_p_level)
1017{
1018 Real value = FE<1,LAGRANGE>::shape( elem->type(), order + add_p_level*elem->p_level(), i, p);
1019 return libMesh::RealGradient( value );
1020}
1021template <> RealGradient FE<1,LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
1022 const unsigned int i, const unsigned int j,
1023 const Point & p,
1024 const bool add_p_level)
1025{
1026 Real value = FE<1,LAGRANGE>::shape_deriv( elem->type(), order + add_p_level*elem->p_level(), i, j, p);
1027 return libMesh::RealGradient( value );
1028}
1029
1030#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1031template <> RealGradient FE<1,LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
1032 const unsigned int i, const unsigned int j,
1033 const Point & p,
1034 const bool add_p_level)
1035{
1036 Real value = FE<1,LAGRANGE>::shape_second_deriv( elem->type(), order + add_p_level*elem->p_level(), i, j, p);
1037 return libMesh::RealGradient( value );
1038}
1039
1040#endif
1041
1042template <> RealGradient FE<1,L2_LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
1043 const unsigned int i, const Point & p,
1044 const bool add_p_level)
1045{
1046 return FE<1,LAGRANGE_VEC>::shape(elem, order, i, p, add_p_level);
1047}
1048template <> RealGradient FE<1,L2_LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
1049 const unsigned int i, const unsigned int j,
1050 const Point & p,
1051 const bool add_p_level)
1052{
1053 return FE<1,LAGRANGE_VEC>::shape_deriv(elem, order, i, j, p, add_p_level);
1054}
1055
1056#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1057
1058template <> RealGradient FE<1,L2_LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
1059 const unsigned int i, const unsigned int j,
1060 const Point & p,
1061 const bool add_p_level)
1062{
1063 return FE<1,LAGRANGE_VEC>::shape_second_deriv(elem, order, i, j, p, add_p_level);
1064}
1065
1066#endif
1067
1068// 2-D
1069template <> RealGradient FE<2,LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
1070 const unsigned int i, const Point & p,
1071 const bool add_p_level)
1072{
1073 Real value = FE<2,LAGRANGE>::shape( elem->type(), order + add_p_level*elem->p_level(), i/2, p );
1074
1075 switch( i%2 )
1076 {
1077 case 0:
1078 return libMesh::RealGradient( value );
1079
1080 case 1:
1081 return libMesh::RealGradient( Real(0), value );
1082
1083 default:
1084 libmesh_error_msg("i%2 must be either 0 or 1!");
1085 }
1086
1087 //dummy
1088 return libMesh::RealGradient();
1089}
1090template <> RealGradient FE<2,LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
1091 const unsigned int i, const unsigned int j,
1092 const Point & p,
1093 const bool add_p_level)
1094{
1095 Real value = FE<2,LAGRANGE>::shape_deriv( elem->type(), order + add_p_level*elem->p_level(), i/2, j, p );
1096
1097 switch( i%2 )
1098 {
1099 case 0:
1100 return libMesh::RealGradient( value );
1101
1102 case 1:
1103 return libMesh::RealGradient( Real(0), value );
1104
1105 default:
1106 libmesh_error_msg("i%2 must be either 0 or 1!");
1107 }
1108
1109 //dummy
1110 return libMesh::RealGradient();
1111}
1112
1113#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1114template <> RealGradient FE<2,LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
1115 const unsigned int i, const unsigned int j,
1116 const Point & p,
1117 const bool add_p_level)
1118{
1119 Real value = FE<2,LAGRANGE>::shape_second_deriv( elem->type(), order + add_p_level*elem->p_level(), i/2, j, p );
1120
1121 switch( i%2 )
1122 {
1123 case 0:
1124 return libMesh::RealGradient( value );
1125
1126 case 1:
1127 return libMesh::RealGradient( Real(0), value );
1128
1129 default:
1130 libmesh_error_msg("i%2 must be either 0 or 1!");
1131 }
1132
1133 //dummy
1134 return libMesh::RealGradient();
1135}
1136
1137#endif
1138
1139template <> RealGradient FE<2,L2_LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
1140 const unsigned int i, const Point & p,
1141 const bool add_p_level)
1142{
1143 return FE<2,LAGRANGE_VEC>::shape(elem, order, i, p, add_p_level);
1144}
1145template <> RealGradient FE<2,L2_LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
1146 const unsigned int i, const unsigned int j,
1147 const Point & p,
1148 const bool add_p_level)
1149{
1150 return FE<2,LAGRANGE_VEC>::shape_deriv(elem, order, i, j, p, add_p_level);
1151}
1152
1153#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1154
1155template <> RealGradient FE<2,L2_LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
1156 const unsigned int i, const unsigned int j,
1157 const Point & p,
1158 const bool add_p_level)
1159{
1160 return FE<2,LAGRANGE_VEC>::shape_second_deriv(elem, order, i, j, p, add_p_level);
1161}
1162
1163#endif
1164
1165// 3-D
1166template <> RealGradient FE<3,LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
1167 const unsigned int i, const Point & p,
1168 const bool add_p_level)
1169{
1170 Real value = FE<3,LAGRANGE>::shape( elem->type(), order + add_p_level*elem->p_level(), i/3, p );
1171
1172 switch( i%3 )
1173 {
1174 case 0:
1175 return libMesh::RealGradient( value );
1176
1177 case 1:
1178 return libMesh::RealGradient( Real(0), value );
1179
1180 case 2:
1181 return libMesh::RealGradient( Real(0), Real(0), value );
1182
1183 default:
1184 libmesh_error_msg("i%3 must be 0, 1, or 2!");
1185 }
1186
1187 //dummy
1188 return libMesh::RealGradient();
1189}
1190template <> RealGradient FE<3,LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
1191 const unsigned int i, const unsigned int j,
1192 const Point & p,
1193 const bool add_p_level)
1194{
1195 Real value = FE<3,LAGRANGE>::shape_deriv( elem->type(), order + add_p_level*elem->p_level(), i/3, j, p );
1196
1197 switch( i%3 )
1198 {
1199 case 0:
1200 return libMesh::RealGradient( value );
1201
1202 case 1:
1203 return libMesh::RealGradient( Real(0), value );
1204
1205 case 2:
1206 return libMesh::RealGradient( Real(0), Real(0), value );
1207
1208 default:
1209 libmesh_error_msg("i%3 must be 0, 1, or 2!");
1210 }
1211
1212 //dummy
1213 return libMesh::RealGradient();
1214}
1215
1216#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1217
1218template <> RealGradient FE<3,LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
1219 const unsigned int i, const unsigned int j,
1220 const Point & p,
1221 const bool add_p_level)
1222{
1223 Real value = FE<3,LAGRANGE>::shape_second_deriv( elem->type(), order + add_p_level*elem->p_level(), i/3, j, p );
1224
1225 switch( i%3 )
1226 {
1227 case 0:
1228 return libMesh::RealGradient( value );
1229
1230 case 1:
1231 return libMesh::RealGradient( Real(0), value );
1232
1233 case 2:
1234 return libMesh::RealGradient( Real(0), Real(0), value );
1235
1236 default:
1237 libmesh_error_msg("i%3 must be 0, 1, or 2!");
1238 }
1239
1240 //dummy
1241 return libMesh::RealGradient();
1242}
1243
1244#endif
1245
1246template <> RealGradient FE<3,L2_LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
1247 const unsigned int i, const Point & p,
1248 const bool add_p_level)
1249{
1250 return FE<3,LAGRANGE_VEC>::shape(elem, order, i, p, add_p_level);
1251}
1252template <> RealGradient FE<3,L2_LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
1253 const unsigned int i, const unsigned int j,
1254 const Point & p,
1255 const bool add_p_level)
1256{
1257 return FE<3,LAGRANGE_VEC>::shape_deriv(elem, order, i, j, p, add_p_level);
1258}
1259
1260#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1261
1262template <> RealGradient FE<3,L2_LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
1263 const unsigned int i, const unsigned int j,
1264 const Point & p,
1265 const bool add_p_level)
1266{
1267 return FE<3,LAGRANGE_VEC>::shape_second_deriv(elem, order, i, j, p, add_p_level);
1268}
1269
1270#endif
1271
1272// Do full-specialization for every dimension, instead
1273// of explicit instantiation at the end of this function.
1274template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<0,LAGRANGE>::n_dofs(t,o); }
1275template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<1,LAGRANGE>::n_dofs(t,o); }
1276template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 2*FE<2,LAGRANGE>::n_dofs(t,o); }
1277template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 3*FE<3,LAGRANGE>::n_dofs(t,o); }
1278
1279template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return FE<0,LAGRANGE>::n_dofs(e,o); }
1280template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return FE<1,LAGRANGE>::n_dofs(e,o); }
1281template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return 2*FE<2,LAGRANGE>::n_dofs(e,o); }
1282template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return 3*FE<3,LAGRANGE>::n_dofs(e,o); }
1283
1284template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<0,L2_LAGRANGE>::n_dofs(t,o); }
1285template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<1,L2_LAGRANGE>::n_dofs(t,o); }
1286template <> unsigned int FE<2,L2_LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 2*FE<2,L2_LAGRANGE>::n_dofs(t,o); }
1287template <> unsigned int FE<3,L2_LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 3*FE<3,L2_LAGRANGE>::n_dofs(t,o); }
1288
1289template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return FE<0,L2_LAGRANGE>::n_dofs(e,o); }
1290template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return FE<1,L2_LAGRANGE>::n_dofs(e,o); }
1291template <> unsigned int FE<2,L2_LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return 2*FE<2,L2_LAGRANGE>::n_dofs(e,o); }
1292template <> unsigned int FE<3,L2_LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return 3*FE<3,L2_LAGRANGE>::n_dofs(e,o); }
1293
1294
1295// Do full-specialization for every dimension, instead
1296// of explicit instantiation at the end of this function.
1297template <> 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); }
1298template <> 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); }
1299template <> 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); }
1300template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return 3*FE<3,LAGRANGE>::n_dofs_at_node(t,o,n); }
1301
1302template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<0,LAGRANGE>::n_dofs_at_node(e,o,n); }
1303template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<1,LAGRANGE>::n_dofs_at_node(e,o,n); }
1304template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return 2*FE<2,LAGRANGE>::n_dofs_at_node(e,o,n); }
1305template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return 3*FE<3,LAGRANGE>::n_dofs_at_node(e,o,n); }
1306
1307template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
1308template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
1309template <> unsigned int FE<2,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
1310template <> unsigned int FE<3,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
1311
1312template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
1313template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
1314template <> unsigned int FE<2,L2_LAGRANGE_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
1315template <> unsigned int FE<3,L2_LAGRANGE_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
1316
1317
1318// Lagrange elements have no dofs per element
1319// (just at the nodes)
1320template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1321template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1322template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1323template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1324
1325template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1326template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1327template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1328template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1329
1330// L2 Lagrange elements have all their dofs on the element
1331template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<0,L2_LAGRANGE_VEC>::n_dofs(t, o); }
1332template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<1,L2_LAGRANGE_VEC>::n_dofs(t, o); }
1333template <> unsigned int FE<2,L2_LAGRANGE_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<2,L2_LAGRANGE_VEC>::n_dofs(t, o); }
1334template <> unsigned int FE<3,L2_LAGRANGE_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<3,L2_LAGRANGE_VEC>::n_dofs(t, o); }
1335
1336template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<0,L2_LAGRANGE_VEC>::n_dofs(&e, o); }
1337template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<1,L2_LAGRANGE_VEC>::n_dofs(&e, o); }
1338template <> unsigned int FE<2,L2_LAGRANGE_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<2,L2_LAGRANGE_VEC>::n_dofs(&e, o); }
1339template <> unsigned int FE<3,L2_LAGRANGE_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<3,L2_LAGRANGE_VEC>::n_dofs(&e, o); }
1340
1341// Lagrange FEMs are always C^0 continuous
1342template <> FEContinuity FE<0,LAGRANGE_VEC>::get_continuity() const { return C_ZERO; }
1343template <> FEContinuity FE<1,LAGRANGE_VEC>::get_continuity() const { return C_ZERO; }
1344template <> FEContinuity FE<2,LAGRANGE_VEC>::get_continuity() const { return C_ZERO; }
1345template <> FEContinuity FE<3,LAGRANGE_VEC>::get_continuity() const { return C_ZERO; }
1346
1347// L2 Lagrange FEMs are always discontinuous
1348template <> FEContinuity FE<0,L2_LAGRANGE_VEC>::get_continuity() const { return DISCONTINUOUS; }
1349template <> FEContinuity FE<1,L2_LAGRANGE_VEC>::get_continuity() const { return DISCONTINUOUS; }
1350template <> FEContinuity FE<2,L2_LAGRANGE_VEC>::get_continuity() const { return DISCONTINUOUS; }
1351template <> FEContinuity FE<3,L2_LAGRANGE_VEC>::get_continuity() const { return DISCONTINUOUS; }
1352
1353// Lagrange FEMs are not hierarchic
1354template <> bool FE<0,LAGRANGE_VEC>::is_hierarchic() const { return false; }
1355template <> bool FE<1,LAGRANGE_VEC>::is_hierarchic() const { return false; }
1356template <> bool FE<2,LAGRANGE_VEC>::is_hierarchic() const { return false; }
1357template <> bool FE<3,LAGRANGE_VEC>::is_hierarchic() const { return false; }
1358template <> bool FE<0,L2_LAGRANGE_VEC>::is_hierarchic() const { return false; }
1359template <> bool FE<1,L2_LAGRANGE_VEC>::is_hierarchic() const { return false; }
1360template <> bool FE<2,L2_LAGRANGE_VEC>::is_hierarchic() const { return false; }
1361template <> bool FE<3,L2_LAGRANGE_VEC>::is_hierarchic() const { return false; }
1362
1363// Lagrange FEM shapes do not need reinit (is this always true?)
1364template <> bool FE<0,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1365template <> bool FE<1,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1366template <> bool FE<2,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1367template <> bool FE<3,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1368template <> bool FE<0,L2_LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1369template <> bool FE<1,L2_LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1370template <> bool FE<2,L2_LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1371template <> bool FE<3,L2_LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1372
1373// Methods for computing Lagrange constraints. Note: we pass the
1374// dimension as the last argument to the anonymous helper function.
1375// Also note: we only need instantiations of this function for
1376// Dim==2 and 3.
1377#ifdef LIBMESH_ENABLE_AMR
1378template <>
1379void FE<2,LAGRANGE_VEC>::compute_constraints (DofConstraints & constraints,
1380 DofMap & dof_map,
1381 const unsigned int variable_number,
1382 const Elem * elem)
1383{ //libmesh_not_implemented();
1384 FEVectorBase::compute_proj_constraints(constraints, dof_map, variable_number, elem);
1385}
1386
1387template <>
1388void FE<3,LAGRANGE_VEC>::compute_constraints (DofConstraints & constraints,
1389 DofMap & dof_map,
1390 const unsigned int variable_number,
1391 const Elem * elem)
1392{ //libmesh_not_implemented();
1393 FEVectorBase::compute_proj_constraints(constraints, dof_map, variable_number, elem);
1394}
1395
1396template <>
1397void FE<2,L2_LAGRANGE_VEC>::compute_constraints (DofConstraints & constraints,
1398 DofMap & dof_map,
1399 const unsigned int variable_number,
1400 const Elem * elem)
1401{ //libmesh_not_implemented();
1402 FEVectorBase::compute_proj_constraints(constraints, dof_map, variable_number, elem);
1403}
1404
1405template <>
1406void FE<3,L2_LAGRANGE_VEC>::compute_constraints (DofConstraints & constraints,
1407 DofMap & dof_map,
1408 const unsigned int variable_number,
1409 const Elem * elem)
1410{ //libmesh_not_implemented();
1411 FEVectorBase::compute_proj_constraints(constraints, dof_map, variable_number, elem);
1412}
1413#endif // LIBMESH_ENABLE_AMR
1414
1415} // namespace libMesh
dim
unsigned int dim
Definition adaptivity_ex3.C:114
libMesh::DofConstraints
The constraint matrix storage format.
Definition dof_map.h:112
libMesh::DofMap
This class handles the numbering of degrees of freedom on a mesh.
Definition dof_map.h:181
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::Elem::n_nodes
virtual unsigned int n_nodes() const =0
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Elem::p_level
unsigned int p_level() const
Definition elem.h:3125
libMesh::FEGenericBase::compute_proj_constraints
static void compute_proj_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
Computes the constraint matrix contributions (for non-conforming adapted meshes) corresponding to var...
Definition fe_base.C:1534
libMesh::FE
A specific instantiation of the FEBase class.
Definition fe.h:128
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::FE::is_hierarchic
virtual bool is_hierarchic() const override
libMesh::FE::nodal_soln
static void nodal_soln(const Elem *elem, const Order o, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln, bool add_p_level=true, const unsigned vdim=1)
Build the nodal soln from the element soln.
libMesh::FE::n_dofs_per_elem
static unsigned int n_dofs_per_elem(const ElemType t, const Order o)
libMesh::FE::get_continuity
virtual FEContinuity get_continuity() const override
libMesh::FE::compute_constraints
static void compute_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
Computes the constraint matrix contributions (for non-conforming adapted meshes) corresponding to var...
libMesh::FE::n_dofs
static unsigned int n_dofs(const ElemType t, const Order o)
libMesh::FE::n_dofs_at_node
static unsigned int n_dofs_at_node(const ElemType t, const Order o, const unsigned int n)
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::FE::shape_deriv
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::FE::shapes_need_reinit
virtual bool shapes_need_reinit() const override
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::LIBMESH_DEFAULT_VECTORIZED_FE
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
Definition fe_bernstein_shape_0D.C:30
libMesh::LIBMESH_FE_NODAL_SOLN_DIM
LIBMESH_FE_NODAL_SOLN_DIM(LIBMESH_FE_NODAL_SOLN_DIM(HIERARCHIC_VEC,(FE< 0, HIERARCHIC >::nodal_soln), 0)
Definition fe_hierarchic_vec.C:96
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::RealGradient
RealVectorValue RealGradient
Definition hp_coarsentest.h:49
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
libMesh::LIBMESH_FE_SIDE_NODAL_SOLN
LIBMESH_FE_SIDE_NODAL_SOLN(HIERARCHIC_VEC)
Definition fe_hierarchic_vec.C:117
n_nodes
const dof_id_type n_nodes
Definition tecplot_io.C:67
value
static const bool value
Definition xdr_io.C:55
Generated on Fri Aug 21 2026 13:32:28 for libMesh by doxygen 1.9.8