fe_lagrange_vec.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 <>
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 <>
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 <>
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 <>
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  }
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  }
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();
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();
618  }
619 
620  //dummy
621  return libMesh::RealGradient();
622  }
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();
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();
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();
689  }
690 
691  //dummy
692  return libMesh::RealGradient();
693  }
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();
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();
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();
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();
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();
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();
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();
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
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 <>
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 <>
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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo