libMesh
fe_monomial.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/fe.h"
22 #include "libmesh/elem.h"
23 #include "libmesh/fe_interface.h"
24 #include "libmesh/enum_to_string.h"
25 
26 namespace libMesh
27 {
28 
29 // Anonymous namespace for local helper functions
30 namespace {
31 
32 void monomial_nodal_soln(const Elem * elem,
33  const Order order,
34  const std::vector<Number> & elem_soln,
35  std::vector<Number> & nodal_soln,
36  const unsigned Dim)
37 {
38  const unsigned int n_nodes = elem->n_nodes();
39 
40  const ElemType elem_type = elem->type();
41 
42  nodal_soln.resize(n_nodes);
43 
44  const Order totalorder = static_cast<Order>(order+elem->p_level());
45 
46  switch (totalorder)
47  {
48  // Constant shape functions
49  case CONSTANT:
50  {
51  libmesh_assert_equal_to (elem_soln.size(), 1);
52 
53  const Number val = elem_soln[0];
54 
55  for (unsigned int n=0; n<n_nodes; n++)
56  nodal_soln[n] = val;
57 
58  return;
59  }
60 
61 
62  // For other orders, do interpolation at the nodes
63  // explicitly.
64  default:
65  {
66  // FEType object to be passed to various FEInterface functions below.
67  FEType fe_type(totalorder, MONOMIAL);
68 
69  const unsigned int n_sf =
70  // FE<Dim,T>::n_shape_functions(elem_type, totalorder);
71  FEInterface::n_shape_functions(Dim, fe_type, elem_type);
72 
73  std::vector<Point> refspace_nodes;
74  FEBase::get_refspace_nodes(elem_type,refspace_nodes);
75  libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
76 
77  for (unsigned int n=0; n<n_nodes; n++)
78  {
79  libmesh_assert_equal_to (elem_soln.size(), n_sf);
80 
81  // Zero before summation
82  nodal_soln[n] = 0;
83 
84  // u_i = Sum (alpha_i phi_i)
85  for (unsigned int i=0; i<n_sf; i++)
86  nodal_soln[n] += elem_soln[i] *
87  // FE<Dim,T>::shape(elem, order, i, mapped_point);
88  FEInterface::shape(Dim, fe_type, elem, i, refspace_nodes[n]);
89  }
90 
91  return;
92  } // default
93  } // switch
94 } // monomial_nodal_soln()
95 
96 
97 
98 
99 unsigned int monomial_n_dofs(const ElemType t, const Order o)
100 {
101  switch (o)
102  {
103 
104  // constant shape functions
105  // no matter what shape there is only one DOF.
106  case CONSTANT:
107  return (t != INVALID_ELEM) ? 1 : 0;
108 
109 
110  // Discontinuous linear shape functions
111  // expressed in the monomials.
112  case FIRST:
113  {
114  switch (t)
115  {
116  case NODEELEM:
117  return 1;
118 
119  case EDGE2:
120  case EDGE3:
121  case EDGE4:
122  return 2;
123 
124  case TRI3:
125  case TRISHELL3:
126  case TRI6:
127  case QUAD4:
128  case QUADSHELL4:
129  case QUAD8:
130  case QUADSHELL8:
131  case QUAD9:
132  return 3;
133 
134  case TET4:
135  case TET10:
136  case HEX8:
137  case HEX20:
138  case HEX27:
139  case PRISM6:
140  case PRISM15:
141  case PRISM18:
142  case PYRAMID5:
143  case PYRAMID13:
144  case PYRAMID14:
145  return 4;
146 
147  case INVALID_ELEM:
148  return 0;
149 
150  default:
151  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
152  }
153  }
154 
155 
156  // Discontinuous quadratic shape functions
157  // expressed in the monomials.
158  case SECOND:
159  {
160  switch (t)
161  {
162  case NODEELEM:
163  return 1;
164 
165  case EDGE2:
166  case EDGE3:
167  case EDGE4:
168  return 3;
169 
170  case TRI3:
171  case TRISHELL3:
172  case TRI6:
173  case QUAD4:
174  case QUADSHELL4:
175  case QUAD8:
176  case QUADSHELL8:
177  case QUAD9:
178  return 6;
179 
180  case TET4:
181  case TET10:
182  case HEX8:
183  case HEX20:
184  case HEX27:
185  case PRISM6:
186  case PRISM15:
187  case PRISM18:
188  case PYRAMID5:
189  case PYRAMID13:
190  case PYRAMID14:
191  return 10;
192 
193  case INVALID_ELEM:
194  return 0;
195 
196  default:
197  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
198  }
199  }
200 
201 
202  // Discontinuous cubic shape functions
203  // expressed in the monomials.
204  case THIRD:
205  {
206  switch (t)
207  {
208  case NODEELEM:
209  return 1;
210 
211  case EDGE2:
212  case EDGE3:
213  case EDGE4:
214  return 4;
215 
216  case TRI3:
217  case TRISHELL3:
218  case TRI6:
219  case QUAD4:
220  case QUADSHELL4:
221  case QUAD8:
222  case QUADSHELL8:
223  case QUAD9:
224  return 10;
225 
226  case TET4:
227  case TET10:
228  case HEX8:
229  case HEX20:
230  case HEX27:
231  case PRISM6:
232  case PRISM15:
233  case PRISM18:
234  case PYRAMID5:
235  case PYRAMID13:
236  case PYRAMID14:
237  return 20;
238 
239  case INVALID_ELEM:
240  return 0;
241 
242  default:
243  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
244  }
245  }
246 
247 
248  // Discontinuous quartic shape functions
249  // expressed in the monomials.
250  case FOURTH:
251  {
252  switch (t)
253  {
254  case NODEELEM:
255  return 1;
256 
257  case EDGE2:
258  case EDGE3:
259  return 5;
260 
261  case TRI3:
262  case TRISHELL3:
263  case TRI6:
264  case QUAD4:
265  case QUADSHELL4:
266  case QUAD8:
267  case QUADSHELL8:
268  case QUAD9:
269  return 15;
270 
271  case TET4:
272  case TET10:
273  case HEX8:
274  case HEX20:
275  case HEX27:
276  case PRISM6:
277  case PRISM15:
278  case PRISM18:
279  case PYRAMID5:
280  case PYRAMID13:
281  case PYRAMID14:
282  return 35;
283 
284  case INVALID_ELEM:
285  return 0;
286 
287  default:
288  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
289  }
290  }
291 
292 
293  default:
294  {
295  const unsigned int order = static_cast<unsigned int>(o);
296  switch (t)
297  {
298  case NODEELEM:
299  return 1;
300 
301  case EDGE2:
302  case EDGE3:
303  return (order+1);
304 
305  case TRI3:
306  case TRISHELL3:
307  case TRI6:
308  case QUAD4:
309  case QUADSHELL4:
310  case QUAD8:
311  case QUADSHELL8:
312  case QUAD9:
313  return (order+1)*(order+2)/2;
314 
315  case TET4:
316  case TET10:
317  case HEX8:
318  case HEX20:
319  case HEX27:
320  case PRISM6:
321  case PRISM15:
322  case PRISM18:
323  case PYRAMID5:
324  case PYRAMID13:
325  case PYRAMID14:
326  return (order+1)*(order+2)*(order+3)/6;
327 
328  case INVALID_ELEM:
329  return 0;
330 
331  default:
332  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
333  }
334  }
335  }
336 } // monomial_n_dofs()
337 
338 
339 } // anonymous namespace
340 
341 
342 
343 
344 
345  // Do full-specialization for every dimension, instead
346  // of explicit instantiation at the end of this file.
347  // This could be macro-ified so that it fits on one line...
348 template <>
349 void FE<0,MONOMIAL>::nodal_soln(const Elem * elem,
350  const Order order,
351  const std::vector<Number> & elem_soln,
352  std::vector<Number> & nodal_soln)
353 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/0); }
354 
355 template <>
356 void FE<1,MONOMIAL>::nodal_soln(const Elem * elem,
357  const Order order,
358  const std::vector<Number> & elem_soln,
359  std::vector<Number> & nodal_soln)
360 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/1); }
361 
362 template <>
363 void FE<2,MONOMIAL>::nodal_soln(const Elem * elem,
364  const Order order,
365  const std::vector<Number> & elem_soln,
366  std::vector<Number> & nodal_soln)
367 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/2); }
368 
369 template <>
370 void FE<3,MONOMIAL>::nodal_soln(const Elem * elem,
371  const Order order,
372  const std::vector<Number> & elem_soln,
373  std::vector<Number> & nodal_soln)
374 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/3); }
375 
376 
377 // Full specialization of n_dofs() function for every dimension
378 template <> unsigned int FE<0,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
379 template <> unsigned int FE<1,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
380 template <> unsigned int FE<2,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
381 template <> unsigned int FE<3,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
382 
383 // Full specialization of n_dofs_at_node() function for every dimension.
384 // Monomials have no dofs at nodes, only element dofs.
385 template <> unsigned int FE<0,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
386 template <> unsigned int FE<1,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
387 template <> unsigned int FE<2,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
388 template <> unsigned int FE<3,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
389 
390 // Full specialization of n_dofs_per_elem() function for every dimension.
391 template <> unsigned int FE<0,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
392 template <> unsigned int FE<1,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
393 template <> unsigned int FE<2,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
394 template <> unsigned int FE<3,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
395 
396 
397 // Full specialization of get_continuity() function for every dimension.
398 template <> FEContinuity FE<0,MONOMIAL>::get_continuity() const { return DISCONTINUOUS; }
399 template <> FEContinuity FE<1,MONOMIAL>::get_continuity() const { return DISCONTINUOUS; }
400 template <> FEContinuity FE<2,MONOMIAL>::get_continuity() const { return DISCONTINUOUS; }
401 template <> FEContinuity FE<3,MONOMIAL>::get_continuity() const { return DISCONTINUOUS; }
402 
403 // Full specialization of is_hierarchic() function for every dimension.
404 // The monomials are hierarchic!
405 template <> bool FE<0,MONOMIAL>::is_hierarchic() const { return true; }
406 template <> bool FE<1,MONOMIAL>::is_hierarchic() const { return true; }
407 template <> bool FE<2,MONOMIAL>::is_hierarchic() const { return true; }
408 template <> bool FE<3,MONOMIAL>::is_hierarchic() const { return true; }
409 
410 #ifdef LIBMESH_ENABLE_AMR
411 
412 // Full specialization of compute_constraints() function for 2D and
413 // 3D only. There are no constraints for discontinuous elements, so
414 // we do nothing.
415 template <> void FE<2,MONOMIAL>::compute_constraints (DofConstraints &, DofMap &, const unsigned int, const Elem *) {}
416 template <> void FE<3,MONOMIAL>::compute_constraints (DofConstraints &, DofMap &, const unsigned int, const Elem *) {}
417 
418 #endif // #ifdef LIBMESH_ENABLE_AMR
419 
420 // Full specialization of shapes_need_reinit() function for every dimension.
421 template <> bool FE<0,MONOMIAL>::shapes_need_reinit() const { return false; }
422 template <> bool FE<1,MONOMIAL>::shapes_need_reinit() const { return false; }
423 template <> bool FE<2,MONOMIAL>::shapes_need_reinit() const { return false; }
424 template <> bool FE<3,MONOMIAL>::shapes_need_reinit() const { return false; }
425 
426 } // namespace libMesh
libMesh::DofConstraints
The constraint matrix storage format.
Definition: dof_map.h:105
libMesh::HEX20
Definition: enum_elem_type.h:48
libMesh::Number
Real Number
Definition: libmesh_common.h:195
libMesh::PRISM6
Definition: enum_elem_type.h:50
libMesh::FEInterface::n_shape_functions
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:446
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...
libMesh::CONSTANT
Definition: enum_order.h:41
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::MONOMIAL
Definition: enum_fe_family.h:39
libMesh::FOURTH
Definition: enum_order.h:45
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::FEAbstract::get_refspace_nodes
static void get_refspace_nodes(const ElemType t, std::vector< Point > &nodes)
Definition: fe_abstract.C:308
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
libMesh::FEInterface::shape
static Real shape(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const Point &p)
Definition: fe_interface.C:687
Generated on Sat Jan 25 2020 12:06:58 for libMesh by   doxygen 1.8.16