libMesh
Loading...
Searching...
No Matches
fe_hierarchic_shape_3D.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3
4// This library is free software; you can redistribute it and/or
5// modify it under the terms of the GNU Lesser General Public
6// License as published by the Free Software Foundation; either
7// version 2.1 of the License, or (at your option) any later version.
8
9// This library is distributed in the hope that it will be useful,
10// but WITHOUT ANY WARRANTY; without even the implied warranty of
11// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12// Lesser General Public License for more details.
13
14// You should have received a copy of the GNU Lesser General Public
15// License along with this library; if not, write to the Free Software
16// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17
18
19// 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
25#include "libmesh/face_tri3.h" // Faster to construct these on the stack
26#include "libmesh/face_quad4.h"
27
28// Anonymous namespace for functions shared by HIERARCHIC and
29// L2_HIERARCHIC implementations. Implementations appear at the bottom
30// of this file.
31namespace
32{
33using namespace libMesh;
34
35unsigned int cube_side(const Point & p);
36
37Point cube_side_point(unsigned int sidenum, const Point & interior_point);
38
39template <FEFamily T>
40Real fe_hierarchic_3D_shape(const Elem * elem,
41 const Order order,
42 const unsigned int i,
43 const Point & p,
44 const bool add_p_level);
45
46template <FEFamily T>
47Real fe_hierarchic_3D_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
56template <FEFamily T>
57Real fe_hierarchic_3D_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#if LIBMESH_DIM > 2
67// Remap non-face-nodes based on point ordering
68template <unsigned int N_nodes>
69unsigned int remap_node(unsigned int n,
70 const Elem & elem,
71 unsigned int nodebegin)
72{
73 std::array<const Point *, N_nodes> points;
74
75 for (auto i : IntRange<unsigned int>(0, N_nodes))
76 points[i] = &elem.point(nodebegin+i);
77
78 std::sort(points.begin(), points.end(),
79 [](const Point * a, const Point * b)
80 { return *a < *b; });
81
82 const Point * pn = points[n-nodebegin];
83
84 for (auto i : IntRange<unsigned int>(nodebegin, nodebegin+N_nodes))
85 if (pn == &elem.point(i))
86 return i;
87
88 libmesh_assert(false);
89 return libMesh::invalid_uint;
90}
91
92
93void cube_remap(unsigned int & side_i,
94 const Elem & elem,
95 unsigned int sidenum,
96 const Elem & side,
97 unsigned int totalorder,
98 Point & sidep)
99{
100 // "vertex" nodes are now decoupled from vertices, so we have
101 // to order them consistently otherwise
102 if (side_i < 4)
103 side_i = remap_node<4>(side_i, side, 0);
104
105 // And "edge" nodes are decoupled from edges, so we have to
106 // reorder them too!
107 else if (side_i < 4u*totalorder)
108 {
109 unsigned int side_node = (side_i - 4)/(totalorder-1)+4;
110 side_node = remap_node<4>(side_node, side, 4);
111 side_i = ((side_i - 4) % (totalorder - 1)) // old local edge_i
112 + 4 + (side_node-4)*(totalorder-1);
113 }
114
115 // Interior dofs in 2D don't care about where xi/eta point in
116 // physical space, but here we need them to match from both
117 // sides of a face!
118 else
119 {
120 // The side numbers its nodes as the side node map of elem does.
121 // Rotate the least node to the origin, then flip about the
122 // diagonal through it if its lesser neighbor is its predecessor
123 const unsigned int orientation = elem.face_orientation(sidenum);
124 const unsigned int min_side_node = orientation / 2;
125 const bool flip = orientation % 2;
126
127 switch (min_side_node) {
128 case 0:
129 if (flip)
130 std::swap(sidep(0), sidep(1));
131 break;
132 case 1:
133 sidep(0) = -sidep(0);
134 if (!flip)
135 std::swap(sidep(0), sidep(1));
136 break;
137 case 2:
138 sidep(0) = -sidep(0);
139 sidep(1) = -sidep(1);
140 if (flip)
141 std::swap(sidep(0), sidep(1));
142 break;
143 case 3:
144 sidep(1) = -sidep(1);
145 if (!flip)
146 std::swap(sidep(0), sidep(1));
147 break;
148 default:
149 libmesh_error();
150 }
151 }
152}
153
154
155void cube_indices(const Elem * elem,
156 const unsigned int totalorder,
157 const unsigned int i,
158 Real & xi, Real & eta, Real & zeta,
159 unsigned int & i0,
160 unsigned int & i1,
161 unsigned int & i2)
162{
163 // The only way to make any sense of this
164 // is to look at the mgflo/mg2/mgf documentation
165 // and make the cut-out cube!
166 // Example i0 and i1 values for totalorder = 3:
167 // FIXME - these examples are incorrect now that we've got truly
168 // hierarchic basis functions
169 // 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
170 // 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
171 // 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};
172 // 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};
173 // 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};
174
175 // the number of DoFs per edge appears everywhere:
176 const unsigned int e = totalorder - 1u;
177
178 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+1u));
179
180 Real xi_saved = xi, eta_saved = eta, zeta_saved = zeta;
181
182 // Vertices:
183 if (i == 0)
184 {
185 i0 = 0;
186 i1 = 0;
187 i2 = 0;
188 }
189 else if (i == 1)
190 {
191 i0 = 1;
192 i1 = 0;
193 i2 = 0;
194 }
195 else if (i == 2)
196 {
197 i0 = 1;
198 i1 = 1;
199 i2 = 0;
200 }
201 else if (i == 3)
202 {
203 i0 = 0;
204 i1 = 1;
205 i2 = 0;
206 }
207 else if (i == 4)
208 {
209 i0 = 0;
210 i1 = 0;
211 i2 = 1;
212 }
213 else if (i == 5)
214 {
215 i0 = 1;
216 i1 = 0;
217 i2 = 1;
218 }
219 else if (i == 6)
220 {
221 i0 = 1;
222 i1 = 1;
223 i2 = 1;
224 }
225 else if (i == 7)
226 {
227 i0 = 0;
228 i1 = 1;
229 i2 = 1;
230 }
231 // Edge 0
232 else if (i < 8 + e)
233 {
234 i0 = i - 6;
235 i1 = 0;
236 i2 = 0;
237 if (elem->positive_edge_orientation(0))
238 xi = -xi_saved;
239 }
240 // Edge 1
241 else if (i < 8 + 2*e)
242 {
243 i0 = 1;
244 i1 = i - e - 6;
245 i2 = 0;
246 if (elem->positive_edge_orientation(1))
247 eta = -eta_saved;
248 }
249 // Edge 2
250 else if (i < 8 + 3*e)
251 {
252 i0 = i - 2*e - 6;
253 i1 = 1;
254 i2 = 0;
255 if (!elem->positive_edge_orientation(2))
256 xi = -xi_saved;
257 }
258 // Edge 3
259 else if (i < 8 + 4*e)
260 {
261 i0 = 0;
262 i1 = i - 3*e - 6;
263 i2 = 0;
264 if (elem->positive_edge_orientation(3))
265 eta = -eta_saved;
266 }
267 // Edge 4
268 else if (i < 8 + 5*e)
269 {
270 i0 = 0;
271 i1 = 0;
272 i2 = i - 4*e - 6;
273 if (elem->positive_edge_orientation(4))
274 zeta = -zeta_saved;
275 }
276 // Edge 5
277 else if (i < 8 + 6*e)
278 {
279 i0 = 1;
280 i1 = 0;
281 i2 = i - 5*e - 6;
282 if (elem->positive_edge_orientation(5))
283 zeta = -zeta_saved;
284 }
285 // Edge 6
286 else if (i < 8 + 7*e)
287 {
288 i0 = 1;
289 i1 = 1;
290 i2 = i - 6*e - 6;
291 if (elem->positive_edge_orientation(6))
292 zeta = -zeta_saved;
293 }
294 // Edge 7
295 else if (i < 8 + 8*e)
296 {
297 i0 = 0;
298 i1 = 1;
299 i2 = i - 7*e - 6;
300 if (elem->positive_edge_orientation(7))
301 zeta = -zeta_saved;
302 }
303 // Edge 8
304 else if (i < 8 + 9*e)
305 {
306 i0 = i - 8*e - 6;
307 i1 = 0;
308 i2 = 1;
309 if (elem->positive_edge_orientation(8))
310 xi = -xi_saved;
311 }
312 // Edge 9
313 else if (i < 8 + 10*e)
314 {
315 i0 = 1;
316 i1 = i - 9*e - 6;
317 i2 = 1;
318 if (elem->positive_edge_orientation(9))
319 eta = -eta_saved;
320 }
321 // Edge 10
322 else if (i < 8 + 11*e)
323 {
324 i0 = i - 10*e - 6;
325 i1 = 1;
326 i2 = 1;
327 if (!elem->positive_edge_orientation(10))
328 xi = -xi_saved;
329 }
330 // Edge 11
331 else if (i < 8 + 12*e)
332 {
333 i0 = 0;
334 i1 = i - 11*e - 6;
335 i2 = 1;
336 if (elem->positive_edge_orientation(11))
337 eta = -eta_saved;
338 }
339 // Faces
340 else if (i < 8 + 12*e + 6*e*e)
341 {
342 const unsigned int face = (i - 8 - 12*e) / (e*e);
343 const unsigned int basisnum = (i - 8 - 12*e) % (e*e);
344
345 // The face's own frame has its origin at the least vertex and its
346 // first axis toward that vertex's lesser neighbor, so every element
347 // sharing the face agrees on it. The origin sits at -1 along both
348 // axes, so the in-face coordinates are projections onto them.
349 const std::vector<unsigned int> v =
350 elem->oriented_face_vertices(face);
351 const Point origin = elem->master_point(v[0]);
352 const Point first = (elem->master_point(v[1]) - origin) / 2;
353 const Point second = (elem->master_point(v[3]) - origin) / 2;
354 const Point p_saved(xi_saved, eta_saved, zeta_saved);
355
356 Real * const coord[3] = {&xi, &eta, &zeta};
357 unsigned int * const index[3] = {&i0, &i1, &i2};
358
359 for (const auto d : make_range(3u))
360 if (first(d) != 0)
361 {
362 *index[d] = square_number_row[basisnum] + 2;
363 *coord[d] = p_saved * first;
364 }
365 else if (second(d) != 0)
366 {
367 *index[d] = square_number_column[basisnum] + 2;
368 *coord[d] = p_saved * second;
369 }
370 else
371 // The vertex mode that is 1 on this face
372 *index[d] = (origin(d) > 0);
373 }
374
375 // Internal DoFs
376 else
377 {
378 unsigned int basisnum = i - 8 - 12*e - 6*e*e;
379 i0 = cube_number_column[basisnum] + 2;
380 i1 = cube_number_row[basisnum] + 2;
381 i2 = cube_number_page[basisnum] + 2;
382 }
383}
384
385
386// Order the barycentric coordinates of triangular prism face \p face_num
387// by the face's vertices. The triangle interior basis is not symmetric
388// in them, so both elements sharing the face must order them alike.
389void orient_triangle_coords(const Elem & elem,
390 const unsigned int face_num,
391 const Point & xi_eta_saved,
392 Point & xi_eta)
393{
394 const std::vector<unsigned int> face_vertex =
395 elem.oriented_face_vertices(face_num);
396
397 const Real barycentric[3] = {1 - xi_eta_saved(0) - xi_eta_saved(1),
398 xi_eta_saved(0),
399 xi_eta_saved(1)};
400
401 // The face's least vertex takes the place of the triangle's vertex 0, so xi and eta are the
402 // barycentric coordinates of the next two vertices in the face's order. Face 4's vertices
403 // 3-5 lie above vertices 0-2, so % 3 maps the vertices of either face to the triangle's.
404 xi_eta(0) = barycentric[face_vertex[1] % 3];
405 xi_eta(1) = barycentric[face_vertex[2] % 3];
406}
407
408
409void prism_indices(const Elem * elem,
410 const unsigned int totalorder,
411 const unsigned int i,
412 Point & xi_eta, Real & zeta,
413 unsigned int & i01,
414 unsigned int & i2)
415{
416 // the number of DoFs per edge appears everywhere:
417 const unsigned int e = totalorder - 1u;
418
419 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+2u)/2u);
420
421 Point xi_eta_saved = xi_eta;
422 Real zeta_saved = zeta;
423
424 // Vertices:
425 if (i == 0)
426 {
427 i01 = 0;
428 i2 = 0;
429 }
430 else if (i == 1)
431 {
432 i01 = 1;
433 i2 = 0;
434 }
435 else if (i == 2)
436 {
437 i01 = 2;
438 i2 = 0;
439 }
440 else if (i == 3)
441 {
442 i01 = 0;
443 i2 = 1;
444 }
445 else if (i == 4)
446 {
447 i01 = 1;
448 i2 = 1;
449 }
450 else if (i == 5)
451 {
452 i01 = 2;
453 i2 = 1;
454 }
455 // Edges 0,1,2 (vertices 6,7,8)
456 else if (i < 6 + 3*e)
457 {
458 // The TRI code will handle any flips here
459 i01 = i - 3;
460 i2 = 0;
461 }
462 // Edge 3,4,5 (vertices 9,10,11)
463 else if (i < 6 + 6*e)
464 {
465 i01 = (i - 6 - 3*e)/e; // which tri DoF are we?
466 i2 = (i - 6 - 3*e)%e+2; // edge DoF? +2 to skip endpoints
467 // EDGE evaluations don't flip, so handle that here
468 if (elem->positive_edge_orientation(i01+3))
469 zeta = -zeta;
470 }
471 // Edge 6,7,8 (vertices 12,13,14)
472 else if (i < 6 + 9*e)
473 {
474 // The TRI code will handle any flips here
475 i01 = i - 3 - 6*e;
476 i2 = 1;
477 }
478 // Faces 1, 2, 3 (nodes 15, 16, 17), the quadrilaterals
479 else if (i < 6 + 9*e + 3*e*e)
480 {
481 const unsigned int face = (i - 6 - 9*e) / (e*e) + 1;
482 const unsigned int basisnum = (i - 6 - 9*e) % (e*e);
483
484 // indexes in edge numbering
485 const unsigned int s0 = square_number_row[basisnum] + 2;
486 const unsigned int s1 = square_number_column[basisnum] + 2;
487
488 // The face's side node map runs along triangle edge face-1 from
489 // n0 to n1, then up to the top triangle and back. Its own frame has
490 // its origin at the least vertex and its first axis toward that
491 // vertex's lesser neighbor, so every element sharing the face
492 // agrees on it.
493 const unsigned int n0 = elem->local_side_node(face, 0),
494 n1 = elem->local_side_node(face, 1);
495 const std::vector<unsigned int> v =
496 elem->oriented_face_vertices(face);
497
498 // The first axis runs along zeta: the edge takes the column mode
499 const bool swap = (v[0] < 3) != (v[1] < 3);
500 i01 = (swap ? s1 : s0) + 1 + (face-1)*e; // edge to triangle side numbering
501 i2 = swap ? s0 : s1;
502
503 // The origin is on the top triangle
504 if (v[0] >= 3)
505 zeta = -zeta_saved;
506
507 // The origin is at the far end of the edge: reflect the triangle by
508 // exchanging the barycentric coordinates of the edge's vertices
509 if (v[0] % 3 == n1 % 3)
510 {
511 Real barycentric[3] = {1 - xi_eta_saved(0) - xi_eta_saved(1),
512 xi_eta_saved(0),
513 xi_eta_saved(1)};
514 std::swap(barycentric[n0 % 3], barycentric[n1 % 3]);
515 xi_eta(0) = barycentric[1];
516 xi_eta(1) = barycentric[2];
517 }
518 }
519 // Face 0, node 18 - node order due to hierarchic numbering
520 else if (i < 6 + 9*e + 3*e*e + e*(e-1)/2)
521 {
522 i01 = i - 3 - 6*e - 3*e*e;
523 i2 = 0;
524 orient_triangle_coords(*elem, 0, xi_eta_saved, xi_eta);
525 }
526 // Face 4
527 else if (i < 6 + 9*e + 3*e*e + e*(e-1))
528 {
529 i01 = i - 3 - 6*e - 3*e*e - e*(e-1)/2;
530 i2 = 1;
531 orient_triangle_coords(*elem, 4, xi_eta_saved, xi_eta);
532 }
533 // Internal DoFs
534 else
535 {
536 // We won't bother with any internal DoF reordering / flipping;
537 // that's fine unless we ever get to 4D.
538 unsigned int basisnum = i - 6 - 9*e - 3*e*e - e*(e-1);
539 i01 = prism_number_triangle[basisnum] + 3 + 3*e;
540 i2 = prism_number_page[basisnum] + 2;
541 }
542}
543
544#endif // LIBMESH_DIM > 2
545
546} // end anonymous namespace
547
548
549
550namespace libMesh
551{
552
553
554LIBMESH_DEFAULT_VECTORIZED_FE(3,HIERARCHIC)
555LIBMESH_DEFAULT_VECTORIZED_FE(3,L2_HIERARCHIC)
556LIBMESH_DEFAULT_VECTORIZED_FE(3,SIDE_HIERARCHIC)
557
558
559template <>
560Real FE<3,HIERARCHIC>::shape(const ElemType,
561 const Order,
562 const unsigned int,
563 const Point &)
564{
565 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
566 return 0.;
567}
568
569
570
571template <>
572Real FE<3,L2_HIERARCHIC>::shape(const ElemType,
573 const Order,
574 const unsigned int,
575 const Point &)
576{
577 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
578 return 0.;
579}
580
581
582
583template <>
584Real FE<3,SIDE_HIERARCHIC>::shape(const ElemType,
585 const Order,
586 const unsigned int,
587 const Point &)
588{
589 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
590 return 0.;
591}
592
593
594
595template <>
596Real FE<3,HIERARCHIC>::shape(const Elem * elem,
597 const Order order,
598 const unsigned int i,
599 const Point & p,
600 const bool add_p_level)
601{
602 return fe_hierarchic_3D_shape<HIERARCHIC>(elem, order, i, p, add_p_level);
603}
604
605
606template <>
607Real FE<3,HIERARCHIC>::shape(const FEType fet,
608 const Elem * elem,
609 const unsigned int i,
610 const Point & p,
611 const bool add_p_level)
612{
613 return fe_hierarchic_3D_shape<HIERARCHIC>(elem, fet.order, i, p, add_p_level);
614}
615
616
617
618
619template <>
620Real FE<3,L2_HIERARCHIC>::shape(const Elem * elem,
621 const Order order,
622 const unsigned int i,
623 const Point & p,
624 const bool add_p_level)
625{
626 return fe_hierarchic_3D_shape<L2_HIERARCHIC>(elem, order, i, p, add_p_level);
627}
628
629
630template <>
631Real FE<3,L2_HIERARCHIC>::shape(const FEType fet,
632 const Elem * elem,
633 const unsigned int i,
634 const Point & p,
635 const bool add_p_level)
636{
637 return fe_hierarchic_3D_shape<L2_HIERARCHIC>(elem, fet.order, i, p, add_p_level);
638}
639
640
641
642template <>
643Real FE<3,SIDE_HIERARCHIC>::shape(const Elem * elem,
644 const Order order,
645 const unsigned int i,
646 const Point & p,
647 const bool add_p_level)
648{
649#if LIBMESH_DIM == 3
650 libmesh_assert(elem);
651 const ElemType type = elem->type();
652
653 const Order totalorder = order + add_p_level*elem->p_level();
654
655 switch (type)
656 {
657 case HEX27:
658 {
659 const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+1u);
660 libmesh_assert_less(i, 6*dofs_per_side);
661
662 const unsigned int sidenum = cube_side(p);
663 if (sidenum > 5)
664 return std::numeric_limits<Real>::quiet_NaN();
665
666 const unsigned int dof_offset = sidenum * dofs_per_side;
667
668 if (i < dof_offset) // i is on a previous side
669 return 0;
670
671 if (i >= dof_offset + dofs_per_side) // i is on a later side
672 return 0;
673
674 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
675 return 1;
676
677 unsigned int side_i = i - dof_offset;
678
679 std::unique_ptr<const Elem> side = elem->build_side_ptr(sidenum);
680
681 Point sidep = cube_side_point(sidenum, p);
682
683 cube_remap(side_i, *elem, sidenum, *side, totalorder, sidep);
684
685 return FE<2,HIERARCHIC>::shape(side.get(), order, side_i, sidep, add_p_level);
686 }
687
688 case TET14:
689 {
690 const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+2u)/2u;
691 libmesh_assert_less(i, 4*dofs_per_side);
692
693 const Real zeta[4] = { Real(1.) - p(0) - p(1) - p(2), p(0), p(1), p(2) };
694
695 unsigned int face_num = 0;
696 if (zeta[0] > zeta[3] &&
697 zeta[1] > zeta[3] &&
698 zeta[2] > zeta[3])
699 {
700 face_num = 0;
701 }
702 else if (zeta[0] > zeta[2] &&
703 zeta[1] > zeta[2] &&
704 zeta[3] > zeta[2])
705 {
706 face_num = 1;
707 }
708 else if (zeta[1] > zeta[0] &&
709 zeta[2] > zeta[0] &&
710 zeta[3] > zeta[0])
711 {
712 face_num = 2;
713 }
714 else
715 {
716 // We'd better not be right between two faces
717 libmesh_assert (zeta[0] > zeta[1] &&
718 zeta[2] > zeta[1] &&
719 zeta[3] > zeta[1]);
720 face_num = 3;
721 }
722
723 if (i < face_num * dofs_per_side ||
724 i >= (face_num+1) * dofs_per_side)
725 return 0;
726
727 if (totalorder == 0)
728 return 1;
729
730 const std::vector<unsigned int> face_vertex =
731 elem->oriented_face_vertices(face_num);
732
733 // We only need a Tri3 to evaluate L2_HIERARCHIC on the affine
734 // master element
735 Tri3 side;
736
737 // We pinky swear not to modify these nodes
738 Elem & e = const_cast<Elem &>(*elem);
739 side.set_node(0, e.node_ptr(face_vertex[0]));
740 side.set_node(1, e.node_ptr(face_vertex[1]));
741 side.set_node(2, e.node_ptr(face_vertex[2]));
742
743 const unsigned int basisnum = i - face_num*dofs_per_side;
744
745 Point sidep {zeta[face_vertex[1]], zeta[face_vertex[2]]};
746
747 return FE<2,L2_HIERARCHIC>::shape(&side, totalorder,
748 basisnum, sidep, false);
749 }
750
751 case PRISM20:
752 case PRISM21:
753 {
754 const unsigned int dofs_per_quad = (totalorder+1u)*(totalorder+1u);
755 const unsigned int dofs_per_tri = (totalorder+1u)*(totalorder+2u)/2u;
756 libmesh_assert_less(i, 3*dofs_per_quad + 2*dofs_per_tri);
757
758 // We only need a Tri3 or Quad4 to evaluate L2_HIERARCHIC on
759 // the affine master element
760 Tri3 tri;
761 Quad4 quad;
762 Elem * side = &quad;
763 unsigned int dofs_on_side = dofs_per_quad;
764
765 // We pinky swear not to modify the nodes we'll point to
766 Elem & e = const_cast<Elem &>(*elem);
767
768 Point sidep;
769
770 // Face number calculation is tricky - the ordering of side
771 // nodes on Prisms does *not* match the ordering of sides!
772 // (the mid-triangle side nodes were added "later")
773 // Here face_num will be the numbering that matches the side
774 // number, but i_offset will have to consider the nodal
775 // ordering.
776 unsigned int face_num = 0;
777 unsigned int i_offset = 0;
778
779 // Triangular coordinates
780 const Real zeta[3] = { Real(1.) - p(0) - p(1), p(0), p(1) };
781
782 // Closeness to midplane
783 const Real zmid = 1 - std::abs(p(2));
784
785 if (zeta[1] > zeta[2] && zeta[0] > zeta[2] &&
786 zmid > 3*zeta[2]) // face 1, quad
787 {
788 face_num = 1;
789 i_offset = 0;
790 }
791 else if (zeta[1] > zeta[0] && zeta[2] > zeta[0] &&
792 zmid > 3*zeta[0]) // face 2, quad
793 {
794 face_num = 2;
795 i_offset = dofs_per_quad;
796 }
797 else if (zeta[0] > zeta[1] && zeta[2] > zeta[1] &&
798 zmid > 3*zeta[1]) // face 3, quad
799 {
800 face_num = 3;
801 i_offset = 2*dofs_per_quad;
802 }
803 else if (p(2) + 1 < 3*zeta[0] &&
804 p(2) + 1 < 3*zeta[1] &&
805 p(2) + 1 < 3*zeta[2]) // face 0, tri
806 {
807 face_num = 0;
808 i_offset = 3*dofs_per_quad;
809 dofs_on_side = dofs_per_tri;
810 side = &tri;
811 }
812 else if (1 - p(2) < 3*zeta[0] &&
813 1 - p(2) < 3*zeta[1] &&
814 1 - p(2) < 3*zeta[2]) // face 4, tri
815 {
816 face_num = 4;
817 i_offset = dofs_per_tri + 3*dofs_per_quad;
818 dofs_on_side = dofs_per_tri;
819 side = &tri;
820 }
821 else
822 {
823 libmesh_error_msg("Evaluating SIDE_HIERARCHIC right between two Prism faces?");
824 }
825
826 if (i < i_offset ||
827 i >= i_offset + dofs_on_side)
828 return 0;
829
830 if (totalorder == 0)
831 return 1;
832
833 const std::vector<unsigned int> face_vertex =
834 elem->oriented_face_vertices(face_num);
835
836 side->set_node(0, e.node_ptr(face_vertex[0]));
837 side->set_node(1, e.node_ptr(face_vertex[1]));
838 side->set_node(2, e.node_ptr(face_vertex[2]));
839 if (face_vertex.size() > 3)
840 side->set_node(3, e.node_ptr(face_vertex[3]));
841
842 if (face_num == 0 || face_num == 4)
843 sidep = {zeta[face_vertex[1]%3], zeta[face_vertex[2]%3]};
844 else
845 {
846 // Transform a coordinate from the master prism to the
847 // master quad, based on two vertex indices defining the
848 // coordinate's direction
849 auto coord_val = [p](int v1, int v2){
850 if (v2-v1 == 3)
851 return p(2);
852 else if (v2-v1 == -3)
853 return -p(2);
854 else if (v1%3 == 0 && v2%3 == 1)
855 return 2*p(0)-1;
856 else if (v2%3 == 0 && v1%3 == 1)
857 return 1-2*p(0);
858 else if (v1%3 == 1 && v2%3 == 2)
859 return p(1)-p(0);
860 else if (v2%3 == 1 && v1%3 == 2)
861 return p(0)-p(1);
862 else if (v1%3 == 2 && v2%3 == 0)
863 return 1-2*p(1);
864 else if (v2%3 == 2 && v1%3 == 0)
865 return 2*p(1)-1;
866 else
867 libmesh_error();
868 };
869
870 sidep = {coord_val(face_vertex[0], face_vertex[1]),
871 coord_val(face_vertex[0], face_vertex[3])};
872 }
873
874 const unsigned int basisnum = i - i_offset;
875
876 return FE<2,L2_HIERARCHIC>::shape(side, totalorder,
877 basisnum, sidep, false);
878 }
879
880
881 default:
882 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
883 }
884
885#else // LIBMESH_DIM != 3
886 libmesh_ignore(elem, order, i, p, add_p_level);
887 libmesh_not_implemented();
888#endif
889}
890
891
892template <>
893Real FE<3,SIDE_HIERARCHIC>::shape(const FEType fet,
894 const Elem * elem,
895 const unsigned int i,
896 const Point & p,
897 const bool add_p_level)
898{
899 return FE<3,SIDE_HIERARCHIC>::shape(elem,fet.order, i, p, add_p_level);
900}
901
902
903template <>
904Real FE<3,HIERARCHIC>::shape_deriv(const ElemType,
905 const Order,
906 const unsigned int,
907 const unsigned int,
908 const Point & )
909{
910 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
911 return 0.;
912}
913
914
915
916template <>
917Real FE<3,L2_HIERARCHIC>::shape_deriv(const ElemType,
918 const Order,
919 const unsigned int,
920 const unsigned int,
921 const Point & )
922{
923 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
924 return 0.;
925}
926
927
928
929template <>
930Real FE<3,SIDE_HIERARCHIC>::shape_deriv(const ElemType,
931 const Order,
932 const unsigned int,
933 const unsigned int,
934 const Point & )
935{
936 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
937 return 0.;
938}
939
940
941
942template <>
943Real FE<3,HIERARCHIC>::shape_deriv(const Elem * elem,
944 const Order order,
945 const unsigned int i,
946 const unsigned int j,
947 const Point & p,
948 const bool add_p_level)
949{
950 return fe_hierarchic_3D_shape_deriv<HIERARCHIC>(elem, order, i, j, p, add_p_level);
951}
952
953
954template <>
955Real FE<3,HIERARCHIC>::shape_deriv(const FEType fet,
956 const Elem * elem,
957 const unsigned int i,
958 const unsigned int j,
959 const Point & p,
960 const bool add_p_level)
961{
962 return fe_hierarchic_3D_shape_deriv<HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
963}
964
965
966
967template <>
968Real FE<3,L2_HIERARCHIC>::shape_deriv(const Elem * elem,
969 const Order order,
970 const unsigned int i,
971 const unsigned int j,
972 const Point & p,
973 const bool add_p_level)
974{
975 return fe_hierarchic_3D_shape_deriv<L2_HIERARCHIC>(elem, order, i, j, p, add_p_level);
976}
977
978
979template <>
980Real FE<3,L2_HIERARCHIC>::shape_deriv(const FEType fet,
981 const Elem * elem,
982 const unsigned int i,
983 const unsigned int j,
984 const Point & p,
985 const bool add_p_level)
986{
987 return fe_hierarchic_3D_shape_deriv<L2_HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
988}
989
990
991
992template <>
993Real FE<3,SIDE_HIERARCHIC>::shape_deriv(const Elem * elem,
994 const Order order,
995 const unsigned int i,
996 const unsigned int j,
997 const Point & p,
998 const bool add_p_level)
999{
1000#if LIBMESH_DIM == 3
1001 libmesh_assert(elem);
1002 const ElemType type = elem->type();
1003
1004 const Order totalorder = order + add_p_level*elem->p_level();
1005
1006 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
1007 return 0; // constants have zero derivative
1008
1009 switch (type)
1010 {
1011 case HEX27:
1012 {
1013 // I need to debug the p>2 case here...
1014 if (totalorder > 2)
1015 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,SIDE_HIERARCHIC>::shape);
1016
1017 const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+1u);
1018 libmesh_assert_less(i, 6*dofs_per_side);
1019
1020 const unsigned int sidenum = cube_side(p);
1021 if (sidenum > 5)
1022 return std::numeric_limits<Real>::quiet_NaN();
1023
1024 const unsigned int dof_offset = sidenum * dofs_per_side;
1025
1026 if (i < dof_offset) // i is on a previous side
1027 return 0;
1028
1029 if (i >= dof_offset + dofs_per_side) // i is on a later side
1030 return 0;
1031
1032 unsigned int side_i = i - dof_offset;
1033
1034 std::unique_ptr<const Elem> side = elem->build_side_ptr(sidenum);
1035
1036 Point sidep = cube_side_point(sidenum, p);
1037
1038 cube_remap(side_i, *elem, sidenum, *side, totalorder, sidep);
1039
1040 // What direction on the side corresponds to the derivative
1041 // direction we want?
1042 unsigned int sidej = 100;
1043
1044 // Do we need a -1 here to flip that direction?
1045 Real f = 1.;
1046
1047 switch (j)
1048 {
1049 case 0: // d()/dxi
1050 {
1051 switch (sidenum)
1052 {
1053 case 0:
1054 sidej = 1;
1055 break;
1056 case 1:
1057 sidej = 0;
1058 break;
1059 case 2:
1060 return 0;
1061 case 3:
1062 sidej = 0;
1063 f = -1;
1064 break;
1065 case 4:
1066 return 0;
1067 case 5:
1068 sidej = 0;
1069 break;
1070 default:
1071 libmesh_error();
1072 }
1073 break;
1074 }
1075 case 1: // d()/deta
1076 {
1077 switch (sidenum)
1078 {
1079 case 0:
1080 sidej = 0;
1081 break;
1082 case 1:
1083 return 0;
1084 case 2:
1085 sidej = 0;
1086 break;
1087 case 3:
1088 return 0;
1089 case 4:
1090 sidej = 0;
1091 f = -1;
1092 break;
1093 case 5:
1094 sidej = 1;
1095 break;
1096 default:
1097 libmesh_error();
1098 }
1099 break;
1100 }
1101 case 2: // d()/dzeta
1102 {
1103 switch (sidenum)
1104 {
1105 case 0:
1106 return 0;
1107 case 1:
1108 case 2:
1109 case 3:
1110 case 4:
1111 sidej = 1;
1112 break;
1113 case 5:
1114 return 0;
1115 default:
1116 libmesh_error();
1117 }
1118 break;
1119 }
1120
1121 default:
1122 libmesh_error_msg("Invalid derivative index j = " << j);
1123 }
1124
1125 return f * FE<2,HIERARCHIC>::shape_deriv(side.get(), order,
1126 side_i, sidej, sidep,
1127 add_p_level);
1128 }
1129
1130 case TET14:
1131 case PRISM20:
1132 case PRISM21:
1133 {
1134 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,SIDE_HIERARCHIC>::shape);
1135 }
1136
1137 default:
1138 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1139 }
1140
1141#else // LIBMESH_DIM != 3
1142 libmesh_ignore(elem, order, i, j, p, add_p_level);
1143 libmesh_not_implemented();
1144#endif
1145}
1146
1147
1148template <>
1149Real FE<3,SIDE_HIERARCHIC>::shape_deriv(const FEType fet,
1150 const Elem * elem,
1151 const unsigned int i,
1152 const unsigned int j,
1153 const Point & p,
1154 const bool add_p_level)
1155{
1156 return FE<3,SIDE_HIERARCHIC>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
1157}
1158
1159
1160#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1161
1162template <>
1163Real FE<3,HIERARCHIC>::shape_second_deriv(const ElemType,
1164 const Order,
1165 const unsigned int,
1166 const unsigned int,
1167 const Point & )
1168{
1169 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1170 return 0.;
1171}
1172
1173
1174
1175template <>
1176Real FE<3,L2_HIERARCHIC>::shape_second_deriv(const ElemType,
1177 const Order,
1178 const unsigned int,
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
1188template <>
1189Real FE<3,SIDE_HIERARCHIC>::shape_second_deriv(const ElemType,
1190 const Order,
1191 const unsigned int,
1192 const unsigned int,
1193 const Point & )
1194{
1195 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1196 return 0.;
1197}
1198
1199
1200
1201template <>
1202Real FE<3,HIERARCHIC>::shape_second_deriv(const Elem * elem,
1203 const Order order,
1204 const unsigned int i,
1205 const unsigned int j,
1206 const Point & p,
1207 const bool add_p_level)
1208{
1209 return fe_hierarchic_3D_shape_second_deriv<HIERARCHIC>(elem, order, i, j, p, add_p_level);
1210}
1211
1212
1213
1214template <>
1215Real FE<3,HIERARCHIC>::shape_second_deriv(const FEType fet,
1216 const Elem * elem,
1217 const unsigned int i,
1218 const unsigned int j,
1219 const Point & p,
1220 const bool add_p_level)
1221{
1222 return fe_hierarchic_3D_shape_second_deriv<HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
1223}
1224
1225
1226
1227template <>
1228Real FE<3,L2_HIERARCHIC>::shape_second_deriv(const Elem * elem,
1229 const Order order,
1230 const unsigned int i,
1231 const unsigned int j,
1232 const Point & p,
1233 const bool add_p_level)
1234{
1235 return fe_hierarchic_3D_shape_second_deriv<L2_HIERARCHIC>(elem, order, i, j, p, add_p_level);
1236}
1237
1238
1239template <>
1240Real FE<3,L2_HIERARCHIC>::shape_second_deriv(const FEType fet,
1241 const Elem * elem,
1242 const unsigned int i,
1243 const unsigned int j,
1244 const Point & p,
1245 const bool add_p_level)
1246{
1247 return fe_hierarchic_3D_shape_second_deriv<L2_HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
1248}
1249
1250
1251template <>
1252Real FE<3,SIDE_HIERARCHIC>::shape_second_deriv(const Elem * elem,
1253 const Order order,
1254 const unsigned int i,
1255 const unsigned int j,
1256 const Point & p,
1257 const bool add_p_level)
1258{
1259#if LIBMESH_DIM == 3
1260 libmesh_assert(elem);
1261 const ElemType type = elem->type();
1262
1263 const Order totalorder = order + add_p_level*elem->p_level();
1264
1265 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
1266 return 0; // constants have zero derivative
1267
1268 switch (type)
1269 {
1270 case HEX27:
1271 {
1272 // I need to debug the p>2 case here...
1273 if (totalorder > 2)
1274 return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
1275 FE<3,SIDE_HIERARCHIC>::shape_deriv);
1276
1277 const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+1u);
1278 libmesh_assert_less(i, 6*dofs_per_side);
1279
1280 const unsigned int sidenum = cube_side(p);
1281 if (sidenum > 5)
1282 return std::numeric_limits<Real>::quiet_NaN();
1283
1284 const unsigned int dof_offset = sidenum * dofs_per_side;
1285
1286 if (i < dof_offset) // i is on a previous side
1287 return 0;
1288
1289 if (i >= dof_offset + dofs_per_side) // i is on a later side
1290 return 0;
1291
1292 unsigned int side_i = i - dof_offset;
1293
1294 std::unique_ptr<const Elem> side = elem->build_side_ptr(sidenum);
1295
1296 Point sidep = cube_side_point(sidenum, p);
1297
1298 cube_remap(side_i, *elem, sidenum, *side, totalorder, sidep);
1299
1300 // What second derivative or mixed derivative on the side
1301 // corresponds to the xi/eta/zeta mix we were asked for?
1302 unsigned int sidej = 100;
1303
1304 // Do we need a -1 here to flip the final derivative value?
1305 Real f = 1.;
1306
1307 switch (j)
1308 {
1309 case 0: // d^2()/dxi^2
1310 {
1311 switch (sidenum)
1312 {
1313 case 0:
1314 sidej = 2;
1315 break;
1316 case 1:
1317 sidej = 0;
1318 break;
1319 case 2:
1320 return 0;
1321 case 3:
1322 sidej = 0;
1323 break;
1324 case 4:
1325 return 0;
1326 case 5:
1327 sidej = 0;
1328 break;
1329 default:
1330 libmesh_error();
1331 }
1332 break;
1333 }
1334 case 1: // d^2()/dxideta
1335 {
1336 switch (sidenum)
1337 {
1338 case 0:
1339 sidej = 1;
1340 break;
1341 case 1:
1342 case 2:
1343 case 3:
1344 case 4:
1345 return 0;
1346 case 5:
1347 sidej = 1;
1348 break;
1349 default:
1350 libmesh_error();
1351 }
1352 break;
1353 }
1354 case 2: // d^2()/deta^2
1355 {
1356 switch (sidenum)
1357 {
1358 case 0:
1359 sidej = 0;
1360 break;
1361 case 1:
1362 return 0;
1363 case 2:
1364 sidej = 0;
1365 break;
1366 case 3:
1367 return 0;
1368 case 4:
1369 sidej = 0;
1370 break;
1371 case 5:
1372 sidej = 2;
1373 break;
1374 default:
1375 libmesh_error();
1376 }
1377 break;
1378 }
1379 case 3: // d^2()/dxidzeta
1380 {
1381 switch (sidenum)
1382 {
1383 case 0:
1384 return 0;
1385 case 1:
1386 sidej = 1;
1387 break;
1388 case 2:
1389 return 0;
1390 case 3:
1391 sidej = 1;
1392 f = -1;
1393 break;
1394 case 4:
1395 case 5:
1396 return 0;
1397 default:
1398 libmesh_error();
1399 }
1400 break;
1401 }
1402 case 4: // d^2()/detadzeta
1403 {
1404 switch (sidenum)
1405 {
1406 case 0:
1407 case 1:
1408 return 0;
1409 case 2:
1410 sidej = 1;
1411 break;
1412 case 3:
1413 return 0;
1414 case 4:
1415 sidej = 1;
1416 f = -1;
1417 break;
1418 case 5:
1419 return 0;
1420 default:
1421 libmesh_error();
1422 }
1423 break;
1424 }
1425 case 5: // d^2()/dzeta^2
1426 {
1427 switch (sidenum)
1428 {
1429 case 0:
1430 return 0;
1431 case 1:
1432 case 2:
1433 case 3:
1434 case 4:
1435 sidej = 2;
1436 break;
1437 case 5:
1438 return 0;
1439 default:
1440 libmesh_error();
1441 }
1442 break;
1443 }
1444
1445 default:
1446 libmesh_error_msg("Invalid derivative index j = " << j);
1447 }
1448
1449 return f * FE<2,HIERARCHIC>::shape_second_deriv(side.get(),
1450 order, side_i,
1451 sidej, sidep,
1452 add_p_level);
1453 }
1454
1455 case TET14:
1456 case PRISM20:
1457 case PRISM21:
1458 {
1459 return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
1460 FE<3,SIDE_HIERARCHIC>::shape_deriv);
1461 }
1462
1463 default:
1464 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1465 }
1466
1467#else // LIBMESH_DIM != 3
1468 libmesh_ignore(elem, order, i, j, p, add_p_level);
1469 libmesh_not_implemented();
1470#endif
1471}
1472
1473
1474template <>
1475Real FE<3,SIDE_HIERARCHIC>::shape_second_deriv(const FEType fet,
1476 const Elem * elem,
1477 const unsigned int i,
1478 const unsigned int j,
1479 const Point & p,
1480 const bool add_p_level)
1481{
1482 return FE<3,SIDE_HIERARCHIC>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
1483}
1484
1485#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1486
1487} // namespace libMesh
1488
1489
1490
1491namespace
1492{
1493using namespace libMesh;
1494
1495
1496unsigned int cube_side (const Point & p)
1497{
1498 const Real xi = p(0), eta = p(1), zeta = p(2);
1499 const Real absxi = std::abs(xi),
1500 abseta = std::abs(eta),
1501 abszeta = std::abs(zeta);
1502 const Real maxabs_xi_eta = std::max(absxi, abseta),
1503 maxabs_xi_zeta = std::max(absxi, abszeta),
1504 maxabs_eta_zeta = std::max(abseta, abszeta);
1505
1506 if (zeta < -maxabs_xi_eta)
1507 return 0;
1508 else if (eta < -maxabs_xi_zeta)
1509 return 1;
1510 else if (xi > maxabs_eta_zeta)
1511 return 2;
1512 else if (eta > maxabs_xi_zeta)
1513 return 3;
1514 else if (xi < -maxabs_eta_zeta)
1515 return 4;
1516 else if (zeta > maxabs_xi_eta)
1517 return 5;
1518
1519 // We need to be able to return invalid values for cases where
1520 // mixed FE are being evaluated together on edges and vertices
1521 return 65535;
1522}
1523
1524
1525
1526Point cube_side_point(unsigned int sidenum, const Point & p)
1527{
1528 Point sidep;
1529
1530 switch (sidenum)
1531 {
1532 case 0:
1533 sidep(0) = p(1);
1534 sidep(1) = p(0);
1535 break;
1536 case 1:
1537 sidep(0) = p(0);
1538 sidep(1) = p(2);
1539 break;
1540 case 2:
1541 sidep(0) = p(1);
1542 sidep(1) = p(2);
1543 break;
1544 case 3:
1545 sidep(0) = -p(0);
1546 sidep(1) = p(2);
1547 break;
1548 case 4:
1549 sidep(0) = -p(1);
1550 sidep(1) = p(2);
1551 break;
1552 case 5:
1553 sidep(0) = p(0);
1554 sidep(1) = p(1);
1555 break;
1556 default:
1557 libmesh_error();
1558 }
1559
1560 return sidep;
1561}
1562
1563
1564
1565template <FEFamily T>
1566Real fe_hierarchic_3D_shape(const Elem * elem,
1567 const Order order,
1568 const unsigned int i,
1569 const Point & p,
1570 const bool add_p_level)
1571{
1572#if LIBMESH_DIM == 3
1573
1574 libmesh_assert(elem);
1575 const ElemType type = elem->type();
1576
1577 const Order totalorder = order + add_p_level*elem->p_level();
1578
1579 switch (type)
1580 {
1581 case HEX8:
1582 case HEX20:
1583 libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
1584 libmesh_fallthrough();
1585 case HEX27:
1586 {
1587 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+1u));
1588
1589 // Compute hex shape functions as a tensor-product
1590 Real xi = p(0);
1591 Real eta = p(1);
1592 Real zeta = p(2);
1593
1594 unsigned int i0, i1, i2;
1595
1596 cube_indices(elem, totalorder, i, xi, eta, zeta, i0, i1, i2);
1597
1598 return (FE<1,T>::shape(EDGE3, totalorder, i0, xi)*
1599 FE<1,T>::shape(EDGE3, totalorder, i1, eta)*
1600 FE<1,T>::shape(EDGE3, totalorder, i2, zeta));
1601 }
1602
1603 case PRISM6:
1604 case PRISM15:
1605 libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
1606 libmesh_fallthrough();
1607 case PRISM18:
1608 libmesh_assert (T == L2_HIERARCHIC || totalorder < 3);
1609 libmesh_fallthrough();
1610 case PRISM20:
1611 case PRISM21:
1612 {
1613 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+2u)/2u);
1614
1615 // Compute prism shape functions as a tensor-product.
1616 // Non-const here, because prism_indices might need to do some
1617 // flips before evaluating edge or face DoFs
1618 Point xi_eta {p(0),p(1)};
1619 Real zeta = p(2);
1620
1621 unsigned int i01, i2;
1622
1623 prism_indices(elem, totalorder, i, xi_eta, zeta, i01, i2);
1624
1625 // We'll use the 2D Tri to handle any basis function flipping
1626 // needed in xi/eta for triangle face+edge bases.
1627 Tri3 tri;
1628
1629 // We pinky swear not to modify these nodes
1630 Elem & e = const_cast<Elem &>(*elem);
1631 if (i2 == 0)
1632 {
1633 tri.set_node(0, e.node_ptr(0));
1634 tri.set_node(1, e.node_ptr(1));
1635 tri.set_node(2, e.node_ptr(2));
1636 }
1637 else if (i2 == 1)
1638 {
1639 tri.set_node(0, e.node_ptr(3));
1640 tri.set_node(1, e.node_ptr(4));
1641 tri.set_node(2, e.node_ptr(5));
1642 }
1643 else
1644 {
1645 // For interior DoFs, no flipping is necessary or done; we
1646 // can just evaluate on any triangle ... but *not* the
1647 // obvious 9,10,11 triangle, because that might not exist
1648 // if we have L2_HIERARCHIC on Prism6.
1649 tri.set_node(0, e.node_ptr(0));
1650 tri.set_node(1, e.node_ptr(1));
1651 tri.set_node(2, e.node_ptr(2));
1652
1653 // For square face DoFs, prism_indices handles flipping,
1654 // and we *can't* override that in the tri shape call.
1655 if (i01 > 2 && i01 < 3u*totalorder)
1656 {
1657 // %(p-1) to find the edge number, %2 for even vs odd
1658 const bool odd_basis = ((i01-3)%(totalorder-1))%2;
1659 if (odd_basis)
1660 {
1661 const int tri_edge = (i01-3)/(totalorder-1);
1662 // Flip nodes now to avoid triggering a shape
1663 // function flip later
1664 if (tri.point(tri_edge) > tri.point((tri_edge+1)%3))
1665 {
1666 Node * n = tri.node_ptr(tri_edge);
1667 tri.set_node(tri_edge, tri.node_ptr((tri_edge+1)%3));
1668 tri.set_node((tri_edge+1)%3, n);
1669 }
1670 }
1671 }
1672 }
1673
1674 return (FE<2,L2_HIERARCHIC>::shape(&tri, totalorder, i01, xi_eta)*
1675 FE<1,L2_HIERARCHIC>::shape(EDGE2, totalorder, i2, zeta));
1676 }
1677
1678 case TET4:
1679 libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
1680 libmesh_fallthrough();
1681 case TET10:
1682 libmesh_assert (T == L2_HIERARCHIC || totalorder < 3);
1683 libmesh_fallthrough();
1684 case TET14:
1685 {
1686 const Real zeta[4] = { 1 - p(0) - p(1) - p(2), p(0), p(1), p(2) };
1687
1688 // Nodal DoFs
1689 if (i < 4)
1690 return zeta[i];
1691
1692 // Edge DoFs
1693 else if (i < 6u*totalorder - 2u)
1694 {
1695 const unsigned int edge_num = (i - 4) / (totalorder - 1u);
1696 // const int edge_node = edge_num + 4;
1697 const unsigned int basisorder = i - 2 - ((totalorder - 1u) * edge_num);
1698
1699 const unsigned int edgevertex0 = Tet4::edge_nodes_map[edge_num][0],
1700 edgevertex1 = Tet4::edge_nodes_map[edge_num][1];
1701
1702 // Get factors to account for edge-flipping
1703 Real flip = 1;
1704 if (basisorder%2 &&
1705 elem->positive_edge_orientation(edge_num))
1706 flip = -1;
1707
1708 const Real crossval = zeta[edgevertex0] + zeta[edgevertex1];
1709 const Real edgenumerator = zeta[edgevertex1] - zeta[edgevertex0];
1710
1711 if (crossval == 0.) // Yes, exact comparison; we seem numerically stable otherwise
1712 {
1713 unsigned int basisfactorial = 1.;
1714 for (unsigned int n=2; n <= basisorder; ++n)
1715 basisfactorial *= n;
1716
1717 return flip * std::pow(edgenumerator, basisorder) / basisfactorial;
1718 }
1719
1720 const Real edgeval = edgenumerator / crossval;
1721 const Real crossfunc = std::pow(crossval, basisorder);
1722
1723 return flip * crossfunc *
1724 FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
1725 basisorder, edgeval);
1726 }
1727
1728 // Face DoFs
1729 else if (i < 2u*totalorder*totalorder + 2u)
1730 {
1731 const int dofs_per_face = (totalorder - 1u) * (totalorder - 2u) / 2;
1732 const int face_num = (i - (6u*totalorder - 2u)) / dofs_per_face;
1733
1734 const std::vector<unsigned int> face_vertex =
1735 elem->oriented_face_vertices(face_num);
1736 const Real zeta0 = zeta[face_vertex[0]],
1737 zeta1 = zeta[face_vertex[1]],
1738 zeta2 = zeta[face_vertex[2]];
1739
1740 const unsigned int basisnum =
1741 i - 4 -
1742 (totalorder - 1u) * /*n_edges*/6 -
1743 (dofs_per_face * face_num);
1744
1745 const unsigned int exp0 = triangular_number_column[basisnum] + 1;
1746 const unsigned int exp1 = triangular_number_row[basisnum] + 1 -
1747 triangular_number_column[basisnum];
1748
1749 Real returnval = 1;
1750 for (unsigned int n = 0; n != exp0; ++n)
1751 returnval *= zeta0;
1752 for (unsigned int n = 0; n != exp1; ++n)
1753 returnval *= zeta1;
1754 returnval *= zeta2;
1755 return returnval;
1756 }
1757
1758 // Interior DoFs
1759 else
1760 {
1761 const unsigned int basisnum = i - 2u*totalorder*totalorder - 2u;
1762 const unsigned int exp0 = tetrahedral_number_column[basisnum] + 1;
1763 const unsigned int exp1 = tetrahedral_number_row[basisnum] + 1 -
1764 tetrahedral_number_column[basisnum] -
1765 tetrahedral_number_page[basisnum];
1766 const unsigned int exp2 = tetrahedral_number_page[basisnum] + 1;
1767
1768 Real returnval = 1;
1769 for (unsigned int n = 0; n != exp0; ++n)
1770 returnval *= zeta[0];
1771 for (unsigned int n = 0; n != exp1; ++n)
1772 returnval *= zeta[1];
1773 for (unsigned int n = 0; n != exp2; ++n)
1774 returnval *= zeta[2];
1775 returnval *= zeta[3];
1776 return returnval;
1777 }
1778 }
1779
1780 default:
1781 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1782 }
1783
1784#else // LIBMESH_DIM != 3
1785 libmesh_ignore(elem, order, i, p, add_p_level);
1786 libmesh_not_implemented();
1787#endif
1788}
1789
1790
1791
1792template <FEFamily T>
1793Real fe_hierarchic_3D_shape_deriv(const Elem * elem,
1794 const Order order,
1795 const unsigned int i,
1796 const unsigned int j,
1797 const Point & p,
1798 const bool add_p_level)
1799{
1800 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,T>::shape);
1801}
1802
1803
1804
1805#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1806
1807template <FEFamily T>
1808Real fe_hierarchic_3D_shape_second_deriv(const Elem * elem,
1809 const Order order,
1810 const unsigned int i,
1811 const unsigned int j,
1812 const Point & p,
1813 const bool add_p_level)
1814{
1815 return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
1816 FE<3,T>::shape_deriv);
1817}
1818
1819#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1820
1821
1822} // anonymous namespace
int
void ErrorVector unsigned int
Definition adjoints_ex3.C:360
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::Elem::set_node
virtual Node *& set_node(const unsigned int i)
Definition elem.h:2688
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition elem.h:2583
libMesh::Elem::local_side_node
virtual unsigned int local_side_node(unsigned int side, unsigned int side_node) const =0
libMesh::Elem::positive_edge_orientation
bool positive_edge_orientation(const unsigned int i) const
Definition elem.C:4037
libMesh::Elem::oriented_face_vertices
std::vector< unsigned int > oriented_face_vertices(const unsigned int i) const
Definition elem.C:4078
libMesh::Elem::node_ptr
const Node * node_ptr(const unsigned int i) const
Definition elem.h:2637
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Elem::master_point
virtual Point master_point(const unsigned int i) const =0
libMesh::Elem::p_level
unsigned int p_level() const
Definition elem.h:3246
libMesh::Elem::face_orientation
unsigned int face_orientation(const unsigned int i) const
Definition elem.C:4053
libMesh::Elem::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i)=0
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition fe_type.h:197
libMesh::FEType::order
OrderWrapper order
The approximation order of the element (at 0 p-refinement level).
Definition fe_type.h:203
libMesh::FE
A specific instantiation of the FEBase class.
Definition fe.h:128
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
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)
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::IntRange
The IntRange templated class is intended to make it easy to loop over integers which are indices of a...
Definition int_range.h:54
libMesh::Node
A Node is like a Point, but with more information.
Definition node.h:55
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:41
libMesh::Quad4
The QUAD4 is an element in 2D composed of 4 nodes.
Definition face_quad4.h:54
libMesh::Tet4::edge_nodes_map
static constexpr 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:207
libMesh::Tri3
The Tri3 is an element in 2D composed of 3 nodes.
Definition face_tri3.h:62
b
static const Real b
Definition disjoint_neighbor_test.C:24
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::tetrahedral_number_column
const unsigned char tetrahedral_number_column[]
Definition number_lookups.C:421
libMesh::tetrahedral_number_row
const unsigned char tetrahedral_number_row[]
Definition number_lookups.C:488
libMesh::prism_number_triangle
const unsigned char prism_number_triangle[]
Definition number_lookups.C:621
libMesh::LIBMESH_DEFAULT_VECTORIZED_FE
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
Definition fe_bernstein_shape_0D.C:30
libMesh::prism_number_page
const unsigned char prism_number_page[]
Definition number_lookups.C:726
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't be...
Definition fe.C:911
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition libmesh_common.h:598
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::libmesh_assert
libmesh_assert(ctx)
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:303
libMesh::triangular_number_row
const unsigned char triangular_number_row[]
Definition number_lookups.C:28
libMesh::cube_number_row
const unsigned char cube_number_row[]
Definition number_lookups.C:196
libMesh::cube_number_page
const unsigned char cube_number_page[]
Definition number_lookups.C:308
libMesh::cube_number_column
const unsigned char cube_number_column[]
Definition number_lookups.C:84
libMesh::square_number_row
const unsigned char square_number_row[]
Definition number_lookups.C:69
libMesh::square_number_column
const unsigned char square_number_column[]
Definition number_lookups.C:56
libMesh::triangular_number_column
const unsigned char triangular_number_column[]
Definition number_lookups.C:41
libMesh::tetrahedral_number_page
const unsigned char tetrahedral_number_page[]
Definition number_lookups.C:550
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:148
libMesh::make_range
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition int_range.h:176
Generated on Thu Nov 12 2026 20:15:21 for libMesh by doxygen 1.9.8