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