libMesh
fe_lagrange_vec.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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 
30 namespace libMesh
31 {
32 
33 
42 
43 
44 // ------------------------------------------------------------
45 // Lagrange-specific implementations
46 
47 
48 // Anonymous namespace for local helper functions
49 namespace {
50 void 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...
637 LIBMESH_FE_NODAL_SOLN_DIM(LAGRANGE_VEC, (FE<1, LAGRANGE>::nodal_soln), 1)
638 
639 template <>
640 void 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 
648 template <>
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 
658 
663 
665 
666 
667 // Specialize for shape function routines by leveraging scalar LAGRANGE elements
668 
669 // 0-D
670 template <> 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 }
676 template <> 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 
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 
695 template <> 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 }
700 template <> 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 
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
718 template <> 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 }
724 template <> 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
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 
742 template <> 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 }
747 template <> 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 
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
765 template <> 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 }
785 template <> 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
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 
832 template <> 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 
838 template <> 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
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
856 template <> 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 }
879 template <> 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 
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 
933 template <> 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 
939 template <> 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
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
958 template <> 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 }
965 template <> 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 
976 template <> 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 
987 template <> 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 }
993 template <> 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 
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
1014 template <> 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 }
1021 template <> 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
1031 template <> 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 
1042 template <> 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 }
1048 template <> 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 
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
1069 template <> 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 }
1090 template <> 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
1114 template <> 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 
1139 template <> 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 }
1145 template <> 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 
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
1166 template <> 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 }
1190 template <> 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 
1218 template <> 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 
1246 template <> 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 }
1252 template <> 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 
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.
1274 template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<0,LAGRANGE>::n_dofs(t,o); }
1275 template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<1,LAGRANGE>::n_dofs(t,o); }
1276 template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 2*FE<2,LAGRANGE>::n_dofs(t,o); }
1277 template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 3*FE<3,LAGRANGE>::n_dofs(t,o); }
1278 
1279 template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return FE<0,LAGRANGE>::n_dofs(e,o); }
1280 template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return FE<1,LAGRANGE>::n_dofs(e,o); }
1281 template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return 2*FE<2,LAGRANGE>::n_dofs(e,o); }
1282 template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return 3*FE<3,LAGRANGE>::n_dofs(e,o); }
1283 
1284 template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<0,L2_LAGRANGE>::n_dofs(t,o); }
1285 template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<1,L2_LAGRANGE>::n_dofs(t,o); }
1286 template <> 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); }
1287 template <> 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 
1289 template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return FE<0,L2_LAGRANGE>::n_dofs(e,o); }
1290 template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return FE<1,L2_LAGRANGE>::n_dofs(e,o); }
1291 template <> 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); }
1292 template <> 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.
1297 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); }
1298 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); }
1299 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); }
1300 template <> 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 
1302 template <> 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); }
1303 template <> 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); }
1304 template <> 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); }
1305 template <> 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 
1307 template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
1308 template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
1309 template <> unsigned int FE<2,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
1310 template <> unsigned int FE<3,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
1311 
1312 template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
1313 template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
1314 template <> unsigned int FE<2,L2_LAGRANGE_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
1315 template <> 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)
1320 template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1321 template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1322 template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1323 template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1324 
1325 template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1326 template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1327 template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1328 template <> 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
1331 template <> 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); }
1332 template <> 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); }
1333 template <> 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); }
1334 template <> 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 
1336 template <> 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); }
1337 template <> 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); }
1338 template <> 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); }
1339 template <> 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
1346 
1347 // L2 Lagrange FEMs are always discontinuous
1352 
1353 // Lagrange FEMs are not hierarchic
1354 template <> bool FE<0,LAGRANGE_VEC>::is_hierarchic() const { return false; }
1355 template <> bool FE<1,LAGRANGE_VEC>::is_hierarchic() const { return false; }
1356 template <> bool FE<2,LAGRANGE_VEC>::is_hierarchic() const { return false; }
1357 template <> bool FE<3,LAGRANGE_VEC>::is_hierarchic() const { return false; }
1358 template <> bool FE<0,L2_LAGRANGE_VEC>::is_hierarchic() const { return false; }
1359 template <> bool FE<1,L2_LAGRANGE_VEC>::is_hierarchic() const { return false; }
1360 template <> bool FE<2,L2_LAGRANGE_VEC>::is_hierarchic() const { return false; }
1361 template <> bool FE<3,L2_LAGRANGE_VEC>::is_hierarchic() const { return false; }
1362 
1363 // Lagrange FEM shapes do not need reinit (is this always true?)
1364 template <> bool FE<0,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1365 template <> bool FE<1,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1366 template <> bool FE<2,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1367 template <> bool FE<3,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1368 template <> bool FE<0,L2_LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1369 template <> bool FE<1,L2_LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1370 template <> bool FE<2,L2_LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1371 template <> 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
1378 template <>
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 
1387 template <>
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 
1396 template <>
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 
1405 template <>
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
static unsigned int n_dofs(const ElemType t, const Order o)
ElemType
Defines an enum for geometric element types.
RealVectorValue RealGradient
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
LIBMESH_FE_NODAL_SOLN_DIM(LIBMESH_FE_NODAL_SOLN_DIM(HIERARCHIC_VEC,(FE< 0, HIERARCHIC >::nodal_soln), 0)
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
unsigned int dim
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
static unsigned int n_dofs_at_node(const ElemType t, const Order o, const unsigned int n)
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
unsigned int p_level() const
Definition: elem.h:3108
virtual bool shapes_need_reinit() const override
The libMesh namespace provides an interface to certain functionality in the library.
virtual bool is_hierarchic() const override
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
A specific instantiation of the FEBase class.
Definition: fe.h:127
LIBMESH_FE_SIDE_NODAL_SOLN(HIERARCHIC_VEC)
const dof_id_type n_nodes
Definition: tecplot_io.C:67
virtual unsigned int n_nodes() const =0
libmesh_assert(ctx)
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:1570
static unsigned int n_dofs_per_elem(const ElemType t, const Order o)
virtual FEContinuity get_continuity() const override
std::string enum_to_string(const T e)
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...
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.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
static const bool value
Definition: xdr_io.C:54
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
The constraint matrix storage format.
Definition: dof_map.h:108
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)