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