Line data Source code
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/libmesh_config.h"
22 :
23 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
24 :
25 : #include "libmesh/elem.h"
26 : #include "libmesh/enum_to_string.h"
27 : #include "libmesh/fe.h"
28 : #include "libmesh/fe_interface.h"
29 : #include "libmesh/fe_macro.h"
30 :
31 : namespace libMesh
32 : {
33 :
34 : // ------------------------------------------------------------
35 : // Bernstein-specific implementations, Steffen Petersen 2005
36 :
37 : // Anonymous namespace for local helper functions
38 : namespace {
39 :
40 0 : void bernstein_nodal_soln(const Elem * elem,
41 : const Order order,
42 : const std::vector<Number> & elem_soln,
43 : std::vector<Number> & nodal_soln,
44 : const bool add_p_level)
45 : {
46 0 : const unsigned int n_nodes = elem->n_nodes();
47 :
48 0 : const ElemType elem_type = elem->type();
49 :
50 0 : nodal_soln.resize(n_nodes);
51 :
52 0 : const Order totalorder = order + add_p_level*elem->p_level();
53 :
54 : // FEType object to be passed to various FEInterface functions below.
55 0 : FEType fe_type(order, BERNSTEIN);
56 :
57 0 : switch (totalorder)
58 : {
59 : // Constant shape functions
60 0 : case CONSTANT:
61 : {
62 0 : libmesh_assert_equal_to (elem_soln.size(), 1);
63 :
64 0 : std::fill(nodal_soln.begin(), nodal_soln.end(), elem_soln[0]);
65 :
66 0 : return;
67 : }
68 :
69 :
70 : // For other bases do interpolation at the nodes
71 : // explicitly.
72 0 : case FIRST:
73 : case SECOND:
74 : case THIRD:
75 : case FOURTH:
76 : case FIFTH:
77 : case SIXTH:
78 : {
79 : const unsigned int n_sf =
80 0 : FEInterface::n_shape_functions(fe_type, elem);
81 :
82 0 : std::vector<Point> refspace_nodes;
83 0 : FEBase::get_refspace_nodes(elem_type,refspace_nodes);
84 0 : libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
85 0 : libmesh_assert_equal_to (elem_soln.size(), n_sf);
86 :
87 : // Zero before summation
88 0 : std::fill(nodal_soln.begin(), nodal_soln.end(), 0);
89 :
90 0 : for (unsigned int n=0; n<n_nodes; n++)
91 : // u_i = Sum (alpha_i phi_i)
92 0 : for (unsigned int i=0; i<n_sf; i++)
93 0 : nodal_soln[n] += elem_soln[i] *
94 0 : FEInterface::shape(fe_type, elem, i, refspace_nodes[n]);
95 :
96 0 : return;
97 : }
98 :
99 0 : default:
100 0 : libmesh_error_msg("ERROR: Invalid total order " << totalorder);
101 : }
102 : } // bernstein_nodal_soln()
103 :
104 :
105 :
106 99408003 : unsigned int BERNSTEIN_n_dofs(const ElemType t, const Order o)
107 : {
108 99408003 : switch (t)
109 : {
110 244365 : case NODEELEM:
111 244365 : return 1;
112 94504 : case EDGE2:
113 : case EDGE3:
114 : case EDGE4:
115 1033293 : return (o+1);
116 9371 : case QUAD4:
117 : case QUADSHELL4:
118 9371 : libmesh_assert_less (o, 2);
119 : libmesh_fallthrough();
120 : case QUAD8:
121 : case QUADSHELL8:
122 : {
123 357464 : if (o == 1)
124 9371 : return 4;
125 247884 : else if (o == 2)
126 20658 : return 8;
127 : else
128 0 : libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for BERNSTEIN FE family!");
129 : }
130 5409064 : case QUAD9:
131 : case QUADSHELL9:
132 57744797 : return ((o+1)*(o+1));
133 279891 : case HEX8:
134 279891 : libmesh_assert_less (o, 2);
135 : libmesh_fallthrough();
136 : case HEX20:
137 : {
138 11313840 : if (o == 1)
139 279891 : return 8;
140 7960308 : else if (o == 2)
141 663360 : return 20;
142 : else
143 0 : libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for BERNSTEIN FE family!");
144 : }
145 1993704 : case HEX27:
146 26380808 : return ((o+1)*(o+1)*(o+1));
147 8729 : case TRI3:
148 : case TRISHELL3:
149 8729 : libmesh_assert_less (o, 2);
150 : libmesh_fallthrough();
151 : case TRI6:
152 : case TRI7:
153 900836 : return ((o+1)*(o+2)/2);
154 31866 : case TET4:
155 31866 : libmesh_assert_less (o, 2);
156 : libmesh_fallthrough();
157 : case TET10:
158 : case TET14:
159 40882 : libmesh_assert_less (o, 3);
160 993166 : return ((o+1)*(o+2)*(o+3)/6);
161 0 : case INVALID_ELEM:
162 0 : return 0;
163 0 : default:
164 0 : libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for BERNSTEIN FE family!");
165 : }
166 : } // BERNSTEIN_n_dofs()
167 :
168 :
169 :
170 91057660 : unsigned int BERNSTEIN_n_dofs(const Elem * e, const Order o)
171 : {
172 8827178 : libmesh_assert(e);
173 99408003 : return BERNSTEIN_n_dofs(e->type(), o);
174 : }
175 :
176 :
177 :
178 19476843 : unsigned int BERNSTEIN_n_dofs_at_node(const ElemType t,
179 : const Order o,
180 : const unsigned int n)
181 : {
182 19476843 : switch (t)
183 : {
184 184751 : case NODEELEM:
185 184751 : return 1;
186 :
187 25948 : case EDGE2:
188 : case EDGE3:
189 25948 : switch (n)
190 : {
191 1846 : case 0:
192 : case 1:
193 1846 : return 1;
194 624 : case 2:
195 624 : libmesh_assert (o>1);
196 6916 : return (o-1);
197 0 : default:
198 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for EDGE2/3!");
199 : }
200 6526 : case EDGE4:
201 6526 : switch (n)
202 : {
203 325 : case 0:
204 : case 1:
205 325 : return 1;
206 143 : case 2:
207 143 : libmesh_assert (o>2);
208 1573 : return (o-1)/2; // floor
209 1573 : case 3:
210 143 : libmesh_assert (o>2);
211 1573 : return o/2; // "ceiling", add extra DoFs here first
212 0 : default:
213 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for EDGE4!");
214 : }
215 278213 : case TRI3:
216 2288 : libmesh_assert_less (n, 3);
217 2288 : libmesh_assert_less (o, 2);
218 : libmesh_fallthrough();
219 : case TRI6:
220 : // Internal DoFs are associated with the elem on a Tri6, or node 6 on a Tri7
221 13715 : libmesh_assert_less (n, 6);
222 : libmesh_fallthrough();
223 : case TRI7:
224 302926 : switch (n)
225 : {
226 14222 : case 0:
227 : case 1:
228 : case 2:
229 14222 : return 1;
230 :
231 10920 : case 3:
232 : case 4:
233 : case 5:
234 122928 : return (o-1);
235 :
236 1859 : case 6:
237 20982 : return ((o-1)*(o-2)/2);
238 0 : default:
239 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI!");
240 : }
241 11494516 : case QUAD4:
242 1482 : libmesh_assert_less (n, 4);
243 1482 : libmesh_assert_less (o, 2);
244 : libmesh_fallthrough();
245 : case QUAD8:
246 1482 : libmesh_assert_less (n, 8);
247 1482 : libmesh_assert_less (o, 3);
248 : libmesh_fallthrough();
249 : case QUAD9:
250 : {
251 12848428 : switch (n)
252 : {
253 611330 : case 0:
254 : case 1:
255 : case 2:
256 : case 3:
257 611330 : return 1;
258 :
259 595064 : case 4:
260 : case 5:
261 : case 6:
262 : case 7:
263 5636615 : return (o-1);
264 :
265 149000 : case 8:
266 1411830 : return ((o-1)*(o-1));
267 :
268 0 : default:
269 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD9!");
270 : }
271 : }
272 1631175 : case HEX8:
273 4602 : libmesh_assert_less (n, 8);
274 4602 : libmesh_assert_less (o, 2);
275 : libmesh_fallthrough();
276 : case HEX20:
277 4602 : libmesh_assert_less (n, 20);
278 4602 : libmesh_assert_less (o, 3);
279 : libmesh_fallthrough();
280 : case HEX27:
281 1761919 : switch (n)
282 : {
283 44094 : case 0:
284 : case 1:
285 : case 2:
286 : case 3:
287 : case 4:
288 : case 5:
289 : case 6:
290 : case 7:
291 44094 : return 1;
292 :
293 57572 : case 8:
294 : case 9:
295 : case 10:
296 : case 11:
297 : case 12:
298 : case 13:
299 : case 14:
300 : case 15:
301 : case 16:
302 : case 17:
303 : case 18:
304 : case 19:
305 754081 : return (o-1);
306 :
307 28820 : case 20:
308 : case 21:
309 : case 22:
310 : case 23:
311 : case 24:
312 : case 25:
313 377359 : return ((o-1)*(o-1));
314 :
315 4860 : case 26:
316 63563 : return ((o-1)*(o-1)*(o-1));
317 :
318 0 : default:
319 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for HEX27!");
320 : }
321 2657174 : case TET4:
322 20046 : libmesh_assert_less (n, 4);
323 20046 : libmesh_assert_less (o, 2);
324 : libmesh_fallthrough();
325 : case TET10:
326 66963 : libmesh_assert_less (n, 10);
327 : libmesh_fallthrough();
328 : case TET14:
329 132548 : libmesh_assert_less (o, 3);
330 2637128 : switch (n)
331 : {
332 60684 : case 0:
333 : case 1:
334 : case 2:
335 : case 3:
336 60684 : return 1;
337 :
338 53196 : case 4:
339 : case 5:
340 : case 6:
341 : case 7:
342 : case 8:
343 : case 9:
344 1033708 : return (o-1);
345 :
346 18668 : case 10:
347 : case 11:
348 : case 12:
349 : case 13:
350 18668 : return 0;
351 :
352 0 : default:
353 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TET4/10/14!");
354 : }
355 0 : case INVALID_ELEM:
356 0 : return 0;
357 0 : default:
358 0 : libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for BERNSTEIN FE family!");
359 : }
360 : } // BERNSTEIN_n_dofs_at_node()
361 :
362 :
363 :
364 2490934 : unsigned int BERNSTEIN_n_dofs_at_node(const Elem & e,
365 : const Order o,
366 : const unsigned int n)
367 : {
368 2657937 : return BERNSTEIN_n_dofs_at_node(e.type(), o, n);
369 : }
370 :
371 :
372 :
373 3454691 : unsigned int BERNSTEIN_n_dofs_per_elem(const ElemType t, const Order o)
374 : {
375 3105511 : switch (t)
376 : {
377 180866 : case NODEELEM:
378 : case EDGE2:
379 : case EDGE3:
380 : case EDGE4:
381 180866 : return 0;
382 1014 : case TRI3:
383 : case TRISHELL3:
384 : case QUAD4:
385 : case QUADSHELL4:
386 1014 : return 0;
387 1625 : case TRI6:
388 18252 : return ((o-1)*(o-2)/2);
389 1625 : case TRI7:
390 1625 : return 0;
391 0 : case QUAD8:
392 : case QUADSHELL8:
393 0 : if (o <= 2)
394 0 : return 0;
395 : libmesh_fallthrough();
396 : case QUAD9:
397 : case QUADSHELL9:
398 147852 : return 0;
399 546 : case HEX8:
400 546 : libmesh_assert_less (o, 2);
401 546 : return 0;
402 0 : case HEX20:
403 0 : libmesh_assert_less (o, 3);
404 0 : return 0;
405 3744 : case HEX27:
406 3744 : return 0;
407 4186 : case TET4:
408 4186 : libmesh_assert_less (o, 2);
409 : libmesh_fallthrough();
410 : case TET10:
411 : case TET14:
412 11908 : libmesh_assert_less (o, 3);
413 11908 : return 0;
414 0 : case INVALID_ELEM:
415 0 : return 0;
416 0 : default:
417 0 : libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for BERNSTEIN FE family!");
418 : }
419 : } // BERNSTEIN_n_dofs_per_elem
420 :
421 :
422 :
423 3105511 : unsigned int BERNSTEIN_n_dofs_per_elem(const Elem & e, const Order o)
424 : {
425 3454691 : return BERNSTEIN_n_dofs_per_elem(e.type(), o);
426 : }
427 :
428 : } // anonymous namespace
429 :
430 :
431 : // Instantiate (side_) nodal_soln() function for every dimension
432 0 : LIBMESH_FE_NODAL_SOLN(BERNSTEIN, bernstein_nodal_soln)
433 0 : LIBMESH_FE_SIDE_NODAL_SOLN(BERNSTEIN)
434 :
435 : // Instantiate n_dofs*() functions for every dimension
436 122339537 : LIBMESH_DEFAULT_NDOFS(BERNSTEIN)
437 :
438 : // Bernstein FEMs are C^0 continuous
439 0 : template <> FEContinuity FE<0,BERNSTEIN>::get_continuity() const { return C_ZERO; }
440 23759 : template <> FEContinuity FE<1,BERNSTEIN>::get_continuity() const { return C_ZERO; }
441 62823 : template <> FEContinuity FE<2,BERNSTEIN>::get_continuity() const { return C_ZERO; }
442 77847 : template <> FEContinuity FE<3,BERNSTEIN>::get_continuity() const { return C_ZERO; }
443 :
444 : // Bernstein FEMs are not hierarchic
445 0 : template <> bool FE<0,BERNSTEIN>::is_hierarchic() const { return false; }
446 60 : template <> bool FE<1,BERNSTEIN>::is_hierarchic() const { return false; }
447 202 : template <> bool FE<2,BERNSTEIN>::is_hierarchic() const { return false; }
448 250 : template <> bool FE<3,BERNSTEIN>::is_hierarchic() const { return false; }
449 :
450 : #ifdef LIBMESH_ENABLE_AMR
451 : // compute_constraints() specializations are only needed for 2 and 3D
452 : template <>
453 2808 : void FE<2,BERNSTEIN>::compute_constraints (DofConstraints & constraints,
454 : DofMap & dof_map,
455 : const unsigned int variable_number,
456 : const Elem * elem)
457 2808 : { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
458 :
459 : template <>
460 11856 : void FE<3,BERNSTEIN>::compute_constraints (DofConstraints & constraints,
461 : DofMap & dof_map,
462 : const unsigned int variable_number,
463 : const Elem * elem)
464 11856 : { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
465 : #endif // #ifdef LIBMESH_ENABLE_AMR
466 :
467 : // Bernstein shapes need reinit only for approximation orders >= 3,
468 : // but we might reach that with p refinement
469 0 : template <> bool FE<0,BERNSTEIN>::shapes_need_reinit() const { return true; }
470 804 : template <> bool FE<1,BERNSTEIN>::shapes_need_reinit() const { return true; }
471 7574 : template <> bool FE<2,BERNSTEIN>::shapes_need_reinit() const { return true; }
472 31482 : template <> bool FE<3,BERNSTEIN>::shapes_need_reinit() const { return true; }
473 :
474 : } // namespace libMesh
475 :
476 : #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
|
