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