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 99401773 : unsigned int BERNSTEIN_n_dofs(const ElemType t, const Order o)
107 : {
108 99401773 : switch (t)
109 : {
110 244365 : case NODEELEM:
111 244365 : return 1;
112 94273 : case EDGE2:
113 : case EDGE3:
114 1032149 : return (o+1);
115 9351 : case QUAD4:
116 : case QUADSHELL4:
117 9351 : libmesh_assert_less (o, 2);
118 : libmesh_fallthrough();
119 : case QUAD8:
120 : case QUADSHELL8:
121 : {
122 357444 : if (o == 1)
123 9351 : return 4;
124 247884 : else if (o == 2)
125 20658 : return 8;
126 : else
127 0 : libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for BERNSTEIN FE family!");
128 : }
129 5408959 : case QUAD9:
130 : case QUADSHELL9:
131 57744417 : return ((o+1)*(o+1));
132 279855 : case HEX8:
133 279855 : libmesh_assert_less (o, 2);
134 : libmesh_fallthrough();
135 : case HEX20:
136 : {
137 11313804 : if (o == 1)
138 279855 : return 8;
139 7960308 : else if (o == 2)
140 663360 : return 20;
141 : else
142 0 : libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for BERNSTEIN FE family!");
143 : }
144 1993385 : case HEX27:
145 26379444 : return ((o+1)*(o+1)*(o+1));
146 8689 : case TRI3:
147 : case TRISHELL3:
148 8689 : libmesh_assert_less (o, 2);
149 : libmesh_fallthrough();
150 : case TRI6:
151 : case TRI7:
152 900066 : return ((o+1)*(o+2)/2);
153 31688 : case TET4:
154 31688 : libmesh_assert_less (o, 2);
155 : libmesh_fallthrough();
156 : case TET10:
157 : case TET14:
158 40130 : libmesh_assert_less (o, 3);
159 990526 : 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 91051801 : unsigned int BERNSTEIN_n_dofs(const Elem * e, const Order o)
170 : {
171 8825417 : libmesh_assert(e);
172 99401773 : return BERNSTEIN_n_dofs(e->type(), o);
173 : }
174 :
175 :
176 :
177 19117631 : unsigned int BERNSTEIN_n_dofs_at_node(const ElemType t,
178 : const Order o,
179 : const unsigned int n)
180 : {
181 19117631 : switch (t)
182 : {
183 184751 : case NODEELEM:
184 184751 : return 1;
185 :
186 23952 : case EDGE2:
187 : case EDGE3:
188 23952 : switch (n)
189 : {
190 1704 : case 0:
191 : case 1:
192 1704 : return 1;
193 576 : case 2:
194 576 : libmesh_assert (o>1);
195 6384 : return (o-1);
196 0 : default:
197 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for EDGE2/3!");
198 : }
199 256812 : case TRI3:
200 2112 : libmesh_assert_less (n, 3);
201 2112 : 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 12660 : libmesh_assert_less (n, 6);
206 : libmesh_fallthrough();
207 : case TRI7:
208 279624 : switch (n)
209 : {
210 13128 : case 0:
211 : case 1:
212 : case 2:
213 13128 : return 1;
214 :
215 10080 : case 3:
216 : case 4:
217 : case 5:
218 113472 : return (o-1);
219 :
220 1716 : case 6:
221 19368 : return ((o-1)*(o-2)/2);
222 0 : default:
223 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI!");
224 : }
225 11483068 : case QUAD4:
226 1368 : libmesh_assert_less (n, 4);
227 1368 : libmesh_assert_less (o, 2);
228 : libmesh_fallthrough();
229 : case QUAD8:
230 1368 : libmesh_assert_less (n, 8);
231 1368 : libmesh_assert_less (o, 3);
232 : libmesh_fallthrough();
233 : case QUAD9:
234 : {
235 12836000 : switch (n)
236 : {
237 610764 : case 0:
238 : case 1:
239 : case 2:
240 : case 3:
241 610764 : return 1;
242 :
243 594644 : case 4:
244 : case 5:
245 : case 6:
246 : case 7:
247 5631751 : return (o-1);
248 :
249 148892 : case 8:
250 1410578 : 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 1537555 : case HEX8:
257 4248 : libmesh_assert_less (n, 8);
258 4248 : libmesh_assert_less (o, 2);
259 : libmesh_fallthrough();
260 : case HEX20:
261 4248 : libmesh_assert_less (n, 20);
262 4248 : libmesh_assert_less (o, 3);
263 : libmesh_fallthrough();
264 : case HEX27:
265 1660011 : switch (n)
266 : {
267 41240 : case 0:
268 : case 1:
269 : case 2:
270 : case 3:
271 : case 4:
272 : case 5:
273 : case 6:
274 : case 7:
275 41240 : return 1;
276 :
277 53920 : 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 710785 : return (o-1);
290 :
291 26992 : case 20:
292 : case 21:
293 : case 22:
294 : case 23:
295 : case 24:
296 : case 25:
297 355687 : return ((o-1)*(o-1));
298 :
299 4552 : case 26:
300 59911 : 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 2452776 : case TET4:
306 18504 : libmesh_assert_less (n, 4);
307 18504 : libmesh_assert_less (o, 2);
308 : libmesh_fallthrough();
309 : case TET10:
310 61812 : libmesh_assert_less (n, 10);
311 : libmesh_fallthrough();
312 : case TET14:
313 122352 : libmesh_assert_less (o, 3);
314 2434272 : switch (n)
315 : {
316 56016 : case 0:
317 : case 1:
318 : case 2:
319 : case 3:
320 56016 : return 1;
321 :
322 49104 : case 4:
323 : case 5:
324 : case 6:
325 : case 7:
326 : case 8:
327 : case 9:
328 954192 : return (o-1);
329 :
330 17232 : case 10:
331 : case 11:
332 : case 12:
333 : case 13:
334 17232 : 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 2356001 : unsigned int BERNSTEIN_n_dofs_at_node(const Elem & e,
349 : const Order o,
350 : const unsigned int n)
351 : {
352 2515489 : return BERNSTEIN_n_dofs_at_node(e.type(), o, n);
353 : }
354 :
355 :
356 :
357 3425570 : unsigned int BERNSTEIN_n_dofs_per_elem(const ElemType t, const Order o)
358 : {
359 3078209 : switch (t)
360 : {
361 180661 : case NODEELEM:
362 : case EDGE2:
363 : case EDGE3:
364 180661 : return 0;
365 936 : case TRI3:
366 : case TRISHELL3:
367 : case QUAD4:
368 : case QUADSHELL4:
369 936 : return 0;
370 1500 : case TRI6:
371 16848 : return ((o-1)*(o-2)/2);
372 1500 : case TRI7:
373 1500 : 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 147760 : return 0;
382 504 : case HEX8:
383 504 : libmesh_assert_less (o, 2);
384 504 : return 0;
385 0 : case HEX20:
386 0 : libmesh_assert_less (o, 3);
387 0 : return 0;
388 3508 : case HEX27:
389 3508 : return 0;
390 3864 : case TET4:
391 3864 : libmesh_assert_less (o, 2);
392 : libmesh_fallthrough();
393 : case TET10:
394 : case TET14:
395 10992 : libmesh_assert_less (o, 3);
396 10992 : 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 3078209 : unsigned int BERNSTEIN_n_dofs_per_elem(const Elem & e, const Order o)
407 : {
408 3425570 : 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 121944974 : 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 17520 : template <> FEContinuity FE<1,BERNSTEIN>::get_continuity() const { return C_ZERO; }
424 58006 : template <> FEContinuity FE<2,BERNSTEIN>::get_continuity() const { return C_ZERO; }
425 71878 : 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 2592 : void FE<2,BERNSTEIN>::compute_constraints (DofConstraints & constraints,
437 : DofMap & dof_map,
438 : const unsigned int variable_number,
439 : const Elem * elem)
440 2592 : { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
441 :
442 : template <>
443 10944 : void FE<3,BERNSTEIN>::compute_constraints (DofConstraints & constraints,
444 : DofMap & dof_map,
445 : const unsigned int variable_number,
446 : const Elem * elem)
447 10944 : { 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 576 : template <> bool FE<1,BERNSTEIN>::shapes_need_reinit() const { return true; }
454 7038 : template <> bool FE<2,BERNSTEIN>::shapes_need_reinit() const { return true; }
455 29118 : template <> bool FE<3,BERNSTEIN>::shapes_need_reinit() const { return true; }
456 :
457 : } // namespace libMesh
458 :
459 : #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
|
