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 : #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 :
29 : namespace {
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 106721784 : 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 115615266 : std::min(elem->point(c),elem->point(d)));
61 : }
62 : } // anonymous namespace
63 :
64 :
65 :
66 : namespace libMesh
67 : {
68 :
69 :
70 50733132 : LIBMESH_DEFAULT_VECTORIZED_FE(3,BERNSTEIN)
71 :
72 :
73 : template <>
74 1146271374 : Real 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 88614784 : libmesh_assert(elem);
84 1146271374 : const ElemType type = elem->type();
85 :
86 : const Order totalorder =
87 1234862830 : order + add_p_level*elem->p_level();
88 :
89 215093870 : auto hex_remap = [i, elem] (const Point & p_in,
90 : const unsigned int * hex_i0,
91 : const unsigned int * hex_i1,
92 720521699 : const unsigned int * hex_i2) {
93 : // Compute hex shape functions as a tensor-product
94 258112644 : const Real xi = p_in(0);
95 258112644 : const Real eta = p_in(1);
96 258112644 : const Real zeta = p_in(2);
97 258112644 : Point p_to_remap = p_in;
98 21509387 : Real & xi_mapped = p_to_remap(0);
99 21509387 : Real & eta_mapped = p_to_remap(1);
100 21509387 : Real & zeta_mapped = p_to_remap(2);
101 :
102 : // handle the edge orientation
103 : {
104 : // Edge 0
105 258112644 : if ((hex_i0[i] >= 2) && (hex_i1[i] == 0) && (hex_i2[i] == 0))
106 : {
107 6803052 : if (elem->positive_edge_orientation(0))
108 0 : xi_mapped = -xi;
109 : }
110 : // Edge 1
111 251309592 : else if ((hex_i0[i] == 1) && (hex_i1[i] >= 2) && (hex_i2[i] == 0))
112 : {
113 7369973 : if (elem->positive_edge_orientation(1))
114 6803052 : eta_mapped = -eta;
115 : }
116 : // edge 2
117 244506540 : else if ((hex_i0[i] >= 2) && (hex_i1[i] == 1) && (hex_i2[i] == 0))
118 : {
119 6803052 : if (!elem->positive_edge_orientation(2))
120 0 : xi_mapped = -xi;
121 : }
122 : // edge 3
123 237703488 : else if ((hex_i0[i] == 0) && (hex_i1[i] >= 2) && (hex_i2[i] == 0))
124 : {
125 7369973 : if (elem->positive_edge_orientation(3))
126 6803052 : eta_mapped = -eta;
127 : }
128 : // edge 4
129 230900436 : else if ((hex_i0[i] == 0) && (hex_i1[i] == 0) && (hex_i2[i] >=2 ))
130 : {
131 7369973 : if (elem->positive_edge_orientation(4))
132 6803052 : zeta_mapped = -zeta;
133 : }
134 : // edge 5
135 224097384 : else if ((hex_i0[i] == 1) && (hex_i1[i] == 0) && (hex_i2[i] >=2 ))
136 : {
137 7369973 : if (elem->positive_edge_orientation(5))
138 6803052 : zeta_mapped = -zeta;
139 : }
140 : // edge 6
141 217294332 : else if ((hex_i0[i] == 1) && (hex_i1[i] == 1) && (hex_i2[i] >=2 ))
142 : {
143 7369973 : if (elem->positive_edge_orientation(6))
144 6803052 : zeta_mapped = -zeta;
145 : }
146 : // edge 7
147 210491280 : else if ((hex_i0[i] == 0) && (hex_i1[i] == 1) && (hex_i2[i] >=2 ))
148 : {
149 7369973 : if (elem->positive_edge_orientation(7))
150 6803052 : zeta_mapped = -zeta;
151 : }
152 : // edge 8
153 203688228 : else if ((hex_i0[i] >=2 ) && (hex_i1[i] == 0) && (hex_i2[i] == 1))
154 : {
155 6803052 : if (elem->positive_edge_orientation(8))
156 0 : xi_mapped = -xi;
157 : }
158 : // edge 9
159 196885176 : else if ((hex_i0[i] == 1) && (hex_i1[i] >=2 ) && (hex_i2[i] == 1))
160 : {
161 7369973 : if (elem->positive_edge_orientation(9))
162 6803052 : eta_mapped = -eta;
163 : }
164 : // edge 10
165 190082124 : else if ((hex_i0[i] >=2 ) && (hex_i1[i] == 1) && (hex_i2[i] == 1))
166 : {
167 6803052 : if (!elem->positive_edge_orientation(10))
168 0 : xi_mapped = -xi;
169 : }
170 : // edge 11
171 183279072 : else if ((hex_i0[i] == 0) && (hex_i1[i] >=2 ) && (hex_i2[i] == 1))
172 : {
173 6803052 : if (elem->positive_edge_orientation(11))
174 6803052 : eta_mapped = -eta;
175 : }
176 : }
177 :
178 :
179 : // handle the face orientation
180 : {
181 : // face 0
182 258112644 : if ((hex_i2[i] == 0) && (hex_i0[i] >= 2) && (hex_i1[i] >= 2))
183 : {
184 17786964 : const Point min_point = get_min_point(elem, 1, 2, 0, 3);
185 :
186 2964494 : if (elem->point(0) == min_point)
187 0 : if (elem->positive_face_orientation(0))
188 : {
189 : // case 1
190 0 : xi_mapped = xi;
191 0 : eta_mapped = eta;
192 : }
193 : else
194 : {
195 : // case 2
196 0 : xi_mapped = eta;
197 0 : eta_mapped = xi;
198 : }
199 :
200 1482247 : else if (elem->point(3) == min_point)
201 17786964 : if (elem->positive_face_orientation(0))
202 : {
203 : // case 3
204 0 : xi_mapped = -eta;
205 0 : eta_mapped = xi;
206 : }
207 : else
208 : {
209 : // case 4
210 17786964 : xi_mapped = xi;
211 17786964 : eta_mapped = -eta;
212 : }
213 :
214 0 : else if (elem->point(2) == min_point)
215 0 : if (elem->positive_face_orientation(0))
216 : {
217 : // case 5
218 0 : xi_mapped = -xi;
219 0 : eta_mapped = -eta;
220 : }
221 : else
222 : {
223 : // case 6
224 0 : xi_mapped = -eta;
225 0 : eta_mapped = -xi;
226 : }
227 :
228 0 : else if (elem->point(1) == min_point)
229 : {
230 0 : if (elem->positive_face_orientation(0))
231 : {
232 : // case 7
233 0 : xi_mapped = eta;
234 0 : eta_mapped = -xi;
235 : }
236 : else
237 : {
238 : // Case 8
239 0 : xi_mapped = -xi;
240 0 : eta_mapped = eta;
241 : }
242 1482247 : }
243 14822470 : }
244 :
245 :
246 : // Face 1
247 240325680 : else if ((hex_i1[i] == 0) && (hex_i0[i] >= 2) && (hex_i2[i] >= 2))
248 : {
249 17786964 : const Point min_point = get_min_point(elem, 0, 1, 5, 4);
250 :
251 2964494 : if (elem->point(0) == min_point)
252 0 : if (!elem->positive_face_orientation(1))
253 : {
254 : // Case 1
255 0 : xi_mapped = xi;
256 0 : zeta_mapped = zeta;
257 : }
258 : else
259 : {
260 : // Case 2
261 0 : xi_mapped = zeta;
262 0 : zeta_mapped = xi;
263 : }
264 :
265 1482247 : else if (elem->point(1) == min_point)
266 0 : if (!elem->positive_face_orientation(1))
267 : {
268 : // Case 3
269 0 : xi_mapped = zeta;
270 0 : zeta_mapped = -xi;
271 : }
272 : else
273 : {
274 : // Case 4
275 0 : xi_mapped = -xi;
276 0 : zeta_mapped = zeta;
277 : }
278 :
279 1482247 : else if (elem->point(5) == min_point)
280 0 : if (!elem->positive_face_orientation(1))
281 : {
282 : // Case 5
283 0 : xi_mapped = -xi;
284 0 : zeta_mapped = -zeta;
285 : }
286 : else
287 : {
288 : // Case 6
289 0 : xi_mapped = -zeta;
290 0 : zeta_mapped = -xi;
291 : }
292 :
293 1482247 : else if (elem->point(4) == min_point)
294 : {
295 17786964 : if (!elem->positive_face_orientation(1))
296 : {
297 : // Case 7
298 17786964 : xi_mapped = -xi;
299 17786964 : zeta_mapped = zeta;
300 : }
301 : else
302 : {
303 : // Case 8
304 0 : xi_mapped = xi;
305 0 : zeta_mapped = -zeta;
306 : }
307 1482247 : }
308 14822470 : }
309 :
310 :
311 : // Face 2
312 222538716 : else if ((hex_i0[i] == 1) && (hex_i1[i] >= 2) && (hex_i2[i] >= 2))
313 : {
314 17786964 : const Point min_point = get_min_point(elem, 1, 2, 6, 5);
315 :
316 2964494 : if (elem->point(1) == min_point)
317 0 : if (!elem->positive_face_orientation(2))
318 : {
319 : // Case 1
320 0 : eta_mapped = eta;
321 0 : zeta_mapped = zeta;
322 : }
323 : else
324 : {
325 : // Case 2
326 0 : eta_mapped = zeta;
327 0 : zeta_mapped = eta;
328 : }
329 :
330 1482247 : else if (elem->point(2) == min_point)
331 0 : if (!elem->positive_face_orientation(2))
332 : {
333 : // Case 3
334 0 : eta_mapped = zeta;
335 0 : zeta_mapped = -eta;
336 : }
337 : else
338 : {
339 : // Case 4
340 0 : eta_mapped = -eta;
341 0 : zeta_mapped = zeta;
342 : }
343 :
344 1482247 : else if (elem->point(6) == min_point)
345 17786964 : if (!elem->positive_face_orientation(2))
346 : {
347 : // Case 5
348 0 : eta_mapped = -eta;
349 0 : zeta_mapped = -zeta;
350 : }
351 : else
352 : {
353 : // Case 6
354 17786964 : eta_mapped = -zeta;
355 17786964 : zeta_mapped = -eta;
356 : }
357 :
358 0 : else if (elem->point(5) == min_point)
359 : {
360 0 : if (!elem->positive_face_orientation(2))
361 : {
362 : // Case 7
363 0 : eta_mapped = -zeta;
364 0 : zeta_mapped = eta;
365 : }
366 : else
367 : {
368 : // Case 8
369 0 : eta_mapped = eta;
370 0 : zeta_mapped = -zeta;
371 : }
372 1482247 : }
373 14822470 : }
374 :
375 :
376 : // Face 3
377 204751752 : else if ((hex_i1[i] == 1) && (hex_i0[i] >= 2) && (hex_i2[i] >= 2))
378 : {
379 17786964 : const Point min_point = get_min_point(elem, 2, 3, 7, 6);
380 :
381 2964494 : if (elem->point(3) == min_point)
382 0 : if (elem->positive_face_orientation(3))
383 : {
384 : // Case 1
385 0 : xi_mapped = xi;
386 0 : zeta_mapped = zeta;
387 : }
388 : else
389 : {
390 : // Case 2
391 0 : xi_mapped = zeta;
392 0 : zeta_mapped = xi;
393 : }
394 :
395 1482247 : else if (elem->point(7) == min_point)
396 17786964 : if (elem->positive_face_orientation(3))
397 : {
398 : // Case 3
399 17786964 : xi_mapped = -zeta;
400 17786964 : zeta_mapped = xi;
401 : }
402 : else
403 : {
404 : // Case 4
405 0 : xi_mapped = xi;
406 0 : zeta_mapped = -zeta;
407 : }
408 :
409 0 : else if (elem->point(6) == min_point)
410 0 : if (elem->positive_face_orientation(3))
411 : {
412 : // Case 5
413 0 : xi_mapped = -xi;
414 0 : zeta_mapped = -zeta;
415 : }
416 : else
417 : {
418 : // Case 6
419 0 : xi_mapped = -zeta;
420 0 : zeta_mapped = -xi;
421 : }
422 :
423 0 : else if (elem->point(2) == min_point)
424 : {
425 0 : if (elem->positive_face_orientation(3))
426 : {
427 : // Case 7
428 0 : xi_mapped = zeta;
429 0 : zeta_mapped = -xi;
430 : }
431 : else
432 : {
433 : // Case 8
434 0 : xi_mapped = -xi;
435 0 : zeta_mapped = zeta;
436 : }
437 1482247 : }
438 14822470 : }
439 :
440 :
441 : // Face 4
442 186964788 : else if ((hex_i0[i] == 0) && (hex_i1[i] >= 2) && (hex_i2[i] >= 2))
443 : {
444 17786964 : const Point min_point = get_min_point(elem, 3, 0, 4, 7);
445 :
446 2964494 : if (elem->point(0) == min_point)
447 0 : if (elem->positive_face_orientation(4))
448 : {
449 : // Case 1
450 0 : eta_mapped = eta;
451 0 : zeta_mapped = zeta;
452 : }
453 : else
454 : {
455 : // Case 2
456 0 : eta_mapped = zeta;
457 0 : zeta_mapped = eta;
458 : }
459 :
460 1482247 : else if (elem->point(4) == min_point)
461 0 : if (elem->positive_face_orientation(4))
462 : {
463 : // Case 3
464 0 : eta_mapped = -zeta;
465 0 : zeta_mapped = eta;
466 : }
467 : else
468 : {
469 : // Case 4
470 0 : eta_mapped = eta;
471 0 : zeta_mapped = -zeta;
472 : }
473 :
474 1482247 : else if (elem->point(7) == min_point)
475 17786964 : if (elem->positive_face_orientation(4))
476 : {
477 : // Case 5
478 0 : eta_mapped = -eta;
479 0 : zeta_mapped = -zeta;
480 : }
481 : else
482 : {
483 : // Case 6
484 17786964 : eta_mapped = -zeta;
485 17786964 : zeta_mapped = -eta;
486 : }
487 :
488 0 : else if (elem->point(3) == min_point)
489 : {
490 0 : if (elem->positive_face_orientation(4))
491 : {
492 : // Case 7
493 0 : eta_mapped = zeta;
494 0 : zeta_mapped = -eta;
495 : }
496 : else
497 : {
498 : // Case 8
499 0 : eta_mapped = -eta;
500 0 : zeta_mapped = zeta;
501 : }
502 1482247 : }
503 14822470 : }
504 :
505 :
506 : // Face 5
507 169177824 : else if ((hex_i2[i] == 1) && (hex_i0[i] >= 2) && (hex_i1[i] >= 2))
508 : {
509 17786964 : const Point min_point = get_min_point(elem, 4, 5, 6, 7);
510 :
511 2964494 : if (elem->point(4) == min_point)
512 0 : if (!elem->positive_face_orientation(5))
513 : {
514 : // Case 1
515 0 : xi_mapped = xi;
516 0 : eta_mapped = eta;
517 : }
518 : else
519 : {
520 : // Case 2
521 0 : xi_mapped = eta;
522 0 : eta_mapped = xi;
523 : }
524 :
525 1482247 : else if (elem->point(5) == min_point)
526 0 : if (!elem->positive_face_orientation(5))
527 : {
528 : // Case 3
529 0 : xi_mapped = eta;
530 0 : eta_mapped = -xi;
531 : }
532 : else
533 : {
534 : // Case 4
535 0 : xi_mapped = -xi;
536 0 : eta_mapped = eta;
537 : }
538 :
539 1482247 : else if (elem->point(6) == min_point)
540 0 : if (!elem->positive_face_orientation(5))
541 : {
542 : // Case 5
543 0 : xi_mapped = -xi;
544 0 : eta_mapped = -eta;
545 : }
546 : else
547 : {
548 : // Case 6
549 0 : xi_mapped = -eta;
550 0 : eta_mapped = -xi;
551 : }
552 :
553 1482247 : else if (elem->point(7) == min_point)
554 : {
555 17786964 : if (!elem->positive_face_orientation(5))
556 : {
557 : // Case 7
558 0 : xi_mapped = -eta;
559 0 : eta_mapped = xi;
560 : }
561 : else
562 : {
563 : // Case 8
564 17786964 : xi_mapped = xi;
565 17786964 : eta_mapped = eta;
566 : }
567 : }
568 : }
569 : }
570 :
571 258112644 : return p_to_remap;
572 1146271374 : };
573 :
574 1146271374 : switch (totalorder)
575 : {
576 : // 1st order Bernstein.
577 46520580 : case FIRST:
578 : {
579 43724620 : switch (type)
580 : {
581 :
582 : // Bernstein shape functions on the tetrahedron.
583 22087620 : case TET4:
584 : case TET10:
585 : case TET14:
586 : {
587 759880 : libmesh_assert_less (i, 4);
588 :
589 : // Area coordinates
590 22087620 : const Real zeta1 = p(0);
591 22087620 : const Real zeta2 = p(1);
592 22087620 : const Real zeta3 = p(2);
593 22087620 : const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
594 :
595 22087620 : switch(i)
596 : {
597 189970 : case 0: return zeta0;
598 5521905 : case 1: return zeta1;
599 5521905 : case 2: return zeta2;
600 5521905 : case 3: return zeta3;
601 :
602 0 : default:
603 0 : libmesh_error_msg("Invalid shape function index i = " << i);
604 : }
605 : }
606 :
607 : // Bernstein shape functions on the hexahedral.
608 24432960 : case HEX8:
609 : case HEX20:
610 : case HEX27:
611 : {
612 2036080 : libmesh_assert_less (i, 8);
613 :
614 : // Compute hex shape functions as a tensor-product
615 24432960 : const Real xi = p(0);
616 24432960 : const Real eta = p(1);
617 24432960 : 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 46829840 : return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
628 24432960 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)*
629 24432960 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta));
630 : }
631 :
632 :
633 0 : default:
634 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
635 : }
636 : }
637 :
638 :
639 :
640 :
641 841638150 : case SECOND:
642 : {
643 841638150 : switch (type)
644 : {
645 :
646 : // Bernstein shape functions on the tetrahedron.
647 8803320 : case TET10:
648 : case TET14:
649 : {
650 337220 : libmesh_assert_less (i, 10);
651 :
652 : // Area coordinates
653 8803320 : const Real zeta1 = p(0);
654 8803320 : const Real zeta2 = p(1);
655 8803320 : const Real zeta3 = p(2);
656 8803320 : const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
657 :
658 8803320 : switch(i)
659 : {
660 880332 : case 0: return zeta0*zeta0;
661 880332 : case 1: return zeta1*zeta1;
662 880332 : case 2: return zeta2*zeta2;
663 880332 : case 3: return zeta3*zeta3;
664 880332 : case 4: return 2.*zeta0*zeta1;
665 880332 : case 5: return 2.*zeta1*zeta2;
666 880332 : case 6: return 2.*zeta0*zeta2;
667 880332 : case 7: return 2.*zeta3*zeta0;
668 880332 : case 8: return 2.*zeta1*zeta3;
669 880332 : case 9: return 2.*zeta2*zeta3;
670 :
671 0 : default:
672 0 : libmesh_error_msg("Invalid shape function index i = " << i);
673 : }
674 : }
675 :
676 : // Bernstein shape functions on the 20-noded hexahedral.
677 151622400 : case HEX20:
678 : {
679 12635200 : libmesh_assert_less (i, 20);
680 :
681 : // Compute hex shape functions as a tensor-product
682 151622400 : const Real xi = p(0);
683 151622400 : const Real eta = p(1);
684 151622400 : const Real zeta = p(2);
685 :
686 290609600 : return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[i], xi)*
687 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[i], eta)*
688 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[i], zeta)
689 176892800 : +hex20_scal20[i]*
690 290609600 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[20], xi)*
691 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[20], eta)*
692 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[20], zeta)
693 176892800 : +hex20_scal21[i]*
694 290609600 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[21], xi)*
695 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[21], eta)*
696 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[21], zeta)
697 176892800 : +hex20_scal22[i]*
698 290609600 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[22], xi)*
699 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[22], eta)*
700 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[22], zeta)
701 176892800 : +hex20_scal23[i]*
702 290609600 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[23], xi)*
703 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[23], eta)*
704 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[23], zeta)
705 176892800 : +hex20_scal24[i]*
706 290609600 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[24], xi)*
707 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[24], eta)*
708 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[24], zeta)
709 176892800 : +hex20_scal25[i]*
710 290609600 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[25], xi)*
711 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[25], eta)*
712 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[25], zeta)
713 164257600 : +hex20_scal26[i]*
714 290609600 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[26], xi)*
715 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[26], eta)*
716 151622400 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[26], zeta));
717 : }
718 :
719 : // Bernstein shape functions on the hexahedral.
720 681212430 : case HEX27:
721 : {
722 51337017 : libmesh_assert_less (i, 27);
723 :
724 : // Compute hex shape functions as a tensor-product
725 681212430 : const Real xi = p(0);
726 681212430 : const Real eta = p(1);
727 681212430 : 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 1311087843 : return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
738 681212430 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)*
739 681212430 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta));
740 : }
741 :
742 :
743 0 : default:
744 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
745 : }
746 :
747 : }
748 :
749 :
750 :
751 : // 3rd-order Bernstein.
752 83910144 : case THIRD:
753 : {
754 83910144 : 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 83910144 : case HEX27:
824 : {
825 6992512 : 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 83910144 : const Point p_mapped = hex_remap(p, i0, i1, i2);
837 90902656 : return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], p_mapped(0))*
838 90902656 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], p_mapped(1))*
839 90902656 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], p_mapped(2)));
840 : }
841 :
842 0 : default:
843 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
844 : }
845 :
846 : }//case THIRD
847 :
848 :
849 : // 4th-order Bernstein.
850 174202500 : case FOURTH:
851 : {
852 174202500 : switch (type)
853 : {
854 :
855 : // Bernstein shape functions on the hexahedral.
856 174202500 : case HEX27:
857 : {
858 14516875 : 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 174202500 : const Point p_mapped = hex_remap(p, i0, i1, i2);
870 188719375 : return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], p_mapped(0))*
871 188719375 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], p_mapped(1))*
872 188719375 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], p_mapped(2)));
873 : }
874 :
875 :
876 0 : default:
877 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
878 : }
879 : }
880 :
881 :
882 0 : default:
883 0 : 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 :
893 : template <>
894 0 : Real FE<3,BERNSTEIN>::shape(const ElemType,
895 : const Order,
896 : const unsigned int,
897 : const Point &)
898 : {
899 0 : libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge and face orientation is needed.");
900 : return 0.;
901 : }
902 :
903 :
904 : template <>
905 0 : Real 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 0 : return FE<3,BERNSTEIN>::shape(elem, fet.order, i, p, add_p_level);
912 : }
913 :
914 : template <>
915 1609013682 : Real 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 129336075 : libmesh_assert(elem);
925 1609013682 : const ElemType type = elem->type();
926 :
927 : const Order totalorder =
928 1738349757 : order + add_p_level*elem->p_level();
929 :
930 129336075 : libmesh_assert_less (j, 3);
931 :
932 1609013682 : switch (totalorder)
933 : {
934 : // 1st order Bernstein.
935 96893364 : case FIRST:
936 : {
937 89292012 : switch (type)
938 : {
939 : // Bernstein shape functions on the tetrahedron.
940 9676020 : case TET4:
941 : case TET10:
942 : case TET14:
943 : {
944 9676020 : 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 87217344 : case HEX8:
950 : case HEX20:
951 : case HEX27:
952 : {
953 7268112 : libmesh_assert_less (i, 8);
954 :
955 : // Compute hex shape functions as a tensor-product
956 87217344 : const Real xi = p(0);
957 87217344 : const Real eta = p(1);
958 87217344 : 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 87217344 : switch (j)
969 : {
970 : // d()/dxi
971 3118480 : case 0:
972 71725040 : return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
973 37421760 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
974 37421760 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
975 :
976 : // d()/deta
977 2422704 : case 1:
978 55722192 : return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
979 29072448 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)*
980 29072448 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
981 :
982 : // d()/dzeta
983 1726928 : case 2:
984 39719344 : return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
985 20723136 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
986 20723136 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta));
987 :
988 0 : default:
989 0 : libmesh_error_msg("Invalid derivative index j = " << j);
990 : }
991 : }
992 :
993 0 : default:
994 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
995 : }
996 : }
997 :
998 :
999 :
1000 :
1001 1390225578 : case SECOND:
1002 : {
1003 1390225578 : switch (type)
1004 : {
1005 : // Bernstein shape functions on the tetrahedron.
1006 4003800 : case TET10:
1007 : case TET14:
1008 : {
1009 4003800 : 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 113716800 : case HEX20:
1014 : {
1015 9476400 : libmesh_assert_less (i, 20);
1016 :
1017 : // Compute hex shape functions as a tensor-product
1018 113716800 : const Real xi = p(0);
1019 113716800 : const Real eta = p(1);
1020 113716800 : const Real zeta = p(2);
1021 :
1022 113716800 : switch (j)
1023 : {
1024 : // d()/dxi
1025 3158800 : case 0:
1026 72652400 : return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[i], 0, xi)*
1027 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[i], eta)*
1028 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[i], zeta)
1029 44223200 : +hex20_scal20[i]*
1030 72652400 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[20], 0, xi)*
1031 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[20], eta)*
1032 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[20], zeta)
1033 44223200 : +hex20_scal21[i]*
1034 72652400 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[21], 0, xi)*
1035 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[21], eta)*
1036 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[21], zeta)
1037 44223200 : +hex20_scal22[i]*
1038 72652400 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[22], 0, xi)*
1039 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[22], eta)*
1040 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[22], zeta)
1041 44223200 : +hex20_scal23[i]*
1042 72652400 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[23], 0, xi)*
1043 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[23], eta)*
1044 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[23], zeta)
1045 44223200 : +hex20_scal24[i]*
1046 72652400 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[24], 0, xi)*
1047 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[24], eta)*
1048 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[24], zeta)
1049 44223200 : +hex20_scal25[i]*
1050 72652400 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[25], 0, xi)*
1051 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[25], eta)*
1052 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[25], zeta)
1053 41064400 : +hex20_scal26[i]*
1054 72652400 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[26], 0, xi)*
1055 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[26], eta)*
1056 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[26], zeta));
1057 :
1058 : // d()/deta
1059 3158800 : case 1:
1060 72652400 : return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[i], xi)*
1061 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[i], 0, eta)*
1062 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[i], zeta)
1063 44223200 : +hex20_scal20[i]*
1064 72652400 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[20], xi)*
1065 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[20], 0, eta)*
1066 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[20], zeta)
1067 44223200 : +hex20_scal21[i]*
1068 72652400 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[21], xi)*
1069 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[21], 0, eta)*
1070 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[21], zeta)
1071 44223200 : +hex20_scal22[i]*
1072 72652400 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[22], xi)*
1073 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[22], 0, eta)*
1074 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[22], zeta)
1075 44223200 : +hex20_scal23[i]*
1076 72652400 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[23], xi)*
1077 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[23], 0, eta)*
1078 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[23], zeta)
1079 44223200 : +hex20_scal24[i]*
1080 72652400 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[24], xi)*
1081 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[24], 0, eta)*
1082 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[24], zeta)
1083 44223200 : +hex20_scal25[i]*
1084 72652400 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[25], xi)*
1085 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[25], 0, eta)*
1086 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[25], zeta)
1087 41064400 : +hex20_scal26[i]*
1088 72652400 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[26], xi)*
1089 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[26], 0, eta)*
1090 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[26], zeta));
1091 :
1092 : // d()/dzeta
1093 3158800 : case 2:
1094 72652400 : return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[i], xi)*
1095 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[i], eta)*
1096 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[i], 0, zeta)
1097 44223200 : +hex20_scal20[i]*
1098 72652400 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[20], xi)*
1099 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[20], eta)*
1100 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[20], 0, zeta)
1101 44223200 : +hex20_scal21[i]*
1102 72652400 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[21], xi)*
1103 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[21], eta)*
1104 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[21], 0, zeta)
1105 44223200 : +hex20_scal22[i]*
1106 72652400 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[22], xi)*
1107 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[22], eta)*
1108 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[22], 0, zeta)
1109 44223200 : +hex20_scal23[i]*
1110 72652400 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[23], xi)*
1111 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[23], eta)*
1112 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[23], 0, zeta)
1113 44223200 : +hex20_scal24[i]*
1114 72652400 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[24], xi)*
1115 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[24], eta)*
1116 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[24], 0, zeta)
1117 44223200 : +hex20_scal25[i]*
1118 72652400 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[25], xi)*
1119 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[25], eta)*
1120 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[25], 0, zeta)
1121 41064400 : +hex20_scal26[i]*
1122 72652400 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[26], xi)*
1123 37905600 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[26], eta)*
1124 37905600 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[26], 0, zeta));
1125 :
1126 0 : default:
1127 0 : libmesh_error_msg("Invalid derivative index j = " << j);
1128 : }
1129 : }
1130 :
1131 : // Bernstein shape functions on the hexahedral.
1132 1272504978 : case HEX27:
1133 : {
1134 101953728 : libmesh_assert_less (i, 27);
1135 :
1136 : // Compute hex shape functions as a tensor-product
1137 1272504978 : const Real xi = p(0);
1138 1272504978 : const Real eta = p(1);
1139 1272504978 : 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 1272504978 : switch (j)
1150 : {
1151 : // d()/dxi
1152 41573628 : case 0:
1153 988900272 : return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1154 515236950 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1155 515236950 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
1156 :
1157 : // d()/deta
1158 33984576 : case 1:
1159 814352076 : return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1160 424168326 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)*
1161 424168326 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
1162 :
1163 : // d()/dzeta
1164 26395524 : case 2:
1165 639803880 : return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1166 333099702 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1167 333099702 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta));
1168 :
1169 0 : default:
1170 0 : libmesh_error_msg("Invalid derivative index j = " << j);
1171 : }
1172 : }
1173 :
1174 :
1175 0 : default:
1176 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1177 : }
1178 : }
1179 :
1180 :
1181 :
1182 : // 3rd-order Bernstein.
1183 39882240 : case THIRD:
1184 : {
1185 39882240 : 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 39882240 : case HEX27:
1238 : {
1239 39882240 : 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 0 : default:
1762 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1763 : }
1764 : }
1765 :
1766 : // 4th-order Bernstein.
1767 82012500 : case FOURTH:
1768 : {
1769 82012500 : switch (type)
1770 : {
1771 :
1772 : // Bernstein shape functions derivatives on the hexahedral.
1773 82012500 : case HEX27:
1774 : {
1775 82012500 : 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 0 : default:
2298 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
2299 : }
2300 : }
2301 :
2302 :
2303 0 : default:
2304 0 : 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 :
2314 : template <>
2315 0 : Real FE<3,BERNSTEIN>::shape_deriv(const ElemType,
2316 : const Order,
2317 : const unsigned int,
2318 : const unsigned int,
2319 : const Point & )
2320 : {
2321 0 : libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge and face orientation is needed.");
2322 : return 0.;
2323 : }
2324 :
2325 : template <>
2326 0 : Real 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 0 : 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 :
2341 : template <>
2342 351367344 : Real 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 351367344 : return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
2350 351367344 : FE<3,BERNSTEIN>::shape_deriv);
2351 : }
2352 :
2353 :
2354 : template <>
2355 0 : Real FE<3,BERNSTEIN>::shape_second_deriv(const ElemType,
2356 : const Order,
2357 : const unsigned int,
2358 : const unsigned int,
2359 : const Point &)
2360 : {
2361 0 : libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge and face orientation is needed.");
2362 : return 0.;
2363 : }
2364 :
2365 :
2366 : template <>
2367 0 : Real 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 0 : 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
|
