libMesh
fe_hierarchic_shape_3D.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 #include "libmesh/cell_tet4.h" // We need edge_nodes_map + side_nodes_map
25 #include "libmesh/cell_prism6.h"
26 #include "libmesh/face_tri3.h" // Faster to construct these on the stack
27 #include "libmesh/face_quad4.h"
28 
29 // Anonymous namespace for functions shared by HIERARCHIC and
30 // L2_HIERARCHIC implementations. Implementations appear at the bottom
31 // of this file.
32 namespace
33 {
34 using namespace libMesh;
35 
36 unsigned int cube_side(const Point & p);
37 
38 Point cube_side_point(unsigned int sidenum, const Point & interior_point);
39 
40 std::array<unsigned int, 4> oriented_prism_nodes(const Elem & elem,
41  unsigned int face_num);
42 
43 std::array<unsigned int, 3> oriented_tet_nodes(const Elem & elem,
44  unsigned int face_num);
45 
46 template <FEFamily T>
47 Real fe_hierarchic_3D_shape(const Elem * elem,
48  const Order order,
49  const unsigned int i,
50  const Point & p,
51  const bool add_p_level);
52 
53 template <FEFamily T>
54 Real fe_hierarchic_3D_shape_deriv(const Elem * elem,
55  const Order order,
56  const unsigned int i,
57  const unsigned int j,
58  const Point & p,
59  const bool add_p_level);
60 
61 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
62 
63 template <FEFamily T>
64 Real fe_hierarchic_3D_shape_second_deriv(const Elem * elem,
65  const Order order,
66  const unsigned int i,
67  const unsigned int j,
68  const Point & p,
69  const bool add_p_level);
70 
71 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
72 
73 #if LIBMESH_DIM > 2
74 Point get_min_point(const Elem * elem,
75  unsigned int a,
76  unsigned int b,
77  unsigned int c,
78  unsigned int d)
79 {
80  return std::min(std::min(elem->point(a),elem->point(b)),
81  std::min(elem->point(c),elem->point(d)));
82 }
83 
84 // Remap non-face-nodes based on point ordering
85 template <unsigned int N_nodes>
86 unsigned int remap_node(unsigned int n,
87  const Elem & elem,
88  unsigned int nodebegin)
89 {
90  std::array<const Point *, N_nodes> points;
91 
92  for (auto i : IntRange<unsigned int>(0, N_nodes))
93  points[i] = &elem.point(nodebegin+i);
94 
95  std::sort(points.begin(), points.end(),
96  [](const Point * a, const Point * b)
97  { return *a < *b; });
98 
99  const Point * pn = points[n-nodebegin];
100 
101  for (auto i : IntRange<unsigned int>(nodebegin, nodebegin+N_nodes))
102  if (pn == &elem.point(i))
103  return i;
104 
105  libmesh_assert(false);
106  return libMesh::invalid_uint;
107 }
108 
109 
110 void cube_remap(unsigned int & side_i,
111  const Elem & side,
112  unsigned int totalorder,
113  Point & sidep)
114 {
115  // "vertex" nodes are now decoupled from vertices, so we have
116  // to order them consistently otherwise
117  if (side_i < 4)
118  side_i = remap_node<4>(side_i, side, 0);
119 
120  // And "edge" nodes are decoupled from edges, so we have to
121  // reorder them too!
122  else if (side_i < 4u*totalorder)
123  {
124  unsigned int side_node = (side_i - 4)/(totalorder-1)+4;
125  side_node = remap_node<4>(side_node, side, 4);
126  side_i = ((side_i - 4) % (totalorder - 1)) // old local edge_i
127  + 4 + (side_node-4)*(totalorder-1);
128  }
129 
130  // Interior dofs in 2D don't care about where xi/eta point in
131  // physical space, but here we need them to match from both
132  // sides of a face!
133  else
134  {
135  unsigned int min_side_node = remap_node<4>(0, side, 0);
136  const bool flip = (side.point(min_side_node) < side.point((min_side_node+1)%4));
137 
138  switch (min_side_node) {
139  case 0:
140  if (flip)
141  std::swap(sidep(0), sidep(1));
142  break;
143  case 1:
144  sidep(0) = -sidep(0);
145  if (!flip)
146  std::swap(sidep(0), sidep(1));
147  break;
148  case 2:
149  sidep(0) = -sidep(0);
150  sidep(1) = -sidep(1);
151  if (flip)
152  std::swap(sidep(0), sidep(1));
153  break;
154  case 3:
155  sidep(1) = -sidep(1);
156  if (!flip)
157  std::swap(sidep(0), sidep(1));
158  break;
159  default:
160  libmesh_error();
161  }
162  }
163 }
164 
165 
166 void cube_indices(const Elem * elem,
167  const unsigned int totalorder,
168  const unsigned int i,
169  Real & xi, Real & eta, Real & zeta,
170  unsigned int & i0,
171  unsigned int & i1,
172  unsigned int & i2)
173 {
174  // The only way to make any sense of this
175  // is to look at the mgflo/mg2/mgf documentation
176  // and make the cut-out cube!
177  // Example i0 and i1 values for totalorder = 3:
178  // FIXME - these examples are incorrect now that we've got truly
179  // hierarchic basis functions
180  // Nodes 0 1 2 3 4 5 6 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 20 20 21 21 21 21 22 22 22 22 23 23 23 23 24 24 24 24 25 25 25 25 26 26 26 26 26 26 26 26
181  // DOFS 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 18 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 60 62 63
182  // static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3};
183  // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3, 0, 0, 0, 0, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3};
184  // static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3};
185 
186  // the number of DoFs per edge appears everywhere:
187  const unsigned int e = totalorder - 1u;
188 
189  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+1u));
190 
191  Real xi_saved = xi, eta_saved = eta, zeta_saved = zeta;
192 
193  // Vertices:
194  if (i == 0)
195  {
196  i0 = 0;
197  i1 = 0;
198  i2 = 0;
199  }
200  else if (i == 1)
201  {
202  i0 = 1;
203  i1 = 0;
204  i2 = 0;
205  }
206  else if (i == 2)
207  {
208  i0 = 1;
209  i1 = 1;
210  i2 = 0;
211  }
212  else if (i == 3)
213  {
214  i0 = 0;
215  i1 = 1;
216  i2 = 0;
217  }
218  else if (i == 4)
219  {
220  i0 = 0;
221  i1 = 0;
222  i2 = 1;
223  }
224  else if (i == 5)
225  {
226  i0 = 1;
227  i1 = 0;
228  i2 = 1;
229  }
230  else if (i == 6)
231  {
232  i0 = 1;
233  i1 = 1;
234  i2 = 1;
235  }
236  else if (i == 7)
237  {
238  i0 = 0;
239  i1 = 1;
240  i2 = 1;
241  }
242  // Edge 0
243  else if (i < 8 + e)
244  {
245  i0 = i - 6;
246  i1 = 0;
247  i2 = 0;
248  if (elem->positive_edge_orientation(0))
249  xi = -xi_saved;
250  }
251  // Edge 1
252  else if (i < 8 + 2*e)
253  {
254  i0 = 1;
255  i1 = i - e - 6;
256  i2 = 0;
257  if (elem->positive_edge_orientation(1))
258  eta = -eta_saved;
259  }
260  // Edge 2
261  else if (i < 8 + 3*e)
262  {
263  i0 = i - 2*e - 6;
264  i1 = 1;
265  i2 = 0;
266  if (!elem->positive_edge_orientation(2))
267  xi = -xi_saved;
268  }
269  // Edge 3
270  else if (i < 8 + 4*e)
271  {
272  i0 = 0;
273  i1 = i - 3*e - 6;
274  i2 = 0;
275  if (elem->positive_edge_orientation(3))
276  eta = -eta_saved;
277  }
278  // Edge 4
279  else if (i < 8 + 5*e)
280  {
281  i0 = 0;
282  i1 = 0;
283  i2 = i - 4*e - 6;
284  if (elem->positive_edge_orientation(4))
285  zeta = -zeta_saved;
286  }
287  // Edge 5
288  else if (i < 8 + 6*e)
289  {
290  i0 = 1;
291  i1 = 0;
292  i2 = i - 5*e - 6;
293  if (elem->positive_edge_orientation(5))
294  zeta = -zeta_saved;
295  }
296  // Edge 6
297  else if (i < 8 + 7*e)
298  {
299  i0 = 1;
300  i1 = 1;
301  i2 = i - 6*e - 6;
302  if (elem->positive_edge_orientation(6))
303  zeta = -zeta_saved;
304  }
305  // Edge 7
306  else if (i < 8 + 8*e)
307  {
308  i0 = 0;
309  i1 = 1;
310  i2 = i - 7*e - 6;
311  if (elem->positive_edge_orientation(7))
312  zeta = -zeta_saved;
313  }
314  // Edge 8
315  else if (i < 8 + 9*e)
316  {
317  i0 = i - 8*e - 6;
318  i1 = 0;
319  i2 = 1;
320  if (elem->positive_edge_orientation(8))
321  xi = -xi_saved;
322  }
323  // Edge 9
324  else if (i < 8 + 10*e)
325  {
326  i0 = 1;
327  i1 = i - 9*e - 6;
328  i2 = 1;
329  if (elem->positive_edge_orientation(9))
330  eta = -eta_saved;
331  }
332  // Edge 10
333  else if (i < 8 + 11*e)
334  {
335  i0 = i - 10*e - 6;
336  i1 = 1;
337  i2 = 1;
338  if (!elem->positive_edge_orientation(10))
339  xi = -xi_saved;
340  }
341  // Edge 11
342  else if (i < 8 + 12*e)
343  {
344  i0 = 0;
345  i1 = i - 11*e - 6;
346  i2 = 1;
347  if (elem->positive_edge_orientation(11))
348  eta = -eta_saved;
349  }
350  // Face 0
351  else if (i < 8 + 12*e + e*e)
352  {
353  unsigned int basisnum = i - 8 - 12*e;
354  i0 = square_number_row[basisnum] + 2;
355  i1 = square_number_column[basisnum] + 2;
356  i2 = 0;
357  const Point min_point = get_min_point(elem, 1, 2, 0, 3);
358 
359  if (elem->point(0) == min_point)
360  if (elem->positive_face_orientation(0))
361  {
362  // Case 1
363  xi = xi_saved;
364  eta = eta_saved;
365  }
366  else
367  {
368  // Case 2
369  xi = eta_saved;
370  eta = xi_saved;
371  }
372 
373  else if (elem->point(3) == min_point)
374  if (elem->positive_face_orientation(0))
375  {
376  // Case 3
377  xi = -eta_saved;
378  eta = xi_saved;
379  }
380  else
381  {
382  // Case 4
383  xi = xi_saved;
384  eta = -eta_saved;
385  }
386 
387  else if (elem->point(2) == min_point)
388  if (elem->positive_face_orientation(0))
389  {
390  // Case 5
391  xi = -xi_saved;
392  eta = -eta_saved;
393  }
394  else
395  {
396  // Case 6
397  xi = -eta_saved;
398  eta = -xi_saved;
399  }
400 
401  else if (elem->point(1) == min_point)
402  {
403  if (elem->positive_face_orientation(0))
404  {
405  // Case 7
406  xi = eta_saved;
407  eta = -xi_saved;
408  }
409  else
410  {
411  // Case 8
412  xi = -xi_saved;
413  eta = eta_saved;
414  }
415  }
416  }
417  // Face 1
418  else if (i < 8 + 12*e + 2*e*e)
419  {
420  unsigned int basisnum = i - 8 - 12*e - e*e;
421  i0 = square_number_row[basisnum] + 2;
422  i1 = 0;
423  i2 = square_number_column[basisnum] + 2;
424  const Point min_point = get_min_point(elem, 0, 1, 5, 4);
425 
426  if (elem->point(0) == min_point)
427  if (!elem->positive_face_orientation(1))
428  {
429  // Case 1
430  xi = xi_saved;
431  zeta = zeta_saved;
432  }
433  else
434  {
435  // Case 2
436  xi = zeta_saved;
437  zeta = xi_saved;
438  }
439 
440  else if (elem->point(1) == min_point)
441  if (!elem->positive_face_orientation(1))
442  {
443  // Case 3
444  xi = zeta_saved;
445  zeta = -xi_saved;
446  }
447  else
448  {
449  // Case 4
450  xi = -xi_saved;
451  zeta = zeta_saved;
452  }
453 
454  else if (elem->point(5) == min_point)
455  if (!elem->positive_face_orientation(1))
456  {
457  // Case 5
458  xi = -xi_saved;
459  zeta = -zeta_saved;
460  }
461  else
462  {
463  // Case 6
464  xi = -zeta_saved;
465  zeta = -xi_saved;
466  }
467 
468  else if (elem->point(4) == min_point)
469  {
470  if (!elem->positive_face_orientation(1))
471  {
472  // Case 7
473  xi = -xi_saved;
474  zeta = zeta_saved;
475  }
476  else
477  {
478  // Case 8
479  xi = xi_saved;
480  zeta = -zeta_saved;
481  }
482  }
483  }
484  // Face 2
485  else if (i < 8 + 12*e + 3*e*e)
486  {
487  unsigned int basisnum = i - 8 - 12*e - 2*e*e;
488  i0 = 1;
489  i1 = square_number_row[basisnum] + 2;
490  i2 = square_number_column[basisnum] + 2;
491  const Point min_point = get_min_point(elem, 1, 2, 6, 5);
492 
493  if (elem->point(1) == min_point)
494  if (!elem->positive_face_orientation(2))
495  {
496  // Case 1
497  eta = eta_saved;
498  zeta = zeta_saved;
499  }
500  else
501  {
502  // Case 2
503  eta = zeta_saved;
504  zeta = eta_saved;
505  }
506 
507  else if (elem->point(2) == min_point)
508  if (!elem->positive_face_orientation(2))
509  {
510  // Case 3
511  eta = zeta_saved;
512  zeta = -eta_saved;
513  }
514  else
515  {
516  // Case 4
517  eta = -eta_saved;
518  zeta = zeta_saved;
519  }
520 
521  else if (elem->point(6) == min_point)
522  if (!elem->positive_face_orientation(2))
523  {
524  // Case 5
525  eta = -eta_saved;
526  zeta = -zeta_saved;
527  }
528  else
529  {
530  // Case 6
531  eta = -zeta_saved;
532  zeta = -eta_saved;
533  }
534 
535  else if (elem->point(5) == min_point)
536  {
537  if (!elem->positive_face_orientation(2))
538  {
539  // Case 7
540  eta = -zeta_saved;
541  zeta = eta_saved;
542  }
543  else
544  {
545  // Case 8
546  eta = eta_saved;
547  zeta = -zeta_saved;
548  }
549  }
550  }
551  // Face 3
552  else if (i < 8 + 12*e + 4*e*e)
553  {
554  unsigned int basisnum = i - 8 - 12*e - 3*e*e;
555  i0 = square_number_row[basisnum] + 2;
556  i1 = 1;
557  i2 = square_number_column[basisnum] + 2;
558  const Point min_point = get_min_point(elem, 2, 3, 7, 6);
559 
560  if (elem->point(3) == min_point)
561  if (elem->positive_face_orientation(3))
562  {
563  // Case 1
564  xi = xi_saved;
565  zeta = zeta_saved;
566  }
567  else
568  {
569  // Case 2
570  xi = zeta_saved;
571  zeta = xi_saved;
572  }
573 
574  else if (elem->point(7) == min_point)
575  if (elem->positive_face_orientation(3))
576  {
577  // Case 3
578  xi = -zeta_saved;
579  zeta = xi_saved;
580  }
581  else
582  {
583  // Case 4
584  xi = xi_saved;
585  zeta = -zeta_saved;
586  }
587 
588  else if (elem->point(6) == min_point)
589  if (elem->positive_face_orientation(3))
590  {
591  // Case 5
592  xi = -xi_saved;
593  zeta = -zeta_saved;
594  }
595  else
596  {
597  // Case 6
598  xi = -zeta_saved;
599  zeta = -xi_saved;
600  }
601 
602  else if (elem->point(2) == min_point)
603  {
604  if (elem->positive_face_orientation(3))
605  {
606  // Case 7
607  xi = zeta_saved;
608  zeta = -xi_saved;
609  }
610  else
611  {
612  // Case 8
613  xi = -xi_saved;
614  zeta = zeta_saved;
615  }
616  }
617  }
618  // Face 4
619  else if (i < 8 + 12*e + 5*e*e)
620  {
621  unsigned int basisnum = i - 8 - 12*e - 4*e*e;
622  i0 = 0;
623  i1 = square_number_row[basisnum] + 2;
624  i2 = square_number_column[basisnum] + 2;
625  const Point min_point = get_min_point(elem, 3, 0, 4, 7);
626 
627  if (elem->point(0) == min_point)
628  if (elem->positive_face_orientation(4))
629  {
630  // Case 1
631  eta = eta_saved;
632  zeta = zeta_saved;
633  }
634  else
635  {
636  // Case 2
637  eta = zeta_saved;
638  zeta = eta_saved;
639  }
640 
641  else if (elem->point(4) == min_point)
642  if (elem->positive_face_orientation(4))
643  {
644  // Case 3
645  eta = -zeta_saved;
646  zeta = eta_saved;
647  }
648  else
649  {
650  // Case 4
651  eta = eta_saved;
652  zeta = -zeta_saved;
653  }
654 
655  else if (elem->point(7) == min_point)
656  if (elem->positive_face_orientation(4))
657  {
658  // Case 5
659  eta = -eta_saved;
660  zeta = -zeta_saved;
661  }
662  else
663  {
664  // Case 6
665  eta = -zeta_saved;
666  zeta = -eta_saved;
667  }
668 
669  else if (elem->point(3) == min_point)
670  {
671  if (elem->positive_face_orientation(4))
672  {
673  // Case 7
674  eta = zeta_saved;
675  zeta = -eta_saved;
676  }
677  else
678  {
679  // Case 8
680  eta = -eta_saved;
681  zeta = zeta_saved;
682  }
683  }
684  }
685  // Face 5
686  else if (i < 8 + 12*e + 6*e*e)
687  {
688  unsigned int basisnum = i - 8 - 12*e - 5*e*e;
689  i0 = square_number_row[basisnum] + 2;
690  i1 = square_number_column[basisnum] + 2;
691  i2 = 1;
692  const Point min_point = get_min_point(elem, 4, 5, 6, 7);
693 
694  if (elem->point(4) == min_point)
695  if (!elem->positive_face_orientation(5))
696  {
697  // Case 1
698  xi = xi_saved;
699  eta = eta_saved;
700  }
701  else
702  {
703  // Case 2
704  xi = eta_saved;
705  eta = xi_saved;
706  }
707 
708  else if (elem->point(5) == min_point)
709  if (!elem->positive_face_orientation(5))
710  {
711  // Case 3
712  xi = eta_saved;
713  eta = -xi_saved;
714  }
715  else
716  {
717  // Case 4
718  xi = -xi_saved;
719  eta = eta_saved;
720  }
721 
722  else if (elem->point(6) == min_point)
723  if (!elem->positive_face_orientation(5))
724  {
725  // Case 5
726  xi = -xi_saved;
727  eta = -eta_saved;
728  }
729  else
730  {
731  // Case 6
732  xi = -eta_saved;
733  eta = -xi_saved;
734  }
735 
736  else if (elem->point(7) == min_point)
737  {
738  if (!elem->positive_face_orientation(5))
739  {
740  // Case 7
741  xi = -eta_saved;
742  eta = xi_saved;
743  }
744  else
745  {
746  // Case 8
747  xi = xi_saved;
748  eta = eta_saved;
749  }
750  }
751  }
752 
753  // Internal DoFs
754  else
755  {
756  unsigned int basisnum = i - 8 - 12*e - 6*e*e;
757  i0 = cube_number_column[basisnum] + 2;
758  i1 = cube_number_row[basisnum] + 2;
759  i2 = cube_number_page[basisnum] + 2;
760  }
761 }
762 
763 
764 void prism_indices(const Elem * elem,
765  const unsigned int totalorder,
766  const unsigned int i,
767  Point & xi_eta, Real & zeta,
768  unsigned int & i01,
769  unsigned int & i2)
770 {
771  // the number of DoFs per edge appears everywhere:
772  const unsigned int e = totalorder - 1u;
773 
774  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+2u)/2u);
775 
776  Point xi_eta_saved = xi_eta;
777  Real zeta_saved = zeta;
778 
779  // Vertices:
780  if (i == 0)
781  {
782  i01 = 0;
783  i2 = 0;
784  }
785  else if (i == 1)
786  {
787  i01 = 1;
788  i2 = 0;
789  }
790  else if (i == 2)
791  {
792  i01 = 2;
793  i2 = 0;
794  }
795  else if (i == 3)
796  {
797  i01 = 0;
798  i2 = 1;
799  }
800  else if (i == 4)
801  {
802  i01 = 1;
803  i2 = 1;
804  }
805  else if (i == 5)
806  {
807  i01 = 2;
808  i2 = 1;
809  }
810  // Edges 0,1,2 (vertices 6,7,8)
811  else if (i < 6 + 3*e)
812  {
813  // The TRI code will handle any flips here
814  i01 = i - 3;
815  i2 = 0;
816  }
817  // Edge 3,4,5 (vertices 9,10,11)
818  else if (i < 6 + 6*e)
819  {
820  i01 = (i - 6 - 3*e)/e; // which tri DoF are we?
821  i2 = (i - 6 - 3*e)%e+2; // edge DoF? +2 to skip endpoints
822  // EDGE evaluations don't flip, so handle that here
823  if (elem->positive_edge_orientation(i01+3))
824  zeta = -zeta;
825  }
826  // Edge 6,7,8 (vertices 12,13,14)
827  else if (i < 6 + 9*e)
828  {
829  // The TRI code will handle any flips here
830  i01 = i - 3 - 6*e;
831  i2 = 1;
832  }
833  // Face 1, node 15 (*before* 0, via node 18 on prism20)
834  else if (i < 6 + 9*e + e*e)
835  {
836  unsigned int basisnum = i - 6 - 9*e;
837 
838  // How wide is the stretch from one side to the other of the
839  // line in the xi-eta plane parallel to this face?
840  const Real xe_scale = 1 - xi_eta_saved(1);
841 
842  // What percentage of the way along that stretch are we?
843  const Real xe_fraction = (xe_scale==0) ?
844  0 : xi_eta_saved(0)/xe_scale;
845 
846  // indexes in edge numbering
847  unsigned int s0 = square_number_row[basisnum] + 2;
848  unsigned int s1 = square_number_column[basisnum] + 2;
849  const Point min_point = get_min_point(elem, 0, 1, 3, 4);
850 
851  if (elem->point(0) == min_point)
852  {
853  if (!elem->positive_face_orientation(1))
854  {
855  // Case 1: no flips needed
856  i01 = s0+1; // edge to triangle side 0 numbering
857  i2 = s1;
858  }
859  else
860  {
861  // Case 2: flip about 0-4 diagonal
862  i01 = s1+1;
863  i2 = s0;
864  zeta = 2*xe_fraction-1;
865  xi_eta(0) = (zeta_saved+1)*xe_scale/2;
866  }
867  }
868  else if (elem->point(3) == min_point)
869  {
870  if (!elem->positive_face_orientation(1))
871  {
872  // Case 3: 0->3->4->1->0 rotation
873  i01 = s1+1;
874  i2 = s0;
875  zeta = 1-2*xe_fraction;
876  xi_eta(0) = (zeta_saved+1)*xe_scale/2;
877  }
878  else
879  {
880  // Case 4: flip about 9-10 midline
881  i01 = s0+1;
882  i2 = s1;
883  zeta = -zeta_saved;
884  }
885  }
886  else if (elem->point(1) == min_point)
887  {
888  if (!elem->positive_face_orientation(1))
889  {
890  // Case 5: 0->1->4->3->0 rotation
891  i01 = s1+1;
892  i2 = s0;
893  zeta = 2*xe_fraction-1;
894  xi_eta(0) = (1-zeta_saved)*xe_scale/2;
895  }
896  else
897  {
898  // Case 6: flip about 6-12 midline
899  i01 = s0+1;
900  i2 = s1;
901  xi_eta(0) = (1-xe_fraction)*xe_scale;
902  }
903  }
904  else if (elem->point(4) == min_point)
905  {
906  if (!elem->positive_face_orientation(1))
907  {
908  // Case 7: 180 degree rotation
909  i01 = s0+1;
910  i2 = s1;
911  xi_eta(0) = (1-xe_fraction)*xe_scale;
912  zeta = -zeta_saved;
913  }
914  else
915  {
916  // Case 8: flip about 1-3 diagonal
917  i01 = s1+1;
918  i2 = s0;
919  zeta = 1-2*xe_fraction;
920  xi_eta(0) = (1-zeta_saved)*xe_scale/2;
921  }
922  }
923  }
924  // Face 2, node 16
925  else if (i < 6 + 9*e + 2*e*e)
926  {
927  unsigned int basisnum = i - 6 - 9*e - e*e;
928 
929  // How wide is the stretch from one side to the other of the
930  // line in the xi-eta plane parallel to this face?
931  const Real xe_scale = xi_eta_saved(0) + xi_eta_saved(1);
932 
933  // What percentage of the way along that stretch are we?
934  const Real xe_fraction = (xe_scale==0) ?
935  0 : xi_eta_saved(0)/xe_scale;
936 
937  // indexes in edge numbering
938  unsigned int s0 = square_number_row[basisnum] + 2;
939  unsigned int s1 = square_number_column[basisnum] + 2;
940  const Point min_point = get_min_point(elem, 1, 2, 4, 5);
941 
942  if (elem->point(1) == min_point)
943  {
944  if (!elem->positive_face_orientation(2))
945  {
946  // Case 1: no flips needed
947  i01 = s0+1+3; // edge to triangle side 1 numbering
948  i2 = s1;
949  }
950  else
951  {
952  // Case 2: flip about 1-5 diagonal
953  i01 = s1+1+e;
954  i2 = s0;
955  zeta = 2*xe_fraction-1;
956  const Real xe = (zeta_saved+1)/2;
957  xi_eta(1) = xe*xe_scale;
958  xi_eta(0) = xe_scale - xi_eta(1);
959  }
960  }
961  else if (elem->point(4) == min_point)
962  {
963  if (!elem->positive_face_orientation(2))
964  {
965  // Case 3: 1->4->5->2->1 rotation
966  i01 = s1+1+e;
967  i2 = s0;
968  zeta = 1-2*xe_fraction;
969  const Real xe = (zeta_saved+1)/2;
970  xi_eta(1) = xe*xe_scale;
971  xi_eta(0) = xe_scale - xi_eta(1);
972  }
973  else
974  {
975  // Case 4: flip about 10-11 midline
976  i01 = s0+1+e;
977  i2 = s1;
978  zeta = -zeta_saved;
979  }
980  }
981  else if (elem->point(2) == min_point)
982  {
983  if (!elem->positive_face_orientation(2))
984  {
985  // Case 5: 1->2->5->4->1 rotation
986  i01 = s1+1+e;
987  i2 = s0;
988  zeta = 2*xe_fraction-1;
989  const Real xe = (1-zeta_saved)/2;
990  xi_eta(1) = xe*xe_scale;
991  xi_eta(0) = xe_scale - xi_eta(1);
992  }
993  else
994  {
995  // Case 6: flip about 7-13 midline
996  i01 = s0+1+e;
997  i2 = s1;
998  const Real xe = (1-xe_fraction);
999  xi_eta(1) = xe*xe_scale;
1000  xi_eta(0) = xe_scale - xi_eta(1);
1001  }
1002  }
1003  else if (elem->point(5) == min_point)
1004  {
1005  if (!elem->positive_face_orientation(2))
1006  {
1007  // Case 7: 180 degree rotation
1008  i01 = s0+1+e;
1009  i2 = s1;
1010  zeta = -zeta_saved;
1011  const Real xe = (1-xe_fraction);
1012  xi_eta(1) = xe*xe_scale;
1013  xi_eta(0) = xe_scale - xi_eta(1);
1014  }
1015  else
1016  {
1017  // Case 8: flip about 1-3 diagonal
1018  i01 = s1+1+e;
1019  i2 = s0;
1020  zeta = 1-2*xe_fraction;
1021  const Real xe = (1-zeta_saved)/2;
1022  xi_eta(1) = xe*xe_scale;
1023  xi_eta(0) = xe_scale - xi_eta(1);
1024  }
1025  }
1026  }
1027  // Face 3, node 17
1028  else if (i < 6 + 9*e + 3*e*e)
1029  {
1030  unsigned int basisnum = i - 6 - 9*e - 2*e*e;
1031 
1032  // How wide is the stretch from one side to the other of the
1033  // line in the xi-eta plane parallel to this face?
1034  const Real xe_scale = 1 - xi_eta_saved(0);
1035 
1036  // What percentage of the way along that stretch are we?
1037  const Real xe_fraction = (xe_scale==0) ?
1038  0 : (xe_scale - xi_eta_saved(1))/xe_scale;
1039 
1040  // indexes in edge numbering
1041  unsigned int s0 = square_number_row[basisnum] + 2;
1042  unsigned int s1 = square_number_column[basisnum] + 2;
1043  const Point min_point = get_min_point(elem, 0, 2, 3, 5);
1044 
1045  if (elem->point(2) == min_point)
1046  {
1047  if (!elem->positive_face_orientation(3))
1048  {
1049  // Case 1: no flips needed
1050  i01 = s0+1+2*e; // edge to triangle side 2 numbering
1051  i2 = s1;
1052  }
1053  else
1054  {
1055  // Case 2: flip about 2-3 diagonal
1056  i01 = s1+1+2*e;
1057  i2 = s0;
1058  zeta = 2*xe_fraction-1;
1059  const Real xe = (zeta_saved+1)/2;
1060  xi_eta(1) = xe_scale - xe*xe_scale;
1061  }
1062  }
1063  else if (elem->point(5) == min_point)
1064  {
1065  if (!elem->positive_face_orientation(3))
1066  {
1067  // Case 3: 2->5->3->0->2 rotation
1068  i01 = s1+1+2*e;
1069  i2 = s0;
1070  zeta = 1-2*xe_fraction;
1071  const Real xe = (zeta_saved+1)/2;
1072  xi_eta(1) = xe_scale - xe*xe_scale;
1073  }
1074  else
1075  {
1076  // Case 4: flip about 11-9 midline
1077  i01 = s0+1+2*e;
1078  i2 = s1;
1079  zeta = -zeta_saved;
1080  }
1081  }
1082  else if (elem->point(0) == min_point)
1083  {
1084  if (!elem->positive_face_orientation(3))
1085  {
1086  // Case 5: 2->0->3->5->2 rotation
1087  i01 = s1+1+2*e;
1088  i2 = s0;
1089  zeta = 2*xe_fraction-1;
1090  const Real xe = (1-zeta_saved)/2;
1091  xi_eta(1) = xe_scale - xe*xe_scale;
1092  }
1093  else
1094  {
1095  // Case 6: flip about 8-14 midline
1096  i01 = s0+1+2*e;
1097  i2 = s1;
1098  const Real xe = (1-xe_fraction);
1099  xi_eta(1) = xe_scale - xe*xe_scale;
1100  }
1101  }
1102  else if (elem->point(3) == min_point)
1103  {
1104  if (!elem->positive_face_orientation(3))
1105  {
1106  // Case 7: 180 degree rotation
1107  i01 = s0+1+2*e;
1108  i2 = s1;
1109  zeta = -zeta_saved;
1110  const Real xe = (1-xe_fraction);
1111  xi_eta(1) = xe_scale - xe*xe_scale;
1112  }
1113  else
1114  {
1115  // Case 8: flip about 0-5 diagonal
1116  i01 = s1+1+2*e;
1117  i2 = s0;
1118  zeta = 1-2*xe_fraction;
1119  const Real xe = (1-zeta_saved)/2;
1120  xi_eta(1) = xe_scale - xe*xe_scale;
1121  }
1122  }
1123  }
1124  // Face 0, node 18 - node order due to hierarchic numbering
1125  else if (i < 6 + 9*e + 3*e*e + e*(e-1)/2)
1126  {
1127  // The TRI code will handle any flips here
1128  i01 = i - 3 - 6*e - 3*e*e;
1129  i2 = 0;
1130  }
1131  // Face 4
1132  else if (i < 6 + 9*e + 3*e*e + e*(e-1))
1133  {
1134  // The TRI code will handle any flips here
1135  i01 = i - 3 - 6*e - 3*e*e - e*(e-1)/2;
1136  i2 = 1;
1137  }
1138  // Internal DoFs
1139  else
1140  {
1141  // We won't bother with any internal DoF reordering / flipping;
1142  // that's fine unless we ever get to 4D.
1143  unsigned int basisnum = i - 6 - 9*e - 3*e*e - e*(e-1);
1144  i01 = prism_number_triangle[basisnum] + 3 + 3*e;
1145  i2 = prism_number_page[basisnum] + 2;
1146  }
1147 }
1148 
1149 #endif // LIBMESH_DIM > 2
1150 
1151 } // end anonymous namespace
1152 
1153 
1154 
1155 namespace libMesh
1156 {
1157 
1158 
1159 LIBMESH_DEFAULT_VECTORIZED_FE(3,HIERARCHIC)
1160 LIBMESH_DEFAULT_VECTORIZED_FE(3,L2_HIERARCHIC)
1161 LIBMESH_DEFAULT_VECTORIZED_FE(3,SIDE_HIERARCHIC)
1162 
1163 
1164 template <>
1165 Real FE<3,HIERARCHIC>::shape(const ElemType,
1166  const Order,
1167  const unsigned int,
1168  const Point &)
1169 {
1170  libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1171  return 0.;
1172 }
1173 
1174 
1175 
1176 template <>
1177 Real FE<3,L2_HIERARCHIC>::shape(const ElemType,
1178  const Order,
1179  const unsigned int,
1180  const Point &)
1181 {
1182  libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1183  return 0.;
1184 }
1185 
1186 
1187 
1188 template <>
1189 Real FE<3,SIDE_HIERARCHIC>::shape(const ElemType,
1190  const Order,
1191  const unsigned int,
1192  const Point &)
1193 {
1194  libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1195  return 0.;
1196 }
1197 
1198 
1199 
1200 template <>
1201 Real FE<3,HIERARCHIC>::shape(const Elem * elem,
1202  const Order order,
1203  const unsigned int i,
1204  const Point & p,
1205  const bool add_p_level)
1206 {
1207  return fe_hierarchic_3D_shape<HIERARCHIC>(elem, order, i, p, add_p_level);
1208 }
1209 
1210 
1211 template <>
1212 Real FE<3,HIERARCHIC>::shape(const FEType fet,
1213  const Elem * elem,
1214  const unsigned int i,
1215  const Point & p,
1216  const bool add_p_level)
1217 {
1218  return fe_hierarchic_3D_shape<HIERARCHIC>(elem, fet.order, i, p, add_p_level);
1219 }
1220 
1221 
1222 
1223 
1224 template <>
1225 Real FE<3,L2_HIERARCHIC>::shape(const Elem * elem,
1226  const Order order,
1227  const unsigned int i,
1228  const Point & p,
1229  const bool add_p_level)
1230 {
1231  return fe_hierarchic_3D_shape<L2_HIERARCHIC>(elem, order, i, p, add_p_level);
1232 }
1233 
1234 
1235 template <>
1236 Real FE<3,L2_HIERARCHIC>::shape(const FEType fet,
1237  const Elem * elem,
1238  const unsigned int i,
1239  const Point & p,
1240  const bool add_p_level)
1241 {
1242  return fe_hierarchic_3D_shape<L2_HIERARCHIC>(elem, fet.order, i, p, add_p_level);
1243 }
1244 
1245 
1246 
1247 template <>
1248 Real FE<3,SIDE_HIERARCHIC>::shape(const Elem * elem,
1249  const Order order,
1250  const unsigned int i,
1251  const Point & p,
1252  const bool add_p_level)
1253 {
1254 #if LIBMESH_DIM == 3
1255  libmesh_assert(elem);
1256  const ElemType type = elem->type();
1257 
1258  const Order totalorder = order + add_p_level*elem->p_level();
1259 
1260  switch (type)
1261  {
1262  case HEX27:
1263  {
1264  const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+1u);
1265  libmesh_assert_less(i, 6*dofs_per_side);
1266 
1267  const unsigned int sidenum = cube_side(p);
1268  if (sidenum > 5)
1269  return std::numeric_limits<Real>::quiet_NaN();
1270 
1271  const unsigned int dof_offset = sidenum * dofs_per_side;
1272 
1273  if (i < dof_offset) // i is on a previous side
1274  return 0;
1275 
1276  if (i >= dof_offset + dofs_per_side) // i is on a later side
1277  return 0;
1278 
1279  if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
1280  return 1;
1281 
1282  unsigned int side_i = i - dof_offset;
1283 
1284  std::unique_ptr<const Elem> side = elem->build_side_ptr(sidenum);
1285 
1286  Point sidep = cube_side_point(sidenum, p);
1287 
1288  cube_remap(side_i, *side, totalorder, sidep);
1289 
1290  return FE<2,HIERARCHIC>::shape(side.get(), order, side_i, sidep, add_p_level);
1291  }
1292 
1293  case TET14:
1294  {
1295  const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+2u)/2u;
1296  libmesh_assert_less(i, 4*dofs_per_side);
1297 
1298  const Real zeta[4] = { Real(1.) - p(0) - p(1) - p(2), p(0), p(1), p(2) };
1299 
1300  unsigned int face_num = 0;
1301  if (zeta[0] > zeta[3] &&
1302  zeta[1] > zeta[3] &&
1303  zeta[2] > zeta[3])
1304  {
1305  face_num = 0;
1306  }
1307  else if (zeta[0] > zeta[2] &&
1308  zeta[1] > zeta[2] &&
1309  zeta[3] > zeta[2])
1310  {
1311  face_num = 1;
1312  }
1313  else if (zeta[1] > zeta[0] &&
1314  zeta[2] > zeta[0] &&
1315  zeta[3] > zeta[0])
1316  {
1317  face_num = 2;
1318  }
1319  else
1320  {
1321  // We'd better not be right between two faces
1322  libmesh_assert (zeta[0] > zeta[1] &&
1323  zeta[2] > zeta[1] &&
1324  zeta[3] > zeta[1]);
1325  face_num = 3;
1326  }
1327 
1328  if (i < face_num * dofs_per_side ||
1329  i >= (face_num+1) * dofs_per_side)
1330  return 0;
1331 
1332  if (totalorder == 0)
1333  return 1;
1334 
1335  const std::array<unsigned int, 3> face_vertex =
1336  oriented_tet_nodes(*elem, face_num);
1337 
1338  // We only need a Tri3 to evaluate L2_HIERARCHIC on the affine
1339  // master element
1340  Tri3 side;
1341 
1342  // We pinky swear not to modify these nodes
1343  Elem & e = const_cast<Elem &>(*elem);
1344  side.set_node(0, e.node_ptr(face_vertex[0]));
1345  side.set_node(1, e.node_ptr(face_vertex[1]));
1346  side.set_node(2, e.node_ptr(face_vertex[2]));
1347 
1348  const unsigned int basisnum = i - face_num*dofs_per_side;
1349 
1350  Point sidep {zeta[face_vertex[1]], zeta[face_vertex[2]]};
1351 
1352  return FE<2,L2_HIERARCHIC>::shape(&side, totalorder,
1353  basisnum, sidep, false);
1354  }
1355 
1356  case PRISM20:
1357  case PRISM21:
1358  {
1359  const unsigned int dofs_per_quad = (totalorder+1u)*(totalorder+1u);
1360  const unsigned int dofs_per_tri = (totalorder+1u)*(totalorder+2u)/2u;
1361  libmesh_assert_less(i, 3*dofs_per_quad + 2*dofs_per_tri);
1362 
1363  // We only need a Tri3 or Quad4 to evaluate L2_HIERARCHIC on
1364  // the affine master element
1365  Tri3 tri;
1366  Quad4 quad;
1367  Elem * side = &quad;
1368  unsigned int dofs_on_side = dofs_per_quad;
1369 
1370  // We pinky swear not to modify the nodes we'll point to
1371  Elem & e = const_cast<Elem &>(*elem);
1372 
1373  Point sidep;
1374 
1375  // Face number calculation is tricky - the ordering of side
1376  // nodes on Prisms does *not* match the ordering of sides!
1377  // (the mid-triangle side nodes were added "later")
1378  // Here face_num will be the numbering that matches the side
1379  // number, but i_offset will have to consider the nodal
1380  // ordering.
1381  unsigned int face_num = 0;
1382  unsigned int i_offset = 0;
1383 
1384  // Triangular coordinates
1385  const Real zeta[3] = { Real(1.) - p(0) - p(1), p(0), p(1) };
1386 
1387  // Closeness to midplane
1388  const Real zmid = 1 - std::abs(p(2));
1389 
1390  if (zeta[1] > zeta[2] && zeta[0] > zeta[2] &&
1391  zmid > 3*zeta[2]) // face 1, quad
1392  {
1393  face_num = 1;
1394  i_offset = 0;
1395  }
1396  else if (zeta[1] > zeta[0] && zeta[2] > zeta[0] &&
1397  zmid > 3*zeta[0]) // face 2, quad
1398  {
1399  face_num = 2;
1400  i_offset = dofs_per_quad;
1401  }
1402  else if (zeta[0] > zeta[1] && zeta[2] > zeta[1] &&
1403  zmid > 3*zeta[1]) // face 3, quad
1404  {
1405  face_num = 3;
1406  i_offset = 2*dofs_per_quad;
1407  }
1408  else if (p(2) + 1 < 3*zeta[0] &&
1409  p(2) + 1 < 3*zeta[1] &&
1410  p(2) + 1 < 3*zeta[2]) // face 0, tri
1411  {
1412  face_num = 0;
1413  i_offset = 3*dofs_per_quad;
1414  dofs_on_side = dofs_per_tri;
1415  side = &tri;
1416  }
1417  else if (1 - p(2) < 3*zeta[0] &&
1418  1 - p(2) < 3*zeta[1] &&
1419  1 - p(2) < 3*zeta[2]) // face 4, tri
1420  {
1421  face_num = 4;
1422  i_offset = dofs_per_tri + 3*dofs_per_quad;
1423  dofs_on_side = dofs_per_tri;
1424  side = &tri;
1425  }
1426  else
1427  {
1428  libmesh_error_msg("Evaluating SIDE_HIERARCHIC right between two Prism faces?");
1429  }
1430 
1431  if (i < i_offset ||
1432  i >= i_offset + dofs_on_side)
1433  return 0;
1434 
1435  if (totalorder == 0)
1436  return 1;
1437 
1438  const std::array<unsigned int, 4> face_vertex =
1439  oriented_prism_nodes(*elem, face_num);
1440 
1441  side->set_node(0, e.node_ptr(face_vertex[0]));
1442  side->set_node(1, e.node_ptr(face_vertex[1]));
1443  side->set_node(2, e.node_ptr(face_vertex[2]));
1444  if (face_vertex[3] < 21)
1445  side->set_node(3, e.node_ptr(face_vertex[3]));
1446 
1447  if (face_num == 0 || face_num == 4)
1448  sidep = {zeta[face_vertex[1]%3], zeta[face_vertex[2]%3]};
1449  else
1450  {
1451  // Transform a coordinate from the master prism to the
1452  // master quad, based on two vertex indices defining the
1453  // coordinate's direction
1454  auto coord_val = [p](int v1, int v2){
1455  if (v2-v1 == 3)
1456  return p(2);
1457  else if (v2-v1 == -3)
1458  return -p(2);
1459  else if (v1%3 == 0 && v2%3 == 1)
1460  return 2*p(0)-1;
1461  else if (v2%3 == 0 && v1%3 == 1)
1462  return 1-2*p(0);
1463  else if (v1%3 == 1 && v2%3 == 2)
1464  return p(1)-p(0);
1465  else if (v2%3 == 1 && v1%3 == 2)
1466  return p(0)-p(1);
1467  else if (v1%3 == 2 && v2%3 == 0)
1468  return 1-2*p(1);
1469  else if (v2%3 == 2 && v1%3 == 0)
1470  return 2*p(1)-1;
1471  else
1472  libmesh_error();
1473  };
1474 
1475  sidep = {coord_val(face_vertex[0], face_vertex[1]),
1476  coord_val(face_vertex[0], face_vertex[3])};
1477  }
1478 
1479  const unsigned int basisnum = i - i_offset;
1480 
1481  return FE<2,L2_HIERARCHIC>::shape(side, totalorder,
1482  basisnum, sidep, false);
1483  }
1484 
1485 
1486  default:
1487  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1488  }
1489 
1490 #else // LIBMESH_DIM != 3
1491  libmesh_ignore(elem, order, i, p, add_p_level);
1492  libmesh_not_implemented();
1493 #endif
1494 }
1495 
1496 
1497 template <>
1498 Real FE<3,SIDE_HIERARCHIC>::shape(const FEType fet,
1499  const Elem * elem,
1500  const unsigned int i,
1501  const Point & p,
1502  const bool add_p_level)
1503 {
1504  return FE<3,SIDE_HIERARCHIC>::shape(elem,fet.order, i, p, add_p_level);
1505 }
1506 
1507 
1508 template <>
1509 Real FE<3,HIERARCHIC>::shape_deriv(const ElemType,
1510  const Order,
1511  const unsigned int,
1512  const unsigned int,
1513  const Point & )
1514 {
1515  libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1516  return 0.;
1517 }
1518 
1519 
1520 
1521 template <>
1522 Real FE<3,L2_HIERARCHIC>::shape_deriv(const ElemType,
1523  const Order,
1524  const unsigned int,
1525  const unsigned int,
1526  const Point & )
1527 {
1528  libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1529  return 0.;
1530 }
1531 
1532 
1533 
1534 template <>
1535 Real FE<3,SIDE_HIERARCHIC>::shape_deriv(const ElemType,
1536  const Order,
1537  const unsigned int,
1538  const unsigned int,
1539  const Point & )
1540 {
1541  libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1542  return 0.;
1543 }
1544 
1545 
1546 
1547 template <>
1548 Real FE<3,HIERARCHIC>::shape_deriv(const Elem * elem,
1549  const Order order,
1550  const unsigned int i,
1551  const unsigned int j,
1552  const Point & p,
1553  const bool add_p_level)
1554 {
1555  return fe_hierarchic_3D_shape_deriv<HIERARCHIC>(elem, order, i, j, p, add_p_level);
1556 }
1557 
1558 
1559 template <>
1560 Real FE<3,HIERARCHIC>::shape_deriv(const FEType fet,
1561  const Elem * elem,
1562  const unsigned int i,
1563  const unsigned int j,
1564  const Point & p,
1565  const bool add_p_level)
1566 {
1567  return fe_hierarchic_3D_shape_deriv<HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
1568 }
1569 
1570 
1571 
1572 template <>
1573 Real FE<3,L2_HIERARCHIC>::shape_deriv(const Elem * elem,
1574  const Order order,
1575  const unsigned int i,
1576  const unsigned int j,
1577  const Point & p,
1578  const bool add_p_level)
1579 {
1580  return fe_hierarchic_3D_shape_deriv<L2_HIERARCHIC>(elem, order, i, j, p, add_p_level);
1581 }
1582 
1583 
1584 template <>
1585 Real FE<3,L2_HIERARCHIC>::shape_deriv(const FEType fet,
1586  const Elem * elem,
1587  const unsigned int i,
1588  const unsigned int j,
1589  const Point & p,
1590  const bool add_p_level)
1591 {
1592  return fe_hierarchic_3D_shape_deriv<L2_HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
1593 }
1594 
1595 
1596 
1597 template <>
1598 Real FE<3,SIDE_HIERARCHIC>::shape_deriv(const Elem * elem,
1599  const Order order,
1600  const unsigned int i,
1601  const unsigned int j,
1602  const Point & p,
1603  const bool add_p_level)
1604 {
1605 #if LIBMESH_DIM == 3
1606  libmesh_assert(elem);
1607  const ElemType type = elem->type();
1608 
1609  const Order totalorder = order + add_p_level*elem->p_level();
1610 
1611  if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
1612  return 0; // constants have zero derivative
1613 
1614  switch (type)
1615  {
1616  case HEX27:
1617  {
1618  // I need to debug the p>2 case here...
1619  if (totalorder > 2)
1620  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,SIDE_HIERARCHIC>::shape);
1621 
1622  const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+1u);
1623  libmesh_assert_less(i, 6*dofs_per_side);
1624 
1625  const unsigned int sidenum = cube_side(p);
1626  if (sidenum > 5)
1627  return std::numeric_limits<Real>::quiet_NaN();
1628 
1629  const unsigned int dof_offset = sidenum * dofs_per_side;
1630 
1631  if (i < dof_offset) // i is on a previous side
1632  return 0;
1633 
1634  if (i >= dof_offset + dofs_per_side) // i is on a later side
1635  return 0;
1636 
1637  unsigned int side_i = i - dof_offset;
1638 
1639  std::unique_ptr<const Elem> side = elem->build_side_ptr(sidenum);
1640 
1641  Point sidep = cube_side_point(sidenum, p);
1642 
1643  cube_remap(side_i, *side, totalorder, sidep);
1644 
1645  // What direction on the side corresponds to the derivative
1646  // direction we want?
1647  unsigned int sidej = 100;
1648 
1649  // Do we need a -1 here to flip that direction?
1650  Real f = 1.;
1651 
1652  switch (j)
1653  {
1654  case 0: // d()/dxi
1655  {
1656  switch (sidenum)
1657  {
1658  case 0:
1659  sidej = 1;
1660  break;
1661  case 1:
1662  sidej = 0;
1663  break;
1664  case 2:
1665  return 0;
1666  case 3:
1667  sidej = 0;
1668  f = -1;
1669  break;
1670  case 4:
1671  return 0;
1672  case 5:
1673  sidej = 0;
1674  break;
1675  default:
1676  libmesh_error();
1677  }
1678  break;
1679  }
1680  case 1: // d()/deta
1681  {
1682  switch (sidenum)
1683  {
1684  case 0:
1685  sidej = 0;
1686  break;
1687  case 1:
1688  return 0;
1689  case 2:
1690  sidej = 0;
1691  break;
1692  case 3:
1693  return 0;
1694  case 4:
1695  sidej = 0;
1696  f = -1;
1697  break;
1698  case 5:
1699  sidej = 1;
1700  break;
1701  default:
1702  libmesh_error();
1703  }
1704  break;
1705  }
1706  case 2: // d()/dzeta
1707  {
1708  switch (sidenum)
1709  {
1710  case 0:
1711  return 0;
1712  case 1:
1713  case 2:
1714  case 3:
1715  case 4:
1716  sidej = 1;
1717  break;
1718  case 5:
1719  return 0;
1720  default:
1721  libmesh_error();
1722  }
1723  break;
1724  }
1725 
1726  default:
1727  libmesh_error_msg("Invalid derivative index j = " << j);
1728  }
1729 
1730  return f * FE<2,HIERARCHIC>::shape_deriv(side.get(), order,
1731  side_i, sidej, sidep,
1732  add_p_level);
1733  }
1734 
1735  case TET14:
1736  case PRISM20:
1737  case PRISM21:
1738  {
1739  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,SIDE_HIERARCHIC>::shape);
1740  }
1741 
1742  default:
1743  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1744  }
1745 
1746 #else // LIBMESH_DIM != 3
1747  libmesh_ignore(elem, order, i, j, p, add_p_level);
1748  libmesh_not_implemented();
1749 #endif
1750 }
1751 
1752 
1753 template <>
1754 Real FE<3,SIDE_HIERARCHIC>::shape_deriv(const FEType fet,
1755  const Elem * elem,
1756  const unsigned int i,
1757  const unsigned int j,
1758  const Point & p,
1759  const bool add_p_level)
1760 {
1761  return FE<3,SIDE_HIERARCHIC>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
1762 }
1763 
1764 
1765 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1766 
1767 template <>
1768 Real FE<3,HIERARCHIC>::shape_second_deriv(const ElemType,
1769  const Order,
1770  const unsigned int,
1771  const unsigned int,
1772  const Point & )
1773 {
1774  libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1775  return 0.;
1776 }
1777 
1778 
1779 
1780 template <>
1781 Real FE<3,L2_HIERARCHIC>::shape_second_deriv(const ElemType,
1782  const Order,
1783  const unsigned int,
1784  const unsigned int,
1785  const Point & )
1786 {
1787  libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1788  return 0.;
1789 }
1790 
1791 
1792 
1793 template <>
1794 Real FE<3,SIDE_HIERARCHIC>::shape_second_deriv(const ElemType,
1795  const Order,
1796  const unsigned int,
1797  const unsigned int,
1798  const Point & )
1799 {
1800  libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1801  return 0.;
1802 }
1803 
1804 
1805 
1806 template <>
1807 Real FE<3,HIERARCHIC>::shape_second_deriv(const Elem * elem,
1808  const Order order,
1809  const unsigned int i,
1810  const unsigned int j,
1811  const Point & p,
1812  const bool add_p_level)
1813 {
1814  return fe_hierarchic_3D_shape_second_deriv<HIERARCHIC>(elem, order, i, j, p, add_p_level);
1815 }
1816 
1817 
1818 
1819 template <>
1820 Real FE<3,HIERARCHIC>::shape_second_deriv(const FEType fet,
1821  const Elem * elem,
1822  const unsigned int i,
1823  const unsigned int j,
1824  const Point & p,
1825  const bool add_p_level)
1826 {
1827  return fe_hierarchic_3D_shape_second_deriv<HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
1828 }
1829 
1830 
1831 
1832 template <>
1833 Real FE<3,L2_HIERARCHIC>::shape_second_deriv(const Elem * elem,
1834  const Order order,
1835  const unsigned int i,
1836  const unsigned int j,
1837  const Point & p,
1838  const bool add_p_level)
1839 {
1840  return fe_hierarchic_3D_shape_second_deriv<L2_HIERARCHIC>(elem, order, i, j, p, add_p_level);
1841 }
1842 
1843 
1844 template <>
1845 Real FE<3,L2_HIERARCHIC>::shape_second_deriv(const FEType fet,
1846  const Elem * elem,
1847  const unsigned int i,
1848  const unsigned int j,
1849  const Point & p,
1850  const bool add_p_level)
1851 {
1852  return fe_hierarchic_3D_shape_second_deriv<L2_HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
1853 }
1854 
1855 
1856 template <>
1857 Real FE<3,SIDE_HIERARCHIC>::shape_second_deriv(const Elem * elem,
1858  const Order order,
1859  const unsigned int i,
1860  const unsigned int j,
1861  const Point & p,
1862  const bool add_p_level)
1863 {
1864 #if LIBMESH_DIM == 3
1865  libmesh_assert(elem);
1866  const ElemType type = elem->type();
1867 
1868  const Order totalorder = order + add_p_level*elem->p_level();
1869 
1870  if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
1871  return 0; // constants have zero derivative
1872 
1873  switch (type)
1874  {
1875  case HEX27:
1876  {
1877  // I need to debug the p>2 case here...
1878  if (totalorder > 2)
1879  return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
1880  FE<3,SIDE_HIERARCHIC>::shape_deriv);
1881 
1882  const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+1u);
1883  libmesh_assert_less(i, 6*dofs_per_side);
1884 
1885  const unsigned int sidenum = cube_side(p);
1886  if (sidenum > 5)
1887  return std::numeric_limits<Real>::quiet_NaN();
1888 
1889  const unsigned int dof_offset = sidenum * dofs_per_side;
1890 
1891  if (i < dof_offset) // i is on a previous side
1892  return 0;
1893 
1894  if (i >= dof_offset + dofs_per_side) // i is on a later side
1895  return 0;
1896 
1897  unsigned int side_i = i - dof_offset;
1898 
1899  std::unique_ptr<const Elem> side = elem->build_side_ptr(sidenum);
1900 
1901  Point sidep = cube_side_point(sidenum, p);
1902 
1903  cube_remap(side_i, *side, totalorder, sidep);
1904 
1905  // What second derivative or mixed derivative on the side
1906  // corresponds to the xi/eta/zeta mix we were asked for?
1907  unsigned int sidej = 100;
1908 
1909  // Do we need a -1 here to flip the final derivative value?
1910  Real f = 1.;
1911 
1912  switch (j)
1913  {
1914  case 0: // d^2()/dxi^2
1915  {
1916  switch (sidenum)
1917  {
1918  case 0:
1919  sidej = 2;
1920  break;
1921  case 1:
1922  sidej = 0;
1923  break;
1924  case 2:
1925  return 0;
1926  case 3:
1927  sidej = 0;
1928  break;
1929  case 4:
1930  return 0;
1931  case 5:
1932  sidej = 0;
1933  break;
1934  default:
1935  libmesh_error();
1936  }
1937  break;
1938  }
1939  case 1: // d^2()/dxideta
1940  {
1941  switch (sidenum)
1942  {
1943  case 0:
1944  sidej = 1;
1945  break;
1946  case 1:
1947  case 2:
1948  case 3:
1949  case 4:
1950  return 0;
1951  case 5:
1952  sidej = 1;
1953  break;
1954  default:
1955  libmesh_error();
1956  }
1957  break;
1958  }
1959  case 2: // d^2()/deta^2
1960  {
1961  switch (sidenum)
1962  {
1963  case 0:
1964  sidej = 0;
1965  break;
1966  case 1:
1967  return 0;
1968  case 2:
1969  sidej = 0;
1970  break;
1971  case 3:
1972  return 0;
1973  case 4:
1974  sidej = 0;
1975  break;
1976  case 5:
1977  sidej = 2;
1978  break;
1979  default:
1980  libmesh_error();
1981  }
1982  break;
1983  }
1984  case 3: // d^2()/dxidzeta
1985  {
1986  switch (sidenum)
1987  {
1988  case 0:
1989  return 0;
1990  case 1:
1991  sidej = 1;
1992  break;
1993  case 2:
1994  return 0;
1995  case 3:
1996  sidej = 1;
1997  f = -1;
1998  break;
1999  case 4:
2000  case 5:
2001  return 0;
2002  default:
2003  libmesh_error();
2004  }
2005  break;
2006  }
2007  case 4: // d^2()/detadzeta
2008  {
2009  switch (sidenum)
2010  {
2011  case 0:
2012  case 1:
2013  return 0;
2014  case 2:
2015  sidej = 1;
2016  break;
2017  case 3:
2018  return 0;
2019  case 4:
2020  sidej = 1;
2021  f = -1;
2022  break;
2023  case 5:
2024  return 0;
2025  default:
2026  libmesh_error();
2027  }
2028  break;
2029  }
2030  case 5: // d^2()/dzeta^2
2031  {
2032  switch (sidenum)
2033  {
2034  case 0:
2035  return 0;
2036  case 1:
2037  case 2:
2038  case 3:
2039  case 4:
2040  sidej = 2;
2041  break;
2042  case 5:
2043  return 0;
2044  default:
2045  libmesh_error();
2046  }
2047  break;
2048  }
2049 
2050  default:
2051  libmesh_error_msg("Invalid derivative index j = " << j);
2052  }
2053 
2054  return f * FE<2,HIERARCHIC>::shape_second_deriv(side.get(),
2055  order, side_i,
2056  sidej, sidep,
2057  add_p_level);
2058  }
2059 
2060  case TET14:
2061  case PRISM20:
2062  case PRISM21:
2063  {
2064  return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
2065  FE<3,SIDE_HIERARCHIC>::shape_deriv);
2066  }
2067 
2068  default:
2069  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
2070  }
2071 
2072 #else // LIBMESH_DIM != 3
2073  libmesh_ignore(elem, order, i, j, p, add_p_level);
2074  libmesh_not_implemented();
2075 #endif
2076 }
2077 
2078 
2079 template <>
2080 Real FE<3,SIDE_HIERARCHIC>::shape_second_deriv(const FEType fet,
2081  const Elem * elem,
2082  const unsigned int i,
2083  const unsigned int j,
2084  const Point & p,
2085  const bool add_p_level)
2086 {
2087  return FE<3,SIDE_HIERARCHIC>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
2088 }
2089 
2090 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
2091 
2092 } // namespace libMesh
2093 
2094 
2095 
2096 namespace
2097 {
2098 using namespace libMesh;
2099 
2100 
2101 unsigned int cube_side (const Point & p)
2102 {
2103  const Real xi = p(0), eta = p(1), zeta = p(2);
2104  const Real absxi = std::abs(xi),
2105  abseta = std::abs(eta),
2106  abszeta = std::abs(zeta);
2107  const Real maxabs_xi_eta = std::max(absxi, abseta),
2108  maxabs_xi_zeta = std::max(absxi, abszeta),
2109  maxabs_eta_zeta = std::max(abseta, abszeta);
2110 
2111  if (zeta < -maxabs_xi_eta)
2112  return 0;
2113  else if (eta < -maxabs_xi_zeta)
2114  return 1;
2115  else if (xi > maxabs_eta_zeta)
2116  return 2;
2117  else if (eta > maxabs_xi_zeta)
2118  return 3;
2119  else if (xi < -maxabs_eta_zeta)
2120  return 4;
2121  else if (zeta > maxabs_xi_eta)
2122  return 5;
2123 
2124  // We need to be able to return invalid values for cases where
2125  // mixed FE are being evaluated together on edges and vertices
2126  return 65535;
2127 }
2128 
2129 
2130 
2131 Point cube_side_point(unsigned int sidenum, const Point & p)
2132 {
2133  Point sidep;
2134 
2135  switch (sidenum)
2136  {
2137  case 0:
2138  sidep(0) = p(1);
2139  sidep(1) = p(0);
2140  break;
2141  case 1:
2142  sidep(0) = p(0);
2143  sidep(1) = p(2);
2144  break;
2145  case 2:
2146  sidep(0) = p(1);
2147  sidep(1) = p(2);
2148  break;
2149  case 3:
2150  sidep(0) = -p(0);
2151  sidep(1) = p(2);
2152  break;
2153  case 4:
2154  sidep(0) = -p(1);
2155  sidep(1) = p(2);
2156  break;
2157  case 5:
2158  sidep(0) = p(0);
2159  sidep(1) = p(1);
2160  break;
2161  default:
2162  libmesh_error();
2163  }
2164 
2165  return sidep;
2166 }
2167 
2168 
2169 void orient_quad(const Elem & elem,
2170  std::array<unsigned int, 4> & face_vertex)
2171 {
2172  // Sort the minimum point into face_vertex[0], the minimum of its
2173  // neighbors into face_vertex[1]. Keep the other two consistent; we
2174  // want to rotate or flip the quad but not to twist it.
2175 
2176  const unsigned int min_pt =
2177  std::min_element(face_vertex.begin(), face_vertex.end(),
2178  [&elem](auto v1, auto v2)
2179  {return elem.point(v1)<elem.point(v2);}) -
2180  face_vertex.begin();
2181 
2182  // Do we flip the quad?
2183  if (elem.point(face_vertex[(min_pt+3)%4]) <
2184  elem.point(face_vertex[(min_pt+1)%4]))
2185  face_vertex = { face_vertex[min_pt], face_vertex[(min_pt+3)%4],
2186  face_vertex[(min_pt+2)%4], face_vertex[(min_pt+1)%4] };
2187  else
2188  face_vertex = { face_vertex[min_pt], face_vertex[(min_pt+1)%4],
2189  face_vertex[(min_pt+2)%4], face_vertex[(min_pt+3)%4] };
2190 }
2191 
2192 
2193 void orient_triangle(const Elem & elem,
2194  unsigned int * face_vertex)
2195 {
2196  // Reorient nodes to account for flipping and rotation.
2197  // We could try to identify indices with symmetric shape
2198  // functions, to skip this in those cases, if we really
2199  // need to optimize later.
2200  //
2201  // With only 3 items, we should bubble sort!
2202  // Programming-for-MechE's class pays off!
2203  bool lastcheck = true;
2204  if (elem.point(face_vertex[0]) > elem.point(face_vertex[1]))
2205  {
2206  std::swap(face_vertex[0], face_vertex[1]);
2207  lastcheck = true;
2208  }
2209  if (elem.point(face_vertex[1]) > elem.point(face_vertex[2]))
2210  std::swap(face_vertex[1], face_vertex[2]);
2211  if (lastcheck && elem.point(face_vertex[0]) > elem.point(face_vertex[1]))
2212  std::swap(face_vertex[0], face_vertex[1]);
2213 }
2214 
2215 
2216 std::array<unsigned int, 4> oriented_prism_nodes(const Elem & elem,
2217  unsigned int face_num)
2218 {
2219  std::array<unsigned int, 4> face_vertex
2220  { Prism6::side_nodes_map[face_num][0],
2221  Prism6::side_nodes_map[face_num][1],
2222  Prism6::side_nodes_map[face_num][2],
2223  Prism6::side_nodes_map[face_num][3] };
2224 
2225  if (face_num > 0 && face_num < 4)
2226  orient_quad(elem, face_vertex);
2227  else
2228  orient_triangle(elem, face_vertex.data());
2229 
2230  return face_vertex;
2231 }
2232 
2233 
2234 std::array<unsigned int, 3> oriented_tet_nodes(const Elem & elem,
2235  unsigned int face_num)
2236 {
2237  std::array<unsigned int, 3> face_vertex
2238  { Tet4::side_nodes_map[face_num][0],
2239  Tet4::side_nodes_map[face_num][1],
2240  Tet4::side_nodes_map[face_num][2] };
2241 
2242  orient_triangle(elem, face_vertex.data());
2243 
2244  return face_vertex;
2245 }
2246 
2247 
2248 template <FEFamily T>
2249 Real fe_hierarchic_3D_shape(const Elem * elem,
2250  const Order order,
2251  const unsigned int i,
2252  const Point & p,
2253  const bool add_p_level)
2254 {
2255 #if LIBMESH_DIM == 3
2256 
2257  libmesh_assert(elem);
2258  const ElemType type = elem->type();
2259 
2260  const Order totalorder = order + add_p_level*elem->p_level();
2261 
2262  switch (type)
2263  {
2264  case HEX8:
2265  case HEX20:
2266  libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
2267  libmesh_fallthrough();
2268  case HEX27:
2269  {
2270  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+1u));
2271 
2272  // Compute hex shape functions as a tensor-product
2273  Real xi = p(0);
2274  Real eta = p(1);
2275  Real zeta = p(2);
2276 
2277  unsigned int i0, i1, i2;
2278 
2279  cube_indices(elem, totalorder, i, xi, eta, zeta, i0, i1, i2);
2280 
2281  return (FE<1,T>::shape(EDGE3, totalorder, i0, xi)*
2282  FE<1,T>::shape(EDGE3, totalorder, i1, eta)*
2283  FE<1,T>::shape(EDGE3, totalorder, i2, zeta));
2284  }
2285 
2286  case PRISM6:
2287  case PRISM15:
2288  libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
2289  libmesh_fallthrough();
2290  case PRISM18:
2291  libmesh_assert (T == L2_HIERARCHIC || totalorder < 3);
2292  libmesh_fallthrough();
2293  case PRISM20:
2294  case PRISM21:
2295  {
2296  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+2u)/2u);
2297 
2298  // Compute prism shape functions as a tensor-product.
2299  // Non-const here, because prism_indices might need to do some
2300  // flips before evaluating edge or face DoFs
2301  Point xi_eta {p(0),p(1)};
2302  Real zeta = p(2);
2303 
2304  unsigned int i01, i2;
2305 
2306  prism_indices(elem, totalorder, i, xi_eta, zeta, i01, i2);
2307 
2308  // We'll use the 2D Tri to handle any basis function flipping
2309  // needed in xi/eta for triangle face+edge bases.
2310  Tri3 tri;
2311 
2312  // We pinky swear not to modify these nodes
2313  Elem & e = const_cast<Elem &>(*elem);
2314  if (i2 == 0)
2315  {
2316  tri.set_node(0, e.node_ptr(0));
2317  tri.set_node(1, e.node_ptr(1));
2318  tri.set_node(2, e.node_ptr(2));
2319  }
2320  else if (i2 == 1)
2321  {
2322  tri.set_node(0, e.node_ptr(3));
2323  tri.set_node(1, e.node_ptr(4));
2324  tri.set_node(2, e.node_ptr(5));
2325  }
2326  else
2327  {
2328  // For interior DoFs, no flipping is necessary or done; we
2329  // can just evaluate on any triangle ... but *not* the
2330  // obvious 9,10,11 triangle, because that might not exist
2331  // if we have L2_HIERARCHIC on Prism6.
2332  tri.set_node(0, e.node_ptr(0));
2333  tri.set_node(1, e.node_ptr(1));
2334  tri.set_node(2, e.node_ptr(2));
2335 
2336  // For square face DoFs, prism_indices handles flipping,
2337  // and we *can't* override that in the tri shape call.
2338  if (i01 > 2 && i01 < 3u*totalorder)
2339  {
2340  // %(p-1) to find the edge number, %2 for even vs odd
2341  const bool odd_basis = ((i01-1)%(totalorder-1))%2;
2342  if (odd_basis)
2343  {
2344  const int tri_edge = (i01-3)/(totalorder-1);
2345  // Flip nodes now to avoid triggering a shape
2346  // function flip later
2347  if (tri.point(tri_edge) > tri.point((tri_edge+1)%3))
2348  {
2349  Node * n = tri.node_ptr(tri_edge);
2350  tri.set_node(tri_edge, tri.node_ptr((tri_edge+1)%3));
2351  tri.set_node((tri_edge+1)%3, n);
2352  }
2353  }
2354  }
2355  }
2356 
2357  return (FE<2,L2_HIERARCHIC>::shape(&tri, totalorder, i01, xi_eta)*
2358  FE<1,L2_HIERARCHIC>::shape(EDGE2, totalorder, i2, zeta));
2359  }
2360 
2361  case TET4:
2362  libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
2363  libmesh_fallthrough();
2364  case TET10:
2365  libmesh_assert (T == L2_HIERARCHIC || totalorder < 3);
2366  libmesh_fallthrough();
2367  case TET14:
2368  {
2369  const Real zeta[4] = { 1 - p(0) - p(1) - p(2), p(0), p(1), p(2) };
2370 
2371  // Nodal DoFs
2372  if (i < 4)
2373  return zeta[i];
2374 
2375  // Edge DoFs
2376  else if (i < 6u*totalorder - 2u)
2377  {
2378  const unsigned int edge_num = (i - 4) / (totalorder - 1u);
2379  // const int edge_node = edge_num + 4;
2380  const unsigned int basisorder = i - 2 - ((totalorder - 1u) * edge_num);
2381 
2382  const unsigned int edgevertex0 = Tet4::edge_nodes_map[edge_num][0],
2383  edgevertex1 = Tet4::edge_nodes_map[edge_num][1];
2384 
2385  // Get factors to account for edge-flipping
2386  Real flip = 1;
2387  if (basisorder%2 &&
2388  elem->positive_edge_orientation(edge_num))
2389  flip = -1;
2390 
2391  const Real crossval = zeta[edgevertex0] + zeta[edgevertex1];
2392  const Real edgenumerator = zeta[edgevertex1] - zeta[edgevertex0];
2393 
2394  if (crossval == 0.) // Yes, exact comparison; we seem numerically stable otherwise
2395  {
2396  unsigned int basisfactorial = 1.;
2397  for (unsigned int n=2; n <= basisorder; ++n)
2398  basisfactorial *= n;
2399 
2400  return std::pow(edgenumerator, basisorder) / basisfactorial;
2401  }
2402 
2403  const Real edgeval = edgenumerator / crossval;
2404  const Real crossfunc = std::pow(crossval, basisorder);
2405 
2406  return flip * crossfunc *
2407  FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
2408  basisorder, edgeval);
2409  }
2410 
2411  // Face DoFs
2412  else if (i < 2u*totalorder*totalorder + 2u)
2413  {
2414  const int dofs_per_face = (totalorder - 1u) * (totalorder - 2u) / 2;
2415  const int face_num = (i - (6u*totalorder - 2u)) / dofs_per_face;
2416 
2417  const std::array<unsigned int, 3> face_vertex =
2418  oriented_tet_nodes(*elem, face_num);
2419  const Real zeta0 = zeta[face_vertex[0]],
2420  zeta1 = zeta[face_vertex[1]],
2421  zeta2 = zeta[face_vertex[2]];
2422 
2423  const unsigned int basisnum =
2424  i - 4 -
2425  (totalorder - 1u) * /*n_edges*/6 -
2426  (dofs_per_face * face_num);
2427 
2428  const unsigned int exp0 = triangular_number_column[basisnum] + 1;
2429  const unsigned int exp1 = triangular_number_row[basisnum] + 1 -
2430  triangular_number_column[basisnum];
2431 
2432  Real returnval = 1;
2433  for (unsigned int n = 0; n != exp0; ++n)
2434  returnval *= zeta0;
2435  for (unsigned int n = 0; n != exp1; ++n)
2436  returnval *= zeta1;
2437  returnval *= zeta2;
2438  return returnval;
2439  }
2440 
2441  // Interior DoFs
2442  else
2443  {
2444  const unsigned int basisnum = i - 2u*totalorder*totalorder - 2u;
2445  const unsigned int exp0 = tetrahedral_number_column[basisnum] + 1;
2446  const unsigned int exp1 = tetrahedral_number_row[basisnum] + 1 -
2447  tetrahedral_number_column[basisnum] -
2448  tetrahedral_number_page[basisnum];
2449  const unsigned int exp2 = tetrahedral_number_page[basisnum] + 1;
2450 
2451  Real returnval = 1;
2452  for (unsigned int n = 0; n != exp0; ++n)
2453  returnval *= zeta[0];
2454  for (unsigned int n = 0; n != exp1; ++n)
2455  returnval *= zeta[1];
2456  for (unsigned int n = 0; n != exp2; ++n)
2457  returnval *= zeta[2];
2458  returnval *= zeta[3];
2459  return returnval;
2460  }
2461  }
2462 
2463  default:
2464  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
2465  }
2466 
2467 #else // LIBMESH_DIM != 3
2468  libmesh_ignore(elem, order, i, p, add_p_level);
2469  libmesh_not_implemented();
2470 #endif
2471 }
2472 
2473 
2474 
2475 template <FEFamily T>
2476 Real fe_hierarchic_3D_shape_deriv(const Elem * elem,
2477  const Order order,
2478  const unsigned int i,
2479  const unsigned int j,
2480  const Point & p,
2481  const bool add_p_level)
2482 {
2483  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,T>::shape);
2484 }
2485 
2486 
2487 
2488 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2489 
2490 template <FEFamily T>
2491 Real fe_hierarchic_3D_shape_second_deriv(const Elem * elem,
2492  const Order order,
2493  const unsigned int i,
2494  const unsigned int j,
2495  const Point & p,
2496  const bool add_p_level)
2497 {
2498  return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
2499  FE<3,T>::shape_deriv);
2500 }
2501 
2502 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
2503 
2504 
2505 } // anonymous namespace
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
libMesh::Tri3
The Tri3 is an element in 2D composed of 3 nodes.
Definition: face_tri3.h:61
libMesh::prism_number_page
const unsigned char prism_number_page[]
Definition: number_lookups.C:726
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::Prism6::side_nodes_map
static const unsigned int side_nodes_map[num_sides][nodes_per_side]
This maps the node of the side to element node numbers.
Definition: cell_prism6.h:172
libMesh::Elem::set_node
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2558
libMesh::Node
A Node is like a Point, but with more information.
Definition: node.h:52
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Quad4
The QUAD4 is an element in 2D composed of 4 nodes.
Definition: face_quad4.h:53
libMesh::IntRange
The IntRange templated class is intended to make it easy to loop over integers which are indices of a...
Definition: int_range.h:53
libMesh::Elem::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i)=0
libMesh::Tet4::side_nodes_map
static const unsigned int side_nodes_map[num_sides][nodes_per_side]
This maps the node of the side to element node numbers.
Definition: cell_tet4.h:195
libMesh::prism_number_triangle
const unsigned char prism_number_triangle[]
Definition: number_lookups.C:621
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::tetrahedral_number_column
const unsigned char tetrahedral_number_column[]
Definition: number_lookups.C:421
libMesh::cube_number_column
const unsigned char cube_number_column[]
Definition: number_lookups.C:84
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:3108
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::tetrahedral_number_row
const unsigned char tetrahedral_number_row[]
Definition: number_lookups.C:488
libMesh::cube_number_row
const unsigned char cube_number_row[]
Definition: number_lookups.C:196
libMesh::tetrahedral_number_page
const unsigned char tetrahedral_number_page[]
Definition: number_lookups.C:550
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::cube_number_page
const unsigned char cube_number_page[]
Definition: number_lookups.C:308
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::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:515
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:3589
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::Elem::node_ptr
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:2507
libMesh::square_number_row
const unsigned char square_number_row[]
Definition: number_lookups.C:69
libMesh::Elem::positive_face_orientation
bool positive_face_orientation(const unsigned int i) const
Definition: elem.C:3598
libMesh::Tet4::edge_nodes_map
static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]
This maps the node of the edge to element node numbers.
Definition: cell_tet4.h:201
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:2453
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 Jul 17 2025 01:29:06 for libMesh by   doxygen 1.8.14