libMesh
Loading...
Searching...
No Matches
fe_bernstein_shape_3D.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3
4// This library is free software; you can redistribute it and/or
5// modify it under the terms of the GNU Lesser General Public
6// License as published by the Free Software Foundation; either
7// version 2.1 of the License, or (at your option) any later version.
8
9// This library is distributed in the hope that it will be useful,
10// but WITHOUT ANY WARRANTY; without even the implied warranty of
11// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12// Lesser General Public License for more details.
13
14// You should have received a copy of the GNU Lesser General Public
15// License along with this library; if not, write to the Free Software
16// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17
18
19#include "libmesh/libmesh_config.h"
20
21#ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
22
23// libmesh includes
24#include "libmesh/fe.h"
25#include "libmesh/elem.h"
26#include "libmesh/enum_to_string.h"
27
28
29namespace {
30 using namespace libMesh;
31
32 // Indices and coefficients for the HEX20
33 //
34 // The only way to make any sense of this
35 // is to look at the mgflo/mg2/mgf documentation
36 // and make the cut-out cube!
37 // 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
38 static const unsigned int hex20_i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
39 static const unsigned int hex20_i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
40 static const unsigned int hex20_i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2};
41 //To compute the hex20 shape functions the original shape functions for hex27 are used.
42 //hex20_scalx[i] tells how often the original x-th shape function has to be added to the original i-th shape function
43 //to compute the new i-th shape function for hex20
44 //example: B_0^HEX20 = B_0^HEX27 - 0.25*B_20^HEX27 - 0.25*B_21^HEX27 + 0*B_22^HEX27 + 0*B_23^HEX27 - 0.25*B_24^HEX27 + 0*B_25^HEX27 - 0.25*B_26^HEX27
45 // B_0^HEX20 = B_0^HEX27 + hex20_scal20[0]*B_20^HEX27 + hex20_scal21[0]*B_21^HEX27 + ...
46
47 static const Real hex20_scal20[] = {-0.25, -0.25, -0.25, -0.25, 0, 0, 0, 0, 0.5, 0.5, 0.5, 0.5, 0, 0, 0, 0, 0, 0, 0, 0};
48 static const Real hex20_scal21[] = {-0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0, 0, 0};
49 static const Real hex20_scal22[] = {0, -0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0, 0};
50 static const Real hex20_scal23[] = {0, 0, -0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0};
51 static const Real hex20_scal24[] = {-0.25, 0, 0, -0.25, -0.25, 0, 0, -0.25, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0, 0, 0, 0.5};
52 static const Real hex20_scal25[] = {0, 0, 0, 0, -0.25, -0.25, -0.25, -0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0.5, 0.5};
53 static const Real hex20_scal26[] = {-0.25, -0.25, -0.25, -0.25, -0.25, -0.25, -0.25, -0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25};
54
55 Point get_min_point(const Elem * elem,
56 unsigned int a, unsigned int b,
57 unsigned int c, unsigned int d)
58 {
59 return std::min(std::min(elem->point(a),elem->point(b)),
60 std::min(elem->point(c),elem->point(d)));
61 }
62} // anonymous namespace
63
64
65
66namespace libMesh
67{
68
69
70LIBMESH_DEFAULT_VECTORIZED_FE(3,BERNSTEIN)
71
72
73template <>
74Real FE<3,BERNSTEIN>::shape(const Elem * elem,
75 const Order order,
76 const unsigned int i,
77 const Point & p,
78 const bool add_p_level)
79{
80
81#if LIBMESH_DIM == 3
82
83 libmesh_assert(elem);
84 const ElemType type = elem->type();
85
86 const Order totalorder =
87 order + add_p_level*elem->p_level();
88
89 auto hex_remap = [i, elem] (const Point & p_in,
90 const unsigned int * hex_i0,
91 const unsigned int * hex_i1,
92 const unsigned int * hex_i2) {
93 // Compute hex shape functions as a tensor-product
94 const Real xi = p_in(0);
95 const Real eta = p_in(1);
96 const Real zeta = p_in(2);
97 Point p_to_remap = p_in;
98 Real & xi_mapped = p_to_remap(0);
99 Real & eta_mapped = p_to_remap(1);
100 Real & zeta_mapped = p_to_remap(2);
101
102 // handle the edge orientation
103 {
104 // Edge 0
105 if ((hex_i0[i] >= 2) && (hex_i1[i] == 0) && (hex_i2[i] == 0))
106 {
107 if (elem->positive_edge_orientation(0))
108 xi_mapped = -xi;
109 }
110 // Edge 1
111 else if ((hex_i0[i] == 1) && (hex_i1[i] >= 2) && (hex_i2[i] == 0))
112 {
113 if (elem->positive_edge_orientation(1))
114 eta_mapped = -eta;
115 }
116 // edge 2
117 else if ((hex_i0[i] >= 2) && (hex_i1[i] == 1) && (hex_i2[i] == 0))
118 {
119 if (!elem->positive_edge_orientation(2))
120 xi_mapped = -xi;
121 }
122 // edge 3
123 else if ((hex_i0[i] == 0) && (hex_i1[i] >= 2) && (hex_i2[i] == 0))
124 {
125 if (elem->positive_edge_orientation(3))
126 eta_mapped = -eta;
127 }
128 // edge 4
129 else if ((hex_i0[i] == 0) && (hex_i1[i] == 0) && (hex_i2[i] >=2 ))
130 {
131 if (elem->positive_edge_orientation(4))
132 zeta_mapped = -zeta;
133 }
134 // edge 5
135 else if ((hex_i0[i] == 1) && (hex_i1[i] == 0) && (hex_i2[i] >=2 ))
136 {
137 if (elem->positive_edge_orientation(5))
138 zeta_mapped = -zeta;
139 }
140 // edge 6
141 else if ((hex_i0[i] == 1) && (hex_i1[i] == 1) && (hex_i2[i] >=2 ))
142 {
143 if (elem->positive_edge_orientation(6))
144 zeta_mapped = -zeta;
145 }
146 // edge 7
147 else if ((hex_i0[i] == 0) && (hex_i1[i] == 1) && (hex_i2[i] >=2 ))
148 {
149 if (elem->positive_edge_orientation(7))
150 zeta_mapped = -zeta;
151 }
152 // edge 8
153 else if ((hex_i0[i] >=2 ) && (hex_i1[i] == 0) && (hex_i2[i] == 1))
154 {
155 if (elem->positive_edge_orientation(8))
156 xi_mapped = -xi;
157 }
158 // edge 9
159 else if ((hex_i0[i] == 1) && (hex_i1[i] >=2 ) && (hex_i2[i] == 1))
160 {
161 if (elem->positive_edge_orientation(9))
162 eta_mapped = -eta;
163 }
164 // edge 10
165 else if ((hex_i0[i] >=2 ) && (hex_i1[i] == 1) && (hex_i2[i] == 1))
166 {
167 if (!elem->positive_edge_orientation(10))
168 xi_mapped = -xi;
169 }
170 // edge 11
171 else if ((hex_i0[i] == 0) && (hex_i1[i] >=2 ) && (hex_i2[i] == 1))
172 {
173 if (elem->positive_edge_orientation(11))
174 eta_mapped = -eta;
175 }
176 }
177
178
179 // handle the face orientation
180 {
181 // face 0
182 if ((hex_i2[i] == 0) && (hex_i0[i] >= 2) && (hex_i1[i] >= 2))
183 {
184 const Point min_point = get_min_point(elem, 1, 2, 0, 3);
185
186 if (elem->point(0) == min_point)
187 if (elem->positive_face_orientation(0))
188 {
189 // case 1
190 xi_mapped = xi;
191 eta_mapped = eta;
192 }
193 else
194 {
195 // case 2
196 xi_mapped = eta;
197 eta_mapped = xi;
198 }
199
200 else if (elem->point(3) == min_point)
201 if (elem->positive_face_orientation(0))
202 {
203 // case 3
204 xi_mapped = -eta;
205 eta_mapped = xi;
206 }
207 else
208 {
209 // case 4
210 xi_mapped = xi;
211 eta_mapped = -eta;
212 }
213
214 else if (elem->point(2) == min_point)
215 if (elem->positive_face_orientation(0))
216 {
217 // case 5
218 xi_mapped = -xi;
219 eta_mapped = -eta;
220 }
221 else
222 {
223 // case 6
224 xi_mapped = -eta;
225 eta_mapped = -xi;
226 }
227
228 else if (elem->point(1) == min_point)
229 {
230 if (elem->positive_face_orientation(0))
231 {
232 // case 7
233 xi_mapped = eta;
234 eta_mapped = -xi;
235 }
236 else
237 {
238 // Case 8
239 xi_mapped = -xi;
240 eta_mapped = eta;
241 }
242 }
243 }
244
245
246 // Face 1
247 else if ((hex_i1[i] == 0) && (hex_i0[i] >= 2) && (hex_i2[i] >= 2))
248 {
249 const Point min_point = get_min_point(elem, 0, 1, 5, 4);
250
251 if (elem->point(0) == min_point)
252 if (!elem->positive_face_orientation(1))
253 {
254 // Case 1
255 xi_mapped = xi;
256 zeta_mapped = zeta;
257 }
258 else
259 {
260 // Case 2
261 xi_mapped = zeta;
262 zeta_mapped = xi;
263 }
264
265 else if (elem->point(1) == min_point)
266 if (!elem->positive_face_orientation(1))
267 {
268 // Case 3
269 xi_mapped = zeta;
270 zeta_mapped = -xi;
271 }
272 else
273 {
274 // Case 4
275 xi_mapped = -xi;
276 zeta_mapped = zeta;
277 }
278
279 else if (elem->point(5) == min_point)
280 if (!elem->positive_face_orientation(1))
281 {
282 // Case 5
283 xi_mapped = -xi;
284 zeta_mapped = -zeta;
285 }
286 else
287 {
288 // Case 6
289 xi_mapped = -zeta;
290 zeta_mapped = -xi;
291 }
292
293 else if (elem->point(4) == min_point)
294 {
295 if (!elem->positive_face_orientation(1))
296 {
297 // Case 7
298 xi_mapped = -xi;
299 zeta_mapped = zeta;
300 }
301 else
302 {
303 // Case 8
304 xi_mapped = xi;
305 zeta_mapped = -zeta;
306 }
307 }
308 }
309
310
311 // Face 2
312 else if ((hex_i0[i] == 1) && (hex_i1[i] >= 2) && (hex_i2[i] >= 2))
313 {
314 const Point min_point = get_min_point(elem, 1, 2, 6, 5);
315
316 if (elem->point(1) == min_point)
317 if (!elem->positive_face_orientation(2))
318 {
319 // Case 1
320 eta_mapped = eta;
321 zeta_mapped = zeta;
322 }
323 else
324 {
325 // Case 2
326 eta_mapped = zeta;
327 zeta_mapped = eta;
328 }
329
330 else if (elem->point(2) == min_point)
331 if (!elem->positive_face_orientation(2))
332 {
333 // Case 3
334 eta_mapped = zeta;
335 zeta_mapped = -eta;
336 }
337 else
338 {
339 // Case 4
340 eta_mapped = -eta;
341 zeta_mapped = zeta;
342 }
343
344 else if (elem->point(6) == min_point)
345 if (!elem->positive_face_orientation(2))
346 {
347 // Case 5
348 eta_mapped = -eta;
349 zeta_mapped = -zeta;
350 }
351 else
352 {
353 // Case 6
354 eta_mapped = -zeta;
355 zeta_mapped = -eta;
356 }
357
358 else if (elem->point(5) == min_point)
359 {
360 if (!elem->positive_face_orientation(2))
361 {
362 // Case 7
363 eta_mapped = -zeta;
364 zeta_mapped = eta;
365 }
366 else
367 {
368 // Case 8
369 eta_mapped = eta;
370 zeta_mapped = -zeta;
371 }
372 }
373 }
374
375
376 // Face 3
377 else if ((hex_i1[i] == 1) && (hex_i0[i] >= 2) && (hex_i2[i] >= 2))
378 {
379 const Point min_point = get_min_point(elem, 2, 3, 7, 6);
380
381 if (elem->point(3) == min_point)
382 if (elem->positive_face_orientation(3))
383 {
384 // Case 1
385 xi_mapped = xi;
386 zeta_mapped = zeta;
387 }
388 else
389 {
390 // Case 2
391 xi_mapped = zeta;
392 zeta_mapped = xi;
393 }
394
395 else if (elem->point(7) == min_point)
396 if (elem->positive_face_orientation(3))
397 {
398 // Case 3
399 xi_mapped = -zeta;
400 zeta_mapped = xi;
401 }
402 else
403 {
404 // Case 4
405 xi_mapped = xi;
406 zeta_mapped = -zeta;
407 }
408
409 else if (elem->point(6) == min_point)
410 if (elem->positive_face_orientation(3))
411 {
412 // Case 5
413 xi_mapped = -xi;
414 zeta_mapped = -zeta;
415 }
416 else
417 {
418 // Case 6
419 xi_mapped = -zeta;
420 zeta_mapped = -xi;
421 }
422
423 else if (elem->point(2) == min_point)
424 {
425 if (elem->positive_face_orientation(3))
426 {
427 // Case 7
428 xi_mapped = zeta;
429 zeta_mapped = -xi;
430 }
431 else
432 {
433 // Case 8
434 xi_mapped = -xi;
435 zeta_mapped = zeta;
436 }
437 }
438 }
439
440
441 // Face 4
442 else if ((hex_i0[i] == 0) && (hex_i1[i] >= 2) && (hex_i2[i] >= 2))
443 {
444 const Point min_point = get_min_point(elem, 3, 0, 4, 7);
445
446 if (elem->point(0) == min_point)
447 if (elem->positive_face_orientation(4))
448 {
449 // Case 1
450 eta_mapped = eta;
451 zeta_mapped = zeta;
452 }
453 else
454 {
455 // Case 2
456 eta_mapped = zeta;
457 zeta_mapped = eta;
458 }
459
460 else if (elem->point(4) == min_point)
461 if (elem->positive_face_orientation(4))
462 {
463 // Case 3
464 eta_mapped = -zeta;
465 zeta_mapped = eta;
466 }
467 else
468 {
469 // Case 4
470 eta_mapped = eta;
471 zeta_mapped = -zeta;
472 }
473
474 else if (elem->point(7) == min_point)
475 if (elem->positive_face_orientation(4))
476 {
477 // Case 5
478 eta_mapped = -eta;
479 zeta_mapped = -zeta;
480 }
481 else
482 {
483 // Case 6
484 eta_mapped = -zeta;
485 zeta_mapped = -eta;
486 }
487
488 else if (elem->point(3) == min_point)
489 {
490 if (elem->positive_face_orientation(4))
491 {
492 // Case 7
493 eta_mapped = zeta;
494 zeta_mapped = -eta;
495 }
496 else
497 {
498 // Case 8
499 eta_mapped = -eta;
500 zeta_mapped = zeta;
501 }
502 }
503 }
504
505
506 // Face 5
507 else if ((hex_i2[i] == 1) && (hex_i0[i] >= 2) && (hex_i1[i] >= 2))
508 {
509 const Point min_point = get_min_point(elem, 4, 5, 6, 7);
510
511 if (elem->point(4) == min_point)
512 if (!elem->positive_face_orientation(5))
513 {
514 // Case 1
515 xi_mapped = xi;
516 eta_mapped = eta;
517 }
518 else
519 {
520 // Case 2
521 xi_mapped = eta;
522 eta_mapped = xi;
523 }
524
525 else if (elem->point(5) == min_point)
526 if (!elem->positive_face_orientation(5))
527 {
528 // Case 3
529 xi_mapped = eta;
530 eta_mapped = -xi;
531 }
532 else
533 {
534 // Case 4
535 xi_mapped = -xi;
536 eta_mapped = eta;
537 }
538
539 else if (elem->point(6) == min_point)
540 if (!elem->positive_face_orientation(5))
541 {
542 // Case 5
543 xi_mapped = -xi;
544 eta_mapped = -eta;
545 }
546 else
547 {
548 // Case 6
549 xi_mapped = -eta;
550 eta_mapped = -xi;
551 }
552
553 else if (elem->point(7) == min_point)
554 {
555 if (!elem->positive_face_orientation(5))
556 {
557 // Case 7
558 xi_mapped = -eta;
559 eta_mapped = xi;
560 }
561 else
562 {
563 // Case 8
564 xi_mapped = xi;
565 eta_mapped = eta;
566 }
567 }
568 }
569 }
570
571 return p_to_remap;
572 };
573
574 switch (totalorder)
575 {
576 // 1st order Bernstein.
577 case FIRST:
578 {
579 switch (type)
580 {
581
582 // Bernstein shape functions on the tetrahedron.
583 case TET4:
584 case TET10:
585 case TET14:
586 {
587 libmesh_assert_less (i, 4);
588
589 // Area coordinates
590 const Real zeta1 = p(0);
591 const Real zeta2 = p(1);
592 const Real zeta3 = p(2);
593 const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
594
595 switch(i)
596 {
597 case 0: return zeta0;
598 case 1: return zeta1;
599 case 2: return zeta2;
600 case 3: return zeta3;
601
602 default:
603 libmesh_error_msg("Invalid shape function index i = " << i);
604 }
605 }
606
607 // Bernstein shape functions on the hexahedral.
608 case HEX8:
609 case HEX20:
610 case HEX27:
611 {
612 libmesh_assert_less (i, 8);
613
614 // Compute hex shape functions as a tensor-product
615 const Real xi = p(0);
616 const Real eta = p(1);
617 const Real zeta = p(2);
618
619 // The only way to make any sense of this
620 // is to look at the mgflo/mg2/mgf documentation
621 // and make the cut-out cube!
622 // 0 1 2 3 4 5 6 7
623 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
624 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
625 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
626
627 return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
628 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)*
629 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta));
630 }
631
632
633 default:
634 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
635 }
636 }
637
638
639
640
641 case SECOND:
642 {
643 switch (type)
644 {
645
646 // Bernstein shape functions on the tetrahedron.
647 case TET10:
648 case TET14:
649 {
650 libmesh_assert_less (i, 10);
651
652 // Area coordinates
653 const Real zeta1 = p(0);
654 const Real zeta2 = p(1);
655 const Real zeta3 = p(2);
656 const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
657
658 switch(i)
659 {
660 case 0: return zeta0*zeta0;
661 case 1: return zeta1*zeta1;
662 case 2: return zeta2*zeta2;
663 case 3: return zeta3*zeta3;
664 case 4: return 2.*zeta0*zeta1;
665 case 5: return 2.*zeta1*zeta2;
666 case 6: return 2.*zeta0*zeta2;
667 case 7: return 2.*zeta3*zeta0;
668 case 8: return 2.*zeta1*zeta3;
669 case 9: return 2.*zeta2*zeta3;
670
671 default:
672 libmesh_error_msg("Invalid shape function index i = " << i);
673 }
674 }
675
676 // Bernstein shape functions on the 20-noded hexahedral.
677 case HEX20:
678 {
679 libmesh_assert_less (i, 20);
680
681 // Compute hex shape functions as a tensor-product
682 const Real xi = p(0);
683 const Real eta = p(1);
684 const Real zeta = p(2);
685
686 return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[i], xi)*
687 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[i], eta)*
688 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[i], zeta)
689 +hex20_scal20[i]*
690 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[20], xi)*
691 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[20], eta)*
692 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[20], zeta)
693 +hex20_scal21[i]*
694 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[21], xi)*
695 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[21], eta)*
696 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[21], zeta)
697 +hex20_scal22[i]*
698 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[22], xi)*
699 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[22], eta)*
700 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[22], zeta)
701 +hex20_scal23[i]*
702 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[23], xi)*
703 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[23], eta)*
704 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[23], zeta)
705 +hex20_scal24[i]*
706 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[24], xi)*
707 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[24], eta)*
708 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[24], zeta)
709 +hex20_scal25[i]*
710 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[25], xi)*
711 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[25], eta)*
712 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[25], zeta)
713 +hex20_scal26[i]*
714 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[26], xi)*
715 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[26], eta)*
716 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[26], zeta));
717 }
718
719 // Bernstein shape functions on the hexahedral.
720 case HEX27:
721 {
722 libmesh_assert_less (i, 27);
723
724 // Compute hex shape functions as a tensor-product
725 const Real xi = p(0);
726 const Real eta = p(1);
727 const Real zeta = p(2);
728
729 // The only way to make any sense of this
730 // is to look at the mgflo/mg2/mgf documentation
731 // and make the cut-out cube!
732 // 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
733 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
734 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
735 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2};
736
737 return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
738 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)*
739 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta));
740 }
741
742
743 default:
744 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
745 }
746
747 }
748
749
750
751 // 3rd-order Bernstein.
752 case THIRD:
753 {
754 switch (type)
755 {
756
757 // // Bernstein shape functions on the tetrahedron.
758 // case TET10:
759 // {
760 // libmesh_assert_less (i, 20);
761
762 // // Area coordinates
763 // const Real zeta1 = p(0);
764 // const Real zeta2 = p(1);
765 // const Real zeta3 = p(2);
766 // const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
767
768
769 // unsigned int shape=i;
770
771 // // handle the edge orientation
772
773 // if ((i== 4||i== 5) && elem->node_id(0) > elem->node_id(1))shape= 9-i; //Edge 0
774 // if ((i== 6||i== 7) && elem->node_id(1) > elem->node_id(2))shape=13-i; //Edge 1
775 // if ((i== 8||i== 9) && elem->node_id(0) > elem->node_id(2))shape=17-i; //Edge 2
776 // if ((i==10||i==11) && elem->node_id(0) > elem->node_id(3))shape=21-i; //Edge 3
777 // if ((i==12||i==13) && elem->node_id(1) > elem->node_id(3))shape=25-i; //Edge 4
778 // if ((i==14||i==15) && elem->node_id(2) > elem->node_id(3))shape=29-i; //Edge 5
779
780 // // No need to handle face orientation in 3rd order.
781
782
783 // switch(shape)
784 // {
785 // //point function
786 // case 0: return zeta0*zeta0*zeta0;
787 // case 1: return zeta1*zeta1*zeta1;
788 // case 2: return zeta2*zeta2*zeta2;
789 // case 3: return zeta3*zeta3*zeta3;
790
791 // //edge functions
792 // case 4: return 3.*zeta0*zeta0*zeta1;
793 // case 5: return 3.*zeta1*zeta1*zeta0;
794
795 // case 6: return 3.*zeta1*zeta1*zeta2;
796 // case 7: return 3.*zeta2*zeta2*zeta1;
797
798 // case 8: return 3.*zeta0*zeta0*zeta2;
799 // case 9: return 3.*zeta2*zeta2*zeta0;
800
801 // case 10: return 3.*zeta0*zeta0*zeta3;
802 // case 11: return 3.*zeta3*zeta3*zeta0;
803
804 // case 12: return 3.*zeta1*zeta1*zeta3;
805 // case 13: return 3.*zeta3*zeta3*zeta1;
806
807 // case 14: return 3.*zeta2*zeta2*zeta3;
808 // case 15: return 3.*zeta3*zeta3*zeta2;
809
810 // //face functions
811 // case 16: return 6.*zeta0*zeta1*zeta2;
812 // case 17: return 6.*zeta0*zeta1*zeta3;
813 // case 18: return 6.*zeta1*zeta2*zeta3;
814 // case 19: return 6.*zeta2*zeta0*zeta3;
815
816 // default:
817 // libmesh_error_msg("Invalid shape function index i = " << i);
818 // }
819 // }
820
821
822 // Bernstein shape functions on the hexahedral.
823 case HEX27:
824 {
825 libmesh_assert_less (i, 64);
826
827 // The only way to make any sense of this
828 // is to look at the mgflo/mg2/mgf documentation
829 // and make the cut-out cube!
830 // 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
831 // 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
832 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};
833 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};
834 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};
835
836 const Point p_mapped = hex_remap(p, i0, i1, i2);
837 return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], p_mapped(0))*
838 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], p_mapped(1))*
839 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], p_mapped(2)));
840 }
841
842 default:
843 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
844 }
845
846 }//case THIRD
847
848
849 // 4th-order Bernstein.
850 case FOURTH:
851 {
852 switch (type)
853 {
854
855 // Bernstein shape functions on the hexahedral.
856 case HEX27:
857 {
858 libmesh_assert_less (i, 125);
859
860 // The only way to make any sense of this
861 // is to look at the mgflo/mg2/mgf documentation
862 // and make the cut-out cube!
863 // 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
864 // 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 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 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
865 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 2, 3, 4, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4};
866 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4};
867 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4};
868
869 const Point p_mapped = hex_remap(p, i0, i1, i2);
870 return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], p_mapped(0))*
871 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], p_mapped(1))*
872 FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], p_mapped(2)));
873 }
874
875
876 default:
877 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
878 }
879 }
880
881
882 default:
883 libmesh_error_msg("Invalid totalorder = " << totalorder);
884 }
885#else // LIBMESH_DIM != 3
886 libmesh_ignore(elem, order, i, p, add_p_level);
887 libmesh_not_implemented();
888#endif
889}
890
891
892
893template <>
894Real FE<3,BERNSTEIN>::shape(const ElemType,
895 const Order,
896 const unsigned int,
897 const Point &)
898{
899 libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge and face orientation is needed.");
900 return 0.;
901}
902
903
904template <>
905Real FE<3,BERNSTEIN>::shape(const FEType fet,
906 const Elem * elem,
907 const unsigned int i,
908 const Point & p,
909 const bool add_p_level)
910{
911 return FE<3,BERNSTEIN>::shape(elem, fet.order, i, p, add_p_level);
912}
913
914template <>
915Real FE<3,BERNSTEIN>::shape_deriv(const Elem * elem,
916 const Order order,
917 const unsigned int i,
918 const unsigned int j,
919 const Point & p,
920 const bool add_p_level)
921{
922
923#if LIBMESH_DIM == 3
924 libmesh_assert(elem);
925 const ElemType type = elem->type();
926
927 const Order totalorder =
928 order + add_p_level*elem->p_level();
929
930 libmesh_assert_less (j, 3);
931
932 switch (totalorder)
933 {
934 // 1st order Bernstein.
935 case FIRST:
936 {
937 switch (type)
938 {
939 // Bernstein shape functions on the tetrahedron.
940 case TET4:
941 case TET10:
942 case TET14:
943 {
944 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,BERNSTEIN>::shape);
945 }
946
947
948 // Bernstein shape functions on the hexahedral.
949 case HEX8:
950 case HEX20:
951 case HEX27:
952 {
953 libmesh_assert_less (i, 8);
954
955 // Compute hex shape functions as a tensor-product
956 const Real xi = p(0);
957 const Real eta = p(1);
958 const Real zeta = p(2);
959
960 // The only way to make any sense of this
961 // is to look at the mgflo/mg2/mgf documentation
962 // and make the cut-out cube!
963 // 0 1 2 3 4 5 6 7
964 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
965 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
966 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
967
968 switch (j)
969 {
970 // d()/dxi
971 case 0:
972 return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
973 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
974 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
975
976 // d()/deta
977 case 1:
978 return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
979 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)*
980 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
981
982 // d()/dzeta
983 case 2:
984 return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
985 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
986 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta));
987
988 default:
989 libmesh_error_msg("Invalid derivative index j = " << j);
990 }
991 }
992
993 default:
994 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
995 }
996 }
997
998
999
1000
1001 case SECOND:
1002 {
1003 switch (type)
1004 {
1005 // Bernstein shape functions on the tetrahedron.
1006 case TET10:
1007 case TET14:
1008 {
1009 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,BERNSTEIN>::shape);
1010 }
1011
1012 // Bernstein shape functions on the hexahedral.
1013 case HEX20:
1014 {
1015 libmesh_assert_less (i, 20);
1016
1017 // Compute hex shape functions as a tensor-product
1018 const Real xi = p(0);
1019 const Real eta = p(1);
1020 const Real zeta = p(2);
1021
1022 switch (j)
1023 {
1024 // d()/dxi
1025 case 0:
1026 return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[i], 0, xi)*
1027 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[i], eta)*
1028 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[i], zeta)
1029 +hex20_scal20[i]*
1030 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[20], 0, xi)*
1031 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[20], eta)*
1032 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[20], zeta)
1033 +hex20_scal21[i]*
1034 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[21], 0, xi)*
1035 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[21], eta)*
1036 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[21], zeta)
1037 +hex20_scal22[i]*
1038 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[22], 0, xi)*
1039 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[22], eta)*
1040 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[22], zeta)
1041 +hex20_scal23[i]*
1042 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[23], 0, xi)*
1043 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[23], eta)*
1044 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[23], zeta)
1045 +hex20_scal24[i]*
1046 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[24], 0, xi)*
1047 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[24], eta)*
1048 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[24], zeta)
1049 +hex20_scal25[i]*
1050 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[25], 0, xi)*
1051 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[25], eta)*
1052 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[25], zeta)
1053 +hex20_scal26[i]*
1054 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[26], 0, xi)*
1055 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[26], eta)*
1056 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[26], zeta));
1057
1058 // d()/deta
1059 case 1:
1060 return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[i], xi)*
1061 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[i], 0, eta)*
1062 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[i], zeta)
1063 +hex20_scal20[i]*
1064 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[20], xi)*
1065 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[20], 0, eta)*
1066 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[20], zeta)
1067 +hex20_scal21[i]*
1068 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[21], xi)*
1069 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[21], 0, eta)*
1070 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[21], zeta)
1071 +hex20_scal22[i]*
1072 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[22], xi)*
1073 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[22], 0, eta)*
1074 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[22], zeta)
1075 +hex20_scal23[i]*
1076 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[23], xi)*
1077 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[23], 0, eta)*
1078 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[23], zeta)
1079 +hex20_scal24[i]*
1080 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[24], xi)*
1081 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[24], 0, eta)*
1082 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[24], zeta)
1083 +hex20_scal25[i]*
1084 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[25], xi)*
1085 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[25], 0, eta)*
1086 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[25], zeta)
1087 +hex20_scal26[i]*
1088 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[26], xi)*
1089 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[26], 0, eta)*
1090 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[26], zeta));
1091
1092 // d()/dzeta
1093 case 2:
1094 return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[i], xi)*
1095 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[i], eta)*
1096 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[i], 0, zeta)
1097 +hex20_scal20[i]*
1098 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[20], xi)*
1099 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[20], eta)*
1100 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[20], 0, zeta)
1101 +hex20_scal21[i]*
1102 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[21], xi)*
1103 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[21], eta)*
1104 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[21], 0, zeta)
1105 +hex20_scal22[i]*
1106 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[22], xi)*
1107 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[22], eta)*
1108 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[22], 0, zeta)
1109 +hex20_scal23[i]*
1110 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[23], xi)*
1111 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[23], eta)*
1112 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[23], 0, zeta)
1113 +hex20_scal24[i]*
1114 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[24], xi)*
1115 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[24], eta)*
1116 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[24], 0, zeta)
1117 +hex20_scal25[i]*
1118 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[25], xi)*
1119 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[25], eta)*
1120 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[25], 0, zeta)
1121 +hex20_scal26[i]*
1122 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[26], xi)*
1123 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[26], eta)*
1124 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[26], 0, zeta));
1125
1126 default:
1127 libmesh_error_msg("Invalid derivative index j = " << j);
1128 }
1129 }
1130
1131 // Bernstein shape functions on the hexahedral.
1132 case HEX27:
1133 {
1134 libmesh_assert_less (i, 27);
1135
1136 // Compute hex shape functions as a tensor-product
1137 const Real xi = p(0);
1138 const Real eta = p(1);
1139 const Real zeta = p(2);
1140
1141 // The only way to make any sense of this
1142 // is to look at the mgflo/mg2/mgf documentation
1143 // and make the cut-out cube!
1144 // 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
1145 static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
1146 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
1147 static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2};
1148
1149 switch (j)
1150 {
1151 // d()/dxi
1152 case 0:
1153 return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1154 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1155 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
1156
1157 // d()/deta
1158 case 1:
1159 return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1160 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)*
1161 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
1162
1163 // d()/dzeta
1164 case 2:
1165 return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1166 FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1167 FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta));
1168
1169 default:
1170 libmesh_error_msg("Invalid derivative index j = " << j);
1171 }
1172 }
1173
1174
1175 default:
1176 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1177 }
1178 }
1179
1180
1181
1182 // 3rd-order Bernstein.
1183 case THIRD:
1184 {
1185 switch (type)
1186 {
1187
1188 // // Bernstein shape functions derivatives.
1189 // case TET10:
1190 // {
1191 // // I have been lazy here and am using finite differences
1192 // // to compute the derivatives!
1193 // const Real eps = 1.e-4;
1194
1195 // libmesh_assert_less (i, 20);
1196 // libmesh_assert_less (j, 3);
1197
1198 // switch (j)
1199 // {
1200 // // d()/dxi
1201 // case 0:
1202 // {
1203 // const Point pp(p(0)+eps, p(1), p(2));
1204 // const Point pm(p(0)-eps, p(1), p(2));
1205
1206 // return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1207 // FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1208 // }
1209
1210 // // d()/deta
1211 // case 1:
1212 // {
1213 // const Point pp(p(0), p(1)+eps, p(2));
1214 // const Point pm(p(0), p(1)-eps, p(2));
1215
1216 // return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1217 // FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1218 // }
1219 // // d()/dzeta
1220 // case 2:
1221 // {
1222 // const Point pp(p(0), p(1), p(2)+eps);
1223 // const Point pm(p(0), p(1), p(2)-eps);
1224
1225 // return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1226 // FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1227 // }
1228 // default:
1229 // libmesh_error_msg("Invalid derivative index j = " << j);
1230 // }
1231
1232
1233 // }
1234
1235
1236 // Bernstein shape functions on the hexahedral.
1237 case HEX27:
1238 {
1239 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,BERNSTEIN>::shape);
1240 }
1241
1242 // // Compute hex shape functions as a tensor-product
1243 // const Real xi = p(0);
1244 // const Real eta = p(1);
1245 // const Real zeta = p(2);
1246 // Real xi_mapped = p(0);
1247 // Real eta_mapped = p(1);
1248 // Real zeta_mapped = p(2);
1249
1250 // // The only way to make any sense of this
1251 // // is to look at the mgflo/mg2/mgf documentation
1252 // // and make the cut-out cube!
1253 // // 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
1254 // // 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
1255 // 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};
1256 // 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};
1257 // 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};
1258
1259
1260
1261 // // handle the edge orientation
1262 // {
1263 // // Edge 0
1264 // if ((i1[i] == 0) && (i2[i] == 0))
1265 // {
1266 // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(1)))
1267 // xi_mapped = -xi;
1268 // }
1269 // // Edge 1
1270 // else if ((i0[i] == 1) && (i2[i] == 0))
1271 // {
1272 // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(2)))
1273 // eta_mapped = -eta;
1274 // }
1275 // // Edge 2
1276 // else if ((i1[i] == 1) && (i2[i] == 0))
1277 // {
1278 // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(2)))
1279 // xi_mapped = -xi;
1280 // }
1281 // // Edge 3
1282 // else if ((i0[i] == 0) && (i2[i] == 0))
1283 // {
1284 // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(3)))
1285 // eta_mapped = -eta;
1286 // }
1287 // // Edge 4
1288 // else if ((i0[i] == 0) && (i1[i] == 0))
1289 // {
1290 // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(4)))
1291 // zeta_mapped = -zeta;
1292 // }
1293 // // Edge 5
1294 // else if ((i0[i] == 1) && (i1[i] == 0))
1295 // {
1296 // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(5)))
1297 // zeta_mapped = -zeta;
1298 // }
1299 // // Edge 6
1300 // else if ((i0[i] == 1) && (i1[i] == 1))
1301 // {
1302 // if (elem->node_id(2) != std::min(elem->node_id(2), elem->node_id(6)))
1303 // zeta_mapped = -zeta;
1304 // }
1305 // // Edge 7
1306 // else if ((i0[i] == 0) && (i1[i] == 1))
1307 // {
1308 // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(7)))
1309 // zeta_mapped = -zeta;
1310 // }
1311 // // Edge 8
1312 // else if ((i1[i] == 0) && (i2[i] == 1))
1313 // {
1314 // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(5)))
1315 // xi_mapped = -xi;
1316 // }
1317 // // Edge 9
1318 // else if ((i0[i] == 1) && (i2[i] == 1))
1319 // {
1320 // if (elem->node_id(5) != std::min(elem->node_id(5), elem->node_id(6)))
1321 // eta_mapped = -eta;
1322 // }
1323 // // Edge 10
1324 // else if ((i1[i] == 1) && (i2[i] == 1))
1325 // {
1326 // if (elem->node_id(7) != std::min(elem->node_id(7), elem->node_id(6)))
1327 // xi_mapped = -xi;
1328 // }
1329 // // Edge 11
1330 // else if ((i0[i] == 0) && (i2[i] == 1))
1331 // {
1332 // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(7)))
1333 // eta_mapped = -eta;
1334 // }
1335 // }
1336
1337
1338 // // handle the face orientation
1339 // {
1340 // // Face 0
1341 // if ((i2[i] == 0) && (i0[i] >= 2) && (i1[i] >= 2))
1342 // {
1343 // const unsigned int min_node = std::min(elem->node_id(1),
1344 // std::min(elem->node_id(2),
1345 // std::min(elem->node_id(0),
1346 // elem->node_id(3))));
1347 // if (elem->node_id(0) == min_node)
1348 // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(3)))
1349 // {
1350 // // Case 1
1351 // xi_mapped = xi;
1352 // eta_mapped = eta;
1353 // }
1354 // else
1355 // {
1356 // // Case 2
1357 // xi_mapped = eta;
1358 // eta_mapped = xi;
1359 // }
1360
1361 // else if (elem->node_id(3) == min_node)
1362 // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(2)))
1363 // {
1364 // // Case 3
1365 // xi_mapped = -eta;
1366 // eta_mapped = xi;
1367 // }
1368 // else
1369 // {
1370 // // Case 4
1371 // xi_mapped = xi;
1372 // eta_mapped = -eta;
1373 // }
1374
1375 // else if (elem->node_id(2) == min_node)
1376 // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(1)))
1377 // {
1378 // // Case 5
1379 // xi_mapped = -xi;
1380 // eta_mapped = -eta;
1381 // }
1382 // else
1383 // {
1384 // // Case 6
1385 // xi_mapped = -eta;
1386 // eta_mapped = -xi;
1387 // }
1388
1389 // else if (elem->node_id(1) == min_node)
1390 // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(0)))
1391 // {
1392 // // Case 7
1393 // xi_mapped = eta;
1394 // eta_mapped = -xi;
1395 // }
1396 // else
1397 // {
1398 // // Case 8
1399 // xi_mapped = -xi;
1400 // eta_mapped = eta;
1401 // }
1402 // }
1403
1404
1405 // // Face 1
1406 // else if ((i1[i] == 0) && (i0[i] >= 2) && (i2[i] >= 2))
1407 // {
1408 // const unsigned int min_node = std::min(elem->node_id(0),
1409 // std::min(elem->node_id(1),
1410 // std::min(elem->node_id(5),
1411 // elem->node_id(4))));
1412 // if (elem->node_id(0) == min_node)
1413 // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(4)))
1414 // {
1415 // // Case 1
1416 // xi_mapped = xi;
1417 // zeta_mapped = zeta;
1418 // }
1419 // else
1420 // {
1421 // // Case 2
1422 // xi_mapped = zeta;
1423 // zeta_mapped = xi;
1424 // }
1425
1426 // else if (elem->node_id(1) == min_node)
1427 // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(0)))
1428 // {
1429 // // Case 3
1430 // xi_mapped = zeta;
1431 // zeta_mapped = -xi;
1432 // }
1433 // else
1434 // {
1435 // // Case 4
1436 // xi_mapped = -xi;
1437 // zeta_mapped = zeta;
1438 // }
1439
1440 // else if (elem->node_id(5) == min_node)
1441 // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(1)))
1442 // {
1443 // // Case 5
1444 // xi_mapped = -xi;
1445 // zeta_mapped = -zeta;
1446 // }
1447 // else
1448 // {
1449 // // Case 6
1450 // xi_mapped = -zeta;
1451 // zeta_mapped = -xi;
1452 // }
1453
1454 // else if (elem->node_id(4) == min_node)
1455 // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(5)))
1456 // {
1457 // // Case 7
1458 // xi_mapped = -xi;
1459 // zeta_mapped = zeta;
1460 // }
1461 // else
1462 // {
1463 // // Case 8
1464 // xi_mapped = xi;
1465 // zeta_mapped = -zeta;
1466 // }
1467 // }
1468
1469
1470 // // Face 2
1471 // else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] >= 2))
1472 // {
1473 // const unsigned int min_node = std::min(elem->node_id(1),
1474 // std::min(elem->node_id(2),
1475 // std::min(elem->node_id(6),
1476 // elem->node_id(5))));
1477 // if (elem->node_id(1) == min_node)
1478 // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(5)))
1479 // {
1480 // // Case 1
1481 // eta_mapped = eta;
1482 // zeta_mapped = zeta;
1483 // }
1484 // else
1485 // {
1486 // // Case 2
1487 // eta_mapped = zeta;
1488 // zeta_mapped = eta;
1489 // }
1490
1491 // else if (elem->node_id(2) == min_node)
1492 // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(1)))
1493 // {
1494 // // Case 3
1495 // eta_mapped = zeta;
1496 // zeta_mapped = -eta;
1497 // }
1498 // else
1499 // {
1500 // // Case 4
1501 // eta_mapped = -eta;
1502 // zeta_mapped = zeta;
1503 // }
1504
1505 // else if (elem->node_id(6) == min_node)
1506 // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(2)))
1507 // {
1508 // // Case 5
1509 // eta_mapped = -eta;
1510 // zeta_mapped = -zeta;
1511 // }
1512 // else
1513 // {
1514 // // Case 6
1515 // eta_mapped = -zeta;
1516 // zeta_mapped = -eta;
1517 // }
1518
1519 // else if (elem->node_id(5) == min_node)
1520 // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(6)))
1521 // {
1522 // // Case 7
1523 // eta_mapped = -zeta;
1524 // zeta_mapped = eta;
1525 // }
1526 // else
1527 // {
1528 // // Case 8
1529 // eta_mapped = eta;
1530 // zeta_mapped = -zeta;
1531 // }
1532 // }
1533
1534
1535 // // Face 3
1536 // else if ((i1[i] == 1) && (i0[i] >= 2) && (i2[i] >= 2))
1537 // {
1538 // const unsigned int min_node = std::min(elem->node_id(2),
1539 // std::min(elem->node_id(3),
1540 // std::min(elem->node_id(7),
1541 // elem->node_id(6))));
1542 // if (elem->node_id(3) == min_node)
1543 // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(7)))
1544 // {
1545 // // Case 1
1546 // xi_mapped = xi;
1547 // zeta_mapped = zeta;
1548 // }
1549 // else
1550 // {
1551 // // Case 2
1552 // xi_mapped = zeta;
1553 // zeta_mapped = xi;
1554 // }
1555
1556 // else if (elem->node_id(7) == min_node)
1557 // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(6)))
1558 // {
1559 // // Case 3
1560 // xi_mapped = -zeta;
1561 // zeta_mapped = xi;
1562 // }
1563 // else
1564 // {
1565 // // Case 4
1566 // xi_mapped = xi;
1567 // zeta_mapped = -zeta;
1568 // }
1569
1570 // else if (elem->node_id(6) == min_node)
1571 // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(2)))
1572 // {
1573 // // Case 5
1574 // xi_mapped = -xi;
1575 // zeta_mapped = -zeta;
1576 // }
1577 // else
1578 // {
1579 // // Case 6
1580 // xi_mapped = -zeta;
1581 // zeta_mapped = -xi;
1582 // }
1583
1584 // else if (elem->node_id(2) == min_node)
1585 // if (elem->node_id(6) == std::min(elem->node_id(3), elem->node_id(6)))
1586 // {
1587 // // Case 7
1588 // xi_mapped = zeta;
1589 // zeta_mapped = -xi;
1590 // }
1591 // else
1592 // {
1593 // // Case 8
1594 // xi_mapped = -xi;
1595 // zeta_mapped = zeta;
1596 // }
1597 // }
1598
1599
1600 // // Face 4
1601 // else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] >= 2))
1602 // {
1603 // const unsigned int min_node = std::min(elem->node_id(3),
1604 // std::min(elem->node_id(0),
1605 // std::min(elem->node_id(4),
1606 // elem->node_id(7))));
1607 // if (elem->node_id(0) == min_node)
1608 // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(4)))
1609 // {
1610 // // Case 1
1611 // eta_mapped = eta;
1612 // zeta_mapped = zeta;
1613 // }
1614 // else
1615 // {
1616 // // Case 2
1617 // eta_mapped = zeta;
1618 // zeta_mapped = eta;
1619 // }
1620
1621 // else if (elem->node_id(4) == min_node)
1622 // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(7)))
1623 // {
1624 // // Case 3
1625 // eta_mapped = -zeta;
1626 // zeta_mapped = eta;
1627 // }
1628 // else
1629 // {
1630 // // Case 4
1631 // eta_mapped = eta;
1632 // zeta_mapped = -zeta;
1633 // }
1634
1635 // else if (elem->node_id(7) == min_node)
1636 // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(3)))
1637 // {
1638 // // Case 5
1639 // eta_mapped = -eta;
1640 // zeta_mapped = -zeta;
1641 // }
1642 // else
1643 // {
1644 // // Case 6
1645 // eta_mapped = -zeta;
1646 // zeta_mapped = -eta;
1647 // }
1648
1649 // else if (elem->node_id(3) == min_node)
1650 // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(0)))
1651 // {
1652 // // Case 7
1653 // eta_mapped = zeta;
1654 // zeta_mapped = -eta;
1655 // }
1656 // else
1657 // {
1658 // // Case 8
1659 // eta_mapped = -eta;
1660 // zeta_mapped = zeta;
1661 // }
1662 // }
1663
1664
1665 // // Face 5
1666 // else if ((i2[i] == 1) && (i0[i] >= 2) && (i1[i] >= 2))
1667 // {
1668 // const unsigned int min_node = std::min(elem->node_id(4),
1669 // std::min(elem->node_id(5),
1670 // std::min(elem->node_id(6),
1671 // elem->node_id(7))));
1672 // if (elem->node_id(4) == min_node)
1673 // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(7)))
1674 // {
1675 // // Case 1
1676 // xi_mapped = xi;
1677 // eta_mapped = eta;
1678 // }
1679 // else
1680 // {
1681 // // Case 2
1682 // xi_mapped = eta;
1683 // eta_mapped = xi;
1684 // }
1685
1686 // else if (elem->node_id(5) == min_node)
1687 // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(4)))
1688 // {
1689 // // Case 3
1690 // xi_mapped = eta;
1691 // eta_mapped = -xi;
1692 // }
1693 // else
1694 // {
1695 // // Case 4
1696 // xi_mapped = -xi;
1697 // eta_mapped = eta;
1698 // }
1699
1700 // else if (elem->node_id(6) == min_node)
1701 // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(5)))
1702 // {
1703 // // Case 5
1704 // xi_mapped = -xi;
1705 // eta_mapped = -eta;
1706 // }
1707 // else
1708 // {
1709 // // Case 6
1710 // xi_mapped = -eta;
1711 // eta_mapped = -xi;
1712 // }
1713
1714 // else if (elem->node_id(7) == min_node)
1715 // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(6)))
1716 // {
1717 // // Case 7
1718 // xi_mapped = -eta;
1719 // eta_mapped = xi;
1720 // }
1721 // else
1722 // {
1723 // // Case 8
1724 // xi_mapped = xi;
1725 // eta_mapped = eta;
1726 // }
1727 // }
1728
1729
1730 // }
1731
1732
1733
1734 // libmesh_assert_less (j, 3);
1735
1736 // switch (j)
1737 // {
1738 // // d()/dxi
1739 // case 0:
1740 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi_mapped)*
1741 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
1742 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
1743
1744 // // d()/deta
1745 // case 1:
1746 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
1747 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta_mapped)*
1748 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
1749
1750 // // d()/dzeta
1751 // case 2:
1752 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
1753 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
1754 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta_mapped));
1755
1756 // default:
1757 // libmesh_error_msg("Invalid derivative index j = " << j);
1758 // }
1759
1760
1761 default:
1762 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1763 }
1764 }
1765
1766 // 4th-order Bernstein.
1767 case FOURTH:
1768 {
1769 switch (type)
1770 {
1771
1772 // Bernstein shape functions derivatives on the hexahedral.
1773 case HEX27:
1774 {
1775 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,BERNSTEIN>::shape);
1776 }
1777
1778 // // Compute hex shape functions as a tensor-product
1779 // const Real xi = p(0);
1780 // const Real eta = p(1);
1781 // const Real zeta = p(2);
1782 // Real xi_mapped = p(0);
1783 // Real eta_mapped = p(1);
1784 // Real zeta_mapped = p(2);
1785
1786 // // The only way to make any sense of this
1787 // // is to look at the mgflo/mg2/mgf documentation
1788 // // and make the cut-out cube!
1789 // // 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
1790 // // 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 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
1791 // static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 2, 3, 4, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4};
1792 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4};
1793 // static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4};
1794
1795
1796
1797 // // handle the edge orientation
1798 // {
1799 // // Edge 0
1800 // if ((i1[i] == 0) && (i2[i] == 0))
1801 // {
1802 // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(1)))
1803 // xi_mapped = -xi;
1804 // }
1805 // // Edge 1
1806 // else if ((i0[i] == 1) && (i2[i] == 0))
1807 // {
1808 // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(2)))
1809 // eta_mapped = -eta;
1810 // }
1811 // // Edge 2
1812 // else if ((i1[i] == 1) && (i2[i] == 0))
1813 // {
1814 // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(2)))
1815 // xi_mapped = -xi;
1816 // }
1817 // // Edge 3
1818 // else if ((i0[i] == 0) && (i2[i] == 0))
1819 // {
1820 // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(3)))
1821 // eta_mapped = -eta;
1822 // }
1823 // // Edge 4
1824 // else if ((i0[i] == 0) && (i1[i] == 0))
1825 // {
1826 // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(4)))
1827 // zeta_mapped = -zeta;
1828 // }
1829 // // Edge 5
1830 // else if ((i0[i] == 1) && (i1[i] == 0))
1831 // {
1832 // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(5)))
1833 // zeta_mapped = -zeta;
1834 // }
1835 // // Edge 6
1836 // else if ((i0[i] == 1) && (i1[i] == 1))
1837 // {
1838 // if (elem->node_id(2) != std::min(elem->node_id(2), elem->node_id(6)))
1839 // zeta_mapped = -zeta;
1840 // }
1841 // // Edge 7
1842 // else if ((i0[i] == 0) && (i1[i] == 1))
1843 // {
1844 // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(7)))
1845 // zeta_mapped = -zeta;
1846 // }
1847 // // Edge 8
1848 // else if ((i1[i] == 0) && (i2[i] == 1))
1849 // {
1850 // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(5)))
1851 // xi_mapped = -xi;
1852 // }
1853 // // Edge 9
1854 // else if ((i0[i] == 1) && (i2[i] == 1))
1855 // {
1856 // if (elem->node_id(5) != std::min(elem->node_id(5), elem->node_id(6)))
1857 // eta_mapped = -eta;
1858 // }
1859 // // Edge 10
1860 // else if ((i1[i] == 1) && (i2[i] == 1))
1861 // {
1862 // if (elem->node_id(7) != std::min(elem->node_id(7), elem->node_id(6)))
1863 // xi_mapped = -xi;
1864 // }
1865 // // Edge 11
1866 // else if ((i0[i] == 0) && (i2[i] == 1))
1867 // {
1868 // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(7)))
1869 // eta_mapped = -eta;
1870 // }
1871 // }
1872
1873
1874 // // handle the face orientation
1875 // {
1876 // // Face 0
1877 // if ((i2[i] == 0) && (i0[i] >= 2) && (i1[i] >= 2))
1878 // {
1879 // const unsigned int min_node = std::min(elem->node_id(1),
1880 // std::min(elem->node_id(2),
1881 // std::min(elem->node_id(0),
1882 // elem->node_id(3))));
1883 // if (elem->node_id(0) == min_node)
1884 // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(3)))
1885 // {
1886 // // Case 1
1887 // xi_mapped = xi;
1888 // eta_mapped = eta;
1889 // }
1890 // else
1891 // {
1892 // // Case 2
1893 // xi_mapped = eta;
1894 // eta_mapped = xi;
1895 // }
1896
1897 // else if (elem->node_id(3) == min_node)
1898 // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(2)))
1899 // {
1900 // // Case 3
1901 // xi_mapped = -eta;
1902 // eta_mapped = xi;
1903 // }
1904 // else
1905 // {
1906 // // Case 4
1907 // xi_mapped = xi;
1908 // eta_mapped = -eta;
1909 // }
1910
1911 // else if (elem->node_id(2) == min_node)
1912 // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(1)))
1913 // {
1914 // // Case 5
1915 // xi_mapped = -xi;
1916 // eta_mapped = -eta;
1917 // }
1918 // else
1919 // {
1920 // // Case 6
1921 // xi_mapped = -eta;
1922 // eta_mapped = -xi;
1923 // }
1924
1925 // else if (elem->node_id(1) == min_node)
1926 // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(0)))
1927 // {
1928 // // Case 7
1929 // xi_mapped = eta;
1930 // eta_mapped = -xi;
1931 // }
1932 // else
1933 // {
1934 // // Case 8
1935 // xi_mapped = -xi;
1936 // eta_mapped = eta;
1937 // }
1938 // }
1939
1940
1941 // // Face 1
1942 // else if ((i1[i] == 0) && (i0[i] >= 2) && (i2[i] >= 2))
1943 // {
1944 // const unsigned int min_node = std::min(elem->node_id(0),
1945 // std::min(elem->node_id(1),
1946 // std::min(elem->node_id(5),
1947 // elem->node_id(4))));
1948 // if (elem->node_id(0) == min_node)
1949 // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(4)))
1950 // {
1951 // // Case 1
1952 // xi_mapped = xi;
1953 // zeta_mapped = zeta;
1954 // }
1955 // else
1956 // {
1957 // // Case 2
1958 // xi_mapped = zeta;
1959 // zeta_mapped = xi;
1960 // }
1961
1962 // else if (elem->node_id(1) == min_node)
1963 // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(0)))
1964 // {
1965 // // Case 3
1966 // xi_mapped = zeta;
1967 // zeta_mapped = -xi;
1968 // }
1969 // else
1970 // {
1971 // // Case 4
1972 // xi_mapped = -xi;
1973 // zeta_mapped = zeta;
1974 // }
1975
1976 // else if (elem->node_id(5) == min_node)
1977 // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(1)))
1978 // {
1979 // // Case 5
1980 // xi_mapped = -xi;
1981 // zeta_mapped = -zeta;
1982 // }
1983 // else
1984 // {
1985 // // Case 6
1986 // xi_mapped = -zeta;
1987 // zeta_mapped = -xi;
1988 // }
1989
1990 // else if (elem->node_id(4) == min_node)
1991 // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(5)))
1992 // {
1993 // // Case 7
1994 // xi_mapped = -xi;
1995 // zeta_mapped = zeta;
1996 // }
1997 // else
1998 // {
1999 // // Case 8
2000 // xi_mapped = xi;
2001 // zeta_mapped = -zeta;
2002 // }
2003 // }
2004
2005
2006 // // Face 2
2007 // else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] >= 2))
2008 // {
2009 // const unsigned int min_node = std::min(elem->node_id(1),
2010 // std::min(elem->node_id(2),
2011 // std::min(elem->node_id(6),
2012 // elem->node_id(5))));
2013 // if (elem->node_id(1) == min_node)
2014 // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(5)))
2015 // {
2016 // // Case 1
2017 // eta_mapped = eta;
2018 // zeta_mapped = zeta;
2019 // }
2020 // else
2021 // {
2022 // // Case 2
2023 // eta_mapped = zeta;
2024 // zeta_mapped = eta;
2025 // }
2026
2027 // else if (elem->node_id(2) == min_node)
2028 // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(1)))
2029 // {
2030 // // Case 3
2031 // eta_mapped = zeta;
2032 // zeta_mapped = -eta;
2033 // }
2034 // else
2035 // {
2036 // // Case 4
2037 // eta_mapped = -eta;
2038 // zeta_mapped = zeta;
2039 // }
2040
2041 // else if (elem->node_id(6) == min_node)
2042 // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(2)))
2043 // {
2044 // // Case 5
2045 // eta_mapped = -eta;
2046 // zeta_mapped = -zeta;
2047 // }
2048 // else
2049 // {
2050 // // Case 6
2051 // eta_mapped = -zeta;
2052 // zeta_mapped = -eta;
2053 // }
2054
2055 // else if (elem->node_id(5) == min_node)
2056 // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(6)))
2057 // {
2058 // // Case 7
2059 // eta_mapped = -zeta;
2060 // zeta_mapped = eta;
2061 // }
2062 // else
2063 // {
2064 // // Case 8
2065 // eta_mapped = eta;
2066 // zeta_mapped = -zeta;
2067 // }
2068 // }
2069
2070
2071 // // Face 3
2072 // else if ((i1[i] == 1) && (i0[i] >= 2) && (i2[i] >= 2))
2073 // {
2074 // const unsigned int min_node = std::min(elem->node_id(2),
2075 // std::min(elem->node_id(3),
2076 // std::min(elem->node_id(7),
2077 // elem->node_id(6))));
2078 // if (elem->node_id(3) == min_node)
2079 // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(7)))
2080 // {
2081 // // Case 1
2082 // xi_mapped = xi;
2083 // zeta_mapped = zeta;
2084 // }
2085 // else
2086 // {
2087 // // Case 2
2088 // xi_mapped = zeta;
2089 // zeta_mapped = xi;
2090 // }
2091
2092 // else if (elem->node_id(7) == min_node)
2093 // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(6)))
2094 // {
2095 // // Case 3
2096 // xi_mapped = -zeta;
2097 // zeta_mapped = xi;
2098 // }
2099 // else
2100 // {
2101 // // Case 4
2102 // xi_mapped = xi;
2103 // zeta_mapped = -zeta;
2104 // }
2105
2106 // else if (elem->node_id(6) == min_node)
2107 // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(2)))
2108 // {
2109 // // Case 5
2110 // xi_mapped = -xi;
2111 // zeta_mapped = -zeta;
2112 // }
2113 // else
2114 // {
2115 // // Case 6
2116 // xi_mapped = -zeta;
2117 // zeta_mapped = -xi;
2118 // }
2119
2120 // else if (elem->node_id(2) == min_node)
2121 // if (elem->node_id(6) == std::min(elem->node_id(3), elem->node_id(6)))
2122 // {
2123 // // Case 7
2124 // xi_mapped = zeta;
2125 // zeta_mapped = -xi;
2126 // }
2127 // else
2128 // {
2129 // // Case 8
2130 // xi_mapped = -xi;
2131 // zeta_mapped = zeta;
2132 // }
2133 // }
2134
2135
2136 // // Face 4
2137 // else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] >= 2))
2138 // {
2139 // const unsigned int min_node = std::min(elem->node_id(3),
2140 // std::min(elem->node_id(0),
2141 // std::min(elem->node_id(4),
2142 // elem->node_id(7))));
2143 // if (elem->node_id(0) == min_node)
2144 // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(4)))
2145 // {
2146 // // Case 1
2147 // eta_mapped = eta;
2148 // zeta_mapped = zeta;
2149 // }
2150 // else
2151 // {
2152 // // Case 2
2153 // eta_mapped = zeta;
2154 // zeta_mapped = eta;
2155 // }
2156
2157 // else if (elem->node_id(4) == min_node)
2158 // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(7)))
2159 // {
2160 // // Case 3
2161 // eta_mapped = -zeta;
2162 // zeta_mapped = eta;
2163 // }
2164 // else
2165 // {
2166 // // Case 4
2167 // eta_mapped = eta;
2168 // zeta_mapped = -zeta;
2169 // }
2170
2171 // else if (elem->node_id(7) == min_node)
2172 // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(3)))
2173 // {
2174 // // Case 5
2175 // eta_mapped = -eta;
2176 // zeta_mapped = -zeta;
2177 // }
2178 // else
2179 // {
2180 // // Case 6
2181 // eta_mapped = -zeta;
2182 // zeta_mapped = -eta;
2183 // }
2184
2185 // else if (elem->node_id(3) == min_node)
2186 // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(0)))
2187 // {
2188 // // Case 7
2189 // eta_mapped = zeta;
2190 // zeta_mapped = -eta;
2191 // }
2192 // else
2193 // {
2194 // // Case 8
2195 // eta_mapped = -eta;
2196 // zeta_mapped = zeta;
2197 // }
2198 // }
2199
2200
2201 // // Face 5
2202 // else if ((i2[i] == 1) && (i0[i] >= 2) && (i1[i] >= 2))
2203 // {
2204 // const unsigned int min_node = std::min(elem->node_id(4),
2205 // std::min(elem->node_id(5),
2206 // std::min(elem->node_id(6),
2207 // elem->node_id(7))));
2208 // if (elem->node_id(4) == min_node)
2209 // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(7)))
2210 // {
2211 // // Case 1
2212 // xi_mapped = xi;
2213 // eta_mapped = eta;
2214 // }
2215 // else
2216 // {
2217 // // Case 2
2218 // xi_mapped = eta;
2219 // eta_mapped = xi;
2220 // }
2221
2222 // else if (elem->node_id(5) == min_node)
2223 // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(4)))
2224 // {
2225 // // Case 3
2226 // xi_mapped = eta;
2227 // eta_mapped = -xi;
2228 // }
2229 // else
2230 // {
2231 // // Case 4
2232 // xi_mapped = -xi;
2233 // eta_mapped = eta;
2234 // }
2235
2236 // else if (elem->node_id(6) == min_node)
2237 // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(5)))
2238 // {
2239 // // Case 5
2240 // xi_mapped = -xi;
2241 // eta_mapped = -eta;
2242 // }
2243 // else
2244 // {
2245 // // Case 6
2246 // xi_mapped = -eta;
2247 // eta_mapped = -xi;
2248 // }
2249
2250 // else if (elem->node_id(7) == min_node)
2251 // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(6)))
2252 // {
2253 // // Case 7
2254 // xi_mapped = -eta;
2255 // eta_mapped = xi;
2256 // }
2257 // else
2258 // {
2259 // // Case 8
2260 // xi_mapped = xi;
2261 // eta_mapped = eta;
2262 // }
2263 // }
2264
2265
2266 // }
2267
2268
2269
2270 // libmesh_assert_less (j, 3);
2271
2272 // switch (j)
2273 // {
2274 // // d()/dxi
2275 // case 0:
2276 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi_mapped)*
2277 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
2278 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
2279
2280 // // d()/deta
2281 // case 1:
2282 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
2283 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta_mapped)*
2284 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
2285
2286 // // d()/dzeta
2287 // case 2:
2288 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
2289 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
2290 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta_mapped));
2291
2292 // default:
2293 // libmesh_error_msg("Invalid derivative index j = " << j);
2294 // }
2295
2296
2297 default:
2298 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
2299 }
2300 }
2301
2302
2303 default:
2304 libmesh_error_msg("Invalid totalorder = " << totalorder);
2305 }
2306
2307#else // LIBMESH_DIM != 3
2308 libmesh_ignore(elem, order, i, j, p, add_p_level);
2309 libmesh_not_implemented();
2310#endif
2311}
2312
2313
2314template <>
2315Real FE<3,BERNSTEIN>::shape_deriv(const ElemType,
2316 const Order,
2317 const unsigned int,
2318 const unsigned int,
2319 const Point & )
2320{
2321 libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge and face orientation is needed.");
2322 return 0.;
2323}
2324
2325template <>
2326Real FE<3,BERNSTEIN>::shape_deriv(const FEType fet,
2327 const Elem * elem,
2328 const unsigned int i,
2329 const unsigned int j,
2330 const Point & p,
2331 const bool add_p_level)
2332{
2333 return FE<3,BERNSTEIN>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
2334}
2335
2336
2337
2338
2339#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2340
2341template <>
2342Real FE<3,BERNSTEIN>::shape_second_deriv(const Elem * elem,
2343 const Order order,
2344 const unsigned int i,
2345 const unsigned int j,
2346 const Point & p,
2347 const bool add_p_level)
2348{
2349 return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
2350 FE<3,BERNSTEIN>::shape_deriv);
2351}
2352
2353
2354template <>
2355Real FE<3,BERNSTEIN>::shape_second_deriv(const ElemType,
2356 const Order,
2357 const unsigned int,
2358 const unsigned int,
2359 const Point &)
2360{
2361 libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge and face orientation is needed.");
2362 return 0.;
2363}
2364
2365
2366template <>
2367Real FE<3,BERNSTEIN>::shape_second_deriv(const FEType fet,
2368 const Elem * elem,
2369 const unsigned int i,
2370 const unsigned int j,
2371 const Point & p,
2372 const bool add_p_level)
2373{
2374 return FE<3,BERNSTEIN>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
2375}
2376
2377
2378
2379
2380#endif
2381
2382} // namespace libMesh
2383
2384
2385
2386#endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition elem.h:2462
libMesh::Elem::positive_edge_orientation
bool positive_edge_orientation(const unsigned int i) const
Definition elem.C:3634
libMesh::Elem::positive_face_orientation
bool positive_face_orientation(const unsigned int i) const
Definition elem.C:3643
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Elem::p_level
unsigned int p_level() const
Definition elem.h:3125
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition fe_type.h:197
libMesh::FEType::order
OrderWrapper order
The approximation order of the element (at 0 p-refinement level).
Definition fe_type.h:203
libMesh::FE
A specific instantiation of the FEBase class.
Definition fe.h:128
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::FE::shape_deriv
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
b
static const Real b
Definition disjoint_neighbor_test.C:24
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::LIBMESH_DEFAULT_VECTORIZED_FE
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
Definition fe_bernstein_shape_0D.C:30
libMesh::fe_fdm_deriv
OutputShape fe_fdm_deriv(const Elem *elem, const Order order, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level, OutputShape(*shape_func)(const Elem *, const Order, const unsigned int, const Point &, const bool))
Helper functions for finite differenced derivatives in cases where analytical calculations haven't be...
Definition fe.C:911
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition libmesh_common.h:569
libMesh::fe_fdm_second_deriv
OutputShape fe_fdm_second_deriv(const Elem *elem, const Order order, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level, OutputShape(*deriv_func)(const Elem *, const Order, const unsigned int, const unsigned int, const Point &, const bool))
Definition fe.C:969
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
Generated on Fri Aug 21 2026 13:32:28 for libMesh by doxygen 1.9.8