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/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 95935276 : unsigned int BERNSTEIN_n_dofs(const ElemType t, const Order o)
107 : {
108 95935276 : switch (t)
109 : {
110 244365 : case NODEELEM:
111 244365 : return 1;
112 92790 : case EDGE2:
113 : case EDGE3:
114 1016777 : return (o+1);
115 8255 : case QUAD4:
116 : case QUADSHELL4:
117 8255 : libmesh_assert_less (o, 2);
118 : libmesh_fallthrough();
119 : case QUAD8:
120 : case QUADSHELL8:
121 : {
122 344832 : if (o == 1)
123 8255 : return 4;
124 247752 : else if (o == 2)
125 20646 : return 8;
126 : else
127 0 : libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for BERNSTEIN FE family!");
128 : }
129 5319331 : case QUAD9:
130 : case QUADSHELL9:
131 56326627 : return ((o+1)*(o+1));
132 259767 : case HEX8:
133 259767 : libmesh_assert_less (o, 2);
134 : libmesh_fallthrough();
135 : case HEX20:
136 : {
137 11073816 : if (o == 1)
138 259767 : return 8;
139 7960176 : else if (o == 2)
140 663348 : return 20;
141 : else
142 0 : libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for BERNSTEIN FE family!");
143 : }
144 1927294 : case HEX27:
145 25589168 : return ((o+1)*(o+1)*(o+1));
146 7633 : case TRI3:
147 : case TRISHELL3:
148 7633 : libmesh_assert_less (o, 2);
149 : libmesh_fallthrough();
150 : case TRI6:
151 : case TRI7:
152 115578 : return ((o+1)*(o+2)/2);
153 25778 : case TET4:
154 25778 : libmesh_assert_less (o, 2);
155 : libmesh_fallthrough();
156 : case TET10:
157 : case TET14:
158 31826 : libmesh_assert_less (o, 3);
159 778135 : return ((o+1)*(o+2)*(o+3)/6);
160 0 : case INVALID_ELEM:
161 0 : return 0;
162 0 : default:
163 0 : libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for BERNSTEIN FE family!");
164 : }
165 : } // BERNSTEIN_n_dofs()
166 :
167 :
168 :
169 87826436 : unsigned int BERNSTEIN_n_dofs(const Elem * e, const Order o)
170 : {
171 8580151 : libmesh_assert(e);
172 95935276 : return BERNSTEIN_n_dofs(e->type(), o);
173 : }
174 :
175 :
176 :
177 18025679 : unsigned int BERNSTEIN_n_dofs_at_node(const ElemType t,
178 : const Order o,
179 : const unsigned int n)
180 : {
181 18025679 : switch (t)
182 : {
183 184751 : case NODEELEM:
184 184751 : return 1;
185 :
186 17964 : case EDGE2:
187 : case EDGE3:
188 17964 : switch (n)
189 : {
190 1278 : case 0:
191 : case 1:
192 1278 : return 1;
193 432 : case 2:
194 432 : libmesh_assert (o>1);
195 4788 : return (o-1);
196 0 : default:
197 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for EDGE2/3!");
198 : }
199 192609 : case TRI3:
200 1584 : libmesh_assert_less (n, 3);
201 1584 : libmesh_assert_less (o, 2);
202 : libmesh_fallthrough();
203 : case TRI6:
204 : // Internal DoFs are associated with the elem on a Tri6, or node 6 on a Tri7
205 9495 : libmesh_assert_less (n, 6);
206 : libmesh_fallthrough();
207 : case TRI7:
208 209718 : switch (n)
209 : {
210 9846 : case 0:
211 : case 1:
212 : case 2:
213 9846 : return 1;
214 :
215 7560 : case 3:
216 : case 4:
217 : case 5:
218 85104 : return (o-1);
219 :
220 1287 : case 6:
221 14526 : return ((o-1)*(o-2)/2);
222 0 : default:
223 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI!");
224 : }
225 11422786 : case QUAD4:
226 1026 : libmesh_assert_less (n, 4);
227 1026 : libmesh_assert_less (o, 2);
228 : libmesh_fallthrough();
229 : case QUAD8:
230 1026 : libmesh_assert_less (n, 8);
231 1026 : libmesh_assert_less (o, 3);
232 : libmesh_fallthrough();
233 : case QUAD9:
234 : {
235 12770438 : switch (n)
236 : {
237 608026 : case 0:
238 : case 1:
239 : case 2:
240 : case 3:
241 608026 : return 1;
242 :
243 592344 : case 4:
244 : case 5:
245 : case 6:
246 : case 7:
247 5604591 : return (o-1);
248 :
249 148308 : case 8:
250 1403680 : return ((o-1)*(o-1));
251 :
252 0 : default:
253 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD9!");
254 : }
255 : }
256 1251511 : case HEX8:
257 3186 : libmesh_assert_less (n, 8);
258 3186 : libmesh_assert_less (o, 2);
259 : libmesh_fallthrough();
260 : case HEX20:
261 3186 : libmesh_assert_less (n, 20);
262 3186 : libmesh_assert_less (o, 3);
263 : libmesh_fallthrough();
264 : case HEX27:
265 1348671 : switch (n)
266 : {
267 32550 : case 0:
268 : case 1:
269 : case 2:
270 : case 3:
271 : case 4:
272 : case 5:
273 : case 6:
274 : case 7:
275 32550 : return 1;
276 :
277 42772 : case 8:
278 : case 9:
279 : case 10:
280 : case 11:
281 : case 12:
282 : case 13:
283 : case 14:
284 : case 15:
285 : case 16:
286 : case 17:
287 : case 18:
288 : case 19:
289 578401 : return (o-1);
290 :
291 21412 : case 20:
292 : case 21:
293 : case 22:
294 : case 23:
295 : case 24:
296 : case 25:
297 289423 : return ((o-1)*(o-1));
298 :
299 3612 : case 26:
300 48747 : return ((o-1)*(o-1)*(o-1));
301 :
302 0 : default:
303 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for HEX27!");
304 : }
305 1839582 : case TET4:
306 13878 : libmesh_assert_less (n, 4);
307 13878 : libmesh_assert_less (o, 2);
308 : libmesh_fallthrough();
309 : case TET10:
310 46359 : libmesh_assert_less (n, 10);
311 : libmesh_fallthrough();
312 : case TET14:
313 91764 : libmesh_assert_less (o, 3);
314 1825704 : switch (n)
315 : {
316 42012 : case 0:
317 : case 1:
318 : case 2:
319 : case 3:
320 42012 : return 1;
321 :
322 36828 : case 4:
323 : case 5:
324 : case 6:
325 : case 7:
326 : case 8:
327 : case 9:
328 715644 : return (o-1);
329 :
330 12924 : case 10:
331 : case 11:
332 : case 12:
333 : case 13:
334 12924 : return 0;
335 :
336 0 : default:
337 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TET4/10/14!");
338 : }
339 0 : case INVALID_ELEM:
340 0 : return 0;
341 0 : default:
342 0 : libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for BERNSTEIN FE family!");
343 : }
344 : } // BERNSTEIN_n_dofs_at_node()
345 :
346 :
347 :
348 1955531 : unsigned int BERNSTEIN_n_dofs_at_node(const Elem & e,
349 : const Order o,
350 : const unsigned int n)
351 : {
352 2093059 : return BERNSTEIN_n_dofs_at_node(e.type(), o, n);
353 : }
354 :
355 :
356 :
357 3339576 : unsigned int BERNSTEIN_n_dofs_per_elem(const ElemType t, const Order o)
358 : {
359 2997519 : switch (t)
360 : {
361 180475 : case NODEELEM:
362 : case EDGE2:
363 : case EDGE3:
364 180475 : return 0;
365 702 : case TRI3:
366 : case TRISHELL3:
367 : case QUAD4:
368 : case QUADSHELL4:
369 702 : return 0;
370 1125 : case TRI6:
371 12636 : return ((o-1)*(o-2)/2);
372 1125 : case TRI7:
373 1125 : return 0;
374 0 : case QUAD8:
375 : case QUADSHELL8:
376 0 : if (o <= 2)
377 0 : return 0;
378 : libmesh_fallthrough();
379 : case QUAD9:
380 : case QUADSHELL9:
381 147224 : return 0;
382 378 : case HEX8:
383 378 : libmesh_assert_less (o, 2);
384 378 : return 0;
385 0 : case HEX20:
386 0 : libmesh_assert_less (o, 3);
387 0 : return 0;
388 2784 : case HEX27:
389 2784 : return 0;
390 2898 : case TET4:
391 2898 : libmesh_assert_less (o, 2);
392 : libmesh_fallthrough();
393 : case TET10:
394 : case TET14:
395 8244 : libmesh_assert_less (o, 3);
396 8244 : return 0;
397 0 : case INVALID_ELEM:
398 0 : return 0;
399 0 : default:
400 0 : libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for BERNSTEIN FE family!");
401 : }
402 : } // BERNSTEIN_n_dofs_per_elem
403 :
404 :
405 :
406 2997519 : unsigned int BERNSTEIN_n_dofs_per_elem(const Elem & e, const Order o)
407 : {
408 3339576 : return BERNSTEIN_n_dofs_per_elem(e.type(), o);
409 : }
410 :
411 : } // anonymous namespace
412 :
413 :
414 : // Instantiate (side_) nodal_soln() function for every dimension
415 0 : LIBMESH_FE_NODAL_SOLN(BERNSTEIN, bernstein_nodal_soln)
416 0 : LIBMESH_FE_SIDE_NODAL_SOLN(BERNSTEIN)
417 :
418 : // Instantiate n_dofs*() functions for every dimension
419 117300531 : LIBMESH_DEFAULT_NDOFS(BERNSTEIN)
420 :
421 : // Bernstein FEMs are C^0 continuous
422 0 : template <> FEContinuity FE<0,BERNSTEIN>::get_continuity() const { return C_ZERO; }
423 13260 : template <> FEContinuity FE<1,BERNSTEIN>::get_continuity() const { return C_ZERO; }
424 44446 : template <> FEContinuity FE<2,BERNSTEIN>::get_continuity() const { return C_ZERO; }
425 54673 : template <> FEContinuity FE<3,BERNSTEIN>::get_continuity() const { return C_ZERO; }
426 :
427 : // Bernstein FEMs are not hierarchic
428 0 : template <> bool FE<0,BERNSTEIN>::is_hierarchic() const { return false; }
429 48 : template <> bool FE<1,BERNSTEIN>::is_hierarchic() const { return false; }
430 202 : template <> bool FE<2,BERNSTEIN>::is_hierarchic() const { return false; }
431 250 : template <> bool FE<3,BERNSTEIN>::is_hierarchic() const { return false; }
432 :
433 : #ifdef LIBMESH_ENABLE_AMR
434 : // compute_constraints() specializations are only needed for 2 and 3D
435 : template <>
436 1944 : void FE<2,BERNSTEIN>::compute_constraints (DofConstraints & constraints,
437 : DofMap & dof_map,
438 : const unsigned int variable_number,
439 : const Elem * elem)
440 1944 : { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
441 :
442 : template <>
443 8208 : void FE<3,BERNSTEIN>::compute_constraints (DofConstraints & constraints,
444 : DofMap & dof_map,
445 : const unsigned int variable_number,
446 : const Elem * elem)
447 8208 : { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
448 : #endif // #ifdef LIBMESH_ENABLE_AMR
449 :
450 : // Bernstein shapes need reinit only for approximation orders >= 3,
451 : // but we might reach that with p refinement
452 0 : template <> bool FE<0,BERNSTEIN>::shapes_need_reinit() const { return true; }
453 468 : template <> bool FE<1,BERNSTEIN>::shapes_need_reinit() const { return true; }
454 5025 : template <> bool FE<2,BERNSTEIN>::shapes_need_reinit() const { return true; }
455 21756 : template <> bool FE<3,BERNSTEIN>::shapes_need_reinit() const { return true; }
456 :
457 : } // namespace libMesh
458 :
459 : #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
|
