libMesh
Loading...
Searching...
No Matches
fe_lagrange.C
Go to the documentation of this file.
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/dof_map.h"
22#include "libmesh/elem.h"
23#include "libmesh/enum_to_string.h"
24#include "libmesh/fe.h"
25#include "libmesh/fe_interface.h"
26#include "libmesh/fe_macro.h"
27#include "libmesh/remote_elem.h"
28#include "libmesh/threads.h"
29
30
31#ifdef LIBMESH_ENABLE_AMR
32#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
33// to have the macro 'inf_fe_switch' available.
34#include "libmesh/fe_interface_macros.h"
35#include "libmesh/inf_fe.h"
36#endif
37#endif
38
39namespace libMesh
40{
41
42// global helper function
43void lagrange_nodal_soln(const Elem * elem,
44 const Order order,
45 const std::vector<Number> & elem_soln,
46 std::vector<Number> & nodal_soln,
47 const bool add_p_level)
48{
49 const unsigned int n_nodes = elem->n_nodes();
50 const ElemType type = elem->type();
51
52 const Order totalorder = order + add_p_level*elem->p_level();
53
54 nodal_soln.resize(n_nodes);
55
56
57
58 switch (totalorder)
59 {
60 // linear Lagrange shape functions
61 case FIRST:
62 {
63 switch (type)
64 {
65 case EDGE3:
66 {
67 libmesh_assert_equal_to (elem_soln.size(), 2);
68 libmesh_assert_equal_to (nodal_soln.size(), 3);
69
70 nodal_soln[0] = elem_soln[0];
71 nodal_soln[1] = elem_soln[1];
72 nodal_soln[2] = .5*(elem_soln[0] + elem_soln[1]);
73
74 return;
75 }
76
77 case EDGE4:
78 {
79 libmesh_assert_equal_to (elem_soln.size(), 2);
80 libmesh_assert_equal_to (nodal_soln.size(), 4);
81
82 nodal_soln[0] = elem_soln[0];
83 nodal_soln[1] = elem_soln[1];
84 nodal_soln[2] = (2.*elem_soln[0] + elem_soln[1])/3.;
85 nodal_soln[3] = (elem_soln[0] + 2.*elem_soln[1])/3.;
86
87 return;
88 }
89
90
91 case TRI7:
92 libmesh_assert_equal_to (nodal_soln.size(), 7);
93 nodal_soln[6] = (elem_soln[0] + elem_soln[1] + elem_soln[2])/3.;
94 libmesh_fallthrough();
95 case TRI6:
96 {
97 libmesh_assert (type == TRI7 || nodal_soln.size() == 6);
98 libmesh_assert_equal_to (elem_soln.size(), 3);
99
100 nodal_soln[0] = elem_soln[0];
101 nodal_soln[1] = elem_soln[1];
102 nodal_soln[2] = elem_soln[2];
103 nodal_soln[3] = .5*(elem_soln[0] + elem_soln[1]);
104 nodal_soln[4] = .5*(elem_soln[1] + elem_soln[2]);
105 nodal_soln[5] = .5*(elem_soln[2] + elem_soln[0]);
106
107 return;
108 }
109
110
111 case QUAD8:
112 case QUAD9:
113 {
114 libmesh_assert_equal_to (elem_soln.size(), 4);
115
116 if (type == QUAD8)
117 libmesh_assert_equal_to (nodal_soln.size(), 8);
118 else
119 libmesh_assert_equal_to (nodal_soln.size(), 9);
120
121
122 nodal_soln[0] = elem_soln[0];
123 nodal_soln[1] = elem_soln[1];
124 nodal_soln[2] = elem_soln[2];
125 nodal_soln[3] = elem_soln[3];
126 nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]);
127 nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]);
128 nodal_soln[6] = .5*(elem_soln[2] + elem_soln[3]);
129 nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]);
130
131 if (type == QUAD9)
132 nodal_soln[8] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
133
134 return;
135 }
136
137
138 case TET14:
139 libmesh_assert_equal_to (nodal_soln.size(), 14);
140 nodal_soln[10] = (elem_soln[0] + elem_soln[1] + elem_soln[2])/3.;
141 nodal_soln[11] = (elem_soln[0] + elem_soln[1] + elem_soln[3])/3.;
142 nodal_soln[12] = (elem_soln[1] + elem_soln[2] + elem_soln[3])/3.;
143 nodal_soln[13] = (elem_soln[0] + elem_soln[2] + elem_soln[3])/3.;
144 libmesh_fallthrough();
145 case TET10:
146 {
147 libmesh_assert_equal_to (elem_soln.size(), 4);
148 libmesh_assert (type == TET14 || nodal_soln.size() == 10);
149
150 nodal_soln[0] = elem_soln[0];
151 nodal_soln[1] = elem_soln[1];
152 nodal_soln[2] = elem_soln[2];
153 nodal_soln[3] = elem_soln[3];
154 nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]);
155 nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]);
156 nodal_soln[6] = .5*(elem_soln[2] + elem_soln[0]);
157 nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]);
158 nodal_soln[8] = .5*(elem_soln[3] + elem_soln[1]);
159 nodal_soln[9] = .5*(elem_soln[3] + elem_soln[2]);
160
161 return;
162 }
163
164
165 case HEX20:
166 case HEX27:
167 {
168 libmesh_assert_equal_to (elem_soln.size(), 8);
169
170 if (type == HEX20)
171 libmesh_assert_equal_to (nodal_soln.size(), 20);
172 else
173 libmesh_assert_equal_to (nodal_soln.size(), 27);
174
175 nodal_soln[0] = elem_soln[0];
176 nodal_soln[1] = elem_soln[1];
177 nodal_soln[2] = elem_soln[2];
178 nodal_soln[3] = elem_soln[3];
179 nodal_soln[4] = elem_soln[4];
180 nodal_soln[5] = elem_soln[5];
181 nodal_soln[6] = elem_soln[6];
182 nodal_soln[7] = elem_soln[7];
183 nodal_soln[8] = .5*(elem_soln[0] + elem_soln[1]);
184 nodal_soln[9] = .5*(elem_soln[1] + elem_soln[2]);
185 nodal_soln[10] = .5*(elem_soln[2] + elem_soln[3]);
186 nodal_soln[11] = .5*(elem_soln[3] + elem_soln[0]);
187 nodal_soln[12] = .5*(elem_soln[0] + elem_soln[4]);
188 nodal_soln[13] = .5*(elem_soln[1] + elem_soln[5]);
189 nodal_soln[14] = .5*(elem_soln[2] + elem_soln[6]);
190 nodal_soln[15] = .5*(elem_soln[3] + elem_soln[7]);
191 nodal_soln[16] = .5*(elem_soln[4] + elem_soln[5]);
192 nodal_soln[17] = .5*(elem_soln[5] + elem_soln[6]);
193 nodal_soln[18] = .5*(elem_soln[6] + elem_soln[7]);
194 nodal_soln[19] = .5*(elem_soln[4] + elem_soln[7]);
195
196 if (type == HEX27)
197 {
198 nodal_soln[20] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
199 nodal_soln[21] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[5]);
200 nodal_soln[22] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[6]);
201 nodal_soln[23] = .25*(elem_soln[2] + elem_soln[3] + elem_soln[6] + elem_soln[7]);
202 nodal_soln[24] = .25*(elem_soln[3] + elem_soln[0] + elem_soln[7] + elem_soln[4]);
203 nodal_soln[25] = .25*(elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]);
204
205 nodal_soln[26] = .125*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3] +
206 elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]);
207 }
208
209 return;
210 }
211
212
213 case PRISM21:
214 nodal_soln[20] = (elem_soln[9] + elem_soln[10] + elem_soln[11])/Real(3);
215 libmesh_fallthrough();
216 case PRISM20:
217 if (type == PRISM20)
218 libmesh_assert_equal_to (nodal_soln.size(), 20);
219 nodal_soln[18] = (elem_soln[0] + elem_soln[1] + elem_soln[2])/Real(3);
220 nodal_soln[19] = (elem_soln[3] + elem_soln[4] + elem_soln[5])/Real(3);
221 libmesh_fallthrough();
222 case PRISM18:
223 if (type == PRISM18)
224 libmesh_assert_equal_to (nodal_soln.size(), 18);
225 nodal_soln[15] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[3]);
226 nodal_soln[16] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[4]);
227 nodal_soln[17] = .25*(elem_soln[2] + elem_soln[0] + elem_soln[3] + elem_soln[5]);
228 libmesh_fallthrough();
229 case PRISM15:
230 {
231 libmesh_assert_equal_to (elem_soln.size(), 6);
232
233 if (type == PRISM15)
234 libmesh_assert_equal_to (nodal_soln.size(), 15);
235
236 nodal_soln[0] = elem_soln[0];
237 nodal_soln[1] = elem_soln[1];
238 nodal_soln[2] = elem_soln[2];
239 nodal_soln[3] = elem_soln[3];
240 nodal_soln[4] = elem_soln[4];
241 nodal_soln[5] = elem_soln[5];
242 nodal_soln[6] = .5*(elem_soln[0] + elem_soln[1]);
243 nodal_soln[7] = .5*(elem_soln[1] + elem_soln[2]);
244 nodal_soln[8] = .5*(elem_soln[0] + elem_soln[2]);
245 nodal_soln[9] = .5*(elem_soln[0] + elem_soln[3]);
246 nodal_soln[10] = .5*(elem_soln[1] + elem_soln[4]);
247 nodal_soln[11] = .5*(elem_soln[2] + elem_soln[5]);
248 nodal_soln[12] = .5*(elem_soln[3] + elem_soln[4]);
249 nodal_soln[13] = .5*(elem_soln[4] + elem_soln[5]);
250 nodal_soln[14] = .5*(elem_soln[3] + elem_soln[5]);
251
252 return;
253 }
254
255 case PYRAMID18:
256 {
257 libmesh_assert_equal_to (nodal_soln.size(), 18);
258
259 nodal_soln[14] = (elem_soln[0] + elem_soln[1] + elem_soln[4])/Real(3);
260 nodal_soln[15] = (elem_soln[1] + elem_soln[2] + elem_soln[4])/Real(3);
261 nodal_soln[16] = (elem_soln[2] + elem_soln[3] + elem_soln[4])/Real(3);
262 nodal_soln[17] = (elem_soln[0] + elem_soln[3] + elem_soln[4])/Real(3);
263
264 libmesh_fallthrough();
265 }
266
267 case PYRAMID14:
268 {
269 if (type == PYRAMID14)
270 libmesh_assert_equal_to (nodal_soln.size(), 14);
271
272 nodal_soln[13] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
273
274 libmesh_fallthrough();
275 }
276
277 case PYRAMID13:
278 {
279 libmesh_assert_equal_to (elem_soln.size(), 5);
280
281 if (type == PYRAMID13)
282 libmesh_assert_equal_to (nodal_soln.size(), 13);
283
284 nodal_soln[0] = elem_soln[0];
285 nodal_soln[1] = elem_soln[1];
286 nodal_soln[2] = elem_soln[2];
287 nodal_soln[3] = elem_soln[3];
288 nodal_soln[4] = elem_soln[4];
289 nodal_soln[5] = .5*(elem_soln[0] + elem_soln[1]);
290 nodal_soln[6] = .5*(elem_soln[1] + elem_soln[2]);
291 nodal_soln[7] = .5*(elem_soln[2] + elem_soln[3]);
292 nodal_soln[8] = .5*(elem_soln[3] + elem_soln[0]);
293 nodal_soln[9] = .5*(elem_soln[0] + elem_soln[4]);
294 nodal_soln[10] = .5*(elem_soln[1] + elem_soln[4]);
295 nodal_soln[11] = .5*(elem_soln[2] + elem_soln[4]);
296 nodal_soln[12] = .5*(elem_soln[3] + elem_soln[4]);
297
298 return;
299 }
300 default:
301 {
302 // By default the element solution _is_ nodal,
303 // so just copy it.
304 nodal_soln = elem_soln;
305
306 return;
307 }
308 }
309 }
310
311 case SECOND:
312 {
313 switch (type)
314 {
315 case EDGE4:
316 {
317 libmesh_assert_equal_to (elem_soln.size(), 3);
318 libmesh_assert_equal_to (nodal_soln.size(), 4);
319
320 // Project quadratic solution onto cubic element nodes
321 nodal_soln[0] = elem_soln[0];
322 nodal_soln[1] = elem_soln[1];
323 nodal_soln[2] = (2.*elem_soln[0] - elem_soln[1] +
324 8.*elem_soln[2])/9.;
325 nodal_soln[3] = (-elem_soln[0] + 2.*elem_soln[1] +
326 8.*elem_soln[2])/9.;
327 return;
328 }
329
330 case TRI7:
331 {
332 libmesh_assert_equal_to (elem_soln.size(), 6);
333 libmesh_assert_equal_to (nodal_soln.size(), 7);
334
335 for (int i=0; i != 6; ++i)
336 nodal_soln[i] = elem_soln[i];
337
338 nodal_soln[6] = -1./9. * (elem_soln[0] + elem_soln[1] + elem_soln[2])
339 +4./9. * (elem_soln[3] + elem_soln[4] + elem_soln[5]);
340
341 return;
342 }
343
344 case TET14:
345 {
346 libmesh_assert_equal_to (elem_soln.size(), 10);
347 libmesh_assert_equal_to (nodal_soln.size(), 14);
348
349 for (int i=0; i != 10; ++i)
350 nodal_soln[i] = elem_soln[i];
351
352 nodal_soln[10] = -1./9. * (elem_soln[0] + elem_soln[1] + elem_soln[2])
353 +4./9. * (elem_soln[4] + elem_soln[5] + elem_soln[6]);
354 nodal_soln[11] = -1./9. * (elem_soln[0] + elem_soln[1] + elem_soln[3])
355 +4./9. * (elem_soln[4] + elem_soln[7] + elem_soln[8]);
356 nodal_soln[12] = -1./9. * (elem_soln[1] + elem_soln[2] + elem_soln[3])
357 +4./9. * (elem_soln[5] + elem_soln[8] + elem_soln[9]);
358 nodal_soln[13] = -1./9. * (elem_soln[0] + elem_soln[2] + elem_soln[3])
359 +4./9. * (elem_soln[6] + elem_soln[7] + elem_soln[9]);
360
361 return;
362 }
363
364 case PRISM21:
365 {
366 nodal_soln[20] = (elem_soln[9] + elem_soln[10] + elem_soln[11])/Real(3);
367 libmesh_fallthrough();
368 }
369 case PRISM20:
370 {
371 if (type == PRISM20)
372 libmesh_assert_equal_to (nodal_soln.size(), 20);
373
374 for (int i=0; i != 18; ++i)
375 nodal_soln[i] = elem_soln[i];
376
377 nodal_soln[18] = (elem_soln[0] + elem_soln[1] + elem_soln[2])/Real(3);
378 nodal_soln[19] = (elem_soln[3] + elem_soln[4] + elem_soln[5])/Real(3);
379 return;
380 }
381
382 case PYRAMID18:
383 {
384 libmesh_assert_equal_to (nodal_soln.size(), 18);
385
386 for (int i=0; i != 14; ++i)
387 nodal_soln[i] = elem_soln[i];
388
389 nodal_soln[14] = (elem_soln[0] + elem_soln[1] + elem_soln[4])/Real(3);
390 nodal_soln[15] = (elem_soln[1] + elem_soln[2] + elem_soln[4])/Real(3);
391 nodal_soln[16] = (elem_soln[2] + elem_soln[3] + elem_soln[4])/Real(3);
392 nodal_soln[17] = (elem_soln[0] + elem_soln[3] + elem_soln[4])/Real(3);
393
394 libmesh_fallthrough();
395 }
396
397 default:
398 {
399 // By default the element solution _is_ nodal, so just
400 // copy the portion relevant to the nodal solution.
401 // (this is the whole nodal solution for true Lagrange
402 // elements, but a smaller part for "L2 Lagrange"
403 libmesh_assert_less_equal(nodal_soln.size(), elem_soln.size());
404 for (const auto i : index_range(nodal_soln))
405 nodal_soln[i] = elem_soln[i];
406
407 return;
408 }
409 }
410 }
411
412
413
414
415 default:
416 {
417 // By default the element solution _is_ nodal, so just copy
418 // the portion relevant to the nodal solution. (this is the
419 // whole nodal solution for true Lagrange elements, but a
420 // smaller part for "L2 Lagrange"
421 libmesh_assert_less_equal(nodal_soln.size(), elem_soln.size());
422 for (const auto i : index_range(nodal_soln))
423 nodal_soln[i] = elem_soln[i];
424
425 return;
426 }
427 }
428}
429
430
431// Anonymous namespace for local helper functions
432namespace {
433unsigned int lagrange_n_dofs(const ElemType t, const Elem * e, const Order o)
434{
435 libmesh_assert(!e || e->type() == t);
436
437 switch (o)
438 {
439 // lagrange can only be constant on a single node
440 case CONSTANT:
441 {
442 switch (t)
443 {
444 case NODEELEM:
445 return 1;
446
447 default:
448 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
449 }
450 }
451
452 // linear Lagrange shape functions
453 case FIRST:
454 {
455 switch (t)
456 {
457 case NODEELEM:
458 return 1;
459
460 case EDGE2:
461 case EDGE3:
462 case EDGE4:
463 return 2;
464
465 case TRI3:
466 case TRISHELL3:
467 case TRI3SUBDIVISION:
468 case TRI6:
469 case TRI7:
470 return 3;
471
472 case QUAD4:
473 case QUADSHELL4:
474 case QUAD8:
475 case QUADSHELL8:
476 case QUAD9:
477 case QUADSHELL9:
478 return 4;
479
480 case TET4:
481 case TET10:
482 case TET14:
483 return 4;
484
485 case HEX8:
486 case HEX20:
487 case HEX27:
488 return 8;
489
490 case PRISM6:
491 case PRISM15:
492 case PRISM18:
493 case PRISM20:
494 case PRISM21:
495 return 6;
496
497 case PYRAMID5:
498 case PYRAMID13:
499 case PYRAMID14:
500 case PYRAMID18:
501 return 5;
502
503 case INVALID_ELEM:
504 return 0;
505
506 case C0POLYGON:
507 case C0POLYHEDRON:
508 // Polygons and polyhedra require using newer FE APIs
509 if (!e)
510 libmesh_error_msg("Code (see stack trace) used an outdated FE function overload.\n"
511 "n_dofs() on a polygon or polyhedron is not defined by ElemType alone.");
512 return e->n_nodes();
513
514 default:
515 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
516 }
517 }
518
519
520 // quadratic Lagrange shape functions
521 case SECOND:
522 {
523 switch (t)
524 {
525 case NODEELEM:
526 return 1;
527
528 case EDGE3:
529 return 3;
530
531 case TRI6:
532 case TRI7:
533 return 6;
534
535 case QUAD8:
536 case QUADSHELL8:
537 return 8;
538
539 case QUAD9:
540 case QUADSHELL9:
541 return 9;
542
543 case TET10:
544 case TET14:
545 return 10;
546
547 case HEX20:
548 return 20;
549
550 case HEX27:
551 return 27;
552
553 case PRISM15:
554 return 15;
555
556 case PRISM18:
557 case PRISM20:
558 case PRISM21:
559 return 18;
560
561 case PYRAMID13:
562 return 13;
563
564 case PYRAMID14:
565 case PYRAMID18:
566 return 14;
567
568 case INVALID_ELEM:
569 return 0;
570
571 default:
572 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
573 }
574 }
575
576 case THIRD:
577 {
578 switch (t)
579 {
580 case NODEELEM:
581 return 1;
582
583 case EDGE4:
584 return 4;
585
586 case PRISM20:
587 return 20;
588
589 case PRISM21:
590 return 21;
591
592 case PYRAMID18:
593 return 18;
594
595 case TRI7:
596 return 7;
597
598 case TET14:
599 return 14;
600
601 case INVALID_ELEM:
602 return 0;
603
604 default:
605 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
606 }
607 }
608
609 default:
610 libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for LAGRANGE FE family!");
611 }
612}
613
614
615
616unsigned int lagrange_n_dofs_at_node(const ElemType t,
617 const Order o,
618 const unsigned int n)
619{
620 switch (o)
621 {
622 // lagrange can only be constant on a single node
623 case CONSTANT:
624 {
625 switch (t)
626 {
627 case NODEELEM:
628 return 1;
629
630 default:
631 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
632 }
633 }
634
635 // linear Lagrange shape functions
636 case FIRST:
637 {
638 switch (t)
639 {
640 case NODEELEM:
641 return 1;
642
643 case EDGE2:
644 case EDGE3:
645 case EDGE4:
646 {
647 switch (n)
648 {
649 case 0:
650 case 1:
651 return 1;
652
653 default:
654 return 0;
655 }
656 }
657
658 case TRI3:
659 case TRISHELL3:
660 case TRI3SUBDIVISION:
661 case TRI6:
662 case TRI7:
663 {
664 switch (n)
665 {
666 case 0:
667 case 1:
668 case 2:
669 return 1;
670
671 default:
672 return 0;
673 }
674 }
675
676 case QUAD4:
677 case QUADSHELL4:
678 case QUAD8:
679 case QUADSHELL8:
680 case QUAD9:
681 case QUADSHELL9:
682 {
683 switch (n)
684 {
685 case 0:
686 case 1:
687 case 2:
688 case 3:
689 return 1;
690
691 default:
692 return 0;
693 }
694 }
695
696 case C0POLYGON:
697 case C0POLYHEDRON:
698 return 1;
699
700
701 case TET4:
702 case TET10:
703 case TET14:
704 {
705 switch (n)
706 {
707 case 0:
708 case 1:
709 case 2:
710 case 3:
711 return 1;
712
713 default:
714 return 0;
715 }
716 }
717
718 case HEX8:
719 case HEX20:
720 case HEX27:
721 {
722 switch (n)
723 {
724 case 0:
725 case 1:
726 case 2:
727 case 3:
728 case 4:
729 case 5:
730 case 6:
731 case 7:
732 return 1;
733
734 default:
735 return 0;
736 }
737 }
738
739 case PRISM6:
740 case PRISM15:
741 case PRISM18:
742 case PRISM20:
743 case PRISM21:
744 {
745 switch (n)
746 {
747 case 0:
748 case 1:
749 case 2:
750 case 3:
751 case 4:
752 case 5:
753 return 1;
754
755 default:
756 return 0;
757 }
758 }
759
760 case PYRAMID5:
761 case PYRAMID13:
762 case PYRAMID14:
763 case PYRAMID18:
764 {
765 switch (n)
766 {
767 case 0:
768 case 1:
769 case 2:
770 case 3:
771 case 4:
772 return 1;
773
774 default:
775 return 0;
776 }
777 }
778
779 case INVALID_ELEM:
780 return 0;
781
782 default:
783 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
784 }
785 }
786
787 // quadratic Lagrange shape functions
788 case SECOND:
789 {
790 switch (t)
791 {
792 // quadratic lagrange has one dof at each node
793 case NODEELEM:
794 case EDGE3:
795 case TRI6:
796 case QUAD8:
797 case QUADSHELL8:
798 case QUAD9:
799 case QUADSHELL9:
800 case TET10:
801 case HEX20:
802 case HEX27:
803 case PRISM15:
804 case PRISM18:
805 case PYRAMID13:
806 case PYRAMID14:
807 return 1;
808
809 case PRISM20:
810 case PRISM21:
811 return (n < 18);
812
813 case PYRAMID18:
814 return (n < 14);
815
816 case TRI7:
817 return (n < 6);
818
819 case TET14:
820 return (n < 10);
821
822 case INVALID_ELEM:
823 return 0;
824
825 default:
826 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
827 }
828 }
829
830 case THIRD:
831 {
832 switch (t)
833 {
834 case NODEELEM:
835 case EDGE4:
836 case PRISM20:
837 case PRISM21:
838 case PYRAMID18:
839 case TRI7:
840 case TET14:
841 return 1;
842
843 case INVALID_ELEM:
844 return 0;
845
846 default:
847 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
848 }
849 }
850
851 default:
852 libmesh_error_msg("Unsupported order: " << o );
853 }
854}
855
856
857
858unsigned int lagrange_n_dofs_at_node(const Elem & e,
859 const Order o,
860 const unsigned int n)
861{
862 return lagrange_n_dofs_at_node(e.type(), o, n);
863}
864
865
866#ifdef LIBMESH_ENABLE_AMR
867void lagrange_compute_constraints (DofConstraints & constraints,
868 DofMap & dof_map,
869 const unsigned int variable_number,
870 const Elem * elem,
871 const unsigned Dim)
872{
873 // Only constrain elements in 2,3D.
874 if (Dim == 1)
875 return;
876
877 libmesh_assert(elem);
878
879 // Only constrain active and ancestor elements
880 if (elem->subactive())
881 return;
882
883#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
884 if (elem->infinite())
885 {
886 const FEType fe_t = dof_map.variable_type(variable_number);
887
888 // expand the infinite_compute_constraint in its template-arguments.
889 switch(Dim)
890 {
891 case 2:
892 {
893 inf_fe_family_mapping_switch(2, inf_compute_constraints (constraints, dof_map, variable_number, elem) , ,; break;);
894 break;
895 }
896 case 3:
897 {
898 inf_fe_family_mapping_switch(3, inf_compute_constraints (constraints, dof_map, variable_number, elem) , ,; break;);
899 break;
900 }
901 default:
902 libmesh_error_msg("Invalid dim = " << Dim);
903 }
904 return;
905 }
906#endif
907
908 const Variable & var = dof_map.variable(variable_number);
909 const FEType fe_type = var.type();
910 const bool add_p_level = fe_type.p_refinement;
911
912 // Pull objects out of the loop to reduce heap operations
913 std::vector<dof_id_type> my_dof_indices, parent_dof_indices;
914 std::unique_ptr<const Elem> my_side, parent_side;
915
916 // Look at the element faces. Check to see if we need to
917 // build constraints.
918 for (auto s : elem->side_index_range())
919 if (const Elem * const neigh = elem->neighbor_ptr(s);
920 neigh && neigh != remote_elem)
921 {
922 // constrain dofs shared between this element and ones coarser
923 if (neigh->level() < elem->level() &&
924 var.active_on_subdomain(neigh->subdomain_id()))
925 {
926 // than this element.
927 // Get pointers to the elements of interest and its parent.
928 const Elem * parent = elem->parent();
929
930 // This can't happen... Only level-0 elements have nullptr
931 // parents, and no level-0 elements can be at a higher
932 // level than their neighbors!
933 libmesh_assert(parent);
934
935 elem->build_side_ptr(my_side, s);
936 parent->build_side_ptr(parent_side, s);
937
938 // We have some elements with interior bubble function bases
939 // and lower-order sides. We need to only ask for the lower
940 // order in those cases.
941 FEType side_fe_type = fe_type;
942 int extra_order = 0;
943 if (my_side->default_order() < fe_type.order)
944 {
945 side_fe_type.order = my_side->default_order();
946 extra_order = (int)(side_fe_type.order) -
947 (int)(fe_type.order);
948 }
949
950 // This function gets called element-by-element, so there
951 // will be a lot of memory allocation going on. We can
952 // at least minimize this for the case of the dof indices
953 // by efficiently preallocating the requisite storage.
954 my_dof_indices.reserve (my_side->n_nodes());
955 parent_dof_indices.reserve (parent_side->n_nodes());
956
957 dof_map.dof_indices (my_side.get(), my_dof_indices,
958 variable_number, extra_order);
959 dof_map.dof_indices (parent_side.get(), parent_dof_indices,
960 variable_number, extra_order);
961
962 const unsigned int n_side_dofs =
963 FEInterface::n_dofs(side_fe_type, my_side.get());
964 const unsigned int n_parent_side_dofs =
965 FEInterface::n_dofs(side_fe_type, parent_side.get());
966 for (unsigned int my_dof=0; my_dof != n_side_dofs; my_dof++)
967 {
968 libmesh_assert_less (my_dof, my_side->n_nodes());
969
970 // My global dof index.
971 const dof_id_type my_dof_g = my_dof_indices[my_dof];
972
973 // Hunt for "constraining against myself" cases before
974 // we bother creating a constraint row
975 bool self_constraint = false;
976 for (unsigned int their_dof=0;
977 their_dof != n_parent_side_dofs; their_dof++)
978 {
979 libmesh_assert_less (their_dof, parent_side->n_nodes());
980
981 // Their global dof index.
982 const dof_id_type their_dof_g =
983 parent_dof_indices[their_dof];
984
985 if (their_dof_g == my_dof_g)
986 {
987 self_constraint = true;
988 break;
989 }
990 }
991
992 if (self_constraint)
993 continue;
994
995 DofConstraintRow * constraint_row;
996
997 // we may be running constraint methods concurrently
998 // on multiple threads, so we need a lock to
999 // ensure that this constraint is "ours"
1000 {
1001 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1002
1003 if (dof_map.is_constrained_dof(my_dof_g))
1004 continue;
1005
1006 constraint_row = &(constraints[my_dof_g]);
1007 libmesh_assert(constraint_row->empty());
1008 }
1009
1010 // The support point of the DOF
1011 const Point & support_point = my_side->point(my_dof);
1012
1013 // Figure out where my node lies on their reference element.
1014 const Point mapped_point = FEMap::inverse_map(Dim-1,
1015 parent_side.get(),
1016 support_point);
1017
1018 // Compute the parent's side shape function values.
1019 for (unsigned int their_dof=0;
1020 their_dof != n_parent_side_dofs; their_dof++)
1021 {
1022 libmesh_assert_less (their_dof, parent_side->n_nodes());
1023
1024 // Their global dof index.
1025 const dof_id_type their_dof_g =
1026 parent_dof_indices[their_dof];
1027
1028 const Real their_dof_value = FEInterface::shape(side_fe_type,
1029 parent_side.get(),
1030 their_dof,
1031 mapped_point);
1032
1033 // Only add non-zero and non-identity values
1034 // for Lagrange basis functions.
1035 if ((std::abs(their_dof_value) > 1.e-5) &&
1036 (std::abs(their_dof_value) < .999))
1037 {
1038 constraint_row->emplace(their_dof_g, their_dof_value);
1039 }
1040 #ifdef DEBUG
1041 // Protect for the case u_i = 0.999 u_j,
1042 // in which case i better equal j.
1043 else if (their_dof_value >= .999)
1044 libmesh_assert_equal_to (my_dof_g, their_dof_g);
1045 #endif
1046 }
1047 }
1048 }
1049
1050 if (elem->active() && add_p_level)
1051 {
1052 // p refinement constraints:
1053 // constrain dofs shared between
1054 // active elements and neighbors with
1055 // lower polynomial degrees
1056 const unsigned int min_p_level =
1057 neigh->min_p_level_by_neighbor(elem, elem->p_level());
1058 if (min_p_level < elem->p_level())
1059 libmesh_error_msg(
1060 "Mismatched p-refinement levels for LAGRANGE is not currently supported");
1061 }
1062 }
1063} // lagrange_compute_constraints()
1064#endif // #ifdef LIBMESH_ENABLE_AMR
1065
1066} // anonymous namespace
1067
1068
1069// Instantiate (side_) nodal_soln() function for every dimension
1070LIBMESH_FE_NODAL_SOLN(LAGRANGE, lagrange_nodal_soln)
1071LIBMESH_FE_SIDE_NODAL_SOLN(LAGRANGE)
1072
1073
1074// Do full-specialization for every dimension, instead
1075// of explicit instantiation at the end of this function.
1076// This could be macro-ified.
1077template <> unsigned int FE<0,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1078template <> unsigned int FE<1,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1079template <> unsigned int FE<2,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1080template <> unsigned int FE<3,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1081
1082template <> unsigned int FE<0,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1083template <> unsigned int FE<1,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1084template <> unsigned int FE<2,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1085template <> unsigned int FE<3,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1086
1087
1088// Do full-specialization for every dimension, instead
1089// of explicit instantiation at the end of this function.
1090template <> unsigned int FE<0,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); }
1091template <> unsigned int FE<1,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); }
1092template <> unsigned int FE<2,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); }
1093template <> unsigned int FE<3,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); }
1094
1095template <> unsigned int FE<0,LAGRANGE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(e, o, n); }
1096template <> unsigned int FE<1,LAGRANGE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(e, o, n); }
1097template <> unsigned int FE<2,LAGRANGE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(e, o, n); }
1098template <> unsigned int FE<3,LAGRANGE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(e, o, n); }
1099
1100
1101// Lagrange elements have no dofs per element
1102// (just at the nodes)
1103template <> unsigned int FE<0,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1104template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1105template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1106template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1107
1108template <> unsigned int FE<0,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1109template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1110template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1111template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1112
1113// Lagrange FEMs are always C^0 continuous
1114template <> FEContinuity FE<0,LAGRANGE>::get_continuity() const { return C_ZERO; }
1115template <> FEContinuity FE<1,LAGRANGE>::get_continuity() const { return C_ZERO; }
1116template <> FEContinuity FE<2,LAGRANGE>::get_continuity() const { return C_ZERO; }
1117template <> FEContinuity FE<3,LAGRANGE>::get_continuity() const { return C_ZERO; }
1118
1119// Lagrange FEMs are not hierarchic
1120template <> bool FE<0,LAGRANGE>::is_hierarchic() const { return false; }
1121template <> bool FE<1,LAGRANGE>::is_hierarchic() const { return false; }
1122template <> bool FE<2,LAGRANGE>::is_hierarchic() const { return false; }
1123template <> bool FE<3,LAGRANGE>::is_hierarchic() const { return false; }
1124
1125// Lagrange FEM shapes do not need reinit (is this always true?)
1126template <> bool FE<0,LAGRANGE>::shapes_need_reinit() const { return false; }
1127template <> bool FE<1,LAGRANGE>::shapes_need_reinit() const { return false; }
1128template <> bool FE<2,LAGRANGE>::shapes_need_reinit() const { return false; }
1129template <> bool FE<3,LAGRANGE>::shapes_need_reinit() const { return false; }
1130
1131// Methods for computing Lagrange constraints. Note: we pass the
1132// dimension as the last argument to the anonymous helper function.
1133// Also note: we only need instantiations of this function for
1134// Dim==2 and 3.
1135#ifdef LIBMESH_ENABLE_AMR
1136template <>
1137void FE<2,LAGRANGE>::compute_constraints (DofConstraints & constraints,
1138 DofMap & dof_map,
1139 const unsigned int variable_number,
1140 const Elem * elem)
1141{ lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/2); }
1142
1143template <>
1144void FE<3,LAGRANGE>::compute_constraints (DofConstraints & constraints,
1145 DofMap & dof_map,
1146 const unsigned int variable_number,
1147 const Elem * elem)
1148{ lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/3); }
1149#endif // LIBMESH_ENABLE_AMR
1150
1151} // namespace libMesh
int
void ErrorVector unsigned int
Definition adjoints_ex3.C:360
libMesh::DofConstraints
The constraint matrix storage format.
Definition dof_map.h:112
libMesh::DofMap
This class handles the numbering of degrees of freedom on a mesh.
Definition dof_map.h:181
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::Elem::n_nodes
virtual unsigned int n_nodes() const =0
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Elem::p_level
unsigned int p_level() const
Definition elem.h:3125
libMesh::FEInterface::shape
static Real shape(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const Point &p)
Definition fe_interface.C:761
libMesh::FEInterface::n_dofs
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition fe_interface.C:353
libMesh::FEMap::inverse_map
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true, const bool extra_checks=true)
Definition fe_map.C:1512
libMesh::FE::is_hierarchic
virtual bool is_hierarchic() const override
libMesh::FE::n_dofs_per_elem
static unsigned int n_dofs_per_elem(const ElemType t, const Order o)
libMesh::FE::get_continuity
virtual FEContinuity get_continuity() const override
libMesh::FE::compute_constraints
static void compute_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
Computes the constraint matrix contributions (for non-conforming adapted meshes) corresponding to var...
libMesh::FE::n_dofs
static unsigned int n_dofs(const ElemType t, const Order o)
libMesh::FE::n_dofs_at_node
static unsigned int n_dofs_at_node(const ElemType t, const Order o, const unsigned int n)
libMesh::FE::shapes_need_reinit
virtual bool shapes_need_reinit() const override
libMesh::Threads::spin_mtx
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition threads.C:30
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::LIBMESH_FE_NODAL_SOLN
LIBMESH_FE_NODAL_SOLN(BERNSTEIN, bernstein_nodal_soln)
Definition fe_bernstein.C:415
libMesh::index_range
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition int_range.h:153
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::lagrange_nodal_soln
void lagrange_nodal_soln(const Elem *elem, const Order order, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln, bool add_p_level=true)
Helper functions for Lagrange-based basis functions.
Definition fe_lagrange.C:43
libMesh::remote_elem
const RemoteElem * remote_elem
Definition remote_elem.C:57
libMesh::DofConstraintRow
std::map< dof_id_type, Real, std::less< dof_id_type >, Threads::scalable_allocator< std::pair< const dof_id_type, Real > > > DofConstraintRow
A row of the Dof constraint matrix.
Definition dof_map.h:100
libMesh::dof_id_type
uint8_t dof_id_type
Definition id_types.h:67
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
libMesh::LIBMESH_FE_SIDE_NODAL_SOLN
LIBMESH_FE_SIDE_NODAL_SOLN(HIERARCHIC_VEC)
Definition fe_hierarchic_vec.C:117
n_nodes
const dof_id_type n_nodes
Definition tecplot_io.C:67
Generated on Fri Aug 21 2026 13:32:28 for libMesh by doxygen 1.9.8