libMesh
fe_l2_lagrange.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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/string_to_enum.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 l2_lagrange_nodal_soln(const Elem * elem,
38  const Order order,
39  const std::vector<Number> & elem_soln,
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(n_nodes);
48 
49 
50 
51  switch (totalorder)
52  {
53  // linear Lagrange shape functions
54  case FIRST:
55  {
56  switch (type)
57  {
58  case EDGE3:
59  {
60  libmesh_assert_equal_to (elem_soln.size(), 2);
61  libmesh_assert_equal_to (nodal_soln.size(), 3);
62 
63  nodal_soln[0] = elem_soln[0];
64  nodal_soln[1] = elem_soln[1];
65  nodal_soln[2] = .5*(elem_soln[0] + elem_soln[1]);
66 
67  return;
68  }
69 
70  case EDGE4:
71  {
72  libmesh_assert_equal_to (elem_soln.size(), 2);
73  libmesh_assert_equal_to (nodal_soln.size(), 4);
74 
75  nodal_soln[0] = elem_soln[0];
76  nodal_soln[1] = elem_soln[1];
77  nodal_soln[2] = (2.*elem_soln[0] + elem_soln[1])/3.;
78  nodal_soln[3] = (elem_soln[0] + 2.*elem_soln[1])/3.;
79 
80  return;
81  }
82 
83 
84  case TRI6:
85  {
86  libmesh_assert_equal_to (elem_soln.size(), 3);
87  libmesh_assert_equal_to (nodal_soln.size(), 6);
88 
89  nodal_soln[0] = elem_soln[0];
90  nodal_soln[1] = elem_soln[1];
91  nodal_soln[2] = elem_soln[2];
92  nodal_soln[3] = .5*(elem_soln[0] + elem_soln[1]);
93  nodal_soln[4] = .5*(elem_soln[1] + elem_soln[2]);
94  nodal_soln[5] = .5*(elem_soln[2] + elem_soln[0]);
95 
96  return;
97  }
98 
99 
100  case QUAD8:
101  case QUAD9:
102  {
103  libmesh_assert_equal_to (elem_soln.size(), 4);
104 
105  if (type == QUAD8)
106  libmesh_assert_equal_to (nodal_soln.size(), 8);
107  else
108  libmesh_assert_equal_to (nodal_soln.size(), 9);
109 
110 
111  nodal_soln[0] = elem_soln[0];
112  nodal_soln[1] = elem_soln[1];
113  nodal_soln[2] = elem_soln[2];
114  nodal_soln[3] = elem_soln[3];
115  nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]);
116  nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]);
117  nodal_soln[6] = .5*(elem_soln[2] + elem_soln[3]);
118  nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]);
119 
120  if (type == QUAD9)
121  nodal_soln[8] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
122 
123  return;
124  }
125 
126 
127  case TET10:
128  {
129  libmesh_assert_equal_to (elem_soln.size(), 4);
130  libmesh_assert_equal_to (nodal_soln.size(), 10);
131 
132  nodal_soln[0] = elem_soln[0];
133  nodal_soln[1] = elem_soln[1];
134  nodal_soln[2] = elem_soln[2];
135  nodal_soln[3] = elem_soln[3];
136  nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]);
137  nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]);
138  nodal_soln[6] = .5*(elem_soln[2] + elem_soln[0]);
139  nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]);
140  nodal_soln[8] = .5*(elem_soln[3] + elem_soln[1]);
141  nodal_soln[9] = .5*(elem_soln[3] + elem_soln[2]);
142 
143  return;
144  }
145 
146 
147  case HEX20:
148  case HEX27:
149  {
150  libmesh_assert_equal_to (elem_soln.size(), 8);
151 
152  if (type == HEX20)
153  libmesh_assert_equal_to (nodal_soln.size(), 20);
154  else
155  libmesh_assert_equal_to (nodal_soln.size(), 27);
156 
157  nodal_soln[0] = elem_soln[0];
158  nodal_soln[1] = elem_soln[1];
159  nodal_soln[2] = elem_soln[2];
160  nodal_soln[3] = elem_soln[3];
161  nodal_soln[4] = elem_soln[4];
162  nodal_soln[5] = elem_soln[5];
163  nodal_soln[6] = elem_soln[6];
164  nodal_soln[7] = elem_soln[7];
165  nodal_soln[8] = .5*(elem_soln[0] + elem_soln[1]);
166  nodal_soln[9] = .5*(elem_soln[1] + elem_soln[2]);
167  nodal_soln[10] = .5*(elem_soln[2] + elem_soln[3]);
168  nodal_soln[11] = .5*(elem_soln[3] + elem_soln[0]);
169  nodal_soln[12] = .5*(elem_soln[0] + elem_soln[4]);
170  nodal_soln[13] = .5*(elem_soln[1] + elem_soln[5]);
171  nodal_soln[14] = .5*(elem_soln[2] + elem_soln[6]);
172  nodal_soln[15] = .5*(elem_soln[3] + elem_soln[7]);
173  nodal_soln[16] = .5*(elem_soln[4] + elem_soln[5]);
174  nodal_soln[17] = .5*(elem_soln[5] + elem_soln[6]);
175  nodal_soln[18] = .5*(elem_soln[6] + elem_soln[7]);
176  nodal_soln[19] = .5*(elem_soln[4] + elem_soln[7]);
177 
178  if (type == HEX27)
179  {
180  nodal_soln[20] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
181  nodal_soln[21] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[5]);
182  nodal_soln[22] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[6]);
183  nodal_soln[23] = .25*(elem_soln[2] + elem_soln[3] + elem_soln[6] + elem_soln[7]);
184  nodal_soln[24] = .25*(elem_soln[3] + elem_soln[0] + elem_soln[7] + elem_soln[4]);
185  nodal_soln[25] = .25*(elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]);
186 
187  nodal_soln[26] = .125*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3] +
188  elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]);
189  }
190 
191  return;
192  }
193 
194 
195  case PRISM15:
196  case PRISM18:
197  {
198  libmesh_assert_equal_to (elem_soln.size(), 6);
199 
200  if (type == PRISM15)
201  libmesh_assert_equal_to (nodal_soln.size(), 15);
202  else
203  libmesh_assert_equal_to (nodal_soln.size(), 18);
204 
205  nodal_soln[0] = elem_soln[0];
206  nodal_soln[1] = elem_soln[1];
207  nodal_soln[2] = elem_soln[2];
208  nodal_soln[3] = elem_soln[3];
209  nodal_soln[4] = elem_soln[4];
210  nodal_soln[5] = elem_soln[5];
211  nodal_soln[6] = .5*(elem_soln[0] + elem_soln[1]);
212  nodal_soln[7] = .5*(elem_soln[1] + elem_soln[2]);
213  nodal_soln[8] = .5*(elem_soln[0] + elem_soln[2]);
214  nodal_soln[9] = .5*(elem_soln[0] + elem_soln[3]);
215  nodal_soln[10] = .5*(elem_soln[1] + elem_soln[4]);
216  nodal_soln[11] = .5*(elem_soln[2] + elem_soln[5]);
217  nodal_soln[12] = .5*(elem_soln[3] + elem_soln[4]);
218  nodal_soln[13] = .5*(elem_soln[4] + elem_soln[5]);
219  nodal_soln[14] = .5*(elem_soln[3] + elem_soln[5]);
220 
221  if (type == PRISM18)
222  {
223  nodal_soln[15] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[3]);
224  nodal_soln[16] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[4]);
225  nodal_soln[17] = .25*(elem_soln[2] + elem_soln[0] + elem_soln[3] + elem_soln[5]);
226  }
227 
228  return;
229  }
230 
231 
232 
233  default:
234  {
235  // By default the element solution _is_ nodal,
236  // so just copy it.
237  nodal_soln = elem_soln;
238 
239  return;
240  }
241  }
242  }
243 
244  case SECOND:
245  {
246  switch (type)
247  {
248  case EDGE4:
249  {
250  libmesh_assert_equal_to (elem_soln.size(), 3);
251  libmesh_assert_equal_to (nodal_soln.size(), 4);
252 
253  // Project quadratic solution onto cubic element nodes
254  nodal_soln[0] = elem_soln[0];
255  nodal_soln[1] = elem_soln[1];
256  nodal_soln[2] = (2.*elem_soln[0] - elem_soln[1] +
257  8.*elem_soln[2])/9.;
258  nodal_soln[3] = (-elem_soln[0] + 2.*elem_soln[1] +
259  8.*elem_soln[2])/9.;
260  return;
261  }
262 
263  default:
264  {
265  // By default the element solution _is_ nodal,
266  // so just copy it.
267  nodal_soln = elem_soln;
268 
269  return;
270  }
271  }
272  }
273 
274 
275 
276 
277  default:
278  {
279  // By default the element solution _is_ nodal,
280  // so just copy it.
281  nodal_soln = elem_soln;
282 
283  return;
284  }
285  }
286 }
287 
288 
289 // TODO: We should make this work, for example, for SECOND on a TRI3
290 // (this is valid with L2_LAGRANGE, but not with LAGRANGE)
291 unsigned int l2_lagrange_n_dofs(const ElemType t, const Order o)
292 {
293  switch (o)
294  {
295 
296  // linear Lagrange shape functions
297  case FIRST:
298  {
299  switch (t)
300  {
301  case NODEELEM:
302  return 1;
303 
304  case EDGE2:
305  case EDGE3:
306  case EDGE4:
307  return 2;
308 
309  case TRI3:
310  case TRISHELL3:
311  case TRI6:
312  return 3;
313 
314  case QUAD4:
315  case QUADSHELL4:
316  case QUAD8:
317  case QUADSHELL8:
318  case QUAD9:
319  return 4;
320 
321  case TET4:
322  case TET10:
323  return 4;
324 
325  case HEX8:
326  case HEX20:
327  case HEX27:
328  return 8;
329 
330  case PRISM6:
331  case PRISM15:
332  case PRISM18:
333  return 6;
334 
335  case PYRAMID5:
336  case PYRAMID13:
337  case PYRAMID14:
338  return 5;
339 
340  case INVALID_ELEM:
341  return 0;
342 
343  default:
344  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
345  }
346  }
347 
348 
349  // quadratic Lagrange shape functions
350  case SECOND:
351  {
352  switch (t)
353  {
354  case NODEELEM:
355  return 1;
356 
357  case EDGE3:
358  return 3;
359 
360  case TRI6:
361  return 6;
362 
363  case QUAD8:
364  case QUADSHELL8:
365  return 8;
366 
367  case QUAD9:
368  return 9;
369 
370  case TET10:
371  return 10;
372 
373  case HEX20:
374  return 20;
375 
376  case HEX27:
377  return 27;
378 
379  case PRISM15:
380  return 15;
381 
382  case PRISM18:
383  return 18;
384 
385  case INVALID_ELEM:
386  return 0;
387 
388  default:
389  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
390  }
391  }
392 
393  case THIRD:
394  {
395  switch (t)
396  {
397  case NODEELEM:
398  return 1;
399 
400  case EDGE4:
401  return 4;
402 
403  case INVALID_ELEM:
404  return 0;
405 
406  default:
407  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
408  }
409  }
410 
411  default:
412  libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for L2_LAGRANGE FE family!");
413  }
414 
415  libmesh_error_msg("We'll never get here!");
416  return 0;
417 }
418 
419 } // anonymous namespace
420 
421 
422  // Do full-specialization for every dimension, instead
423  // of explicit instantiation at the end of this file.
424  // This could be macro-ified so that it fits on one line...
425 template <>
427  const Order order,
428  const std::vector<Number> & elem_soln,
429  std::vector<Number> & nodal_soln)
430 { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
431 
432 template <>
434  const Order order,
435  const std::vector<Number> & elem_soln,
436  std::vector<Number> & nodal_soln)
437 { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
438 
439 template <>
441  const Order order,
442  const std::vector<Number> & elem_soln,
443  std::vector<Number> & nodal_soln)
444 { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
445 
446 template <>
448  const Order order,
449  const std::vector<Number> & elem_soln,
450  std::vector<Number> & nodal_soln)
451 { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
452 
453 
454 // Do full-specialization for every dimension, instead
455 // of explicit instantiation at the end of this function.
456 // This could be macro-ified.
457 template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
458 template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
459 template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
460 template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
461 
462 
463 // L2 Lagrange elements have all dofs on elements (hence none on nodes)
464 template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
465 template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
466 template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
467 template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
468 
469 
470 // L2 Lagrange elements have all dofs on elements
471 template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
472 template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
473 template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
474 template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
475 
476 // L2 Lagrange FEMs are DISCONTINUOUS
481 
482 // Lagrange FEMs are not hierarchic
483 template <> bool FE<0,L2_LAGRANGE>::is_hierarchic() const { return false; }
484 template <> bool FE<1,L2_LAGRANGE>::is_hierarchic() const { return false; }
485 template <> bool FE<2,L2_LAGRANGE>::is_hierarchic() const { return false; }
486 template <> bool FE<3,L2_LAGRANGE>::is_hierarchic() const { return false; }
487 
488 // Lagrange FEM shapes do not need reinit (is this always true?)
489 template <> bool FE<0,L2_LAGRANGE>::shapes_need_reinit() const { return false; }
490 template <> bool FE<1,L2_LAGRANGE>::shapes_need_reinit() const { return false; }
491 template <> bool FE<2,L2_LAGRANGE>::shapes_need_reinit() const { return false; }
492 template <> bool FE<3,L2_LAGRANGE>::shapes_need_reinit() const { return false; }
493 
494 // We don't need any constraints for this DISCONTINUOUS basis, hence these
495 // are NOOPs
496 #ifdef LIBMESH_ENABLE_AMR
497 template <>
499  DofMap &,
500  const unsigned int,
501  const Elem *)
502 { }
503 
504 template <>
506  DofMap &,
507  const unsigned int,
508  const Elem *)
509 { }
510 #endif // LIBMESH_ENABLE_AMR
511 
512 } // namespace libMesh
static unsigned int n_dofs(const ElemType t, const Order o)
virtual bool shapes_need_reinit() const libmesh_override
ElemType
Defines an enum for geometric element types.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
static unsigned int n_dofs_at_node(const ElemType t, const Order o, const unsigned int n)
The libMesh namespace provides an interface to certain functionality in the library.
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
const dof_id_type n_nodes
Definition: tecplot_io.C:67
virtual FEContinuity get_continuity() const libmesh_override
virtual bool is_hierarchic() const libmesh_override
static unsigned int n_dofs_per_elem(const ElemType t, const Order o)
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
std::string enum_to_string(const T e)
static void compute_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
Computes the constraint matrix contributions (for non-conforming adapted meshes) corresponding to var...
static void nodal_soln(const Elem *elem, const Order o, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)
Build the nodal soln from the element soln.
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
The constraint matrix storage format.
Definition: dof_map.h:96