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