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