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 :
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 293257527 : quad_indices(const Elem * elem,
69 : const unsigned int totalorder,
70 : const unsigned int i)
71 : {
72 24088245 : 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 48177277 : if (i == 0)
83 1382150 : { i0 = 0; i1 = 0; }
84 45412934 : else if (i == 1)
85 1382150 : { i0 = 1; i1 = 0; }
86 42648591 : else if (i == 2)
87 1382150 : { i0 = 1; i1 = 1; }
88 39884248 : else if (i == 3)
89 1382150 : { i0 = 0; i1 = 1; }
90 : // Edge DoFs
91 225721487 : else if (i < totalorder + 3u)
92 34444787 : { i0 = i - 2; i1 = 0; }
93 191276700 : else if (i < 2u*totalorder + 2)
94 34444787 : { i0 = 1; i1 = i - totalorder - 1; }
95 156831913 : else if (i < 3u*totalorder + 1)
96 34444787 : { i0 = i - 2u*totalorder; i1 = 1; }
97 122387126 : else if (i < 4u*totalorder)
98 34444787 : { i0 = 0; i1 = i - 3u*totalorder + 1; }
99 : // Interior DoFs
100 : else
101 : {
102 87942339 : unsigned int basisnum = i - 4*totalorder;
103 87942339 : i0 = square_number_column[basisnum] + 2;
104 87942339 : i1 = square_number_row[basisnum] + 2;
105 : }
106 :
107 : // Flip odd degree of freedom values if necessary
108 : // to keep continuity on sides
109 24088245 : Real f = 1.;
110 :
111 293257527 : if ((i0%2) && (i0 > 2) && (i1 == 0))
112 12074514 : f = elem->positive_edge_orientation(0)?-1.:1.;
113 281183013 : else if ((i0%2) && (i0>2) && (i1 == 1))
114 12074514 : f = !elem->positive_edge_orientation(2)?-1.:1.;
115 269108499 : else if ((i0 == 0) && (i1%2) && (i1>2))
116 12074514 : f = !elem->positive_edge_orientation(3)?-1.:1.;
117 257033985 : else if ((i0 == 1) && (i1%2) && (i1>2))
118 12074514 : f = elem->positive_edge_orientation(1)?-1.:1.;
119 :
120 317345772 : return {i0, i1, f};
121 : }
122 :
123 : } // anonymous namespace
124 :
125 :
126 :
127 : namespace libMesh
128 : {
129 :
130 :
131 6299475 : LIBMESH_DEFAULT_VECTORIZED_FE(2,HIERARCHIC)
132 15451578 : LIBMESH_DEFAULT_VECTORIZED_FE(2,L2_HIERARCHIC)
133 14645806 : 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 211781835 : 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 211781835 : 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 2047163678 : 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 2047163678 : 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 42832020 : 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 3569614 : libmesh_assert(elem);
225 42832020 : const ElemType type = elem->type();
226 :
227 46401634 : const Order totalorder = order + add_p_level*elem->p_level();
228 :
229 42832020 : const unsigned int dofs_per_side = totalorder+1u;
230 :
231 39262406 : switch (type)
232 : {
233 41184108 : case TRI6:
234 : case TRI7:
235 : {
236 3432222 : 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 3432222 : Real f = 1.;
243 :
244 41184108 : const Real zeta1 = p(0);
245 41184108 : const Real zeta2 = p(1);
246 41184108 : const Real zeta0 = 1. - zeta1 - zeta2;
247 :
248 41184108 : if (zeta1 > zeta2 && zeta0 > zeta2) // side 0
249 : {
250 13719948 : if (i >= dofs_per_side)
251 761988 : return 0;
252 :
253 4573316 : if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
254 8240 : return 1;
255 :
256 8045560 : if ((i < 2 || i % 2) &&
257 3571194 : elem->positive_edge_orientation(0))
258 130248 : f = -1;
259 :
260 4847120 : return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i, f*(zeta1-zeta0));
261 : }
262 27464160 : else if (zeta1 > zeta0 && zeta2 > zeta0) // side 1
263 : {
264 13741794 : if (i < dofs_per_side ||
265 9161196 : i >= 2*dofs_per_side)
266 763940 : return 0;
267 :
268 4580598 : if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
269 8250 : return 1;
270 :
271 4481538 : const unsigned int side_i = i - dofs_per_side;
272 :
273 8057546 : if ((side_i < 2 || side_i % 2) &&
274 3576008 : elem->positive_edge_orientation(1))
275 162178 : f = -1;
276 :
277 4855258 : return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, side_i, f*(zeta2-zeta1));
278 : }
279 : else // side 2
280 : {
281 1143330 : libmesh_assert (zeta2 >= zeta1 && zeta0 >= zeta1); // On a corner???
282 :
283 13722366 : if (i < 2*dofs_per_side)
284 762220 : return 0;
285 :
286 4574122 : if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
287 8230 : return 1;
288 :
289 4475202 : const unsigned int side_i = i - 2*dofs_per_side;
290 :
291 8047102 : if ((side_i < 2 || side_i % 2) &&
292 3571900 : elem->positive_edge_orientation(2))
293 162178 : f = -1;
294 :
295 4848082 : return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, side_i, f*(zeta0-zeta2));
296 : }
297 : }
298 1647912 : case QUAD8:
299 : case QUADSHELL8:
300 : case QUAD9:
301 : case QUADSHELL9:
302 : {
303 137392 : 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 137392 : Real f = 1.;
310 :
311 1647912 : const Real xi = p(0), eta = p(1);
312 1647912 : if (eta < xi)
313 : {
314 815232 : if (eta < -xi) // side 0
315 : {
316 401168 : if (i >= dofs_per_side)
317 24444 : return 0;
318 :
319 100292 : if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
320 4120 : return 1;
321 :
322 88358 : if ((i < 2 || i % 2) &&
323 36206 : elem->positive_edge_orientation(0))
324 1480 : f = -1;
325 :
326 56180 : return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i, f*xi);
327 : }
328 : else // side 1
329 : {
330 414064 : if (i < dofs_per_side ||
331 310548 : i >= 2*dofs_per_side)
332 26838 : return 0;
333 :
334 103516 : if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
335 4130 : return 1;
336 :
337 55336 : const unsigned int side_i = i - dofs_per_side;
338 :
339 93866 : if ((side_i < 2 || side_i % 2) &&
340 38530 : elem->positive_edge_orientation(1))
341 740 : f = -1;
342 :
343 60152 : return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, side_i, f*eta);
344 : }
345 : }
346 : else // xi < eta
347 : {
348 832680 : if (eta > -xi) // side 2
349 : {
350 411344 : if (i < 2*dofs_per_side ||
351 205672 : i >= 3*dofs_per_side)
352 25188 : return 0;
353 :
354 102836 : if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
355 4120 : return 1;
356 :
357 54676 : const unsigned int side_i = i - 2*dofs_per_side;
358 :
359 92802 : if ((side_i < 2 || side_i % 2) &&
360 38126 : !elem->positive_edge_orientation(2))
361 1110 : f = -1;
362 :
363 58952 : return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, side_i, f*xi);
364 : }
365 : else // side 3
366 : {
367 421336 : if (i < 3*dofs_per_side)
368 26574 : return 0;
369 :
370 105334 : if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
371 4130 : return 1;
372 :
373 57124 : const unsigned int side_i = i - 3*dofs_per_side;
374 :
375 96968 : if ((side_i < 2 || side_i % 2) &&
376 39844 : !elem->positive_edge_orientation(3))
377 740 : f = -1;
378 :
379 61852 : 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 42878854 : 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 42878854 : 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 1088110698 : 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 93504117 : Real flip = 1;
924 1088110698 : if (basisorder%2 && (elem.point(noden) > elem.point((noden+1)%3)))
925 13122225 : 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 1088110698 : if (crossval == 0.)
932 : {
933 35982 : unsigned int basisfactorial = 1.;
934 1190704 : for (unsigned int n=2; n <= basisorder; ++n)
935 727972 : basisfactorial *= n;
936 :
937 462732 : return std::pow(edgenumerator, basisorder) / basisfactorial;
938 : }
939 : // Experimentally, as c -> 0, n propto c, I'm still seeing good
940 : // behavior from the default implementation below:
941 :
942 1087647966 : const Real edgeval = edgenumerator / crossval;
943 93468135 : const Real crossfunc = std::pow(crossval, basisorder);
944 :
945 1087647966 : return flip * crossfunc *
946 1087647966 : FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
947 1087647966 : basisorder, edgeval);
948 : }
949 :
950 : template <FEFamily T>
951 2258945513 : Real fe_hierarchic_2D_shape(const Elem * elem,
952 : const Order order,
953 : const unsigned int i,
954 : const Point & p,
955 : const bool add_p_level)
956 : {
957 193094111 : libmesh_assert(elem);
958 :
959 2452041971 : const Order totalorder = order + add_p_level*elem->p_level();
960 193094111 : libmesh_assert_greater (totalorder, 0);
961 :
962 2258945513 : switch (elem->type())
963 : {
964 1982945368 : case TRI3:
965 : case TRISHELL3:
966 : case TRI6:
967 : case TRI7:
968 : {
969 1982945368 : const Real zeta1 = p(0);
970 1982945368 : const Real zeta2 = p(1);
971 1982945368 : const Real zeta0 = 1. - zeta1 - zeta2;
972 :
973 170233916 : libmesh_assert_less (i, (totalorder+1u)*(totalorder+2u)/2);
974 9735327 : libmesh_assert (T == L2_HIERARCHIC || elem->type() == TRI6 ||
975 : elem->type() == TRI7 || totalorder < 2);
976 :
977 : // Vertex DoFs
978 1982945368 : if (i == 0)
979 17972469 : return zeta0;
980 1773040271 : else if (i == 1)
981 17972469 : return zeta1;
982 1563135174 : else if (i == 2)
983 17972469 : return zeta2;
984 : // Edge DoFs
985 1353230077 : else if (i < totalorder + 2u)
986 : {
987 362703566 : const unsigned int basisorder = i - 1;
988 :
989 362703566 : const Real crossval = zeta0 + zeta1;
990 362703566 : const Real edgenumerator = zeta1 - zeta0;
991 :
992 362703566 : return fe_triangle_helper(*elem, edgenumerator, crossval,
993 362703566 : basisorder, totalorder, 0);
994 : }
995 990526511 : else if (i < 2u*totalorder + 1)
996 : {
997 362703566 : const unsigned int basisorder = i - totalorder;
998 :
999 362703566 : const Real crossval = zeta2 + zeta1;
1000 362703566 : const Real edgenumerator = zeta2 - zeta1;
1001 :
1002 362703566 : return fe_triangle_helper(*elem, edgenumerator, crossval,
1003 362703566 : basisorder, totalorder, 1);
1004 : }
1005 627822945 : else if (i < 3u*totalorder)
1006 : {
1007 362703566 : const unsigned int basisorder = i - (2u*totalorder) + 1;
1008 :
1009 362703566 : const Real crossval = zeta0 + zeta2;
1010 362703566 : const Real edgenumerator = zeta0 - zeta2;
1011 :
1012 362703566 : return fe_triangle_helper(*elem, edgenumerator, crossval,
1013 362703566 : basisorder, totalorder, 2);
1014 : }
1015 : // Interior DoFs
1016 : else
1017 : {
1018 265119379 : const unsigned int basisnum = i - (3u*totalorder);
1019 265119379 : unsigned int exp0 = triangular_number_column[basisnum] + 1;
1020 265119379 : unsigned int exp1 = triangular_number_row[basisnum] + 1 -
1021 : triangular_number_column[basisnum];
1022 :
1023 22812392 : Real returnval = 1;
1024 600204608 : for (unsigned int n = 0; n != exp0; ++n)
1025 335085229 : returnval *= zeta0;
1026 600204608 : for (unsigned int n = 0; n != exp1; ++n)
1027 335085229 : returnval *= zeta1;
1028 265119379 : returnval *= zeta2;
1029 265119379 : return returnval;
1030 : }
1031 : }
1032 :
1033 : // Hierarchic shape functions on the quadrilateral.
1034 268052965 : case QUAD4:
1035 : case QUADSHELL4:
1036 14913015 : libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
1037 : libmesh_fallthrough();
1038 : case QUAD8:
1039 : case QUADSHELL8:
1040 : case QUAD9:
1041 : case QUADSHELL9:
1042 : {
1043 : // Compute quad shape functions as a tensor-product
1044 276000145 : auto [i0, i1, f] = quad_indices(elem, totalorder, i);
1045 :
1046 298861127 : return f*(FE<1,T>::shape(EDGE3, totalorder, i0, p(0))*
1047 298861127 : FE<1,T>::shape(EDGE3, totalorder, i1, p(1)));
1048 : }
1049 :
1050 0 : default:
1051 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(elem->type()));
1052 : }
1053 :
1054 : return 0.;
1055 : }
1056 :
1057 :
1058 : template <FEFamily T>
1059 149215830 : Real fe_hierarchic_2D_shape_deriv(const Elem * elem,
1060 : const Order order,
1061 : const unsigned int i,
1062 : const unsigned int j,
1063 : const Point & p,
1064 : const bool add_p_level)
1065 : {
1066 12287378 : libmesh_assert(elem);
1067 :
1068 149215830 : const ElemType type = elem->type();
1069 :
1070 161503208 : const Order totalorder = order + add_p_level*elem->p_level();
1071 :
1072 12287378 : libmesh_assert_greater (totalorder, 0);
1073 :
1074 136928452 : switch (type)
1075 : {
1076 : // 1st & 2nd-order Hierarchics.
1077 131958448 : case TRI3:
1078 : case TRISHELL3:
1079 : case TRI6:
1080 : case TRI7:
1081 : {
1082 131958448 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,T>::shape);
1083 : }
1084 :
1085 16087804 : case QUAD4:
1086 : case QUADSHELL4:
1087 58472 : libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
1088 : libmesh_fallthrough();
1089 : case QUAD8:
1090 : case QUADSHELL8:
1091 : case QUAD9:
1092 : case QUADSHELL9:
1093 : {
1094 : // Compute quad shape functions as a tensor-product
1095 17257382 : auto [i0, i1, f] = quad_indices(elem, totalorder, i);
1096 :
1097 17257382 : switch (j)
1098 : {
1099 : // d()/dxi
1100 10326014 : case 0:
1101 11063965 : return f*(FE<1,T>::shape_deriv(EDGE3, totalorder, i0, 0, p(0))*
1102 11063965 : FE<1,T>::shape (EDGE3, totalorder, i1, p(1)));
1103 :
1104 : // d()/deta
1105 6931368 : case 1:
1106 7421467 : return f*(FE<1,T>::shape (EDGE3, totalorder, i0, p(0))*
1107 7421467 : FE<1,T>::shape_deriv(EDGE3, totalorder, i1, 0, p(1)));
1108 :
1109 0 : default:
1110 0 : libmesh_error_msg("Invalid derivative index j = " << j);
1111 : }
1112 : }
1113 :
1114 0 : default:
1115 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
1116 : }
1117 :
1118 : return 0.;
1119 : }
1120 :
1121 :
1122 :
1123 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1124 :
1125 : template <FEFamily T>
1126 3215832 : Real fe_hierarchic_2D_shape_second_deriv(const Elem * elem,
1127 : const Order order,
1128 : const unsigned int i,
1129 : const unsigned int j,
1130 : const Point & p,
1131 : const bool add_p_level)
1132 : {
1133 38972673 : return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
1134 3215832 : FE<2,T>::shape_deriv);
1135 : }
1136 :
1137 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1138 :
1139 : } // anonymous namespace
|
