libMesh
Loading...
Searching...
No Matches
fe_monomial.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 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
27namespace libMesh
28{
29
30unsigned int monomial_n_dofs(const Elem * e, const Order o)
31{
33 return monomial_n_dofs(e->type(), o);
34}
35
36
37unsigned 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
316namespace {
317
318void 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
381LIBMESH_FE_NODAL_SOLN(MONOMIAL, monomial_nodal_soln)
383
384
385// Full specialization of n_dofs() function for every dimension
386template <> unsigned int FE<0,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
387template <> unsigned int FE<1,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
388template <> unsigned int FE<2,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
389template <> unsigned int FE<3,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
390
391template <> unsigned int FE<0,MONOMIAL>::n_dofs(const Elem * e, const Order o) { return monomial_n_dofs(e, o); }
392template <> unsigned int FE<1,MONOMIAL>::n_dofs(const Elem * e, const Order o) { return monomial_n_dofs(e, o); }
393template <> unsigned int FE<2,MONOMIAL>::n_dofs(const Elem * e, const Order o) { return monomial_n_dofs(e, o); }
394template <> 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.
398template <> unsigned int FE<0,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
399template <> unsigned int FE<1,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
400template <> unsigned int FE<2,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
401template <> unsigned int FE<3,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
402
403template <> unsigned int FE<0,MONOMIAL>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
404template <> unsigned int FE<1,MONOMIAL>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
405template <> unsigned int FE<2,MONOMIAL>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
406template <> 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.
409template <> unsigned int FE<0,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
410template <> unsigned int FE<1,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
411template <> unsigned int FE<2,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
412template <> unsigned int FE<3,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
413
414template <> unsigned int FE<0,MONOMIAL>::n_dofs_per_elem(const Elem & e, const Order o) { return monomial_n_dofs(&e, o); }
415template <> unsigned int FE<1,MONOMIAL>::n_dofs_per_elem(const Elem & e, const Order o) { return monomial_n_dofs(&e, o); }
416template <> unsigned int FE<2,MONOMIAL>::n_dofs_per_elem(const Elem & e, const Order o) { return monomial_n_dofs(&e, o); }
417template <> 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.
425
426// Full specialization of is_hierarchic() function for every dimension.
427// The monomials are hierarchic!
428template <> bool FE<0,MONOMIAL>::is_hierarchic() const { return true; }
429template <> bool FE<1,MONOMIAL>::is_hierarchic() const { return true; }
430template <> bool FE<2,MONOMIAL>::is_hierarchic() const { return true; }
431template <> 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.
438template <> void FE<2,MONOMIAL>::compute_constraints (DofConstraints &, DofMap &, const unsigned int, const Elem *) {}
439template <> 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.
444template <> bool FE<0,MONOMIAL>::shapes_need_reinit() const { return false; }
445template <> bool FE<1,MONOMIAL>::shapes_need_reinit() const { return false; }
446template <> bool FE<2,MONOMIAL>::shapes_need_reinit() const { return false; }
447template <> bool FE<3,MONOMIAL>::shapes_need_reinit() const { return false; }
448
449} // namespace libMesh
The constraint matrix storage format.
Definition dof_map.h:112
This class handles the numbering of degrees of freedom on a mesh.
Definition dof_map.h:181
This is the base class from which all geometric element types are derived.
Definition elem.h:96
virtual ElemType type() const =0
static void get_refspace_nodes(const ElemType t, std::vector< Point > &nodes)
static Real shape(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const Point &p)
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
virtual bool is_hierarchic() const override
static unsigned int n_dofs_per_elem(const ElemType t, const Order o)
virtual FEContinuity get_continuity() const override
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 unsigned int n_dofs(const ElemType t, const Order o)
static unsigned int n_dofs_at_node(const ElemType t, const Order o, const unsigned int n)
virtual bool shapes_need_reinit() const override
std::string enum_to_string(const T e)
The libMesh namespace provides an interface to certain functionality in the library.
LIBMESH_FE_NODAL_SOLN(BERNSTEIN, bernstein_nodal_soln)
ElemType
Defines an enum for geometric element types.
libmesh_assert(ctx)
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_SIDE_NODAL_SOLN(HIERARCHIC_VEC)
const dof_id_type n_nodes
Definition tecplot_io.C:67