libMesh
fe_lagrange.C
Go to the documentation of this file.
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/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 
39 namespace libMesh
40 {
41 
42 // global helper function
43 void 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
432 namespace {
433 unsigned 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 
616 unsigned 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 
858 unsigned 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
867 void 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 = dof_map.should_p_refine_var(variable_number);
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
1070 LIBMESH_FE_NODAL_SOLN(LAGRANGE, lagrange_nodal_soln)
1071 LIBMESH_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.
1077 template <> unsigned int FE<0,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1078 template <> unsigned int FE<1,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1079 template <> unsigned int FE<2,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1080 template <> unsigned int FE<3,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1081 
1082 template <> unsigned int FE<0,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1083 template <> unsigned int FE<1,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1084 template <> unsigned int FE<2,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1085 template <> 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.
1090 template <> 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); }
1091 template <> 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); }
1092 template <> 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); }
1093 template <> 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 
1095 template <> 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); }
1096 template <> 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); }
1097 template <> 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); }
1098 template <> 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)
1103 template <> unsigned int FE<0,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1104 template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1105 template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1106 template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1107 
1108 template <> unsigned int FE<0,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1109 template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1110 template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1111 template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1112 
1113 // Lagrange FEMs are always C^0 continuous
1114 template <> FEContinuity FE<0,LAGRANGE>::get_continuity() const { return C_ZERO; }
1115 template <> FEContinuity FE<1,LAGRANGE>::get_continuity() const { return C_ZERO; }
1116 template <> FEContinuity FE<2,LAGRANGE>::get_continuity() const { return C_ZERO; }
1117 template <> FEContinuity FE<3,LAGRANGE>::get_continuity() const { return C_ZERO; }
1118 
1119 // Lagrange FEMs are not hierarchic
1120 template <> bool FE<0,LAGRANGE>::is_hierarchic() const { return false; }
1121 template <> bool FE<1,LAGRANGE>::is_hierarchic() const { return false; }
1122 template <> bool FE<2,LAGRANGE>::is_hierarchic() const { return false; }
1123 template <> bool FE<3,LAGRANGE>::is_hierarchic() const { return false; }
1124 
1125 // Lagrange FEM shapes do not need reinit (is this always true?)
1126 template <> bool FE<0,LAGRANGE>::shapes_need_reinit() const { return false; }
1127 template <> bool FE<1,LAGRANGE>::shapes_need_reinit() const { return false; }
1128 template <> bool FE<2,LAGRANGE>::shapes_need_reinit() const { return false; }
1129 template <> 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
1136 template <>
1137 void 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 
1143 template <>
1144 void 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
libMesh::FE::n_dofs
static unsigned int n_dofs(const ElemType t, const Order o)
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::Order
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
libMesh::FEInterface::n_dofs
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:355
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:1628
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
libMesh::FE::n_dofs_at_node
static unsigned int n_dofs_at_node(const ElemType t, const Order o, const unsigned int n)
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:3108
libMesh::FE::shapes_need_reinit
virtual bool shapes_need_reinit() const override
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::LIBMESH_FE_NODAL_SOLN
LIBMESH_FE_NODAL_SOLN(LIBMESH_FE_SIDE_NODAL_SOLN() LIBMESH_DEFAULT_NDOFS(BERNSTEIN) template<> FEContinuity FE< 0 BERNSTEIN, bernstein_nodal_soln)
Definition: fe_bernstein.C:415
libMesh::FE::is_hierarchic
virtual bool is_hierarchic() const override
libMesh::DofMap
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
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
libMesh::Elem::n_nodes
virtual unsigned int n_nodes() const =0
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:760
libMesh::libmesh_assert
libmesh_assert(ctx)
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::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::Utility::enum_to_string
std::string enum_to_string(const T e)
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::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::FEContinuity
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
Definition: enum_fe_family.h:84
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::Elem::type
virtual ElemType type() const =0
libMesh::DofConstraints
The constraint matrix storage format.
Definition: dof_map.h:108
int
void ErrorVector unsigned int
Definition: adjoints_ex3.C:360
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:117
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67
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::remote_elem
const RemoteElem * remote_elem
Definition: remote_elem.C:57
Generated on Thu Jul 17 2025 01:29:06 for libMesh by   doxygen 1.8.14