Line data Source code
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 17038542 : unsigned int monomial_n_dofs(const Elem * e, const Order o)
31 : {
32 2083160 : libmesh_assert(e);
33 17038542 : return monomial_n_dofs(e->type(), o);
34 : }
35 :
36 :
37 17038542 : unsigned int monomial_n_dofs(const ElemType t, const Order o)
38 : {
39 17038542 : switch (o)
40 : {
41 :
42 : // constant shape functions
43 : // no matter what shape there is only one DOF.
44 8305483 : case CONSTANT:
45 8305483 : return (t != INVALID_ELEM) ? 1 : 0;
46 :
47 :
48 : // Discontinuous linear shape functions
49 : // expressed in the monomials.
50 5377572 : case FIRST:
51 : {
52 5377572 : switch (t)
53 : {
54 0 : case NODEELEM:
55 0 : return 1;
56 :
57 623 : case EDGE2:
58 : case EDGE3:
59 : case EDGE4:
60 623 : return 2;
61 :
62 168292 : 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 168292 : return 3;
74 :
75 861119 : 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 861119 : return 4;
92 :
93 0 : case INVALID_ELEM:
94 0 : return 0;
95 :
96 0 : default:
97 0 : 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 1461190 : case SECOND:
105 : {
106 1461190 : switch (t)
107 : {
108 0 : case NODEELEM:
109 0 : return 1;
110 :
111 588 : case EDGE2:
112 : case EDGE3:
113 : case EDGE4:
114 588 : return 3;
115 :
116 52247 : 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 52247 : return 6;
128 :
129 218942 : 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 218942 : return 10;
146 :
147 0 : case INVALID_ELEM:
148 0 : return 0;
149 :
150 0 : default:
151 0 : 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 1101946 : case THIRD:
159 : {
160 1101946 : switch (t)
161 : {
162 0 : case NODEELEM:
163 0 : return 1;
164 :
165 588 : case EDGE2:
166 : case EDGE3:
167 : case EDGE4:
168 588 : return 4;
169 :
170 9144 : 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 9144 : return 10;
182 :
183 180347 : 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 180347 : return 20;
200 :
201 0 : case INVALID_ELEM:
202 0 : return 0;
203 :
204 0 : default:
205 0 : 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 528906 : case FOURTH:
213 : {
214 528906 : switch (t)
215 : {
216 0 : case NODEELEM:
217 0 : return 1;
218 :
219 588 : case EDGE2:
220 : case EDGE3:
221 588 : return 5;
222 :
223 49880 : 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 49880 : return 15;
235 :
236 476056 : 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 476056 : return 35;
252 :
253 0 : case INVALID_ELEM:
254 0 : return 0;
255 :
256 0 : default:
257 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
258 : }
259 : }
260 :
261 :
262 263445 : default:
263 : {
264 263445 : const unsigned int order = static_cast<unsigned int>(o);
265 263445 : switch (t)
266 : {
267 0 : case NODEELEM:
268 0 : return 1;
269 :
270 1485 : case EDGE2:
271 : case EDGE3:
272 1485 : return (order+1);
273 :
274 24922 : 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 24922 : return (order+1)*(order+2)/2;
286 :
287 237038 : 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 237038 : return (order+1)*(order+2)*(order+3)/6;
303 :
304 0 : case INVALID_ELEM:
305 0 : return 0;
306 :
307 0 : default:
308 0 : 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 1210595 : 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 1210595 : const unsigned int n_nodes = elem->n_nodes();
325 :
326 1210595 : const ElemType elem_type = elem->type();
327 :
328 1210595 : nodal_soln.resize(n_nodes);
329 :
330 1321910 : const Order totalorder = order + add_p_level*elem->p_level();
331 :
332 1210595 : switch (totalorder)
333 : {
334 : // Constant shape functions
335 211474 : case CONSTANT:
336 : {
337 105737 : libmesh_assert_equal_to (elem_soln.size(), 1);
338 :
339 317211 : std::fill(nodal_soln.begin(), nodal_soln.end(), elem_soln[0]);
340 :
341 105737 : return;
342 : }
343 :
344 :
345 : // For other orders, do interpolation at the nodes
346 : // explicitly.
347 5578 : default:
348 : {
349 : // FEType object to be passed to various FEInterface functions below.
350 5578 : FEType fe_type(order, MONOMIAL);
351 :
352 : const unsigned int n_sf =
353 52769 : FEInterface::n_shape_functions(fe_type, elem);
354 :
355 11156 : std::vector<Point> refspace_nodes;
356 52769 : FEBase::get_refspace_nodes(elem_type,refspace_nodes);
357 5578 : libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
358 5578 : libmesh_assert_equal_to (elem_soln.size(), n_sf);
359 :
360 : // Zero before summation
361 5578 : std::fill(nodal_soln.begin(), nodal_soln.end(), 0);
362 :
363 586605 : for (unsigned int n=0; n<n_nodes; n++)
364 : // u_i = Sum (alpha_i phi_i)
365 2553080 : for (unsigned int i=0; i<n_sf; i++)
366 2225293 : nodal_soln[n] += elem_soln[i] *
367 2225293 : FEInterface::shape(fe_type, elem, i, refspace_nodes[n]);
368 :
369 5578 : 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 1210595 : LIBMESH_FE_NODAL_SOLN(MONOMIAL, monomial_nodal_soln)
382 0 : LIBMESH_FE_SIDE_NODAL_SOLN(MONOMIAL)
383 :
384 :
385 : // Full specialization of n_dofs() function for every dimension
386 0 : template <> unsigned int FE<0,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
387 0 : template <> unsigned int FE<1,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
388 0 : template <> unsigned int FE<2,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
389 0 : template <> unsigned int FE<3,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
390 :
391 2 : template <> unsigned int FE<0,MONOMIAL>::n_dofs(const Elem * e, const Order o) { return monomial_n_dofs(e, o); }
392 7532 : template <> unsigned int FE<1,MONOMIAL>::n_dofs(const Elem * e, const Order o) { return monomial_n_dofs(e, o); }
393 524424 : template <> unsigned int FE<2,MONOMIAL>::n_dofs(const Elem * e, const Order o) { return monomial_n_dofs(e, o); }
394 2653659 : 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 31 : template <> unsigned int FE<0,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
399 70581 : template <> unsigned int FE<1,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
400 9374058 : template <> unsigned int FE<2,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
401 84795597 : template <> unsigned int FE<3,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
402 :
403 296 : template <> unsigned int FE<0,MONOMIAL>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
404 81342 : template <> unsigned int FE<1,MONOMIAL>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
405 2962498 : template <> unsigned int FE<2,MONOMIAL>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
406 19704043 : 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 0 : template <> unsigned int FE<0,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
410 0 : template <> unsigned int FE<1,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
411 0 : template <> unsigned int FE<2,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
412 0 : template <> unsigned int FE<3,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
413 :
414 197 : template <> unsigned int FE<0,MONOMIAL>::n_dofs_per_elem(const Elem & e, const Order o) { return monomial_n_dofs(&e, o); }
415 53558 : template <> unsigned int FE<1,MONOMIAL>::n_dofs_per_elem(const Elem & e, const Order o) { return monomial_n_dofs(&e, o); }
416 1768353 : template <> unsigned int FE<2,MONOMIAL>::n_dofs_per_elem(const Elem & e, const Order o) { return monomial_n_dofs(&e, o); }
417 7624503 : 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 361 : template <> FEContinuity FE<0,MONOMIAL>::get_continuity() const { return DISCONTINUOUS; }
422 20819 : template <> FEContinuity FE<1,MONOMIAL>::get_continuity() const { return DISCONTINUOUS; }
423 90741 : template <> FEContinuity FE<2,MONOMIAL>::get_continuity() const { return DISCONTINUOUS; }
424 153928 : 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 0 : template <> bool FE<0,MONOMIAL>::is_hierarchic() const { return true; }
429 60 : template <> bool FE<1,MONOMIAL>::is_hierarchic() const { return true; }
430 318 : template <> bool FE<2,MONOMIAL>::is_hierarchic() const { return true; }
431 790 : 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 0 : template <> void FE<2,MONOMIAL>::compute_constraints (DofConstraints &, DofMap &, const unsigned int, const Elem *) {}
439 0 : 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 0 : template <> bool FE<0,MONOMIAL>::shapes_need_reinit() const { return false; }
445 1314 : template <> bool FE<1,MONOMIAL>::shapes_need_reinit() const { return false; }
446 761897 : template <> bool FE<2,MONOMIAL>::shapes_need_reinit() const { return false; }
447 2741309 : template <> bool FE<3,MONOMIAL>::shapes_need_reinit() const { return false; }
448 :
449 : } // namespace libMesh
|
