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] = -1./9. * (elem_soln[9] + elem_soln[10] + elem_soln[11])
367 +4./9. * (elem_soln[15] + elem_soln[16] + elem_soln[17]);
368 libmesh_fallthrough();
369 }
370 case PRISM20:
371 {
372 libmesh_assert_equal_to (elem_soln.size(), 18);
373 if (type == PRISM20)
374 libmesh_assert_equal_to (nodal_soln.size(), 20);
375
376 for (int i=0; i != 18; ++i)
377 nodal_soln[i] = elem_soln[i];
378
379 nodal_soln[18] = -1./9. * (elem_soln[0] + elem_soln[1] + elem_soln[2])
380 +4./9. * (elem_soln[6] + elem_soln[7] + elem_soln[8]);
381 nodal_soln[19] = -1./9. * (elem_soln[3] + elem_soln[4] + elem_soln[5])
382 +4./9. * (elem_soln[12] + elem_soln[13] + elem_soln[14]);
383 return;
384 }
385
386 case PYRAMID18:
387 {
388 libmesh_assert_equal_to (elem_soln.size(), 14);
389 libmesh_assert_equal_to (nodal_soln.size(), 18);
390
391 for (int i=0; i != 14; ++i)
392 nodal_soln[i] = elem_soln[i];
393
394 nodal_soln[14] = -1./9. * (elem_soln[0] + elem_soln[1] + elem_soln[4])
395 +4./9. * (elem_soln[5] + elem_soln[9] + elem_soln[10]);
396 nodal_soln[15] = -1./9. * (elem_soln[1] + elem_soln[2] + elem_soln[4])
397 +4./9. * (elem_soln[6] + elem_soln[10] + elem_soln[11]);
398 nodal_soln[16] = -1./9. * (elem_soln[2] + elem_soln[3] + elem_soln[4])
399 +4./9. * (elem_soln[7] + elem_soln[11] + elem_soln[12]);
400 nodal_soln[17] = -1./9. * (elem_soln[0] + elem_soln[3] + elem_soln[4])
401 +4./9. * (elem_soln[8] + elem_soln[9] + elem_soln[12]);
402
403 return;
404 }
405
406 default:
407 {
408 // By default the element solution _is_ nodal, so just
409 // copy the portion relevant to the nodal solution.
410 // (this is the whole nodal solution for true Lagrange
411 // elements, but a smaller part for "L2 Lagrange"
412 libmesh_assert_less_equal(nodal_soln.size(), elem_soln.size());
413 for (const auto i : index_range(nodal_soln))
414 nodal_soln[i] = elem_soln[i];
415
416 return;
417 }
418 }
419 }
420
421
422
423
424 default:
425 {
426 // By default the element solution _is_ nodal, so just copy
427 // the portion relevant to the nodal solution. (this is the
428 // whole nodal solution for true Lagrange elements, but a
429 // smaller part for "L2 Lagrange"
430 libmesh_assert_less_equal(nodal_soln.size(), elem_soln.size());
431 for (const auto i : index_range(nodal_soln))
432 nodal_soln[i] = elem_soln[i];
433
434 return;
435 }
436 }
437}
438
439
440// Anonymous namespace for local helper functions
441namespace {
442unsigned int lagrange_n_dofs(const ElemType t, const Elem * e, const Order o)
443{
444 libmesh_assert(!e || e->type() == t);
445
446 switch (o)
447 {
448 // lagrange can only be constant on a single node
449 case CONSTANT:
450 {
451 switch (t)
452 {
453 case NODEELEM:
454 return 1;
455
456 default:
457 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
458 }
459 }
460
461 // linear Lagrange shape functions
462 case FIRST:
463 {
464 switch (t)
465 {
466 case NODEELEM:
467 return 1;
468
469 case EDGE2:
470 case EDGE3:
471 case EDGE4:
472 return 2;
473
474 case TRI3:
475 case TRISHELL3:
476 case TRI3SUBDIVISION:
477 case TRI6:
478 case TRI7:
479 return 3;
480
481 case QUAD4:
482 case QUADSHELL4:
483 case QUAD8:
484 case QUADSHELL8:
485 case QUAD9:
486 case QUADSHELL9:
487 return 4;
488
489 case TET4:
490 case TET10:
491 case TET14:
492 return 4;
493
494 case HEX8:
495 case HEX20:
496 case HEX27:
497 return 8;
498
499 case PRISM6:
500 case PRISM15:
501 case PRISM18:
502 case PRISM20:
503 case PRISM21:
504 return 6;
505
506 case PYRAMID5:
507 case PYRAMID13:
508 case PYRAMID14:
509 case PYRAMID18:
510 return 5;
511
512 case INVALID_ELEM:
513 return 0;
514
515 case C0POLYGON:
516 case C0POLYHEDRON:
517 // Polygons and polyhedra require using newer FE APIs
518 if (!e)
519 libmesh_error_msg("Code (see stack trace) used an outdated FE function overload.\n"
520 "n_dofs() on a polygon or polyhedron is not defined by ElemType alone.");
521 return e->n_nodes();
522
523 default:
524 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
525 }
526 }
527
528
529 // quadratic Lagrange shape functions
530 case SECOND:
531 {
532 switch (t)
533 {
534 case NODEELEM:
535 return 1;
536
537 case EDGE3:
538 return 3;
539
540 case TRI6:
541 case TRI7:
542 return 6;
543
544 case QUAD8:
545 case QUADSHELL8:
546 return 8;
547
548 case QUAD9:
549 case QUADSHELL9:
550 return 9;
551
552 case TET10:
553 case TET14:
554 return 10;
555
556 case HEX20:
557 return 20;
558
559 case HEX27:
560 return 27;
561
562 case PRISM15:
563 return 15;
564
565 case PRISM18:
566 case PRISM20:
567 case PRISM21:
568 return 18;
569
570 case PYRAMID13:
571 return 13;
572
573 case PYRAMID14:
574 case PYRAMID18:
575 return 14;
576
577 case INVALID_ELEM:
578 return 0;
579
580 default:
581 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
582 }
583 }
584
585 case THIRD:
586 {
587 switch (t)
588 {
589 case NODEELEM:
590 return 1;
591
592 case EDGE4:
593 return 4;
594
595 case PRISM20:
596 return 20;
597
598 case PRISM21:
599 return 21;
600
601 case PYRAMID18:
602 return 18;
603
604 case TRI7:
605 return 7;
606
607 case TET14:
608 return 14;
609
610 case INVALID_ELEM:
611 return 0;
612
613 default:
614 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
615 }
616 }
617
618 default:
619 libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for LAGRANGE FE family!");
620 }
621}
622
623
624
625unsigned int lagrange_n_dofs_at_node(const ElemType t,
626 const Order o,
627 const unsigned int n)
628{
629 switch (o)
630 {
631 // lagrange can only be constant on a single node
632 case CONSTANT:
633 {
634 switch (t)
635 {
636 case NODEELEM:
637 return 1;
638
639 default:
640 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
641 }
642 }
643
644 // linear Lagrange shape functions
645 case FIRST:
646 {
647 switch (t)
648 {
649 case NODEELEM:
650 return 1;
651
652 case EDGE2:
653 case EDGE3:
654 case EDGE4:
655 {
656 switch (n)
657 {
658 case 0:
659 case 1:
660 return 1;
661
662 default:
663 return 0;
664 }
665 }
666
667 case TRI3:
668 case TRISHELL3:
669 case TRI3SUBDIVISION:
670 case TRI6:
671 case TRI7:
672 {
673 switch (n)
674 {
675 case 0:
676 case 1:
677 case 2:
678 return 1;
679
680 default:
681 return 0;
682 }
683 }
684
685 case QUAD4:
686 case QUADSHELL4:
687 case QUAD8:
688 case QUADSHELL8:
689 case QUAD9:
690 case QUADSHELL9:
691 {
692 switch (n)
693 {
694 case 0:
695 case 1:
696 case 2:
697 case 3:
698 return 1;
699
700 default:
701 return 0;
702 }
703 }
704
705 case C0POLYGON:
706 case C0POLYHEDRON:
707 return 1;
708
709
710 case TET4:
711 case TET10:
712 case TET14:
713 {
714 switch (n)
715 {
716 case 0:
717 case 1:
718 case 2:
719 case 3:
720 return 1;
721
722 default:
723 return 0;
724 }
725 }
726
727 case HEX8:
728 case HEX20:
729 case HEX27:
730 {
731 switch (n)
732 {
733 case 0:
734 case 1:
735 case 2:
736 case 3:
737 case 4:
738 case 5:
739 case 6:
740 case 7:
741 return 1;
742
743 default:
744 return 0;
745 }
746 }
747
748 case PRISM6:
749 case PRISM15:
750 case PRISM18:
751 case PRISM20:
752 case PRISM21:
753 {
754 switch (n)
755 {
756 case 0:
757 case 1:
758 case 2:
759 case 3:
760 case 4:
761 case 5:
762 return 1;
763
764 default:
765 return 0;
766 }
767 }
768
769 case PYRAMID5:
770 case PYRAMID13:
771 case PYRAMID14:
772 case PYRAMID18:
773 {
774 switch (n)
775 {
776 case 0:
777 case 1:
778 case 2:
779 case 3:
780 case 4:
781 return 1;
782
783 default:
784 return 0;
785 }
786 }
787
788 case INVALID_ELEM:
789 return 0;
790
791 default:
792 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
793 }
794 }
795
796 // quadratic Lagrange shape functions
797 case SECOND:
798 {
799 switch (t)
800 {
801 // quadratic lagrange has one dof at each node
802 case NODEELEM:
803 case EDGE3:
804 case TRI6:
805 case QUAD8:
806 case QUADSHELL8:
807 case QUAD9:
808 case QUADSHELL9:
809 case TET10:
810 case HEX20:
811 case HEX27:
812 case PRISM15:
813 case PRISM18:
814 case PYRAMID13:
815 case PYRAMID14:
816 return 1;
817
818 case PRISM20:
819 case PRISM21:
820 return (n < 18);
821
822 case PYRAMID18:
823 return (n < 14);
824
825 case TRI7:
826 return (n < 6);
827
828 case TET14:
829 return (n < 10);
830
831 case INVALID_ELEM:
832 return 0;
833
834 default:
835 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
836 }
837 }
838
839 case THIRD:
840 {
841 switch (t)
842 {
843 case NODEELEM:
844 case EDGE4:
845 case PRISM20:
846 case PRISM21:
847 case PYRAMID18:
848 case TRI7:
849 case TET14:
850 return 1;
851
852 case INVALID_ELEM:
853 return 0;
854
855 default:
856 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
857 }
858 }
859
860 default:
861 libmesh_error_msg("Unsupported order: " << o );
862 }
863}
864
865
866
867unsigned int lagrange_n_dofs_at_node(const Elem & e,
868 const Order o,
869 const unsigned int n)
870{
871 return lagrange_n_dofs_at_node(e.type(), o, n);
872}
873
874
875#ifdef LIBMESH_ENABLE_AMR
876void lagrange_compute_constraints (DofConstraints & constraints,
877 DofMap & dof_map,
878 const unsigned int variable_number,
879 const Elem * elem,
880 const unsigned Dim)
881{
882 // Only constrain elements in 2,3D.
883 if (Dim == 1)
884 return;
885
886 libmesh_assert(elem);
887
888 // Only constrain active and ancestor elements
889 if (elem->subactive())
890 return;
891
892#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
893 if (elem->infinite())
894 {
895 const FEType fe_t = dof_map.variable_type(variable_number);
896
897 // expand the infinite_compute_constraint in its template-arguments.
898 switch(Dim)
899 {
900 case 2:
901 {
902 inf_fe_family_mapping_switch(2, inf_compute_constraints (constraints, dof_map, variable_number, elem) , ,; break;);
903 break;
904 }
905 case 3:
906 {
907 inf_fe_family_mapping_switch(3, inf_compute_constraints (constraints, dof_map, variable_number, elem) , ,; break;);
908 break;
909 }
910 default:
911 libmesh_error_msg("Invalid dim = " << Dim);
912 }
913 return;
914 }
915#endif
916
917 const Variable & var = dof_map.variable(variable_number);
918 const FEType fe_type = var.type();
919 const bool add_p_level = fe_type.p_refinement;
920
921 // Pull objects out of the loop to reduce heap operations
922 std::vector<dof_id_type> my_dof_indices, parent_dof_indices;
923 std::unique_ptr<const Elem> my_side, parent_side;
924
925 // Look at the element faces. Check to see if we need to
926 // build constraints.
927 for (auto s : elem->side_index_range())
928 if (const Elem * const neigh = elem->neighbor_ptr(s);
929 neigh && neigh != remote_elem)
930 {
931 // constrain dofs shared between this element and ones coarser
932 if (neigh->level() < elem->level() &&
933 var.active_on_subdomain(neigh->subdomain_id()))
934 {
935 // than this element.
936 // Get pointers to the elements of interest and its parent.
937 const Elem * parent = elem->parent();
938
939 // This can't happen... Only level-0 elements have nullptr
940 // parents, and no level-0 elements can be at a higher
941 // level than their neighbors!
942 libmesh_assert(parent);
943
944 elem->build_side_ptr(my_side, s);
945 parent->build_side_ptr(parent_side, s);
946
947 // We have some elements with interior bubble function bases
948 // and lower-order sides. We need to only ask for the lower
949 // order in those cases.
950 FEType side_fe_type = fe_type;
951 int extra_order = 0;
952 if (my_side->default_order() < fe_type.order)
953 {
954 side_fe_type.order = my_side->default_order();
955 extra_order = (int)(side_fe_type.order) -
956 (int)(fe_type.order);
957 }
958
959 // This function gets called element-by-element, so there
960 // will be a lot of memory allocation going on. We can
961 // at least minimize this for the case of the dof indices
962 // by efficiently preallocating the requisite storage.
963 my_dof_indices.reserve (my_side->n_nodes());
964 parent_dof_indices.reserve (parent_side->n_nodes());
965
966 dof_map.dof_indices (my_side.get(), my_dof_indices,
967 variable_number, extra_order);
968 dof_map.dof_indices (parent_side.get(), parent_dof_indices,
969 variable_number, extra_order);
970
971 const unsigned int n_side_dofs =
972 FEInterface::n_dofs(side_fe_type, my_side.get());
973 const unsigned int n_parent_side_dofs =
974 FEInterface::n_dofs(side_fe_type, parent_side.get());
975 for (unsigned int my_dof=0; my_dof != n_side_dofs; my_dof++)
976 {
977 libmesh_assert_less (my_dof, my_side->n_nodes());
978
979 // My global dof index.
980 const dof_id_type my_dof_g = my_dof_indices[my_dof];
981
982 // Hunt for "constraining against myself" cases before
983 // we bother creating a constraint row
984 bool self_constraint = false;
985 for (unsigned int their_dof=0;
986 their_dof != n_parent_side_dofs; their_dof++)
987 {
988 libmesh_assert_less (their_dof, parent_side->n_nodes());
989
990 // Their global dof index.
991 const dof_id_type their_dof_g =
992 parent_dof_indices[their_dof];
993
994 if (their_dof_g == my_dof_g)
995 {
996 self_constraint = true;
997 break;
998 }
999 }
1000
1001 if (self_constraint)
1002 continue;
1003
1004 DofConstraintRow * constraint_row;
1005
1006 // we may be running constraint methods concurrently
1007 // on multiple threads, so we need a lock to
1008 // ensure that this constraint is "ours"
1009 {
1010 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1011
1012 if (dof_map.is_constrained_dof(my_dof_g))
1013 continue;
1014
1015 constraint_row = &(constraints[my_dof_g]);
1016 libmesh_assert(constraint_row->empty());
1017 }
1018
1019 // The support point of the DOF
1020 const Point & support_point = my_side->point(my_dof);
1021
1022 // Figure out where my node lies on their reference element.
1023 const Point mapped_point = FEMap::inverse_map(Dim-1,
1024 parent_side.get(),
1025 support_point);
1026
1027 // Compute the parent's side shape function values.
1028 for (unsigned int their_dof=0;
1029 their_dof != n_parent_side_dofs; their_dof++)
1030 {
1031 libmesh_assert_less (their_dof, parent_side->n_nodes());
1032
1033 // Their global dof index.
1034 const dof_id_type their_dof_g =
1035 parent_dof_indices[their_dof];
1036
1037 const Real their_dof_value = FEInterface::shape(side_fe_type,
1038 parent_side.get(),
1039 their_dof,
1040 mapped_point);
1041
1042 // Only add non-zero and non-identity values
1043 // for Lagrange basis functions.
1044 if ((std::abs(their_dof_value) > 1.e-5) &&
1045 (std::abs(their_dof_value) < .999))
1046 {
1047 constraint_row->emplace(their_dof_g, their_dof_value);
1048 }
1049 #ifdef DEBUG
1050 // Protect for the case u_i = 0.999 u_j,
1051 // in which case i better equal j.
1052 else if (their_dof_value >= .999)
1053 libmesh_assert_equal_to (my_dof_g, their_dof_g);
1054 #endif
1055 }
1056 }
1057 }
1058
1059 if (elem->active() && add_p_level)
1060 {
1061 // p refinement constraints:
1062 // constrain dofs shared between
1063 // active elements and neighbors with
1064 // lower polynomial degrees
1065 const unsigned int min_p_level =
1066 neigh->min_p_level_by_neighbor(elem, elem->p_level());
1067 if (min_p_level < elem->p_level())
1068 libmesh_error_msg(
1069 "Mismatched p-refinement levels for LAGRANGE is not currently supported");
1070 }
1071 }
1072} // lagrange_compute_constraints()
1073#endif // #ifdef LIBMESH_ENABLE_AMR
1074
1075} // anonymous namespace
1076
1077
1078// Instantiate (side_) nodal_soln() function for every dimension
1079LIBMESH_FE_NODAL_SOLN(LAGRANGE, lagrange_nodal_soln)
1080LIBMESH_FE_SIDE_NODAL_SOLN(LAGRANGE)
1081
1082
1083// Do full-specialization for every dimension, instead
1084// of explicit instantiation at the end of this function.
1085// This could be macro-ified.
1086template <> unsigned int FE<0,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1087template <> unsigned int FE<1,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1088template <> unsigned int FE<2,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1089template <> unsigned int FE<3,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1090
1091template <> unsigned int FE<0,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1092template <> unsigned int FE<1,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1093template <> unsigned int FE<2,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1094template <> unsigned int FE<3,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1095
1096
1097// Do full-specialization for every dimension, instead
1098// of explicit instantiation at the end of this function.
1099template <> 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); }
1100template <> 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); }
1101template <> 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); }
1102template <> 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); }
1103
1104template <> 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); }
1105template <> 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); }
1106template <> 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); }
1107template <> 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); }
1108
1109
1110// Lagrange elements have no dofs per element
1111// (just at the nodes)
1112template <> unsigned int FE<0,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1113template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1114template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1115template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1116
1117template <> unsigned int FE<0,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1118template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1119template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1120template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1121
1122// Lagrange FEMs are always C^0 continuous
1123template <> FEContinuity FE<0,LAGRANGE>::get_continuity() const { return C_ZERO; }
1124template <> FEContinuity FE<1,LAGRANGE>::get_continuity() const { return C_ZERO; }
1125template <> FEContinuity FE<2,LAGRANGE>::get_continuity() const { return C_ZERO; }
1126template <> FEContinuity FE<3,LAGRANGE>::get_continuity() const { return C_ZERO; }
1127
1128// Lagrange FEMs are not hierarchic
1129template <> bool FE<0,LAGRANGE>::is_hierarchic() const { return false; }
1130template <> bool FE<1,LAGRANGE>::is_hierarchic() const { return false; }
1131template <> bool FE<2,LAGRANGE>::is_hierarchic() const { return false; }
1132template <> bool FE<3,LAGRANGE>::is_hierarchic() const { return false; }
1133
1134// Lagrange FEM shapes do not need reinit (is this always true?)
1135template <> bool FE<0,LAGRANGE>::shapes_need_reinit() const { return false; }
1136template <> bool FE<1,LAGRANGE>::shapes_need_reinit() const { return false; }
1137template <> bool FE<2,LAGRANGE>::shapes_need_reinit() const { return false; }
1138template <> bool FE<3,LAGRANGE>::shapes_need_reinit() const { return false; }
1139
1140// Methods for computing Lagrange constraints. Note: we pass the
1141// dimension as the last argument to the anonymous helper function.
1142// Also note: we only need instantiations of this function for
1143// Dim==2 and 3.
1144#ifdef LIBMESH_ENABLE_AMR
1145template <>
1146void FE<2,LAGRANGE>::compute_constraints (DofConstraints & constraints,
1147 DofMap & dof_map,
1148 const unsigned int variable_number,
1149 const Elem * elem)
1150{ lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/2); }
1151
1152template <>
1153void FE<3,LAGRANGE>::compute_constraints (DofConstraints & constraints,
1154 DofMap & dof_map,
1155 const unsigned int variable_number,
1156 const Elem * elem)
1157{ lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/3); }
1158#endif // LIBMESH_ENABLE_AMR
1159
1160} // 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:3246
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:764
libMesh::FEInterface::n_dofs
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition fe_interface.C:356
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:148
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 Thu Nov 12 2026 20:15:21 for libMesh by doxygen 1.9.8