fe_l2_lagrange.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/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 
416 } // anonymous namespace
417 
418 
419  // Do full-specialization for every dimension, instead
420  // of explicit instantiation at the end of this file.
421  // This could be macro-ified so that it fits on one line...
422 template <>
423 void FE<0,L2_LAGRANGE>::nodal_soln(const Elem * elem,
424  const Order order,
425  const std::vector<Number> & elem_soln,
426  std::vector<Number> & nodal_soln)
427 { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
428 
429 template <>
430 void FE<1,L2_LAGRANGE>::nodal_soln(const Elem * elem,
431  const Order order,
432  const std::vector<Number> & elem_soln,
433  std::vector<Number> & nodal_soln)
434 { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
435 
436 template <>
437 void FE<2,L2_LAGRANGE>::nodal_soln(const Elem * elem,
438  const Order order,
439  const std::vector<Number> & elem_soln,
440  std::vector<Number> & nodal_soln)
441 { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
442 
443 template <>
444 void FE<3,L2_LAGRANGE>::nodal_soln(const Elem * elem,
445  const Order order,
446  const std::vector<Number> & elem_soln,
447  std::vector<Number> & nodal_soln)
448 { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
449 
450 
451 // Do full-specialization for every dimension, instead
452 // of explicit instantiation at the end of this function.
453 // This could be macro-ified.
454 template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
455 template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
456 template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
457 template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
458 
459 
460 // L2 Lagrange elements have all dofs on elements (hence none on nodes)
461 template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
462 template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
463 template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
464 template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
465 
466 
467 // L2 Lagrange elements have all dofs on elements
468 template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
469 template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
470 template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
471 template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
472 
473 // L2 Lagrange FEMs are DISCONTINUOUS
474 template <> FEContinuity FE<0,L2_LAGRANGE>::get_continuity() const { return DISCONTINUOUS; }
475 template <> FEContinuity FE<1,L2_LAGRANGE>::get_continuity() const { return DISCONTINUOUS; }
476 template <> FEContinuity FE<2,L2_LAGRANGE>::get_continuity() const { return DISCONTINUOUS; }
477 template <> FEContinuity FE<3,L2_LAGRANGE>::get_continuity() const { return DISCONTINUOUS; }
478 
479 // Lagrange FEMs are not hierarchic
480 template <> bool FE<0,L2_LAGRANGE>::is_hierarchic() const { return false; }
481 template <> bool FE<1,L2_LAGRANGE>::is_hierarchic() const { return false; }
482 template <> bool FE<2,L2_LAGRANGE>::is_hierarchic() const { return false; }
483 template <> bool FE<3,L2_LAGRANGE>::is_hierarchic() const { return false; }
484 
485 // Lagrange FEM shapes do not need reinit (is this always true?)
486 template <> bool FE<0,L2_LAGRANGE>::shapes_need_reinit() const { return false; }
487 template <> bool FE<1,L2_LAGRANGE>::shapes_need_reinit() const { return false; }
488 template <> bool FE<2,L2_LAGRANGE>::shapes_need_reinit() const { return false; }
489 template <> bool FE<3,L2_LAGRANGE>::shapes_need_reinit() const { return false; }
490 
491 // We don't need any constraints for this DISCONTINUOUS basis, hence these
492 // are NOOPs
493 #ifdef LIBMESH_ENABLE_AMR
494 template <>
495 void FE<2,L2_LAGRANGE>::compute_constraints (DofConstraints &,
496  DofMap &,
497  const unsigned int,
498  const Elem *)
499 { }
500 
501 template <>
502 void FE<3,L2_LAGRANGE>::compute_constraints (DofConstraints &,
503  DofMap &,
504  const unsigned int,
505  const Elem *)
506 { }
507 #endif // LIBMESH_ENABLE_AMR
508 
509 } // namespace libMesh
libMesh::FE::n_dofs
static unsigned int n_dofs(const ElemType t, const Order o)
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)
fe_interface.h
libMesh::FE::shapes_need_reinit
virtual bool shapes_need_reinit() const override
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::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::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh::FE::compute_constraints
static void compute_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
dof_map.h
string_to_enum.h
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)
generated by doxygen