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