libMesh
fe_hierarchic_shape_2D.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 // Local includes
20 #include "libmesh/fe.h"
21 #include "libmesh/elem.h"
22 #include "libmesh/number_lookups.h"
23 #include "libmesh/enum_to_string.h"
24 
25 // Anonymous namespace for functions shared by HIERARCHIC and
26 // L2_HIERARCHIC implementations. Implementations appear at the bottom
27 // of this file.
28 namespace
29 {
30 using namespace libMesh;
31 
32 Real fe_triangle_helper (const Elem & elem,
33  const Real edgenumerator,
34  const Real crossval,
35  const unsigned int basisorder,
36  const Order totalorder,
37  const unsigned int noden);
38 
39 template <FEFamily T>
40 Real fe_hierarchic_2D_shape(const Elem * elem,
41  const Order order,
42  const unsigned int i,
43  const Point & p,
44  const bool add_p_level);
45 
46 template <FEFamily T>
47 Real fe_hierarchic_2D_shape_deriv(const Elem * elem,
48  const Order order,
49  const unsigned int i,
50  const unsigned int j,
51  const Point & p,
52  const bool add_p_level);
53 
54 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
55 
56 template <FEFamily T>
57 Real fe_hierarchic_2D_shape_second_deriv(const Elem * elem,
58  const Order order,
59  const unsigned int i,
60  const unsigned int j,
61  const Point & p,
62  const bool add_p_level);
63 
64 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
65 
66 
67 std::tuple<unsigned int, unsigned int, Real>
68 quad_indices(const Elem * elem,
69  const unsigned int totalorder,
70  const unsigned int i)
71 {
72  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
73 
74  // Example i, i0, i1 values for totalorder = 5:
75  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
76  // static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 3, 2, 4, 4, 4, 3, 2, 5, 5, 5, 5, 4, 3, 2};
77  // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 3, 3, 2, 3, 4, 4, 4, 2, 3, 4, 5, 5, 5, 5};
78 
79  unsigned int i0, i1;
80 
81  // Vertex DoFs
82  if (i == 0)
83  { i0 = 0; i1 = 0; }
84  else if (i == 1)
85  { i0 = 1; i1 = 0; }
86  else if (i == 2)
87  { i0 = 1; i1 = 1; }
88  else if (i == 3)
89  { i0 = 0; i1 = 1; }
90  // Edge DoFs
91  else if (i < totalorder + 3u)
92  { i0 = i - 2; i1 = 0; }
93  else if (i < 2u*totalorder + 2)
94  { i0 = 1; i1 = i - totalorder - 1; }
95  else if (i < 3u*totalorder + 1)
96  { i0 = i - 2u*totalorder; i1 = 1; }
97  else if (i < 4u*totalorder)
98  { i0 = 0; i1 = i - 3u*totalorder + 1; }
99  // Interior DoFs
100  else
101  {
102  unsigned int basisnum = i - 4*totalorder;
103  i0 = square_number_column[basisnum] + 2;
104  i1 = square_number_row[basisnum] + 2;
105  }
106 
107  // Flip odd degree of freedom values if necessary
108  // to keep continuity on sides
109  Real f = 1.;
110 
111  if ((i0%2) && (i0 > 2) && (i1 == 0))
112  f = elem->positive_edge_orientation(0)?-1.:1.;
113  else if ((i0%2) && (i0>2) && (i1 == 1))
114  f = !elem->positive_edge_orientation(2)?-1.:1.;
115  else if ((i0 == 0) && (i1%2) && (i1>2))
116  f = !elem->positive_edge_orientation(3)?-1.:1.;
117  else if ((i0 == 1) && (i1%2) && (i1>2))
118  f = elem->positive_edge_orientation(1)?-1.:1.;
119 
120  return {i0, i1, f};
121 }
122 
123 } // anonymous namespace
124 
125 
126 
127 namespace libMesh
128 {
129 
130 
131 LIBMESH_DEFAULT_VECTORIZED_FE(2,HIERARCHIC)
132 LIBMESH_DEFAULT_VECTORIZED_FE(2,L2_HIERARCHIC)
133 LIBMESH_DEFAULT_VECTORIZED_FE(2,SIDE_HIERARCHIC)
134 
135 
136 template <>
137 Real FE<2,HIERARCHIC>::shape(const ElemType,
138  const Order,
139  const unsigned int,
140  const Point &)
141 {
142  libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
143  return 0.;
144 }
145 
146 
147 
148 template <>
149 Real FE<2,L2_HIERARCHIC>::shape(const ElemType,
150  const Order,
151  const unsigned int,
152  const Point &)
153 {
154  libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
155  return 0.;
156 }
157 
158 
159 
160 template <>
161 Real FE<2,SIDE_HIERARCHIC>::shape(const ElemType,
162  const Order,
163  const unsigned int,
164  const Point &)
165 {
166  libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
167  return 0.;
168 }
169 
170 
171 
172 template <>
173 Real FE<2,HIERARCHIC>::shape(const Elem * elem,
174  const Order order,
175  const unsigned int i,
176  const Point & p,
177  const bool add_p_level)
178 {
179  return fe_hierarchic_2D_shape<HIERARCHIC>(elem, order, i, p, add_p_level);
180 }
181 
182 
183 
184 template <>
185 Real FE<2,HIERARCHIC>::shape(const FEType fet,
186  const Elem * elem,
187  const unsigned int i,
188  const Point & p,
189  const bool add_p_level)
190 {
191  return fe_hierarchic_2D_shape<HIERARCHIC>(elem, fet.order, i, p, add_p_level);
192 }
193 
194 
195 template <>
196 Real FE<2,L2_HIERARCHIC>::shape(const Elem * elem,
197  const Order order,
198  const unsigned int i,
199  const Point & p,
200  const bool add_p_level)
201 {
202  return fe_hierarchic_2D_shape<L2_HIERARCHIC>(elem, order, i, p, add_p_level);
203 }
204 
205 
206 template <>
207 Real FE<2,L2_HIERARCHIC>::shape(const FEType fet,
208  const Elem * elem,
209  const unsigned int i,
210  const Point & p,
211  const bool add_p_level)
212 {
213  return fe_hierarchic_2D_shape<L2_HIERARCHIC>(elem, fet.order, i, p, add_p_level);
214 }
215 
216 
217 template <>
218 Real FE<2,SIDE_HIERARCHIC>::shape(const Elem * elem,
219  const Order order,
220  const unsigned int i,
221  const Point & p,
222  const bool add_p_level)
223 {
224  libmesh_assert(elem);
225  const ElemType type = elem->type();
226 
227  const Order totalorder = order + add_p_level*elem->p_level();
228 
229  const unsigned int dofs_per_side = totalorder+1u;
230 
231  switch (type)
232  {
233  case TRI6:
234  case TRI7:
235  {
236  libmesh_assert_less(i, 3*dofs_per_side);
237 
238  // Flip odd degree of freedom values if necessary
239  // to keep continuity on sides. We'll flip xi/eta rather than
240  // flipping phi, so that we can use this to handle the "nodal"
241  // degrees of freedom too.
242  Real f = 1.;
243 
244  const Real zeta1 = p(0);
245  const Real zeta2 = p(1);
246  const Real zeta0 = 1. - zeta1 - zeta2;
247 
248  if (zeta1 > zeta2 && zeta0 > zeta2) // side 0
249  {
250  if (i >= dofs_per_side)
251  return 0;
252 
253  if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
254  return 1;
255 
256  if ((i < 2 || i % 2) &&
257  elem->positive_edge_orientation(0))
258  f = -1;
259 
260  return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i, f*(zeta1-zeta0));
261  }
262  else if (zeta1 > zeta0 && zeta2 > zeta0) // side 1
263  {
264  if (i < dofs_per_side ||
265  i >= 2*dofs_per_side)
266  return 0;
267 
268  if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
269  return 1;
270 
271  const unsigned int side_i = i - dofs_per_side;
272 
273  if ((side_i < 2 || side_i % 2) &&
274  elem->positive_edge_orientation(1))
275  f = -1;
276 
277  return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, side_i, f*(zeta2-zeta1));
278  }
279  else // side 2
280  {
281  libmesh_assert (zeta2 >= zeta1 && zeta0 >= zeta1); // On a corner???
282 
283  if (i < 2*dofs_per_side)
284  return 0;
285 
286  if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
287  return 1;
288 
289  const unsigned int side_i = i - 2*dofs_per_side;
290 
291  if ((side_i < 2 || side_i % 2) &&
292  elem->positive_edge_orientation(2))
293  f = -1;
294 
295  return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, side_i, f*(zeta0-zeta2));
296  }
297  }
298  case QUAD8:
299  case QUADSHELL8:
300  case QUAD9:
301  case QUADSHELL9:
302  {
303  libmesh_assert_less(i, 4*dofs_per_side);
304 
305  // Flip odd degree of freedom values if necessary
306  // to keep continuity on sides. We'll flip xi/eta rather than
307  // flipping phi, so that we can use this to handle the "nodal"
308  // degrees of freedom too.
309  Real f = 1.;
310 
311  const Real xi = p(0), eta = p(1);
312  if (eta < xi)
313  {
314  if (eta < -xi) // side 0
315  {
316  if (i >= dofs_per_side)
317  return 0;
318 
319  if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
320  return 1;
321 
322  if ((i < 2 || i % 2) &&
323  elem->positive_edge_orientation(0))
324  f = -1;
325 
326  return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i, f*xi);
327  }
328  else // side 1
329  {
330  if (i < dofs_per_side ||
331  i >= 2*dofs_per_side)
332  return 0;
333 
334  if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
335  return 1;
336 
337  const unsigned int side_i = i - dofs_per_side;
338 
339  if ((side_i < 2 || side_i % 2) &&
340  elem->positive_edge_orientation(1))
341  f = -1;
342 
343  return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, side_i, f*eta);
344  }
345  }
346  else // xi < eta
347  {
348  if (eta > -xi) // side 2
349  {
350  if (i < 2*dofs_per_side ||
351  i >= 3*dofs_per_side)
352  return 0;
353 
354  if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
355  return 1;
356 
357  const unsigned int side_i = i - 2*dofs_per_side;
358 
359  if ((side_i < 2 || side_i % 2) &&
360  !elem->positive_edge_orientation(2))
361  f = -1;
362 
363  return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, side_i, f*xi);
364  }
365  else // side 3
366  {
367  if (i < 3*dofs_per_side)
368  return 0;
369 
370  if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
371  return 1;
372 
373  const unsigned int side_i = i - 3*dofs_per_side;
374 
375  if ((side_i < 2 || side_i % 2) &&
376  !elem->positive_edge_orientation(3))
377  f = -1;
378 
379  return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, side_i, f*eta);
380  }
381  }
382  }
383  default:
384  libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(elem->type()));
385  }
386  return 0;
387 }
388 
389 
390 template <>
391 Real FE<2,SIDE_HIERARCHIC>::shape(const FEType fet,
392  const Elem * elem,
393  const unsigned int i,
394  const Point & p,
395  const bool add_p_level)
396 {
397  return FE<2,SIDE_HIERARCHIC>::shape(elem, fet.order, i, p, add_p_level);
398 }
399 
400 
401 template <>
402 Real FE<2,HIERARCHIC>::shape_deriv(const ElemType,
403  const Order,
404  const unsigned int,
405  const unsigned int,
406  const Point &)
407 {
408  libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
409  return 0.;
410 }
411 
412 
413 
414 template <>
415 Real FE<2,L2_HIERARCHIC>::shape_deriv(const ElemType,
416  const Order,
417  const unsigned int,
418  const unsigned int,
419  const Point &)
420 {
421  libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
422  return 0.;
423 }
424 
425 
426 
427 template <>
428 Real FE<2,SIDE_HIERARCHIC>::shape_deriv(const ElemType,
429  const Order,
430  const unsigned int,
431  const unsigned int,
432  const Point &)
433 {
434  libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
435  return 0.;
436 }
437 
438 
439 
440 template <>
441 Real FE<2,HIERARCHIC>::shape_deriv(const Elem * elem,
442  const Order order,
443  const unsigned int i,
444  const unsigned int j,
445  const Point & p,
446  const bool add_p_level)
447 {
448  return fe_hierarchic_2D_shape_deriv<HIERARCHIC>(elem, order, i, j, p, add_p_level);
449 }
450 
451 
452 template <>
453 Real FE<2,HIERARCHIC>::shape_deriv(const FEType fet,
454  const Elem * elem,
455  const unsigned int i,
456  const unsigned int j,
457  const Point & p,
458  const bool add_p_level)
459 {
460  return fe_hierarchic_2D_shape_deriv<HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
461 }
462 
463 
464 
465 
466 template <>
467 Real FE<2,L2_HIERARCHIC>::shape_deriv(const Elem * elem,
468  const Order order,
469  const unsigned int i,
470  const unsigned int j,
471  const Point & p,
472  const bool add_p_level)
473 {
474  return fe_hierarchic_2D_shape_deriv<L2_HIERARCHIC>(elem, order, i, j, p, add_p_level);
475 }
476 
477 
478 template <>
479 Real FE<2,L2_HIERARCHIC>::shape_deriv(const FEType fet,
480  const Elem * elem,
481  const unsigned int i,
482  const unsigned int j,
483  const Point & p,
484  const bool add_p_level)
485 {
486  return fe_hierarchic_2D_shape_deriv<L2_HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
487 }
488 
489 
490 
491 template <>
492 Real FE<2,SIDE_HIERARCHIC>::shape_deriv(const Elem * elem,
493  const Order order,
494  const unsigned int i,
495  const unsigned int j,
496  const Point & p,
497  const bool add_p_level)
498 {
499  libmesh_assert(elem);
500 
501  const ElemType type = elem->type();
502 
503  const Order totalorder = order + add_p_level*elem->p_level();
504 
505  if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
506  return 0;
507 
508  const unsigned int dofs_per_side = totalorder+1u;
509 
510  switch (type)
511  {
512  case TRI6:
513  case TRI7:
514  {
515  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SIDE_HIERARCHIC>::shape);
516  }
517 #if 0
518  {
519  libmesh_assert_less(i, 3*dofs_per_side);
520  libmesh_assert_less (j, 2);
521 
522  // Flip odd degree of freedom values if necessary
523  // to keep continuity on sides. We'll flip xi/eta rather than
524  // flipping phi, so that we can use this to handle the "nodal"
525  // degrees of freedom too.
526  Real f = 1.;
527 
528  const Real zeta1 = p(0);
529  const Real zeta2 = p(1);
530  const Real zeta0 = 1. - zeta1 - zeta2;
531 
532  if (zeta1 > zeta2 && zeta0 > zeta2) // side 0
533  {
534  if (j == 1) // d/deta is perpendicular here
535  return 0;
536 
537  if (i >= dofs_per_side)
538  return 0;
539 
540  if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
541  return 0;
542 
543  if ((i < 2 || i % 2) &&
544  elem->positive_edge_orientation(0))
545  f = -1;
546 
547  return f*FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, i, 0, f*(zeta1-zeta0));
548  }
549  else if (zeta1 > zeta0 && zeta2 > zeta0) // side 1
550  {
551  if (i < dofs_per_side ||
552  i >= 2*dofs_per_side)
553  return 0;
554 
555  if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
556  return 0;
557 
558  const unsigned int side_i = i - dofs_per_side;
559 
560  if ((side_i < 2 || side_i % 2) &&
561  elem->positive_edge_orientation(1))
562  f = -1;
563 
564  Real g = 1;
565  if (j == 0) // 2D d/dxi is in the opposite direction on this edge
566  g = -1;
567 
568  return f*g*FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, side_i, 0, f*(zeta2-zeta1));
569  }
570  else // side 2
571  {
572  libmesh_assert (zeta2 >= zeta1 && zeta0 >= zeta1); // On a corner???
573 
574  if (j == 0) // d/dxi is perpendicular here
575  return 0;
576 
577  if (i < 2*dofs_per_side)
578  return 0;
579 
580  if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
581  return 0;
582 
583  const unsigned int side_i = i - 2*dofs_per_side;
584 
585  if ((side_i < 2 || side_i % 2) &&
586  elem->positive_edge_orientation(2))
587  f = -1;
588 
589  return -f*FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, side_i, 0, f*(zeta0-zeta2));
590  }
591  }
592 #endif
593  case QUAD8:
594  case QUADSHELL8:
595  case QUAD9:
596  case QUADSHELL9:
597  {
598  libmesh_assert_less(i, 4*dofs_per_side);
599 
600  // Flip odd degree of freedom values if necessary
601  // to keep continuity on sides. We'll flip xi/eta rather than
602  // flipping phi, so that we can use this to handle the "nodal"
603  // degrees of freedom too.
604  Real f = 1.;
605 
606  const Real xi = p(0), eta = p(1);
607  if (eta < xi)
608  {
609  if (eta < -xi) // side 0
610  {
611  if (i >= dofs_per_side)
612  return 0;
613  if (j != 0)
614  return 0;
615  if ((i < 2 || i % 2) &&
616  elem->positive_edge_orientation(0))
617  f = -1;
618 
619  return f*FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, i, 0, f*xi);
620  }
621  else // side 1
622  {
623  if (i < dofs_per_side ||
624  i >= 2*dofs_per_side)
625  return 0;
626  if (j != 1)
627  return 0;
628 
629  const unsigned int side_i = i - dofs_per_side;
630 
631  if ((side_i < 2 || side_i % 2) &&
632  elem->positive_edge_orientation(1))
633  f = -1;
634 
635  return f*FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, side_i, 0, f*eta);
636  }
637  }
638  else // xi < eta
639  {
640  if (eta > -xi) // side 2
641  {
642  if (i < 2*dofs_per_side ||
643  i >= 3*dofs_per_side)
644  return 0;
645  if (j != 0)
646  return 0;
647 
648  const unsigned int side_i = i - 2*dofs_per_side;
649 
650  if ((side_i < 2 || side_i % 2) &&
651  !elem->positive_edge_orientation(2))
652  f = -1;
653 
654  return f*FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, side_i, 0, f*xi);
655  }
656  else // side 3
657  {
658  if (i < 3*dofs_per_side)
659  return 0;
660  if (j != 1)
661  return 0;
662 
663  const unsigned int side_i = i - 3*dofs_per_side;
664 
665  if ((side_i < 2 || side_i % 2) &&
666  !elem->positive_edge_orientation(3))
667  f = -1;
668 
669  return f*FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, side_i, 0, f*eta);
670  }
671  }
672  }
673  default:
674  libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(elem->type()));
675  }
676  return 0;
677 }
678 
679 
680 template <>
681 Real FE<2,SIDE_HIERARCHIC>::shape_deriv(const FEType fet,
682  const Elem * elem,
683  const unsigned int i,
684  const unsigned int j,
685  const Point & p,
686  const bool add_p_level)
687 {
688  return FE<2,SIDE_HIERARCHIC>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
689 }
690 
691 
692 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
693 
694 template <>
695 Real FE<2,HIERARCHIC>::shape_second_deriv(const ElemType,
696  const Order,
697  const unsigned int,
698  const unsigned int,
699  const Point &)
700 {
701  libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
702  return 0.;
703 }
704 
705 
706 
707 template <>
708 Real FE<2,L2_HIERARCHIC>::shape_second_deriv(const ElemType,
709  const Order,
710  const unsigned int,
711  const unsigned int,
712  const Point &)
713 {
714  libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
715  return 0.;
716 }
717 
718 
719 
720 template <>
721 Real FE<2,SIDE_HIERARCHIC>::shape_second_deriv(const ElemType,
722  const Order,
723  const unsigned int,
724  const unsigned int,
725  const Point &)
726 {
727  libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
728  return 0.;
729 }
730 
731 
732 
733 template <>
734 Real FE<2,HIERARCHIC>::shape_second_deriv(const Elem * elem,
735  const Order order,
736  const unsigned int i,
737  const unsigned int j,
738  const Point & p,
739  const bool add_p_level)
740 {
741  return fe_hierarchic_2D_shape_second_deriv<HIERARCHIC>(elem, order, i, j, p, add_p_level);
742 }
743 
744 
745 template <>
746 Real FE<2,HIERARCHIC>::shape_second_deriv(const FEType fet,
747  const Elem * elem,
748  const unsigned int i,
749  const unsigned int j,
750  const Point & p,
751  const bool add_p_level)
752 {
753  return fe_hierarchic_2D_shape_second_deriv<HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
754 }
755 
756 
757 template <>
758 Real FE<2,L2_HIERARCHIC>::shape_second_deriv(const Elem * elem,
759  const Order order,
760  const unsigned int i,
761  const unsigned int j,
762  const Point & p,
763  const bool add_p_level)
764 {
765  return fe_hierarchic_2D_shape_second_deriv<L2_HIERARCHIC>(elem, order, i, j, p, add_p_level);
766 }
767 
768 
769 template <>
770 Real FE<2,L2_HIERARCHIC>::shape_second_deriv(const FEType fet,
771  const Elem * elem,
772  const unsigned int i,
773  const unsigned int j,
774  const Point & p,
775  const bool add_p_level)
776 {
777  return fe_hierarchic_2D_shape_second_deriv<L2_HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
778 }
779 
780 
781 template <>
782 Real FE<2,SIDE_HIERARCHIC>::shape_second_deriv(const Elem * elem,
783  const Order order,
784  const unsigned int i,
785  const unsigned int j,
786  const Point & p,
787  const bool add_p_level)
788 {
789  libmesh_assert(elem);
790  const ElemType type = elem->type();
791 
792  const Order totalorder = order + add_p_level*elem->p_level();
793 
794  if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
795  return 0;
796 
797  const unsigned int dofs_per_side = totalorder+1u;
798 
799  switch (type)
800  {
801  case TRI6:
802  case TRI7:
803  {
804  return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
805  FE<2,SIDE_HIERARCHIC>::shape_deriv);
806  }
807  case QUAD8:
808  case QUADSHELL8:
809  case QUAD9:
810  case QUADSHELL9:
811  {
812  libmesh_assert_less(i, 4*dofs_per_side);
813 
814  // Flip odd degree of freedom values if necessary
815  // to keep continuity on sides. We'll flip xi/eta rather than
816  // flipping phi, so that we can use this to handle the "nodal"
817  // degrees of freedom too.
818  Real f = 1.;
819 
820  const Real xi = p(0), eta = p(1);
821  if (eta < xi)
822  {
823  if (eta < -xi) // side 0
824  {
825  if (i >= dofs_per_side)
826  return 0;
827  if (j != 0)
828  return 0;
829  if ((i < 2 || i % 2) &&
830  elem->positive_edge_orientation(0))
831  f = -1;
832 
833  return FE<1,HIERARCHIC>::shape_second_deriv(EDGE3, totalorder, i, 0, f*xi);
834  }
835  else // side 1
836  {
837  if (i < dofs_per_side ||
838  i >= 2*dofs_per_side)
839  return 0;
840  if (j != 2)
841  return 0;
842 
843  const unsigned int side_i = i - dofs_per_side;
844 
845  if ((side_i < 2 || side_i % 2) &&
846  elem->positive_edge_orientation(1))
847  f = -1;
848 
849  return FE<1,HIERARCHIC>::shape_second_deriv(EDGE3, totalorder, side_i, 0, f*eta);
850  }
851  }
852  else // xi < eta
853  {
854  if (eta > -xi) // side 2
855  {
856  if (i < 2*dofs_per_side ||
857  i >= 3*dofs_per_side)
858  return 0;
859  if (j != 0)
860  return 0;
861 
862  const unsigned int side_i = i - 2*dofs_per_side;
863 
864  if ((side_i < 2 || side_i % 2) &&
865  !elem->positive_edge_orientation(2))
866  f = -1;
867 
868  return FE<1,HIERARCHIC>::shape_second_deriv(EDGE3, totalorder, side_i, 0, f*xi);
869  }
870  else // side 3
871  {
872  if (i < 3*dofs_per_side)
873  return 0;
874  if (j != 2)
875  return 0;
876 
877  const unsigned int side_i = i - 3*dofs_per_side;
878 
879  if ((side_i < 2 || side_i % 2) &&
880  !elem->positive_edge_orientation(3))
881  f = -1;
882 
883  return FE<1,HIERARCHIC>::shape_second_deriv(EDGE3, totalorder, side_i, 0, f*eta);
884  }
885  }
886  }
887  default:
888  libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(elem->type()));
889  }
890  return 0;
891 }
892 
893 
894 template <>
895 Real FE<2,SIDE_HIERARCHIC>::shape_second_deriv(const FEType fet,
896  const Elem * elem,
897  const unsigned int i,
898  const unsigned int j,
899  const Point & p,
900  const bool add_p_level)
901 {
902  return FE<2,SIDE_HIERARCHIC>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
903 }
904 
905 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
906 
907 } // namespace libMesh
908 
909 
910 
911 namespace
912 {
913 using namespace libMesh;
914 
915 Real fe_triangle_helper (const Elem & elem,
916  const Real edgenumerator,
917  const Real crossval,
918  const unsigned int basisorder,
919  const Order totalorder,
920  const unsigned int noden)
921 {
922  // Get factors to account for edge-flipping
923  Real flip = 1;
924  if (basisorder%2 && (elem.point(noden) > elem.point((noden+1)%3)))
925  flip = -1.;
926 
927  // Avoid NaN around vertices ... but we still have to match the true
928  // function, even when we're *outside* the triangle (crossval==0 on
929  // a line, not just at a point!), to handle imprecise queries and
930  // FDM derivatives correctly!
931  if (crossval == 0.)
932  {
933  unsigned int basisfactorial = 1.;
934  for (unsigned int n=2; n <= basisorder; ++n)
935  basisfactorial *= n;
936 
937  return std::pow(edgenumerator, basisorder) / basisfactorial;
938  }
939  // Experimentally, as c -> 0, n propto c, I'm still seeing good
940  // behavior from the default implementation below:
941 
942  const Real edgeval = edgenumerator / crossval;
943  const Real crossfunc = std::pow(crossval, basisorder);
944 
945  return flip * crossfunc *
946  FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
947  basisorder, edgeval);
948 }
949 
950 template <FEFamily T>
951 Real fe_hierarchic_2D_shape(const Elem * elem,
952  const Order order,
953  const unsigned int i,
954  const Point & p,
955  const bool add_p_level)
956 {
957  libmesh_assert(elem);
958 
959  const Order totalorder = order + add_p_level*elem->p_level();
960  libmesh_assert_greater (totalorder, 0);
961 
962  switch (elem->type())
963  {
964  case TRI3:
965  case TRISHELL3:
966  case TRI6:
967  case TRI7:
968  {
969  const Real zeta1 = p(0);
970  const Real zeta2 = p(1);
971  const Real zeta0 = 1. - zeta1 - zeta2;
972 
973  libmesh_assert_less (i, (totalorder+1u)*(totalorder+2u)/2);
974  libmesh_assert (T == L2_HIERARCHIC || elem->type() == TRI6 ||
975  elem->type() == TRI7 || totalorder < 2);
976 
977  // Vertex DoFs
978  if (i == 0)
979  return zeta0;
980  else if (i == 1)
981  return zeta1;
982  else if (i == 2)
983  return zeta2;
984  // Edge DoFs
985  else if (i < totalorder + 2u)
986  {
987  const unsigned int basisorder = i - 1;
988 
989  const Real crossval = zeta0 + zeta1;
990  const Real edgenumerator = zeta1 - zeta0;
991 
992  return fe_triangle_helper(*elem, edgenumerator, crossval,
993  basisorder, totalorder, 0);
994  }
995  else if (i < 2u*totalorder + 1)
996  {
997  const unsigned int basisorder = i - totalorder;
998 
999  const Real crossval = zeta2 + zeta1;
1000  const Real edgenumerator = zeta2 - zeta1;
1001 
1002  return fe_triangle_helper(*elem, edgenumerator, crossval,
1003  basisorder, totalorder, 1);
1004  }
1005  else if (i < 3u*totalorder)
1006  {
1007  const unsigned int basisorder = i - (2u*totalorder) + 1;
1008 
1009  const Real crossval = zeta0 + zeta2;
1010  const Real edgenumerator = zeta0 - zeta2;
1011 
1012  return fe_triangle_helper(*elem, edgenumerator, crossval,
1013  basisorder, totalorder, 2);
1014  }
1015  // Interior DoFs
1016  else
1017  {
1018  const unsigned int basisnum = i - (3u*totalorder);
1019  unsigned int exp0 = triangular_number_column[basisnum] + 1;
1020  unsigned int exp1 = triangular_number_row[basisnum] + 1 -
1021  triangular_number_column[basisnum];
1022 
1023  Real returnval = 1;
1024  for (unsigned int n = 0; n != exp0; ++n)
1025  returnval *= zeta0;
1026  for (unsigned int n = 0; n != exp1; ++n)
1027  returnval *= zeta1;
1028  returnval *= zeta2;
1029  return returnval;
1030  }
1031  }
1032 
1033  // Hierarchic shape functions on the quadrilateral.
1034  case QUAD4:
1035  case QUADSHELL4:
1036  libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
1037  libmesh_fallthrough();
1038  case QUAD8:
1039  case QUADSHELL8:
1040  case QUAD9:
1041  case QUADSHELL9:
1042  {
1043  // Compute quad shape functions as a tensor-product
1044  auto [i0, i1, f] = quad_indices(elem, totalorder, i);
1045 
1046  return f*(FE<1,T>::shape(EDGE3, totalorder, i0, p(0))*
1047  FE<1,T>::shape(EDGE3, totalorder, i1, p(1)));
1048  }
1049 
1050  default:
1051  libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(elem->type()));
1052  }
1053 
1054  return 0.;
1055 }
1056 
1057 
1058 template <FEFamily T>
1059 Real fe_hierarchic_2D_shape_deriv(const Elem * elem,
1060  const Order order,
1061  const unsigned int i,
1062  const unsigned int j,
1063  const Point & p,
1064  const bool add_p_level)
1065 {
1066  libmesh_assert(elem);
1067 
1068  const ElemType type = elem->type();
1069 
1070  const Order totalorder = order + add_p_level*elem->p_level();
1071 
1072  libmesh_assert_greater (totalorder, 0);
1073 
1074  switch (type)
1075  {
1076  // 1st & 2nd-order Hierarchics.
1077  case TRI3:
1078  case TRISHELL3:
1079  case TRI6:
1080  case TRI7:
1081  {
1082  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,T>::shape);
1083  }
1084 
1085  case QUAD4:
1086  case QUADSHELL4:
1087  libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
1088  libmesh_fallthrough();
1089  case QUAD8:
1090  case QUADSHELL8:
1091  case QUAD9:
1092  case QUADSHELL9:
1093  {
1094  // Compute quad shape functions as a tensor-product
1095  auto [i0, i1, f] = quad_indices(elem, totalorder, i);
1096 
1097  switch (j)
1098  {
1099  // d()/dxi
1100  case 0:
1101  return f*(FE<1,T>::shape_deriv(EDGE3, totalorder, i0, 0, p(0))*
1102  FE<1,T>::shape (EDGE3, totalorder, i1, p(1)));
1103 
1104  // d()/deta
1105  case 1:
1106  return f*(FE<1,T>::shape (EDGE3, totalorder, i0, p(0))*
1107  FE<1,T>::shape_deriv(EDGE3, totalorder, i1, 0, p(1)));
1108 
1109  default:
1110  libmesh_error_msg("Invalid derivative index j = " << j);
1111  }
1112  }
1113 
1114  default:
1115  libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
1116  }
1117 
1118  return 0.;
1119 }
1120 
1121 
1122 
1123 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1124 
1125 template <FEFamily T>
1126 Real fe_hierarchic_2D_shape_second_deriv(const Elem * elem,
1127  const Order order,
1128  const unsigned int i,
1129  const unsigned int j,
1130  const Point & p,
1131  const bool add_p_level)
1132 {
1133  return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
1134  FE<2,T>::shape_deriv);
1135 }
1136 
1137 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1138 
1139 } // anonymous namespace
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
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::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
libMesh::triangular_number_row
const unsigned char triangular_number_row[]
Definition: number_lookups.C:28
libMesh::FE::shape_deriv
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:3109
libMesh::FEType::order
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:215
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::fe_fdm_deriv
OutputShape fe_fdm_deriv(const Elem *elem, const Order order, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level, OutputShape(*shape_func)(const Elem *, const Order, const unsigned int, const Point &, const bool))
Helper functions for finite differenced derivatives in cases where analytical calculations haven&#39;t be...
Definition: fe.C:911
libMesh::square_number_column
const unsigned char square_number_column[]
Definition: number_lookups.C:56
libMesh::LIBMESH_DEFAULT_VECTORIZED_FE
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
Definition: fe_bernstein_shape_0D.C:30
libMesh::FE
A specific instantiation of the FEBase class.
Definition: fe.h:127
libMesh::Utility::pow
T pow(const T &x)
Definition: utility.h:328
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Elem::positive_edge_orientation
bool positive_edge_orientation(const unsigned int i) const
Definition: elem.C:3604
libMesh::fe_fdm_second_deriv
OutputShape fe_fdm_second_deriv(const Elem *elem, const Order order, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level, OutputShape(*deriv_func)(const Elem *, const Order, const unsigned int, const unsigned int, const Point &, const bool))
Definition: fe.C:969
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::triangular_number_column
const unsigned char triangular_number_column[]
Definition: number_lookups.C:41
libMesh::square_number_row
const unsigned char square_number_row[]
Definition: number_lookups.C:69
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:2454
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
Generated on Thu Oct 30 2025 20:01:50 for libMesh by   doxygen 1.8.14