libMesh
fe_lagrange_vec.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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 = static_cast<Order>(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...
636 template <>
638  const Order order,
639  const std::vector<Number> & elem_soln,
640  std::vector<Number> & nodal_soln,
641  const bool add_p_level)
642 { FE<0,LAGRANGE>::nodal_soln(elem, order, elem_soln, nodal_soln, add_p_level); }
643 
644 template <>
646  const Order order,
647  const std::vector<Number> & elem_soln,
648  std::vector<Number> & nodal_soln,
649  const bool add_p_level)
650 { FE<1,LAGRANGE>::nodal_soln(elem, order, elem_soln, nodal_soln, add_p_level); }
651 
652 template <>
654  const Order order,
655  const std::vector<Number> & elem_soln,
656  std::vector<Number> & nodal_soln,
657  const bool add_p_level)
658 { lagrange_vec_nodal_soln(elem, order, elem_soln, 2 /*dimension*/, nodal_soln, add_p_level); }
659 
660 template <>
662  const Order order,
663  const std::vector<Number> & elem_soln,
664  std::vector<Number> & nodal_soln,
665  const bool add_p_level)
666 { lagrange_vec_nodal_soln(elem, order, elem_soln, 3 /*dimension*/, nodal_soln, add_p_level); }
667 
669 
670 template <>
671 void FE<0,L2_LAGRANGE_VEC>::nodal_soln(const Elem * elem,
672  const Order order,
673  const std::vector<Number> & elem_soln,
674  std::vector<Number> & nodal_soln,
675  const bool add_p_level)
676 { FE<0,LAGRANGE_VEC>::nodal_soln(elem, order, elem_soln, nodal_soln, add_p_level); }
677 
678 template <>
680  const Order order,
681  const std::vector<Number> & elem_soln,
682  std::vector<Number> & nodal_soln,
683  const bool add_p_level)
684 { FE<1,LAGRANGE_VEC>::nodal_soln(elem, order, elem_soln, nodal_soln, add_p_level); }
685 
686 template <>
688  const Order order,
689  const std::vector<Number> & elem_soln,
690  std::vector<Number> & nodal_soln,
691  const bool add_p_level)
692 { FE<2,LAGRANGE_VEC>::nodal_soln(elem, order, elem_soln, nodal_soln, add_p_level); }
693 
694 template <>
696  const Order order,
697  const std::vector<Number> & elem_soln,
698  std::vector<Number> & nodal_soln,
699  const bool add_p_level)
700 { FE<3,LAGRANGE_VEC>::nodal_soln(elem, order, elem_soln, nodal_soln, add_p_level); }
701 
703 
704 
705 // Specialize for shape function routines by leveraging scalar LAGRANGE elements
706 
707 // 0-D
708 template <> RealGradient FE<0,LAGRANGE_VEC>::shape(const ElemType type, const Order order,
709  const unsigned int i, const Point & p)
710 {
711  Real value = FE<0,LAGRANGE>::shape( type, order, i, p );
712  return libMesh::RealGradient( value );
713 }
714 template <> RealGradient FE<0,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
715  const unsigned int i, const unsigned int j,
716  const Point & p)
717 {
718  Real value = FE<0,LAGRANGE>::shape_deriv( type, order, i, j, p );
719  return libMesh::RealGradient( value );
720 }
721 
722 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
723 
725  const unsigned int i, const unsigned int j,
726  const Point & p)
727 {
728  Real value = FE<0,LAGRANGE>::shape_second_deriv( type, order, i, j, p );
729  return libMesh::RealGradient( value );
730 }
731 #endif
732 
733 template <> RealGradient FE<0,L2_LAGRANGE_VEC>::shape(const ElemType type, const Order order,
734  const unsigned int i, const Point & p)
735 {
736  return FE<0,LAGRANGE_VEC>::shape(type, order, i, p);
737 }
738 template <> RealGradient FE<0,L2_LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
739  const unsigned int i, const unsigned int j,
740  const Point & p)
741 {
742  return FE<0,LAGRANGE_VEC>::shape_deriv(type, order, i, j, p);
743 }
744 
745 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
746 
748  const unsigned int i, const unsigned int j,
749  const Point & p)
750 {
751  return FE<0,LAGRANGE_VEC>::shape_second_deriv(type, order, i, j, p);
752 }
753 #endif
754 
755 // 1-D
756 template <> RealGradient FE<1,LAGRANGE_VEC>::shape(const ElemType type, const Order order,
757  const unsigned int i, const Point & p)
758 {
759  Real value = FE<1,LAGRANGE>::shape( type, order, i, p );
760  return libMesh::RealGradient( value );
761 }
762 template <> RealGradient FE<1,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
763  const unsigned int i, const unsigned int j,
764  const Point & p)
765 {
766  Real value = FE<1,LAGRANGE>::shape_deriv( type, order, i, j, p );
767  return libMesh::RealGradient( value );
768 }
769 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
771  const unsigned int i, const unsigned int j,
772  const Point & p)
773 {
774  Real value = FE<1,LAGRANGE>::shape_second_deriv( type, order, i, j, p );
775  return libMesh::RealGradient( value );
776 }
777 
778 #endif
779 
780 template <> RealGradient FE<1,L2_LAGRANGE_VEC>::shape(const ElemType type, const Order order,
781  const unsigned int i, const Point & p)
782 {
783  return FE<1,LAGRANGE_VEC>::shape(type, order, i, p);
784 }
785 template <> RealGradient FE<1,L2_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  return FE<1,LAGRANGE_VEC>::shape_deriv(type, order, i, j, p);
790 }
791 
792 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
793 
795  const unsigned int i, const unsigned int j,
796  const Point & p)
797 {
798  return FE<1,LAGRANGE_VEC>::shape_second_deriv(type, order, i, j, p);
799 }
800 #endif
801 
802 // 2-D
803 template <> RealGradient FE<2,LAGRANGE_VEC>::shape(const ElemType type, const Order order,
804  const unsigned int i, const Point & p)
805 {
806  Real value = FE<2,LAGRANGE>::shape( type, order, i/2, p );
807 
808  switch( i%2 )
809  {
810  case 0:
811  return libMesh::RealGradient( value );
812 
813  case 1:
814  return libMesh::RealGradient( Real(0), value );
815 
816  default:
817  libmesh_error_msg("i%2 must be either 0 or 1!");
818  }
819 
820  //dummy
821  return libMesh::RealGradient();
822 }
823 template <> RealGradient FE<2,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
824  const unsigned int i, const unsigned int j,
825  const Point & p)
826 {
827  Real value = FE<2,LAGRANGE>::shape_deriv( type, order, i/2, j, p );
828 
829  switch( i%2 )
830  {
831  case 0:
832  return libMesh::RealGradient( value );
833 
834  case 1:
835  return libMesh::RealGradient( Real(0), value );
836 
837  default:
838  libmesh_error_msg("i%2 must be either 0 or 1!");
839  }
840 
841  //dummy
842  return libMesh::RealGradient();
843 }
844 
845 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
847  const unsigned int i, const unsigned int j,
848  const Point & p)
849 {
850  Real value = FE<2,LAGRANGE>::shape_second_deriv( type, order, i/2, j, p );
851 
852  switch( i%2 )
853  {
854  case 0:
855  return libMesh::RealGradient( value );
856 
857  case 1:
858  return libMesh::RealGradient( Real(0), value );
859 
860  default:
861  libmesh_error_msg("i%2 must be either 0 or 1!");
862  }
863 
864  //dummy
865  return libMesh::RealGradient();
866 }
867 
868 #endif
869 
870 template <> RealGradient FE<2,L2_LAGRANGE_VEC>::shape(const ElemType type, const Order order,
871  const unsigned int i, const Point & p)
872 {
873  return FE<2,LAGRANGE_VEC>::shape(type, order, i, p);
874 }
875 
876 template <> RealGradient FE<2,L2_LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
877  const unsigned int i, const unsigned int j,
878  const Point & p)
879 {
880  return FE<2,LAGRANGE_VEC>::shape_deriv(type, order, i, j, p);
881 }
882 
883 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
885  const unsigned int i, const unsigned int j,
886  const Point & p)
887 {
888  return FE<2,LAGRANGE_VEC>::shape_second_deriv(type, order, i, j, p);
889 }
890 
891 #endif
892 
893 // 3-D
894 template <> RealGradient FE<3,LAGRANGE_VEC>::shape(const ElemType type, const Order order,
895  const unsigned int i, const Point & p)
896 {
897  Real value = FE<3,LAGRANGE>::shape( type, order, i/3, p );
898 
899  switch( i%3 )
900  {
901  case 0:
902  return libMesh::RealGradient( value );
903 
904  case 1:
905  return libMesh::RealGradient( Real(0), value );
906 
907  case 2:
908  return libMesh::RealGradient( Real(0), Real(0), value );
909 
910  default:
911  libmesh_error_msg("i%3 must be 0, 1, or 2!");
912  }
913 
914  //dummy
915  return libMesh::RealGradient();
916 }
917 template <> RealGradient FE<3,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
918  const unsigned int i, const unsigned int j,
919  const Point & p)
920 {
921  Real value = FE<3,LAGRANGE>::shape_deriv( type, order, i/3, j, p );
922 
923  switch( i%3 )
924  {
925  case 0:
926  return libMesh::RealGradient( value );
927 
928  case 1:
929  return libMesh::RealGradient( Real(0), value );
930 
931  case 2:
932  return libMesh::RealGradient( Real(0), Real(0), value );
933 
934  default:
935  libmesh_error_msg("i%3 must be 0, 1, or 2!");
936  }
937 
938  //dummy
939  return libMesh::RealGradient();
940 }
941 
942 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
943 
945  const unsigned int i, const unsigned int j,
946  const Point & p)
947 {
948  Real value = FE<3,LAGRANGE>::shape_second_deriv( type, order, i/3, j, p );
949 
950  switch( i%3 )
951  {
952  case 0:
953  return libMesh::RealGradient( value );
954 
955  case 1:
956  return libMesh::RealGradient( Real(0), value );
957 
958  case 2:
959  return libMesh::RealGradient( Real(0), Real(0), value );
960 
961  default:
962  libmesh_error_msg("i%3 must be 0, 1, or 2!");
963  }
964 
965  //dummy
966  return libMesh::RealGradient();
967 }
968 
969 #endif
970 
971 template <> RealGradient FE<3,L2_LAGRANGE_VEC>::shape(const ElemType type, const Order order,
972  const unsigned int i, const Point & p)
973 {
974  return FE<3,LAGRANGE_VEC>::shape(type, order, i, p);
975 }
976 
977 template <> RealGradient FE<3,L2_LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
978  const unsigned int i, const unsigned int j,
979  const Point & p)
980 {
981  return FE<3,LAGRANGE_VEC>::shape_deriv(type, order, i, j, p);
982 }
983 
984 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
986  const unsigned int i, const unsigned int j,
987  const Point & p)
988 {
989  return FE<3,LAGRANGE_VEC>::shape_second_deriv(type, order, i, j, p);
990 }
991 
992 #endif
993 
994 
995 // 0-D
996 template <> RealGradient FE<0,LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
997  const unsigned int i, const Point & p,
998  const bool add_p_level)
999 {
1000  Real value = FE<0,LAGRANGE>::shape( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, p);
1001  return libMesh::RealGradient( value );
1002 }
1003 template <> RealGradient FE<0,LAGRANGE_VEC>::shape_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  Real value = FE<0,LAGRANGE>::shape_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
1009  return libMesh::RealGradient( value );
1010 }
1011 
1012 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1013 
1014 template <> RealGradient FE<0,LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
1015  const unsigned int i, const unsigned int j,
1016  const Point & p,
1017  const bool add_p_level)
1018 {
1019  Real value = FE<0,LAGRANGE>::shape_second_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
1020  return libMesh::RealGradient( value );
1021 }
1022 
1023 #endif
1024 
1025 template <> RealGradient FE<0,L2_LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
1026  const unsigned int i, const Point & p,
1027  const bool add_p_level)
1028 {
1029  return FE<0,LAGRANGE_VEC>::shape(elem, order, i, p, add_p_level);
1030 }
1031 template <> RealGradient FE<0,L2_LAGRANGE_VEC>::shape_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  return FE<0,LAGRANGE_VEC>::shape_deriv(elem, order, i, j, p, add_p_level);
1037 }
1038 
1039 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1040 
1042  const unsigned int i, const unsigned int j,
1043  const Point & p,
1044  const bool add_p_level)
1045 {
1046  return FE<0,LAGRANGE_VEC>::shape_second_deriv(elem, order, i, j, p, add_p_level);
1047 }
1048 
1049 #endif
1050 
1051 // 1-D
1052 template <> RealGradient FE<1,LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
1053  const unsigned int i, const Point & p,
1054  const bool add_p_level)
1055 {
1056  Real value = FE<1,LAGRANGE>::shape( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, p);
1057  return libMesh::RealGradient( value );
1058 }
1059 template <> RealGradient FE<1,LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
1060  const unsigned int i, const unsigned int j,
1061  const Point & p,
1062  const bool add_p_level)
1063 {
1064  Real value = FE<1,LAGRANGE>::shape_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
1065  return libMesh::RealGradient( value );
1066 }
1067 
1068 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1069 template <> RealGradient FE<1,LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
1070  const unsigned int i, const unsigned int j,
1071  const Point & p,
1072  const bool add_p_level)
1073 {
1074  Real value = FE<1,LAGRANGE>::shape_second_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
1075  return libMesh::RealGradient( value );
1076 }
1077 
1078 #endif
1079 
1080 template <> RealGradient FE<1,L2_LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
1081  const unsigned int i, const Point & p,
1082  const bool add_p_level)
1083 {
1084  return FE<1,LAGRANGE_VEC>::shape(elem, order, i, p, add_p_level);
1085 }
1086 template <> RealGradient FE<1,L2_LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
1087  const unsigned int i, const unsigned int j,
1088  const Point & p,
1089  const bool add_p_level)
1090 {
1091  return FE<1,LAGRANGE_VEC>::shape_deriv(elem, order, i, j, p, add_p_level);
1092 }
1093 
1094 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1095 
1097  const unsigned int i, const unsigned int j,
1098  const Point & p,
1099  const bool add_p_level)
1100 {
1101  return FE<1,LAGRANGE_VEC>::shape_second_deriv(elem, order, i, j, p, add_p_level);
1102 }
1103 
1104 #endif
1105 
1106 // 2-D
1107 template <> RealGradient FE<2,LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
1108  const unsigned int i, const Point & p,
1109  const bool add_p_level)
1110 {
1111  Real value = FE<2,LAGRANGE>::shape( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i/2, p );
1112 
1113  switch( i%2 )
1114  {
1115  case 0:
1116  return libMesh::RealGradient( value );
1117 
1118  case 1:
1119  return libMesh::RealGradient( Real(0), value );
1120 
1121  default:
1122  libmesh_error_msg("i%2 must be either 0 or 1!");
1123  }
1124 
1125  //dummy
1126  return libMesh::RealGradient();
1127 }
1128 template <> RealGradient FE<2,LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
1129  const unsigned int i, const unsigned int j,
1130  const Point & p,
1131  const bool add_p_level)
1132 {
1133  Real value = FE<2,LAGRANGE>::shape_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i/2, j, p );
1134 
1135  switch( i%2 )
1136  {
1137  case 0:
1138  return libMesh::RealGradient( value );
1139 
1140  case 1:
1141  return libMesh::RealGradient( Real(0), value );
1142 
1143  default:
1144  libmesh_error_msg("i%2 must be either 0 or 1!");
1145  }
1146 
1147  //dummy
1148  return libMesh::RealGradient();
1149 }
1150 
1151 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1152 template <> RealGradient FE<2,LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
1153  const unsigned int i, const unsigned int j,
1154  const Point & p,
1155  const bool add_p_level)
1156 {
1157  Real value = FE<2,LAGRANGE>::shape_second_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i/2, j, p );
1158 
1159  switch( i%2 )
1160  {
1161  case 0:
1162  return libMesh::RealGradient( value );
1163 
1164  case 1:
1165  return libMesh::RealGradient( Real(0), value );
1166 
1167  default:
1168  libmesh_error_msg("i%2 must be either 0 or 1!");
1169  }
1170 
1171  //dummy
1172  return libMesh::RealGradient();
1173 }
1174 
1175 #endif
1176 
1177 template <> RealGradient FE<2,L2_LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
1178  const unsigned int i, const Point & p,
1179  const bool add_p_level)
1180 {
1181  return FE<2,LAGRANGE_VEC>::shape(elem, order, i, p, add_p_level);
1182 }
1183 template <> RealGradient FE<2,L2_LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
1184  const unsigned int i, const unsigned int j,
1185  const Point & p,
1186  const bool add_p_level)
1187 {
1188  return FE<2,LAGRANGE_VEC>::shape_deriv(elem, order, i, j, p, add_p_level);
1189 }
1190 
1191 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1192 
1194  const unsigned int i, const unsigned int j,
1195  const Point & p,
1196  const bool add_p_level)
1197 {
1198  return FE<2,LAGRANGE_VEC>::shape_second_deriv(elem, order, i, j, p, add_p_level);
1199 }
1200 
1201 #endif
1202 
1203 // 3-D
1204 template <> RealGradient FE<3,LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
1205  const unsigned int i, const Point & p,
1206  const bool add_p_level)
1207 {
1208  Real value = FE<3,LAGRANGE>::shape( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i/3, p );
1209 
1210  switch( i%3 )
1211  {
1212  case 0:
1213  return libMesh::RealGradient( value );
1214 
1215  case 1:
1216  return libMesh::RealGradient( Real(0), value );
1217 
1218  case 2:
1219  return libMesh::RealGradient( Real(0), Real(0), value );
1220 
1221  default:
1222  libmesh_error_msg("i%3 must be 0, 1, or 2!");
1223  }
1224 
1225  //dummy
1226  return libMesh::RealGradient();
1227 }
1228 template <> RealGradient FE<3,LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
1229  const unsigned int i, const unsigned int j,
1230  const Point & p,
1231  const bool add_p_level)
1232 {
1233  Real value = FE<3,LAGRANGE>::shape_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i/3, j, p );
1234 
1235  switch( i%3 )
1236  {
1237  case 0:
1238  return libMesh::RealGradient( value );
1239 
1240  case 1:
1241  return libMesh::RealGradient( Real(0), value );
1242 
1243  case 2:
1244  return libMesh::RealGradient( Real(0), Real(0), value );
1245 
1246  default:
1247  libmesh_error_msg("i%3 must be 0, 1, or 2!");
1248  }
1249 
1250  //dummy
1251  return libMesh::RealGradient();
1252 }
1253 
1254 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1255 
1256 template <> RealGradient FE<3,LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
1257  const unsigned int i, const unsigned int j,
1258  const Point & p,
1259  const bool add_p_level)
1260 {
1261  Real value = FE<3,LAGRANGE>::shape_second_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i/3, j, p );
1262 
1263  switch( i%3 )
1264  {
1265  case 0:
1266  return libMesh::RealGradient( value );
1267 
1268  case 1:
1269  return libMesh::RealGradient( Real(0), value );
1270 
1271  case 2:
1272  return libMesh::RealGradient( Real(0), Real(0), value );
1273 
1274  default:
1275  libmesh_error_msg("i%3 must be 0, 1, or 2!");
1276  }
1277 
1278  //dummy
1279  return libMesh::RealGradient();
1280 }
1281 
1282 #endif
1283 
1284 template <> RealGradient FE<3,L2_LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
1285  const unsigned int i, const Point & p,
1286  const bool add_p_level)
1287 {
1288  return FE<3,LAGRANGE_VEC>::shape(elem, order, i, p, add_p_level);
1289 }
1290 template <> RealGradient FE<3,L2_LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
1291  const unsigned int i, const unsigned int j,
1292  const Point & p,
1293  const bool add_p_level)
1294 {
1295  return FE<3,LAGRANGE_VEC>::shape_deriv(elem, order, i, j, p, add_p_level);
1296 }
1297 
1298 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1299 
1301  const unsigned int i, const unsigned int j,
1302  const Point & p,
1303  const bool add_p_level)
1304 {
1305  return FE<3,LAGRANGE_VEC>::shape_second_deriv(elem, order, i, j, p, add_p_level);
1306 }
1307 
1308 #endif
1309 
1310 // Do full-specialization for every dimension, instead
1311 // of explicit instantiation at the end of this function.
1312 // This could be macro-ified.
1313 template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<0,LAGRANGE>::n_dofs(t,o); }
1314 template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<1,LAGRANGE>::n_dofs(t,o); }
1315 template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 2*FE<2,LAGRANGE>::n_dofs(t,o); }
1316 template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 3*FE<3,LAGRANGE>::n_dofs(t,o); }
1317 
1318 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); }
1319 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); }
1320 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); }
1321 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); }
1322 
1323 
1324 // Do full-specialization for every dimension, instead
1325 // of explicit instantiation at the end of this function.
1326 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); }
1327 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); }
1328 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); }
1329 template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return 3*FE<2,LAGRANGE>::n_dofs_at_node(t,o,n); }
1330 
1331 template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
1332 template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
1333 template <> unsigned int FE<2,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
1334 template <> unsigned int FE<3,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
1335 
1336 
1337 // Lagrange elements have no dofs per element
1338 // (just at the nodes)
1339 template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1340 template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1341 template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1342 template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1343 
1344 // L2 Lagrange elements have all their dofs on the element
1345 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); }
1346 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); }
1347 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); }
1348 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); }
1349 
1350 // Lagrange FEMs are always C^0 continuous
1355 
1356 // L2 Lagrange FEMs are always discontinuous
1361 
1362 // Lagrange FEMs are not hierarchic
1363 template <> bool FE<0,LAGRANGE_VEC>::is_hierarchic() const { return false; }
1364 template <> bool FE<1,LAGRANGE_VEC>::is_hierarchic() const { return false; }
1365 template <> bool FE<2,LAGRANGE_VEC>::is_hierarchic() const { return false; }
1366 template <> bool FE<3,LAGRANGE_VEC>::is_hierarchic() const { return false; }
1367 template <> bool FE<0,L2_LAGRANGE_VEC>::is_hierarchic() const { return false; }
1368 template <> bool FE<1,L2_LAGRANGE_VEC>::is_hierarchic() const { return false; }
1369 template <> bool FE<2,L2_LAGRANGE_VEC>::is_hierarchic() const { return false; }
1370 template <> bool FE<3,L2_LAGRANGE_VEC>::is_hierarchic() const { return false; }
1371 
1372 // Lagrange FEM shapes do not need reinit (is this always true?)
1373 template <> bool FE<0,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1374 template <> bool FE<1,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1375 template <> bool FE<2,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1376 template <> bool FE<3,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1377 template <> bool FE<0,L2_LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1378 template <> bool FE<1,L2_LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1379 template <> bool FE<2,L2_LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1380 template <> bool FE<3,L2_LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
1381 
1382 // Methods for computing Lagrange constraints. Note: we pass the
1383 // dimension as the last argument to the anonymous helper function.
1384 // Also note: we only need instantiations of this function for
1385 // Dim==2 and 3.
1386 #ifdef LIBMESH_ENABLE_AMR
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 
1414 template <>
1416  DofMap & dof_map,
1417  const unsigned int variable_number,
1418  const Elem * elem)
1419 { //libmesh_not_implemented();
1420  FEVectorBase::compute_proj_constraints(constraints, dof_map, variable_number, elem);
1421 }
1422 #endif // LIBMESH_ENABLE_AMR
1423 
1424 } // 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
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:2945
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
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:169
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
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 nodal_soln(const Elem *elem, const Order o, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln, bool add_p_level=true)
Build the nodal soln from the element soln.
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:1575
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...
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:98
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)