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