libMesh
fe_lagrange_vec.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2019 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/fe.h"
23 #include "libmesh/fe_interface.h"
24 #include "libmesh/elem.h"
25 #include "libmesh/tensor_value.h"
26 
27 namespace libMesh
28 {
29 
30 // ------------------------------------------------------------
31 // Lagrange-specific implementations
32 
33 
34 // Anonymous namespace for local helper functions
35 namespace {
36 void lagrange_vec_nodal_soln(const Elem * elem,
37  const Order order,
38  const std::vector<Number> & elem_soln,
39  const int dim,
40  std::vector<Number> & nodal_soln)
41 {
42  const unsigned int n_nodes = elem->n_nodes();
43  const ElemType type = elem->type();
44 
45  const Order totalorder = static_cast<Order>(order+elem->p_level());
46 
47  nodal_soln.resize(dim*n_nodes);
48 
49  switch (totalorder)
50  {
51  // linear Lagrange shape functions
52  case FIRST:
53  {
54  switch (type)
55  {
56  case TRI6:
57  {
58  libmesh_assert_equal_to (elem_soln.size(), 2*3);
59  libmesh_assert_equal_to (nodal_soln.size(), 2*6);
60 
61  // node 0 components
62  nodal_soln[0] = elem_soln[0];
63  nodal_soln[1] = elem_soln[1];
64 
65  // node 1 components
66  nodal_soln[2] = elem_soln[2];
67  nodal_soln[3] = elem_soln[3];
68 
69  // node 2 components
70  nodal_soln[4] = elem_soln[4];
71  nodal_soln[5] = elem_soln[5];
72 
73  // node 3 components
74  nodal_soln[6] = .5*(elem_soln[0] + elem_soln[2]);
75  nodal_soln[7] = .5*(elem_soln[1] + elem_soln[3]);
76 
77  // node 4 components
78  nodal_soln[8] = .5*(elem_soln[2] + elem_soln[4]);
79  nodal_soln[9] = .5*(elem_soln[3] + elem_soln[5]);
80 
81  // node 5 components
82  nodal_soln[10] = .5*(elem_soln[0] + elem_soln[4]);
83  nodal_soln[11] = .5*(elem_soln[1] + elem_soln[5]);
84 
85  return;
86  }
87 
88 
89  case QUAD8:
90  case QUAD9:
91  {
92  libmesh_assert_equal_to (elem_soln.size(), 2*4);
93 
94  if (type == QUAD8)
95  libmesh_assert_equal_to (nodal_soln.size(), 2*8);
96  else
97  libmesh_assert_equal_to (nodal_soln.size(), 2*9);
98 
99  // node 0 components
100  nodal_soln[0] = elem_soln[0];
101  nodal_soln[1] = elem_soln[1];
102 
103  // node 1 components
104  nodal_soln[2] = elem_soln[2];
105  nodal_soln[3] = elem_soln[3];
106 
107  // node 2 components
108  nodal_soln[4] = elem_soln[4];
109  nodal_soln[5] = elem_soln[5];
110 
111  // node 3 components
112  nodal_soln[6] = elem_soln[6];
113  nodal_soln[7] = elem_soln[7];
114 
115  // node 4 components
116  nodal_soln[8] = .5*(elem_soln[0] + elem_soln[2]);
117  nodal_soln[9] = .5*(elem_soln[1] + elem_soln[3]);
118 
119  // node 5 components
120  nodal_soln[10] = .5*(elem_soln[2] + elem_soln[4]);
121  nodal_soln[11] = .5*(elem_soln[3] + elem_soln[5]);
122 
123  // node 6 components
124  nodal_soln[12] = .5*(elem_soln[4] + elem_soln[6]);
125  nodal_soln[13] = .5*(elem_soln[5] + elem_soln[7]);
126 
127  // node 7 components
128  nodal_soln[14] = .5*(elem_soln[6] + elem_soln[0]);
129  nodal_soln[15] = .5*(elem_soln[7] + elem_soln[1]);
130 
131  if (type == QUAD9)
132  {
133  // node 8 components
134  nodal_soln[16] = .25*(elem_soln[0] + elem_soln[2] + elem_soln[4] + elem_soln[6]);
135  nodal_soln[17] = .25*(elem_soln[1] + elem_soln[3] + elem_soln[5] + elem_soln[7]);
136  }
137 
138  return;
139  }
140 
141 
142  case TET10:
143  {
144  libmesh_assert_equal_to (elem_soln.size(), 3*4);
145  libmesh_assert_equal_to (nodal_soln.size(), 3*10);
146 
147  // node 0 components
148  nodal_soln[0] = elem_soln[0];
149  nodal_soln[1] = elem_soln[1];
150  nodal_soln[2] = elem_soln[2];
151 
152  // node 1 components
153  nodal_soln[3] = elem_soln[3];
154  nodal_soln[4] = elem_soln[4];
155  nodal_soln[5] = elem_soln[5];
156 
157  // node 2 components
158  nodal_soln[6] = elem_soln[6];
159  nodal_soln[7] = elem_soln[7];
160  nodal_soln[8] = elem_soln[8];
161 
162  // node 3 components
163  nodal_soln[9] = elem_soln[9];
164  nodal_soln[10] = elem_soln[10];
165  nodal_soln[11] = elem_soln[11];
166 
167  // node 4 components
168  nodal_soln[12] = .5*(elem_soln[0] + elem_soln[3]);
169  nodal_soln[13] = .5*(elem_soln[1] + elem_soln[4]);
170  nodal_soln[14] = .5*(elem_soln[2] + elem_soln[5]);
171 
172  // node 5 components
173  nodal_soln[15] = .5*(elem_soln[3] + elem_soln[6]);
174  nodal_soln[16] = .5*(elem_soln[4] + elem_soln[7]);
175  nodal_soln[17] = .5*(elem_soln[5] + elem_soln[8]);
176 
177  // node 6 components
178  nodal_soln[18] = .5*(elem_soln[6] + elem_soln[0]);
179  nodal_soln[19] = .5*(elem_soln[7] + elem_soln[1]);
180  nodal_soln[20] = .5*(elem_soln[8] + elem_soln[2]);
181 
182  // node 7 components
183  nodal_soln[21] = .5*(elem_soln[9] + elem_soln[0]);
184  nodal_soln[22] = .5*(elem_soln[10] + elem_soln[1]);
185  nodal_soln[23] = .5*(elem_soln[11] + elem_soln[2]);
186 
187  // node 8 components
188  nodal_soln[24] = .5*(elem_soln[9] + elem_soln[3]);
189  nodal_soln[25] = .5*(elem_soln[10] + elem_soln[4]);
190  nodal_soln[26] = .5*(elem_soln[11] + elem_soln[5]);
191 
192  // node 9 components
193  nodal_soln[27] = .5*(elem_soln[9] + elem_soln[6]);
194  nodal_soln[28] = .5*(elem_soln[10] + elem_soln[7]);
195  nodal_soln[29] = .5*(elem_soln[11] + elem_soln[8]);
196 
197  return;
198  }
199 
200 
201  case HEX20:
202  case HEX27:
203  {
204  libmesh_assert_equal_to (elem_soln.size(), 3*8);
205 
206  if (type == HEX20)
207  libmesh_assert_equal_to (nodal_soln.size(), 3*20);
208  else
209  libmesh_assert_equal_to (nodal_soln.size(), 3*27);
210 
211  // node 0 components
212  nodal_soln[0] = elem_soln[0];
213  nodal_soln[1] = elem_soln[1];
214  nodal_soln[2] = elem_soln[2];
215 
216  // node 1 components
217  nodal_soln[3] = elem_soln[3];
218  nodal_soln[4] = elem_soln[4];
219  nodal_soln[5] = elem_soln[5];
220 
221  // node 2 components
222  nodal_soln[6] = elem_soln[6];
223  nodal_soln[7] = elem_soln[7];
224  nodal_soln[8] = elem_soln[8];
225 
226  // node 3 components
227  nodal_soln[9] = elem_soln[9];
228  nodal_soln[10] = elem_soln[10];
229  nodal_soln[11] = elem_soln[11];
230 
231  // node 4 components
232  nodal_soln[12] = elem_soln[12];
233  nodal_soln[13] = elem_soln[13];
234  nodal_soln[14] = elem_soln[14];
235 
236  // node 5 components
237  nodal_soln[15] = elem_soln[15];
238  nodal_soln[16] = elem_soln[16];
239  nodal_soln[17] = elem_soln[17];
240 
241  // node 6 components
242  nodal_soln[18] = elem_soln[18];
243  nodal_soln[19] = elem_soln[19];
244  nodal_soln[20] = elem_soln[20];
245 
246  // node 7 components
247  nodal_soln[21] = elem_soln[21];
248  nodal_soln[22] = elem_soln[22];
249  nodal_soln[23] = elem_soln[23];
250 
251  // node 8 components
252  nodal_soln[24] = .5*(elem_soln[0] + elem_soln[3]);
253  nodal_soln[25] = .5*(elem_soln[1] + elem_soln[4]);
254  nodal_soln[26] = .5*(elem_soln[2] + elem_soln[5]);
255 
256  // node 9 components
257  nodal_soln[27] = .5*(elem_soln[3] + elem_soln[6]);
258  nodal_soln[28] = .5*(elem_soln[4] + elem_soln[7]);
259  nodal_soln[29] = .5*(elem_soln[5] + elem_soln[8]);
260 
261  // node 10 components
262  nodal_soln[30] = .5*(elem_soln[6] + elem_soln[9]);
263  nodal_soln[31] = .5*(elem_soln[7] + elem_soln[10]);
264  nodal_soln[32] = .5*(elem_soln[8] + elem_soln[11]);
265 
266  // node 11 components
267  nodal_soln[33] = .5*(elem_soln[9] + elem_soln[0]);
268  nodal_soln[34] = .5*(elem_soln[10] + elem_soln[1]);
269  nodal_soln[35] = .5*(elem_soln[11] + elem_soln[2]);
270 
271  // node 12 components
272  nodal_soln[36] = .5*(elem_soln[0] + elem_soln[12]);
273  nodal_soln[37] = .5*(elem_soln[1] + elem_soln[13]);
274  nodal_soln[38] = .5*(elem_soln[2] + elem_soln[14]);
275 
276  // node 13 components
277  nodal_soln[39] = .5*(elem_soln[3] + elem_soln[15]);
278  nodal_soln[40] = .5*(elem_soln[4] + elem_soln[16]);
279  nodal_soln[41] = .5*(elem_soln[5] + elem_soln[17]);
280 
281  // node 14 components
282  nodal_soln[42] = .5*(elem_soln[6] + elem_soln[18]);
283  nodal_soln[43] = .5*(elem_soln[7] + elem_soln[19]);
284  nodal_soln[44] = .5*(elem_soln[8] + elem_soln[20]);
285 
286  // node 15 components
287  nodal_soln[45] = .5*(elem_soln[9] + elem_soln[21]);
288  nodal_soln[46] = .5*(elem_soln[10] + elem_soln[22]);
289  nodal_soln[47] = .5*(elem_soln[11] + elem_soln[23]);
290 
291  // node 16 components
292  nodal_soln[48] = .5*(elem_soln[12] + elem_soln[15]);
293  nodal_soln[49] = .5*(elem_soln[13] + elem_soln[16]);
294  nodal_soln[50] = .5*(elem_soln[14] + elem_soln[17]);
295 
296  // node 17 components
297  nodal_soln[51] = .5*(elem_soln[15] + elem_soln[18]);
298  nodal_soln[52] = .5*(elem_soln[16] + elem_soln[19]);
299  nodal_soln[53] = .5*(elem_soln[17] + elem_soln[20]);
300 
301  // node 18 components
302  nodal_soln[54] = .5*(elem_soln[18] + elem_soln[21]);
303  nodal_soln[55] = .5*(elem_soln[19] + elem_soln[22]);
304  nodal_soln[56] = .5*(elem_soln[20] + elem_soln[23]);
305 
306  // node 19 components
307  nodal_soln[57] = .5*(elem_soln[12] + elem_soln[21]);
308  nodal_soln[58] = .5*(elem_soln[13] + elem_soln[22]);
309  nodal_soln[59] = .5*(elem_soln[14] + elem_soln[23]);
310 
311  if (type == HEX27)
312  {
313  // node 20 components
314  nodal_soln[60] = .25*(elem_soln[0] + elem_soln[3] + elem_soln[6] + elem_soln[9]);
315  nodal_soln[61] = .25*(elem_soln[1] + elem_soln[4] + elem_soln[7] + elem_soln[10]);
316  nodal_soln[62] = .25*(elem_soln[2] + elem_soln[5] + elem_soln[8] + elem_soln[11]);
317 
318  // node 21 components
319  nodal_soln[63] = .25*(elem_soln[0] + elem_soln[3] + elem_soln[12] + elem_soln[15]);
320  nodal_soln[64] = .25*(elem_soln[1] + elem_soln[4] + elem_soln[13] + elem_soln[16]);
321  nodal_soln[65] = .25*(elem_soln[2] + elem_soln[5] + elem_soln[14] + elem_soln[17]);
322 
323  // node 22 components
324  nodal_soln[66] = .25*(elem_soln[3] + elem_soln[6] + elem_soln[15] + elem_soln[18]);
325  nodal_soln[67] = .25*(elem_soln[4] + elem_soln[7] + elem_soln[16] + elem_soln[19]);
326  nodal_soln[68] = .25*(elem_soln[5] + elem_soln[8] + elem_soln[17] + elem_soln[20]);
327 
328  // node 23 components
329  nodal_soln[69] = .25*(elem_soln[6] + elem_soln[9] + elem_soln[18] + elem_soln[21]);
330  nodal_soln[70] = .25*(elem_soln[7] + elem_soln[10] + elem_soln[19] + elem_soln[22]);
331  nodal_soln[71] = .25*(elem_soln[8] + elem_soln[11] + elem_soln[20] + elem_soln[23]);
332 
333  // node 24 components
334  nodal_soln[72] = .25*(elem_soln[9] + elem_soln[0] + elem_soln[21] + elem_soln[12]);
335  nodal_soln[73] = .25*(elem_soln[10] + elem_soln[1] + elem_soln[22] + elem_soln[13]);
336  nodal_soln[74] = .25*(elem_soln[11] + elem_soln[2] + elem_soln[23] + elem_soln[14]);
337 
338  // node 25 components
339  nodal_soln[75] = .25*(elem_soln[12] + elem_soln[15] + elem_soln[18] + elem_soln[21]);
340  nodal_soln[76] = .25*(elem_soln[13] + elem_soln[16] + elem_soln[19] + elem_soln[22]);
341  nodal_soln[77] = .25*(elem_soln[14] + elem_soln[17] + elem_soln[20] + elem_soln[23]);
342 
343  // node 26 components
344  nodal_soln[78] = .125*(elem_soln[0] + elem_soln[3] + elem_soln[6] + elem_soln[9] +
345  elem_soln[12] + elem_soln[15] + elem_soln[18] + elem_soln[21]);
346 
347  nodal_soln[79] = .125*(elem_soln[1] + elem_soln[4] + elem_soln[7] + elem_soln[10] +
348  elem_soln[13] + elem_soln[16] + elem_soln[19] + elem_soln[22]);
349 
350  nodal_soln[80] = .125*(elem_soln[2] + elem_soln[5] + elem_soln[8] + elem_soln[11] +
351  elem_soln[14] + elem_soln[17] + elem_soln[20] + elem_soln[23]);
352  }
353 
354  return;
355  }
356 
357 
358  case PRISM15:
359  case PRISM18:
360  {
361  libmesh_assert_equal_to (elem_soln.size(), 3*6);
362 
363  if (type == PRISM15)
364  libmesh_assert_equal_to (nodal_soln.size(), 3*15);
365  else
366  libmesh_assert_equal_to (nodal_soln.size(), 3*18);
367 
368  // node 0 components
369  nodal_soln[0] = elem_soln[0];
370  nodal_soln[1] = elem_soln[1];
371  nodal_soln[2] = elem_soln[2];
372 
373  // node 1 components
374  nodal_soln[3] = elem_soln[3];
375  nodal_soln[4] = elem_soln[4];
376  nodal_soln[5] = elem_soln[5];
377 
378  // node 2 components
379  nodal_soln[6] = elem_soln[6];
380  nodal_soln[7] = elem_soln[7];
381  nodal_soln[8] = elem_soln[8];
382 
383  // node 3 components
384  nodal_soln[9] = elem_soln[9];
385  nodal_soln[10] = elem_soln[10];
386  nodal_soln[11] = elem_soln[11];
387 
388  // node 4 components
389  nodal_soln[12] = elem_soln[12];
390  nodal_soln[13] = elem_soln[13];
391  nodal_soln[14] = elem_soln[14];
392 
393  // node 5 components
394  nodal_soln[15] = elem_soln[15];
395  nodal_soln[16] = elem_soln[16];
396  nodal_soln[17] = elem_soln[17];
397 
398  // node 6 components
399  nodal_soln[18] = .5*(elem_soln[0] + elem_soln[3]);
400  nodal_soln[19] = .5*(elem_soln[1] + elem_soln[4]);
401  nodal_soln[20] = .5*(elem_soln[2] + elem_soln[5]);
402 
403  // node 7 components
404  nodal_soln[21] = .5*(elem_soln[3] + elem_soln[6]);
405  nodal_soln[22] = .5*(elem_soln[4] + elem_soln[7]);
406  nodal_soln[23] = .5*(elem_soln[5] + elem_soln[8]);
407 
408  // node 8 components
409  nodal_soln[24] = .5*(elem_soln[0] + elem_soln[6]);
410  nodal_soln[25] = .5*(elem_soln[1] + elem_soln[7]);
411  nodal_soln[26] = .5*(elem_soln[2] + elem_soln[8]);
412 
413  // node 9 components
414  nodal_soln[27] = .5*(elem_soln[0] + elem_soln[9]);
415  nodal_soln[28] = .5*(elem_soln[1] + elem_soln[10]);
416  nodal_soln[29] = .5*(elem_soln[2] + elem_soln[11]);
417 
418  // node 10 components
419  nodal_soln[30] = .5*(elem_soln[3] + elem_soln[12]);
420  nodal_soln[31] = .5*(elem_soln[4] + elem_soln[13]);
421  nodal_soln[32] = .5*(elem_soln[5] + elem_soln[14]);
422 
423  // node 11 components
424  nodal_soln[33] = .5*(elem_soln[6] + elem_soln[15]);
425  nodal_soln[34] = .5*(elem_soln[7] + elem_soln[16]);
426  nodal_soln[35] = .5*(elem_soln[8] + elem_soln[17]);
427 
428  // node 12 components
429  nodal_soln[36] = .5*(elem_soln[9] + elem_soln[12]);
430  nodal_soln[37] = .5*(elem_soln[10] + elem_soln[13]);
431  nodal_soln[38] = .5*(elem_soln[11] + elem_soln[14]);
432 
433  // node 13 components
434  nodal_soln[39] = .5*(elem_soln[12] + elem_soln[15]);
435  nodal_soln[40] = .5*(elem_soln[13] + elem_soln[16]);
436  nodal_soln[41] = .5*(elem_soln[14] + elem_soln[17]);
437 
438  // node 14 components
439  nodal_soln[42] = .5*(elem_soln[12] + elem_soln[15]);
440  nodal_soln[43] = .5*(elem_soln[13] + elem_soln[16]);
441  nodal_soln[44] = .5*(elem_soln[14] + elem_soln[17]);
442 
443  if (type == PRISM18)
444  {
445  // node 15 components
446  nodal_soln[45] = .25*(elem_soln[0] + elem_soln[3] + elem_soln[12] + elem_soln[9]);
447  nodal_soln[46] = .25*(elem_soln[1] + elem_soln[4] + elem_soln[13] + elem_soln[10]);
448  nodal_soln[47] = .25*(elem_soln[2] + elem_soln[5] + elem_soln[14] + elem_soln[11]);
449 
450  // node 16 components
451  nodal_soln[48] = .25*(elem_soln[3] + elem_soln[6] + elem_soln[15] + elem_soln[12]);
452  nodal_soln[49] = .25*(elem_soln[4] + elem_soln[7] + elem_soln[16] + elem_soln[13]);
453  nodal_soln[50] = .25*(elem_soln[5] + elem_soln[8] + elem_soln[17] + elem_soln[14]);
454 
455  // node 17 components
456  nodal_soln[51] = .25*(elem_soln[6] + elem_soln[0] + elem_soln[9] + elem_soln[15]);
457  nodal_soln[52] = .25*(elem_soln[7] + elem_soln[1] + elem_soln[10] + elem_soln[16]);
458  nodal_soln[53] = .25*(elem_soln[8] + elem_soln[2] + elem_soln[11] + elem_soln[17]);
459  }
460 
461  return;
462  }
463 
464  default:
465  {
466  // By default the element solution _is_ nodal,
467  // so just copy it.
468  nodal_soln = elem_soln;
469 
470  return;
471  }
472  }
473  }
474 
475  case SECOND:
476  {
477  switch (type)
478  {
479  default:
480  {
481  // By default the element solution _is_ nodal,
482  // so just copy it.
483  nodal_soln = elem_soln;
484 
485  return;
486  }
487  }
488  }
489 
490  default:
491  {
492 
493  }
494 
495  } // switch(totalorder)
496 
497 }// void lagrange_vec_nodal_soln
498 
499 } // anonymous namespace
500 
501 
502  // Do full-specialization for every dimension, instead
503  // of explicit instantiation at the end of this file.
504  // This could be macro-ified so that it fits on one line...
505 template <>
506 void FE<0,LAGRANGE_VEC>::nodal_soln(const Elem * elem,
507  const Order order,
508  const std::vector<Number> & elem_soln,
509  std::vector<Number> & nodal_soln)
510 { FE<0,LAGRANGE>::nodal_soln(elem, order, elem_soln, nodal_soln); }
511 
512 template <>
513 void FE<1,LAGRANGE_VEC>::nodal_soln(const Elem * elem,
514  const Order order,
515  const std::vector<Number> & elem_soln,
516  std::vector<Number> & nodal_soln)
517 { FE<1,LAGRANGE>::nodal_soln(elem, order, elem_soln, nodal_soln); }
518 
519 template <>
520 void FE<2,LAGRANGE_VEC>::nodal_soln(const Elem * elem,
521  const Order order,
522  const std::vector<Number> & elem_soln,
523  std::vector<Number> & nodal_soln)
524 { lagrange_vec_nodal_soln(elem, order, elem_soln, 2 /*dimension*/, nodal_soln); }
525 
526 template <>
527 void FE<3,LAGRANGE_VEC>::nodal_soln(const Elem * elem,
528  const Order order,
529  const std::vector<Number> & elem_soln,
530  std::vector<Number> & nodal_soln)
531 { lagrange_vec_nodal_soln(elem, order, elem_soln, 3 /*dimension*/, nodal_soln); }
532 
533 
534 // Specialize for shape function routines by leveraging scalar LAGRANGE elements
535 
536 // 0-D
537 template <> RealGradient FE<0,LAGRANGE_VEC>::shape(const ElemType type, const Order order,
538  const unsigned int i, const Point & p)
539 {
540  Real value = FE<0,LAGRANGE>::shape( type, order, i, p );
541  return libMesh::RealGradient( value );
542 }
543 template <> RealGradient FE<0,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
544  const unsigned int i, const unsigned int j,
545  const Point & p)
546 {
547  Real value = FE<0,LAGRANGE>::shape_deriv( type, order, i, j, p );
548  return libMesh::RealGradient( value );
549 }
550 
551 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
552 
553 template <> RealGradient FE<0,LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order,
554  const unsigned int i, const unsigned int j,
555  const Point & p)
556 {
557  Real value = FE<0,LAGRANGE>::shape_second_deriv( type, order, i, j, p );
558  return libMesh::RealGradient( value );
559 }
560 #endif
561 
562 // 1-D
563 template <> RealGradient FE<1,LAGRANGE_VEC>::shape(const ElemType type, const Order order,
564  const unsigned int i, const Point & p)
565 {
566  Real value = FE<1,LAGRANGE>::shape( type, order, i, p );
567  return libMesh::RealGradient( value );
568 }
569 template <> RealGradient FE<1,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
570  const unsigned int i, const unsigned int j,
571  const Point & p)
572 {
573  Real value = FE<1,LAGRANGE>::shape_deriv( type, order, i, j, p );
574  return libMesh::RealGradient( value );
575 }
576 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
577 template <> RealGradient FE<1,LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order,
578  const unsigned int i, const unsigned int j,
579  const Point & p)
580 {
581  Real value = FE<1,LAGRANGE>::shape_second_deriv( type, order, i, j, p );
582  return libMesh::RealGradient( value );
583 }
584 
585 #endif
586 
587 // 2-D
588 template <> RealGradient FE<2,LAGRANGE_VEC>::shape(const ElemType type, const Order order,
589  const unsigned int i, const Point & p)
590 {
591  Real value = FE<2,LAGRANGE>::shape( type, order, i/2, p );
592 
593  switch( i%2 )
594  {
595  case 0:
596  return libMesh::RealGradient( value );
597 
598  case 1:
599  return libMesh::RealGradient( Real(0), value );
600 
601  default:
602  libmesh_error_msg("i%2 must be either 0 or 1!");
603  }
604 
605  //dummy
606  return libMesh::RealGradient();
607 }
608 template <> RealGradient FE<2,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
609  const unsigned int i, const unsigned int j,
610  const Point & p)
611 {
612  Real value = FE<2,LAGRANGE>::shape_deriv( type, order, i/2, j, p );
613 
614  switch( i%2 )
615  {
616  case 0:
617  return libMesh::RealGradient( value );
618 
619  case 1:
620  return libMesh::RealGradient( Real(0), value );
621 
622  default:
623  libmesh_error_msg("i%2 must be either 0 or 1!");
624  }
625 
626  //dummy
627  return libMesh::RealGradient();
628 }
629 
630 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
631 template <> RealGradient FE<2,LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order,
632  const unsigned int i, const unsigned int j,
633  const Point & p)
634 {
635  Real value = FE<2,LAGRANGE>::shape_second_deriv( type, order, i/2, j, p );
636 
637  switch( i%2 )
638  {
639  case 0:
640  return libMesh::RealGradient( value );
641 
642  case 1:
643  return libMesh::RealGradient( Real(0), value );
644 
645  default:
646  libmesh_error_msg("i%2 must be either 0 or 1!");
647  }
648 
649  //dummy
650  return libMesh::RealGradient();
651 }
652 
653 #endif
654 
655 
656 // 3-D
657 template <> RealGradient FE<3,LAGRANGE_VEC>::shape(const ElemType type, const Order order,
658  const unsigned int i, const Point & p)
659 {
660  Real value = FE<3,LAGRANGE>::shape( type, order, i/3, p );
661 
662  switch( i%3 )
663  {
664  case 0:
665  return libMesh::RealGradient( value );
666 
667  case 1:
668  return libMesh::RealGradient( Real(0), value );
669 
670  case 2:
671  return libMesh::RealGradient( Real(0), Real(0), value );
672 
673  default:
674  libmesh_error_msg("i%3 must be 0, 1, or 2!");
675  }
676 
677  //dummy
678  return libMesh::RealGradient();
679 }
680 template <> RealGradient FE<3,LAGRANGE_VEC>::shape_deriv(const ElemType type, const Order order,
681  const unsigned int i, const unsigned int j,
682  const Point & p)
683 {
684  Real value = FE<3,LAGRANGE>::shape_deriv( type, order, i/3, j, p );
685 
686  switch( i%3 )
687  {
688  case 0:
689  return libMesh::RealGradient( value );
690 
691  case 1:
692  return libMesh::RealGradient( Real(0), value );
693 
694  case 2:
695  return libMesh::RealGradient( Real(0), Real(0), value );
696 
697  default:
698  libmesh_error_msg("i%3 must be 0, 1, or 2!");
699  }
700 
701  //dummy
702  return libMesh::RealGradient();
703 }
704 
705 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
706 
707 template <> RealGradient FE<3,LAGRANGE_VEC>::shape_second_deriv(const ElemType type, const Order order,
708  const unsigned int i, const unsigned int j,
709  const Point & p)
710 {
711  Real value = FE<3,LAGRANGE>::shape_second_deriv( type, order, i/3, j, p );
712 
713  switch( i%3 )
714  {
715  case 0:
716  return libMesh::RealGradient( value );
717 
718  case 1:
719  return libMesh::RealGradient( Real(0), value );
720 
721  case 2:
722  return libMesh::RealGradient( Real(0), Real(0), value );
723 
724  default:
725  libmesh_error_msg("i%3 must be 0, 1, or 2!");
726  }
727 
728  //dummy
729  return libMesh::RealGradient();
730 }
731 
732 #endif
733 
734 
735 // 0-D
736 template <> RealGradient FE<0,LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
737  const unsigned int i, const Point & p,
738  const bool add_p_level)
739 {
740  Real value = FE<0,LAGRANGE>::shape( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, p);
741  return libMesh::RealGradient( value );
742 }
743 template <> RealGradient FE<0,LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
744  const unsigned int i, const unsigned int j,
745  const Point & p,
746  const bool add_p_level)
747 {
748  Real value = FE<0,LAGRANGE>::shape_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
749  return libMesh::RealGradient( value );
750 }
751 
752 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
753 
754 template <> RealGradient FE<0,LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
755  const unsigned int i, const unsigned int j,
756  const Point & p,
757  const bool add_p_level)
758 {
759  Real value = FE<0,LAGRANGE>::shape_second_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
760  return libMesh::RealGradient( value );
761 }
762 
763 #endif
764 
765 // 1-D
766 template <> RealGradient FE<1,LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
767  const unsigned int i, const Point & p,
768  const bool add_p_level)
769 {
770  Real value = FE<1,LAGRANGE>::shape( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, p);
771  return libMesh::RealGradient( value );
772 }
773 template <> RealGradient FE<1,LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
774  const unsigned int i, const unsigned int j,
775  const Point & p,
776  const bool add_p_level)
777 {
778  Real value = FE<1,LAGRANGE>::shape_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
779  return libMesh::RealGradient( value );
780 }
781 
782 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
783 template <> RealGradient FE<1,LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
784  const unsigned int i, const unsigned int j,
785  const Point & p,
786  const bool add_p_level)
787 {
788  Real value = FE<1,LAGRANGE>::shape_second_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
789  return libMesh::RealGradient( value );
790 }
791 
792 #endif
793 
794 // 2-D
795 template <> RealGradient FE<2,LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
796  const unsigned int i, const Point & p,
797  const bool add_p_level)
798 {
799  Real value = FE<2,LAGRANGE>::shape( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i/2, p );
800 
801  switch( i%2 )
802  {
803  case 0:
804  return libMesh::RealGradient( value );
805 
806  case 1:
807  return libMesh::RealGradient( Real(0), value );
808 
809  default:
810  libmesh_error_msg("i%2 must be either 0 or 1!");
811  }
812 
813  //dummy
814  return libMesh::RealGradient();
815 }
816 template <> RealGradient FE<2,LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
817  const unsigned int i, const unsigned int j,
818  const Point & p,
819  const bool add_p_level)
820 {
821  Real value = FE<2,LAGRANGE>::shape_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i/2, j, p );
822 
823  switch( i%2 )
824  {
825  case 0:
826  return libMesh::RealGradient( value );
827 
828  case 1:
829  return libMesh::RealGradient( Real(0), value );
830 
831  default:
832  libmesh_error_msg("i%2 must be either 0 or 1!");
833  }
834 
835  //dummy
836  return libMesh::RealGradient();
837 }
838 
839 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
840 template <> RealGradient FE<2,LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
841  const unsigned int i, const unsigned int j,
842  const Point & p,
843  const bool add_p_level)
844 {
845  Real value = FE<2,LAGRANGE>::shape_second_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i/2, j, p );
846 
847  switch( i%2 )
848  {
849  case 0:
850  return libMesh::RealGradient( value );
851 
852  case 1:
853  return libMesh::RealGradient( Real(0), value );
854 
855  default:
856  libmesh_error_msg("i%2 must be either 0 or 1!");
857  }
858 
859  //dummy
860  return libMesh::RealGradient();
861 }
862 
863 #endif
864 
865 // 3-D
866 template <> RealGradient FE<3,LAGRANGE_VEC>::shape(const Elem * elem, const Order order,
867  const unsigned int i, const Point & p,
868  const bool add_p_level)
869 {
870  Real value = FE<3,LAGRANGE>::shape( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i/3, p );
871 
872  switch( i%3 )
873  {
874  case 0:
875  return libMesh::RealGradient( value );
876 
877  case 1:
878  return libMesh::RealGradient( Real(0), value );
879 
880  case 2:
881  return libMesh::RealGradient( Real(0), Real(0), value );
882 
883  default:
884  libmesh_error_msg("i%3 must be 0, 1, or 2!");
885  }
886 
887  //dummy
888  return libMesh::RealGradient();
889 }
890 template <> RealGradient FE<3,LAGRANGE_VEC>::shape_deriv(const Elem * elem, const Order order,
891  const unsigned int i, const unsigned int j,
892  const Point & p,
893  const bool add_p_level)
894 {
895  Real value = FE<3,LAGRANGE>::shape_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i/3, j, p );
896 
897  switch( i%3 )
898  {
899  case 0:
900  return libMesh::RealGradient( value );
901 
902  case 1:
903  return libMesh::RealGradient( Real(0), value );
904 
905  case 2:
906  return libMesh::RealGradient( Real(0), Real(0), value );
907 
908  default:
909  libmesh_error_msg("i%3 must be 0, 1, or 2!");
910  }
911 
912  //dummy
913  return libMesh::RealGradient();
914 }
915 
916 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
917 
918 template <> RealGradient FE<3,LAGRANGE_VEC>::shape_second_deriv(const Elem * elem, const Order order,
919  const unsigned int i, const unsigned int j,
920  const Point & p,
921  const bool add_p_level)
922 {
923  Real value = FE<3,LAGRANGE>::shape_second_deriv( elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i/3, j, p );
924 
925  switch( i%3 )
926  {
927  case 0:
928  return libMesh::RealGradient( value );
929 
930  case 1:
931  return libMesh::RealGradient( Real(0), value );
932 
933  case 2:
934  return libMesh::RealGradient( Real(0), Real(0), value );
935 
936  default:
937  libmesh_error_msg("i%3 must be 0, 1, or 2!");
938  }
939 
940  //dummy
941  return libMesh::RealGradient();
942 }
943 
944 #endif
945 
946 // Do full-specialization for every dimension, instead
947 // of explicit instantiation at the end of this function.
948 // This could be macro-ified.
949 template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<0,LAGRANGE>::n_dofs(t,o); }
950 template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<1,LAGRANGE>::n_dofs(t,o); }
951 template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 2*FE<2,LAGRANGE>::n_dofs(t,o); }
952 template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 3*FE<3,LAGRANGE>::n_dofs(t,o); }
953 
954 
955 // Do full-specialization for every dimension, instead
956 // of explicit instantiation at the end of this function.
957 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); }
958 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); }
959 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); }
960 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); }
961 
962 
963 // Lagrange elements have no dofs per element
964 // (just at the nodes)
965 template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
966 template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
967 template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
968 template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
969 
970 // Lagrange FEMs are always C^0 continuous
971 template <> FEContinuity FE<0,LAGRANGE_VEC>::get_continuity() const { return C_ZERO; }
972 template <> FEContinuity FE<1,LAGRANGE_VEC>::get_continuity() const { return C_ZERO; }
973 template <> FEContinuity FE<2,LAGRANGE_VEC>::get_continuity() const { return C_ZERO; }
974 template <> FEContinuity FE<3,LAGRANGE_VEC>::get_continuity() const { return C_ZERO; }
975 
976 // Lagrange FEMs are not hierarchic
977 template <> bool FE<0,LAGRANGE_VEC>::is_hierarchic() const { return false; }
978 template <> bool FE<1,LAGRANGE_VEC>::is_hierarchic() const { return false; }
979 template <> bool FE<2,LAGRANGE_VEC>::is_hierarchic() const { return false; }
980 template <> bool FE<3,LAGRANGE_VEC>::is_hierarchic() const { return false; }
981 
982 // Lagrange FEM shapes do not need reinit (is this always true?)
983 template <> bool FE<0,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
984 template <> bool FE<1,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
985 template <> bool FE<2,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
986 template <> bool FE<3,LAGRANGE_VEC>::shapes_need_reinit() const { return false; }
987 
988 // Methods for computing Lagrange constraints. Note: we pass the
989 // dimension as the last argument to the anonymous helper function.
990 // Also note: we only need instantiations of this function for
991 // Dim==2 and 3.
992 #ifdef LIBMESH_ENABLE_AMR
993 template <>
994 void FE<2,LAGRANGE_VEC>::compute_constraints (DofConstraints & constraints,
995  DofMap & dof_map,
996  const unsigned int variable_number,
997  const Elem * elem)
998 { //libmesh_not_implemented();
999  FEVectorBase::compute_proj_constraints(constraints, dof_map, variable_number, elem);
1000 }
1001 
1002 template <>
1003 void FE<3,LAGRANGE_VEC>::compute_constraints (DofConstraints & constraints,
1004  DofMap & dof_map,
1005  const unsigned int variable_number,
1006  const Elem * elem)
1007 { //libmesh_not_implemented();
1008  FEVectorBase::compute_proj_constraints(constraints, dof_map, variable_number, elem);
1009 }
1010 #endif // LIBMESH_ENABLE_AMR
1011 
1012 } // namespace libMesh
libMesh::DofConstraints
The constraint matrix storage format.
Definition: dof_map.h:105
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::HEX20
Definition: enum_elem_type.h:48
libMesh::RealGradient
RealVectorValue RealGradient
Definition: hp_coarsentest.h:49
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
Definition: factoryfunction.C:55
libMesh::TET10
Definition: enum_elem_type.h:46
libMesh::Order
Order
Definition: enum_order.h:40
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2512
libMesh::SECOND
Definition: enum_order.h:43
dim
unsigned int dim
Definition: adaptivity_ex3.C:113
libMesh::VectorValue< Real >
libMesh::PRISM15
Definition: enum_elem_type.h:51
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::HEX27
Definition: enum_elem_type.h:49
libMesh::FE::shape_deriv
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
libMesh::C_ZERO
Definition: enum_fe_family.h:79
libMesh::FE::compute_constraints
static void compute_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
Computes the constraint matrix contributions (for non-conforming adapted meshes) corresponding to var...
n_nodes
const dof_id_type n_nodes
Definition: tecplot_io.C:68
libMesh::TRI6
Definition: enum_elem_type.h:40
libMesh::FE::shapes_need_reinit
virtual bool shapes_need_reinit() const override
libMesh::FE::n_dofs_at_node
static unsigned int n_dofs_at_node(const ElemType t, const Order o, const unsigned int n)
libMesh::FE::nodal_soln
static void nodal_soln(const Elem *elem, const Order o, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)
Build the nodal soln from the element soln.
libMesh::FE::get_continuity
virtual FEContinuity get_continuity() const override
value
static const bool value
Definition: xdr_io.C:56
libMesh::DofMap
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:176
libMesh::FE::n_dofs_per_elem
static unsigned int n_dofs_per_elem(const ElemType t, const Order o)
libMesh::FEGenericBase::compute_proj_constraints
static void compute_proj_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
Computes the constraint matrix contributions (for non-conforming adapted meshes) corresponding to var...
Definition: fe_base.C:1396
libMesh::FE::is_hierarchic
virtual bool is_hierarchic() const override
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:100
libMesh::FE::n_dofs
static unsigned int n_dofs(const ElemType t, const Order o)
libMesh::QUAD9
Definition: enum_elem_type.h:43
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
libMesh::FEContinuity
FEContinuity
Definition: enum_fe_family.h:77
libMesh::PRISM18
Definition: enum_elem_type.h:52
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::FIRST
Definition: enum_order.h:42
libMesh::QUAD8
Definition: enum_elem_type.h:42
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
Generated on Sat Jan 25 2020 12:06:58 for libMesh by   doxygen 1.8.16