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/elem.h"
22 : #include "libmesh/enum_to_string.h"
23 : #include "libmesh/fe.h"
24 : #include "libmesh/fe_interface.h"
25 : #include "libmesh/fe_macro.h"
26 :
27 : namespace libMesh
28 : {
29 :
30 : // ------------------------------------------------------------
31 : // Hierarchic-specific implementations
32 :
33 : // Anonymous namespace for local helper functions
34 : namespace {
35 :
36 4824 : void hierarchic_nodal_soln(const Elem * elem,
37 : const Order order,
38 : const std::vector<Number> & elem_soln,
39 : std::vector<Number> & nodal_soln,
40 : const bool add_p_level)
41 : {
42 4824 : const unsigned int n_nodes = elem->n_nodes();
43 :
44 4824 : const ElemType elem_type = elem->type();
45 :
46 4824 : nodal_soln.resize(n_nodes);
47 :
48 5226 : const Order totalorder = order + add_p_level*elem->p_level();
49 :
50 : // FEType object to be passed to various FEInterface functions below.
51 402 : FEType fe_type(order, HIERARCHIC);
52 :
53 4824 : switch (totalorder)
54 : {
55 : // Constant shape functions
56 0 : case CONSTANT:
57 : {
58 0 : libmesh_assert_equal_to (elem_soln.size(), 1);
59 :
60 0 : std::fill(nodal_soln.begin(), nodal_soln.end(), elem_soln[0]);
61 :
62 0 : return;
63 : }
64 :
65 :
66 : // For other orders do interpolation at the nodes
67 : // explicitly.
68 4824 : default:
69 : {
70 : const unsigned int n_sf =
71 4824 : FEInterface::n_shape_functions(fe_type, elem);
72 :
73 804 : std::vector<Point> refspace_nodes;
74 4824 : FEBase::get_refspace_nodes(elem_type,refspace_nodes);
75 402 : libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
76 402 : libmesh_assert_equal_to (elem_soln.size(), n_sf);
77 :
78 : // Zero before summation
79 402 : std::fill(nodal_soln.begin(), nodal_soln.end(), 0);
80 :
81 61920 : for (unsigned int n=0; n<n_nodes; n++)
82 : // u_i = Sum (alpha_i phi_i)
83 1130328 : for (unsigned int i=0; i<n_sf; i++)
84 1162668 : nodal_soln[n] += elem_soln[i] *
85 1162668 : FEInterface::shape(fe_type, elem, i, refspace_nodes[n]);
86 :
87 402 : return;
88 : }
89 : }
90 : } // hierarchic_nodal_soln()
91 :
92 :
93 :
94 1571689 : unsigned int HIERARCHIC_n_dofs(const ElemType t, const Order o)
95 : {
96 193520 : libmesh_assert_greater (o, 0);
97 1571689 : switch (t)
98 : {
99 0 : case NODEELEM:
100 0 : return 1;
101 4814 : case EDGE2:
102 : case EDGE3:
103 25969 : return (o+1);
104 507 : case QUAD4:
105 : case QUADSHELL4:
106 507 : libmesh_assert_less (o, 2);
107 : libmesh_fallthrough();
108 : case QUAD8:
109 : case QUADSHELL8:
110 : case QUAD9:
111 : case QUADSHELL9:
112 312101 : return ((o+1)*(o+1));
113 3439 : case HEX8:
114 : case HEX20:
115 3439 : libmesh_assert_less (o, 2);
116 : libmesh_fallthrough();
117 : case HEX27:
118 633807 : return ((o+1)*(o+1)*(o+1));
119 7180 : case PRISM6:
120 : case PRISM15:
121 7180 : libmesh_assert_less (o, 2);
122 : libmesh_fallthrough();
123 : case PRISM18:
124 11463 : libmesh_assert_less (o, 3);
125 : libmesh_fallthrough();
126 : case PRISM20:
127 : case PRISM21:
128 390879 : return ((o+1)*(o+1)*(o+2)/2);
129 609 : case TRI3:
130 609 : libmesh_assert_less (o, 2);
131 : libmesh_fallthrough();
132 : case TRI6:
133 : case TRI7:
134 123696 : return ((o+1)*(o+2)/2);
135 1257 : case TET4:
136 1257 : libmesh_assert_less (o, 2);
137 : libmesh_fallthrough();
138 : case TET10:
139 4281 : libmesh_assert_less (o, 3);
140 : libmesh_fallthrough();
141 : case TET14:
142 122986 : return ((o+1)*(o+2)*(o+3)/6);
143 0 : case INVALID_ELEM:
144 0 : return 0;
145 0 : default:
146 0 : libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for HIERARCHIC FE family!");
147 : }
148 : } // HIERARCHIC_n_dofs()
149 :
150 :
151 :
152 1451880 : unsigned int HIERARCHIC_n_dofs(const Elem * e, const Order o)
153 : {
154 193520 : libmesh_assert(e);
155 1571689 : return HIERARCHIC_n_dofs(e->type(), o);
156 : }
157 :
158 :
159 :
160 24552910 : unsigned int HIERARCHIC_n_dofs_at_node(const ElemType t,
161 : const Order o,
162 : const unsigned int n)
163 : {
164 1718859 : libmesh_assert_greater (o, 0);
165 24552910 : switch (t)
166 : {
167 0 : case NODEELEM:
168 0 : return 1;
169 160815 : case EDGE2:
170 : case EDGE3:
171 160815 : switch (n)
172 : {
173 9000 : case 0:
174 : case 1:
175 9000 : return 1;
176 : // Internal DoFs are associated with the elem, not its nodes
177 8329 : case 2:
178 4166 : libmesh_assert_equal_to(t, EDGE3);
179 8329 : return 0;
180 0 : default:
181 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for EDGE2/3!");
182 : }
183 816336 : case TRI3:
184 792 : libmesh_assert_less (n, 3);
185 792 : libmesh_assert_less (o, 2);
186 : libmesh_fallthrough();
187 : case TRI6:
188 : // Internal DoFs are associated with the elem on a Tri6, or node 6 on a Tri7
189 42128 : libmesh_assert_less (n, 6);
190 : libmesh_fallthrough();
191 : case TRI7:
192 879253 : switch (n)
193 : {
194 31798 : case 0:
195 : case 1:
196 : case 2:
197 31798 : return 1;
198 :
199 28819 : case 3:
200 : case 4:
201 : case 5:
202 395972 : return (o-1);
203 :
204 3092 : case 6:
205 39263 : return ((o-1)*(o-2)/2);
206 0 : default:
207 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI!");
208 : }
209 3914603 : case QUAD4:
210 : case QUADSHELL4:
211 513 : libmesh_assert_less (n, 4);
212 513 : libmesh_assert_less (o, 2);
213 : libmesh_fallthrough();
214 : case QUAD8:
215 : case QUADSHELL8:
216 : case QUAD9:
217 : case QUADSHELL9:
218 4245922 : switch (n)
219 : {
220 153896 : case 0:
221 : case 1:
222 : case 2:
223 : case 3:
224 153896 : return 1;
225 :
226 142428 : case 4:
227 : case 5:
228 : case 6:
229 : case 7:
230 1819146 : return (o-1);
231 :
232 : // Internal DoFs are associated with the elem, not its nodes
233 70840 : case 8:
234 70840 : return 0;
235 :
236 0 : default:
237 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD4/8/9!");
238 : }
239 12964603 : case HEX8:
240 1593 : libmesh_assert_less (n, 8);
241 1593 : libmesh_assert_less (o, 2);
242 : libmesh_fallthrough();
243 : case HEX20:
244 1593 : libmesh_assert_less (n, 20);
245 1593 : libmesh_assert_less (o, 2);
246 : libmesh_fallthrough();
247 : case HEX27:
248 13914242 : switch (n)
249 : {
250 289053 : case 0:
251 : case 1:
252 : case 2:
253 : case 3:
254 : case 4:
255 : case 5:
256 : case 6:
257 : case 7:
258 289053 : return 1;
259 :
260 418027 : case 8:
261 : case 9:
262 : case 10:
263 : case 11:
264 : case 12:
265 : case 13:
266 : case 14:
267 : case 15:
268 : case 16:
269 : case 17:
270 : case 18:
271 : case 19:
272 6084822 : return (o-1);
273 :
274 209388 : case 20:
275 : case 21:
276 : case 22:
277 : case 23:
278 : case 24:
279 : case 25:
280 3045864 : return ((o-1)*(o-1));
281 :
282 : // Internal DoFs are associated with the elem, not its nodes
283 68637 : case 26:
284 68637 : return 0;
285 0 : default:
286 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for HEX8/20/27!");
287 : }
288 :
289 1789578 : case PRISM6:
290 2322 : libmesh_assert_less (n, 6);
291 : libmesh_fallthrough();
292 : case PRISM15:
293 7740 : libmesh_assert_less (n, 15);
294 7740 : libmesh_assert_less (o, 2);
295 : libmesh_fallthrough();
296 : case PRISM18:
297 26721 : libmesh_assert_less (n, 18);
298 26721 : libmesh_assert_less (o, 3);
299 : libmesh_fallthrough();
300 : case PRISM20:
301 : // Internal DoFs are associated with the elem on a Prism20, or
302 : // node 20 on a Prism21
303 97794 : libmesh_assert_less (n, 20);
304 : libmesh_fallthrough();
305 : case PRISM21:
306 172422 : libmesh_assert_less (n, 21);
307 1959678 : switch (n)
308 : {
309 54756 : case 0:
310 : case 1:
311 : case 2:
312 : case 3:
313 : case 4:
314 : case 5:
315 54756 : return 1;
316 :
317 75807 : case 6:
318 : case 7:
319 : case 8:
320 : case 9:
321 : case 10:
322 : case 11:
323 : case 12:
324 : case 13:
325 : case 14:
326 861876 : return (o-1);
327 :
328 24300 : case 15:
329 : case 16:
330 : case 17:
331 276732 : return ((o-1)*(o-1));
332 :
333 14004 : case 18:
334 : case 19:
335 159624 : return ((o-1)*(o-2)/2);
336 3555 : case 20:
337 44109 : return ((o-1)*(o-1)*(o-2)/2);
338 0 : default:
339 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI!");
340 : }
341 :
342 3213441 : case TET4:
343 6939 : libmesh_assert_less (o, 2);
344 6939 : libmesh_assert_less (n, 4);
345 : libmesh_fallthrough();
346 : case TET10:
347 39420 : libmesh_assert_less (o, 3);
348 39420 : libmesh_assert_less (n, 10);
349 : libmesh_fallthrough();
350 : case TET14:
351 186498 : libmesh_assert_less (n, 14);
352 3206502 : switch (n)
353 : {
354 65556 : case 0:
355 : case 1:
356 : case 2:
357 : case 3:
358 65556 : return 1;
359 :
360 79866 : case 4:
361 : case 5:
362 : case 6:
363 : case 7:
364 : case 8:
365 : case 9:
366 1402218 : return (o-1);
367 :
368 41076 : case 10:
369 : case 11:
370 : case 12:
371 : case 13:
372 696546 : return ((o-1)*(o-2)/2);
373 :
374 0 : default:
375 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TET!");
376 : }
377 :
378 :
379 0 : case INVALID_ELEM:
380 0 : return 0;
381 :
382 0 : default:
383 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
384 : }
385 : } // HIERARCHIC_n_dofs_at_node()
386 :
387 :
388 :
389 5837126 : unsigned int HIERARCHIC_n_dofs_at_node(const Elem & e,
390 : const Order o,
391 : const unsigned int n)
392 : {
393 6238459 : return HIERARCHIC_n_dofs_at_node(e.type(), o, n);
394 : }
395 :
396 :
397 :
398 1337092 : unsigned int HIERARCHIC_n_dofs_per_elem(const ElemType t,
399 : const Order o)
400 : {
401 94871 : libmesh_assert_greater (o, 0);
402 1337092 : switch (t)
403 : {
404 0 : case NODEELEM:
405 0 : return 0;
406 4135 : case EDGE2:
407 : case EDGE3:
408 50673 : return (o-1);
409 351 : case TRI3:
410 : case TRISHELL3:
411 : case QUAD4:
412 : case QUADSHELL4:
413 351 : return 0;
414 6076 : case TRI6:
415 86812 : return ((o-1)*(o-2)/2);
416 2732 : case TRI7:
417 2732 : return 0;
418 32816 : case QUAD8:
419 : case QUADSHELL8:
420 : case QUAD9:
421 : case QUADSHELL9:
422 419437 : return ((o-1)*(o-1));
423 189 : case HEX8:
424 : case HEX20:
425 189 : libmesh_assert_less (o, 2);
426 189 : return 0;
427 29033 : case HEX27:
428 428755 : return ((o-1)*(o-1)*(o-1));
429 720 : case PRISM6:
430 : case PRISM15:
431 720 : libmesh_assert_less (o, 2);
432 : libmesh_fallthrough();
433 : case PRISM18:
434 1557 : libmesh_assert_less (o, 3);
435 1557 : return 0;
436 2799 : case PRISM20:
437 34551 : return ((o-1)*(o-1)*(o-2)/2);
438 2799 : case PRISM21: // Store interior DoFs on interior node
439 2799 : return 0;
440 1449 : case TET4:
441 1449 : libmesh_assert_less (o, 2);
442 : libmesh_fallthrough();
443 : case TET10:
444 4122 : libmesh_assert_less (o, 3);
445 : libmesh_fallthrough();
446 : case TET14:
447 242019 : return ((o-3)*(o-2)*(o-1)/6);
448 0 : case INVALID_ELEM:
449 0 : return 0;
450 0 : default:
451 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
452 : }
453 : } // HIERARCHIC_n_dofs_per_elem()
454 :
455 :
456 :
457 1185387 : unsigned int HIERARCHIC_n_dofs_per_elem(const Elem & e,
458 : const Order o)
459 : {
460 1273784 : return HIERARCHIC_n_dofs_per_elem(e.type(), o);
461 : }
462 :
463 : } // anonymous namespace
464 :
465 :
466 : // Instantiate (side_) nodal_soln() function for every dimension
467 4824 : LIBMESH_FE_NODAL_SOLN(HIERARCHIC, hierarchic_nodal_soln)
468 1248 : LIBMESH_FE_SIDE_NODAL_SOLN(HIERARCHIC)
469 :
470 :
471 : // Instantiate n_dofs*() functions for every dimension
472 27461691 : LIBMESH_DEFAULT_NDOFS(HIERARCHIC)
473 :
474 :
475 : // Hierarchic FEMs are C^0 continuous
476 0 : template <> FEContinuity FE<0,HIERARCHIC>::get_continuity() const { return C_ZERO; }
477 27744 : template <> FEContinuity FE<1,HIERARCHIC>::get_continuity() const { return C_ZERO; }
478 107090 : template <> FEContinuity FE<2,HIERARCHIC>::get_continuity() const { return C_ZERO; }
479 237044 : template <> FEContinuity FE<3,HIERARCHIC>::get_continuity() const { return C_ZERO; }
480 :
481 : // Hierarchic FEMs are hierarchic (duh!)
482 0 : template <> bool FE<0,HIERARCHIC>::is_hierarchic() const { return true; }
483 48 : template <> bool FE<1,HIERARCHIC>::is_hierarchic() const { return true; }
484 250 : template <> bool FE<2,HIERARCHIC>::is_hierarchic() const { return true; }
485 516 : template <> bool FE<3,HIERARCHIC>::is_hierarchic() const { return true; }
486 :
487 : #ifdef LIBMESH_ENABLE_AMR
488 : // compute_constraints() specializations are only needed for 2 and 3D
489 : template <>
490 17428 : void FE<2,HIERARCHIC>::compute_constraints (DofConstraints & constraints,
491 : DofMap & dof_map,
492 : const unsigned int variable_number,
493 : const Elem * elem)
494 17428 : { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
495 :
496 : template <>
497 18972 : void FE<3,HIERARCHIC>::compute_constraints (DofConstraints & constraints,
498 : DofMap & dof_map,
499 : const unsigned int variable_number,
500 : const Elem * elem)
501 18972 : { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
502 : #endif // #ifdef LIBMESH_ENABLE_AMR
503 :
504 : // Hierarchic FEM shapes need reinit
505 0 : template <> bool FE<0,HIERARCHIC>::shapes_need_reinit() const { return true; }
506 15954 : template <> bool FE<1,HIERARCHIC>::shapes_need_reinit() const { return true; }
507 177631 : template <> bool FE<2,HIERARCHIC>::shapes_need_reinit() const { return true; }
508 160714 : template <> bool FE<3,HIERARCHIC>::shapes_need_reinit() const { return true; }
509 :
510 : } // namespace libMesh
|
