Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : // C++ includes
20 :
21 : // Local includes
22 : #include "libmesh/fe.h"
23 :
24 : #include "libmesh/elem.h"
25 : #include "libmesh/int_range.h"
26 :
27 :
28 : namespace libMesh
29 : {
30 :
31 :
32 2736 : LIBMESH_DEFAULT_VECTORIZED_FE(3,XYZ)
33 :
34 :
35 : template <>
36 25543291 : Real FE<3,XYZ>::shape(const Elem * elem,
37 : const Order libmesh_dbg_var(order),
38 : const unsigned int i,
39 : const Point & point_in,
40 : const bool libmesh_dbg_var(add_p_level))
41 : {
42 : #if LIBMESH_DIM == 3
43 6467651 : libmesh_assert(elem);
44 :
45 25543291 : Point avg = elem->vertex_average();
46 6467651 : Point max_distance = Point(0.,0.,0.);
47 411188081 : for (auto p : make_range(elem->n_nodes()))
48 1542579160 : for (unsigned int d = 0; d < 3; d++)
49 : {
50 1743793572 : const Real distance = std::abs(avg(d) - elem->point(p)(d));
51 1156934370 : max_distance(d) = std::max(distance, max_distance(d));
52 : }
53 :
54 25543291 : const Real x = point_in(0);
55 25543291 : const Real y = point_in(1);
56 25543291 : const Real z = point_in(2);
57 25543291 : const Real xc = avg(0);
58 25543291 : const Real yc = avg(1);
59 25543291 : const Real zc = avg(2);
60 25543291 : const Real distx = max_distance(0);
61 25543291 : const Real disty = max_distance(1);
62 25543291 : const Real distz = max_distance(2);
63 25543291 : const Real dx = (x - xc)/distx;
64 25543291 : const Real dy = (y - yc)/disty;
65 25543291 : const Real dz = (z - zc)/distz;
66 :
67 : #ifndef NDEBUG
68 : // totalorder is only used in the assertion below, so
69 : // we avoid declaring it when asserts are not active.
70 6467651 : const unsigned int totalorder = order + add_p_level * elem->p_level();
71 : #endif
72 6467651 : libmesh_assert_less (i, (totalorder+1) * (totalorder+2) *
73 : (totalorder+3)/6);
74 :
75 : // monomials. since they are hierarchic we only need one case block.
76 25543291 : switch (i)
77 : {
78 : // constant
79 830325 : case 0:
80 830325 : return 1.;
81 :
82 : // linears
83 3306288 : case 1:
84 3306288 : return dx;
85 :
86 3306288 : case 2:
87 3306288 : return dy;
88 :
89 3306288 : case 3:
90 3306288 : return dz;
91 :
92 : // quadratics
93 1080984 : case 4:
94 1080984 : return dx*dx;
95 :
96 1080984 : case 5:
97 1080984 : return dx*dy;
98 :
99 1080984 : case 6:
100 1080984 : return dy*dy;
101 :
102 1080984 : case 7:
103 1080984 : return dx*dz;
104 :
105 1080984 : case 8:
106 1080984 : return dz*dy;
107 :
108 1080984 : case 9:
109 1080984 : return dz*dz;
110 :
111 : // cubics
112 332761 : case 10:
113 332761 : return dx*dx*dx;
114 :
115 332761 : case 11:
116 332761 : return dx*dx*dy;
117 :
118 332761 : case 12:
119 332761 : return dx*dy*dy;
120 :
121 332761 : case 13:
122 332761 : return dy*dy*dy;
123 :
124 332761 : case 14:
125 332761 : return dx*dx*dz;
126 :
127 332761 : case 15:
128 332761 : return dx*dy*dz;
129 :
130 332761 : case 16:
131 332761 : return dy*dy*dz;
132 :
133 332761 : case 17:
134 332761 : return dx*dz*dz;
135 :
136 332761 : case 18:
137 332761 : return dy*dz*dz;
138 :
139 332761 : case 19:
140 332761 : return dz*dz*dz;
141 :
142 : // quartics
143 166975 : case 20:
144 166975 : return dx*dx*dx*dx;
145 :
146 166975 : case 21:
147 166975 : return dx*dx*dx*dy;
148 :
149 166975 : case 22:
150 166975 : return dx*dx*dy*dy;
151 :
152 166975 : case 23:
153 166975 : return dx*dy*dy*dy;
154 :
155 166975 : case 24:
156 166975 : return dy*dy*dy*dy;
157 :
158 166975 : case 25:
159 166975 : return dx*dx*dx*dz;
160 :
161 166975 : case 26:
162 166975 : return dx*dx*dy*dz;
163 :
164 166975 : case 27:
165 166975 : return dx*dy*dy*dz;
166 :
167 166975 : case 28:
168 166975 : return dy*dy*dy*dz;
169 :
170 166975 : case 29:
171 166975 : return dx*dx*dz*dz;
172 :
173 166975 : case 30:
174 166975 : return dx*dy*dz*dz;
175 :
176 166975 : case 31:
177 166975 : return dy*dy*dz*dz;
178 :
179 166975 : case 32:
180 166975 : return dx*dz*dz*dz;
181 :
182 166975 : case 33:
183 166975 : return dy*dz*dz*dz;
184 :
185 166975 : case 34:
186 166975 : return dz*dz*dz*dz;
187 :
188 0 : default:
189 0 : unsigned int o = 0;
190 0 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
191 0 : unsigned int i2 = i - (o*(o+1)*(o+2)/6);
192 0 : unsigned int block=o, nz = 0;
193 0 : for (; block < i2; block += (o-nz+1)) { nz++; }
194 0 : const unsigned int nx = block - i2;
195 0 : const unsigned int ny = o - nx - nz;
196 0 : Real val = 1.;
197 0 : for (unsigned int index=0; index != nx; index++)
198 0 : val *= dx;
199 0 : for (unsigned int index=0; index != ny; index++)
200 0 : val *= dy;
201 0 : for (unsigned int index=0; index != nz; index++)
202 0 : val *= dz;
203 0 : return val;
204 : }
205 :
206 : #else // LIBMESH_DIM != 3
207 : libmesh_assert(true || order || add_p_level);
208 : libmesh_ignore(elem, i, point_in);
209 : libmesh_not_implemented();
210 : #endif
211 : }
212 :
213 :
214 :
215 : template <>
216 0 : Real FE<3,XYZ>::shape(const ElemType,
217 : const Order,
218 : const unsigned int,
219 : const Point &)
220 : {
221 0 : libmesh_error_msg("XYZ polynomials require the element.");
222 : return 0.;
223 : }
224 :
225 :
226 :
227 : template <>
228 0 : Real FE<3,XYZ>::shape(const FEType fet,
229 : const Elem * elem,
230 : const unsigned int i,
231 : const Point & p,
232 : const bool add_p_level)
233 : {
234 0 : return FE<3,XYZ>::shape(elem, fet.order, i, p, add_p_level);
235 : }
236 :
237 :
238 :
239 :
240 :
241 : template <>
242 56587872 : Real FE<3,XYZ>::shape_deriv(const Elem * elem,
243 : const Order libmesh_dbg_var(order),
244 : const unsigned int i,
245 : const unsigned int j,
246 : const Point & point_in,
247 : const bool libmesh_dbg_var(add_p_level))
248 : {
249 : #if LIBMESH_DIM == 3
250 :
251 14379165 : libmesh_assert(elem);
252 14379165 : libmesh_assert_less (j, 3);
253 :
254 56587872 : Point avg = elem->vertex_average();
255 14379165 : Point max_distance = Point(0.,0.,0.);
256 920194230 : for (auto p : make_range(elem->n_nodes()))
257 3454425432 : for (unsigned int d = 0; d < 3; d++)
258 : {
259 3910073148 : const Real distance = std::abs(avg(d) - elem->point(p)(d));
260 2590819074 : max_distance(d) = std::max(distance, max_distance(d));
261 : }
262 :
263 56587872 : const Real x = point_in(0);
264 56587872 : const Real y = point_in(1);
265 56587872 : const Real z = point_in(2);
266 56587872 : const Real xc = avg(0);
267 56587872 : const Real yc = avg(1);
268 56587872 : const Real zc = avg(2);
269 56587872 : const Real distx = max_distance(0);
270 56587872 : const Real disty = max_distance(1);
271 56587872 : const Real distz = max_distance(2);
272 56587872 : const Real dx = (x - xc)/distx;
273 56587872 : const Real dy = (y - yc)/disty;
274 56587872 : const Real dz = (z - zc)/distz;
275 :
276 : #ifndef NDEBUG
277 : // totalorder is only used in the assertion below, so
278 : // we avoid declaring it when asserts are not active.
279 14379165 : const unsigned int totalorder = order + add_p_level*elem->p_level();
280 : #endif
281 14379165 : libmesh_assert_less (i, (totalorder+1) * (totalorder+2) *
282 : (totalorder+3)/6);
283 :
284 56587872 : switch (j)
285 : {
286 : // d()/dx
287 18862624 : case 0:
288 : {
289 18862624 : switch (i)
290 : {
291 : // constant
292 731085 : case 0:
293 731085 : return 0.;
294 :
295 : // linear
296 2909986 : case 1:
297 2909986 : return 1./distx;
298 :
299 731085 : case 2:
300 731085 : return 0.;
301 :
302 731085 : case 3:
303 731085 : return 0.;
304 :
305 : // quadratic
306 779405 : case 4:
307 779405 : return 2.*dx/distx;
308 :
309 779405 : case 5:
310 779405 : return dy/distx;
311 :
312 198060 : case 6:
313 198060 : return 0.;
314 :
315 779405 : case 7:
316 779405 : return dz/distx;
317 :
318 198060 : case 8:
319 198060 : return 0.;
320 :
321 198060 : case 9:
322 198060 : return 0.;
323 :
324 : // cubic
325 144390 : case 10:
326 144390 : return 3.*dx*dx/distx;
327 :
328 144390 : case 11:
329 144390 : return 2.*dx*dy/distx;
330 :
331 144390 : case 12:
332 144390 : return dy*dy/distx;
333 :
334 38583 : case 13:
335 38583 : return 0.;
336 :
337 144390 : case 14:
338 144390 : return 2.*dx*dz/distx;
339 :
340 144390 : case 15:
341 144390 : return dy*dz/distx;
342 :
343 38583 : case 16:
344 38583 : return 0.;
345 :
346 144390 : case 17:
347 144390 : return dz*dz/distx;
348 :
349 38583 : case 18:
350 38583 : return 0.;
351 :
352 38583 : case 19:
353 38583 : return 0.;
354 :
355 : // quartics
356 73490 : case 20:
357 73490 : return 4.*dx*dx*dx/distx;
358 :
359 73490 : case 21:
360 73490 : return 3.*dx*dx*dy/distx;
361 :
362 73490 : case 22:
363 73490 : return 2.*dx*dy*dy/distx;
364 :
365 73490 : case 23:
366 73490 : return dy*dy*dy/distx;
367 :
368 19635 : case 24:
369 19635 : return 0.;
370 :
371 73490 : case 25:
372 73490 : return 3.*dx*dx*dz/distx;
373 :
374 73490 : case 26:
375 73490 : return 2.*dx*dy*dz/distx;
376 :
377 73490 : case 27:
378 73490 : return dy*dy*dz/distx;
379 :
380 19635 : case 28:
381 19635 : return 0.;
382 :
383 73490 : case 29:
384 73490 : return 2.*dx*dz*dz/distx;
385 :
386 73490 : case 30:
387 73490 : return dy*dz*dz/distx;
388 :
389 19635 : case 31:
390 19635 : return 0.;
391 :
392 73490 : case 32:
393 73490 : return dz*dz*dz/distx;
394 :
395 19635 : case 33:
396 19635 : return 0.;
397 :
398 19635 : case 34:
399 19635 : return 0.;
400 :
401 0 : default:
402 0 : unsigned int o = 0;
403 0 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
404 0 : unsigned int i2 = i - (o*(o+1)*(o+2)/6);
405 0 : unsigned int block=o, nz = 0;
406 0 : for (; block < i2; block += (o-nz+1)) { nz++; }
407 0 : const unsigned int nx = block - i2;
408 0 : const unsigned int ny = o - nx - nz;
409 0 : Real val = nx;
410 0 : for (unsigned int index=1; index < nx; index++)
411 0 : val *= dx;
412 0 : for (unsigned int index=0; index != ny; index++)
413 0 : val *= dy;
414 0 : for (unsigned int index=0; index != nz; index++)
415 0 : val *= dz;
416 0 : return val/distx;
417 : }
418 : }
419 :
420 :
421 : // d()/dy
422 18862624 : case 1:
423 : {
424 18862624 : switch (i)
425 : {
426 : // constant
427 731085 : case 0:
428 731085 : return 0.;
429 :
430 : // linear
431 731085 : case 1:
432 731085 : return 0.;
433 :
434 2909986 : case 2:
435 2909986 : return 1./disty;
436 :
437 731085 : case 3:
438 731085 : return 0.;
439 :
440 : // quadratic
441 198060 : case 4:
442 198060 : return 0.;
443 :
444 779405 : case 5:
445 779405 : return dx/disty;
446 :
447 779405 : case 6:
448 779405 : return 2.*dy/disty;
449 :
450 198060 : case 7:
451 198060 : return 0.;
452 :
453 779405 : case 8:
454 779405 : return dz/disty;
455 :
456 198060 : case 9:
457 198060 : return 0.;
458 :
459 : // cubic
460 38583 : case 10:
461 38583 : return 0.;
462 :
463 144390 : case 11:
464 144390 : return dx*dx/disty;
465 :
466 144390 : case 12:
467 144390 : return 2.*dx*dy/disty;
468 :
469 144390 : case 13:
470 144390 : return 3.*dy*dy/disty;
471 :
472 38583 : case 14:
473 38583 : return 0.;
474 :
475 144390 : case 15:
476 144390 : return dx*dz/disty;
477 :
478 144390 : case 16:
479 144390 : return 2.*dy*dz/disty;
480 :
481 38583 : case 17:
482 38583 : return 0.;
483 :
484 144390 : case 18:
485 144390 : return dz*dz/disty;
486 :
487 38583 : case 19:
488 38583 : return 0.;
489 :
490 : // quartics
491 19635 : case 20:
492 19635 : return 0.;
493 :
494 73490 : case 21:
495 73490 : return dx*dx*dx/disty;
496 :
497 73490 : case 22:
498 73490 : return 2.*dx*dx*dy/disty;
499 :
500 73490 : case 23:
501 73490 : return 3.*dx*dy*dy/disty;
502 :
503 73490 : case 24:
504 73490 : return 4.*dy*dy*dy/disty;
505 :
506 19635 : case 25:
507 19635 : return 0.;
508 :
509 73490 : case 26:
510 73490 : return dx*dx*dz/disty;
511 :
512 73490 : case 27:
513 73490 : return 2.*dx*dy*dz/disty;
514 :
515 73490 : case 28:
516 73490 : return 3.*dy*dy*dz/disty;
517 :
518 19635 : case 29:
519 19635 : return 0.;
520 :
521 73490 : case 30:
522 73490 : return dx*dz*dz/disty;
523 :
524 73490 : case 31:
525 73490 : return 2.*dy*dz*dz/disty;
526 :
527 19635 : case 32:
528 19635 : return 0.;
529 :
530 73490 : case 33:
531 73490 : return dz*dz*dz/disty;
532 :
533 19635 : case 34:
534 19635 : return 0.;
535 :
536 0 : default:
537 0 : unsigned int o = 0;
538 0 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
539 0 : unsigned int i2 = i - (o*(o+1)*(o+2)/6);
540 0 : unsigned int block=o, nz = 0;
541 0 : for (; block < i2; block += (o-nz+1)) { nz++; }
542 0 : const unsigned int nx = block - i2;
543 0 : const unsigned int ny = o - nx - nz;
544 0 : Real val = ny;
545 0 : for (unsigned int index=0; index != nx; index++)
546 0 : val *= dx;
547 0 : for (unsigned int index=1; index < ny; index++)
548 0 : val *= dy;
549 0 : for (unsigned int index=0; index != nz; index++)
550 0 : val *= dz;
551 0 : return val/disty;
552 : }
553 : }
554 :
555 :
556 : // d()/dz
557 18862624 : case 2:
558 : {
559 18862624 : switch (i)
560 : {
561 : // constant
562 731085 : case 0:
563 731085 : return 0.;
564 :
565 : // linear
566 731085 : case 1:
567 731085 : return 0.;
568 :
569 731085 : case 2:
570 731085 : return 0.;
571 :
572 2909986 : case 3:
573 2909986 : return 1./distz;
574 :
575 : // quadratic
576 198060 : case 4:
577 198060 : return 0.;
578 :
579 198060 : case 5:
580 198060 : return 0.;
581 :
582 198060 : case 6:
583 198060 : return 0.;
584 :
585 779405 : case 7:
586 779405 : return dx/distz;
587 :
588 779405 : case 8:
589 779405 : return dy/distz;
590 :
591 779405 : case 9:
592 779405 : return 2.*dz/distz;
593 :
594 : // cubic
595 38583 : case 10:
596 38583 : return 0.;
597 :
598 38583 : case 11:
599 38583 : return 0.;
600 :
601 38583 : case 12:
602 38583 : return 0.;
603 :
604 38583 : case 13:
605 38583 : return 0.;
606 :
607 144390 : case 14:
608 144390 : return dx*dx/distz;
609 :
610 144390 : case 15:
611 144390 : return dx*dy/distz;
612 :
613 144390 : case 16:
614 144390 : return dy*dy/distz;
615 :
616 144390 : case 17:
617 144390 : return 2.*dx*dz/distz;
618 :
619 144390 : case 18:
620 144390 : return 2.*dy*dz/distz;
621 :
622 144390 : case 19:
623 144390 : return 3.*dz*dz/distz;
624 :
625 : // quartics
626 19635 : case 20:
627 19635 : return 0.;
628 :
629 19635 : case 21:
630 19635 : return 0.;
631 :
632 19635 : case 22:
633 19635 : return 0.;
634 :
635 19635 : case 23:
636 19635 : return 0.;
637 :
638 19635 : case 24:
639 19635 : return 0.;
640 :
641 73490 : case 25:
642 73490 : return dx*dx*dx/distz;
643 :
644 73490 : case 26:
645 73490 : return dx*dx*dy/distz;
646 :
647 73490 : case 27:
648 73490 : return dx*dy*dy/distz;
649 :
650 73490 : case 28:
651 73490 : return dy*dy*dy/distz;
652 :
653 73490 : case 29:
654 73490 : return 2.*dx*dx*dz/distz;
655 :
656 73490 : case 30:
657 73490 : return 2.*dx*dy*dz/distz;
658 :
659 73490 : case 31:
660 73490 : return 2.*dy*dy*dz/distz;
661 :
662 73490 : case 32:
663 73490 : return 3.*dx*dz*dz/distz;
664 :
665 73490 : case 33:
666 73490 : return 3.*dy*dz*dz/distz;
667 :
668 73490 : case 34:
669 73490 : return 4.*dz*dz*dz/distz;
670 :
671 0 : default:
672 0 : unsigned int o = 0;
673 0 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
674 0 : unsigned int i2 = i - (o*(o+1)*(o+2)/6);
675 0 : unsigned int block=o, nz = 0;
676 0 : for (; block < i2; block += (o-nz+1)) { nz++; }
677 0 : const unsigned int nx = block - i2;
678 0 : const unsigned int ny = o - nx - nz;
679 0 : Real val = nz;
680 0 : for (unsigned int index=0; index != nx; index++)
681 0 : val *= dx;
682 0 : for (unsigned int index=0; index != ny; index++)
683 0 : val *= dy;
684 0 : for (unsigned int index=1; index < nz; index++)
685 0 : val *= dz;
686 0 : return val/distz;
687 : }
688 : }
689 :
690 :
691 0 : default:
692 0 : libmesh_error_msg("Invalid j = " << j);
693 : }
694 :
695 : #else // LIBMESH_DIM != 3
696 : libmesh_assert(true || order || add_p_level);
697 : libmesh_ignore(elem, i, j, point_in);
698 : libmesh_not_implemented();
699 : #endif
700 : }
701 :
702 :
703 : template <>
704 0 : Real FE<3,XYZ>::shape_deriv(const ElemType,
705 : const Order,
706 : const unsigned int,
707 : const unsigned int,
708 : const Point &)
709 : {
710 0 : libmesh_error_msg("XYZ polynomials require the element.");
711 : return 0.;
712 : }
713 :
714 :
715 :
716 : template <>
717 0 : Real FE<3,XYZ>::shape_deriv(const FEType fet,
718 : const Elem * elem,
719 : const unsigned int i,
720 : const unsigned int j,
721 : const Point & p,
722 : const bool add_p_level)
723 : {
724 0 : return FE<3,XYZ>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
725 : }
726 :
727 :
728 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
729 :
730 :
731 : template <>
732 51981504 : Real FE<3,XYZ>::shape_second_deriv(const Elem * elem,
733 : const Order libmesh_dbg_var(order),
734 : const unsigned int i,
735 : const unsigned int j,
736 : const Point & point_in,
737 : const bool libmesh_dbg_var(add_p_level))
738 : {
739 : #if LIBMESH_DIM == 3
740 :
741 13459770 : libmesh_assert(elem);
742 13459770 : libmesh_assert_less (j, 6);
743 :
744 51981504 : Point avg = elem->vertex_average();
745 13459770 : Point max_distance = Point(0.,0.,0.);
746 922474860 : for (const Point & p : elem->node_ref_range())
747 3481973424 : for (unsigned int d = 0; d < 3; d++)
748 : {
749 3288194622 : const Real distance = std::abs(avg(d) - p(d));
750 2611480068 : max_distance(d) = std::max(distance, max_distance(d));
751 : }
752 :
753 51981504 : const Real x = point_in(0);
754 51981504 : const Real y = point_in(1);
755 51981504 : const Real z = point_in(2);
756 51981504 : const Real xc = avg(0);
757 51981504 : const Real yc = avg(1);
758 51981504 : const Real zc = avg(2);
759 51981504 : const Real distx = max_distance(0);
760 51981504 : const Real disty = max_distance(1);
761 51981504 : const Real distz = max_distance(2);
762 51981504 : const Real dx = (x - xc)/distx;
763 51981504 : const Real dy = (y - yc)/disty;
764 51981504 : const Real dz = (z - zc)/distz;
765 51981504 : const Real dist2x = pow(distx,2.);
766 51981504 : const Real dist2y = pow(disty,2.);
767 51981504 : const Real dist2z = pow(distz,2.);
768 51981504 : const Real distxy = distx * disty;
769 51981504 : const Real distxz = distx * distz;
770 51981504 : const Real distyz = disty * distz;
771 :
772 : #ifndef NDEBUG
773 : // totalorder is only used in the assertion below, so
774 : // we avoid declaring it when asserts are not active.
775 13459770 : const unsigned int totalorder = order + add_p_level*elem->p_level();
776 : #endif
777 13459770 : libmesh_assert_less (i, (totalorder+1) * (totalorder+2) *
778 : (totalorder+3)/6);
779 :
780 : // monomials. since they are hierarchic we only need one case block.
781 51981504 : switch (j)
782 : {
783 : // d^2()/dx^2
784 8663584 : case 0:
785 : {
786 8663584 : switch (i)
787 : {
788 : // constant
789 761652 : case 0:
790 :
791 : // linear
792 : case 1:
793 : case 2:
794 : case 3:
795 761652 : return 0.;
796 :
797 : // quadratic
798 521357 : case 4:
799 521357 : return 2./dist2x;
800 :
801 667740 : case 5:
802 : case 6:
803 : case 7:
804 : case 8:
805 : case 9:
806 667740 : return 0.;
807 :
808 : // cubic
809 144390 : case 10:
810 144390 : return 6.*dx/dist2x;
811 :
812 144390 : case 11:
813 144390 : return 2.*dy/dist2x;
814 :
815 77166 : case 12:
816 : case 13:
817 77166 : return 0.;
818 :
819 144390 : case 14:
820 144390 : return 2.*dz/dist2x;
821 :
822 192915 : case 15:
823 : case 16:
824 : case 17:
825 : case 18:
826 : case 19:
827 192915 : return 0.;
828 :
829 : // quartics
830 73490 : case 20:
831 73490 : return 12.*dx*dx/dist2x;
832 :
833 73490 : case 21:
834 73490 : return 6.*dx*dy/dist2x;
835 :
836 73490 : case 22:
837 73490 : return 2.*dy*dy/dist2x;
838 :
839 39270 : case 23:
840 : case 24:
841 39270 : return 0.;
842 :
843 73490 : case 25:
844 73490 : return 6.*dx*dz/dist2x;
845 :
846 73490 : case 26:
847 73490 : return 2.*dy*dz/dist2x;
848 :
849 39270 : case 27:
850 : case 28:
851 39270 : return 0.;
852 :
853 73490 : case 29:
854 73490 : return 2.*dz*dz/dist2x;
855 :
856 98175 : case 30:
857 : case 31:
858 : case 32:
859 : case 33:
860 : case 34:
861 98175 : return 0.;
862 :
863 0 : default:
864 0 : unsigned int o = 0;
865 0 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
866 0 : unsigned int i2 = i - (o*(o+1)*(o+2)/6);
867 0 : unsigned int block=o, nz = 0;
868 0 : for (; block < i2; block += (o-nz+1)) { nz++; }
869 0 : const unsigned int nx = block - i2;
870 0 : const unsigned int ny = o - nx - nz;
871 0 : Real val = nx * (nx - 1);
872 0 : for (unsigned int index=2; index < nx; index++)
873 0 : val *= dx;
874 0 : for (unsigned int index=0; index != ny; index++)
875 0 : val *= dy;
876 0 : for (unsigned int index=0; index != nz; index++)
877 0 : val *= dz;
878 0 : return val/dist2x;
879 : }
880 : }
881 :
882 :
883 : // d^2()/dxdy
884 8663584 : case 1:
885 : {
886 8663584 : switch (i)
887 : {
888 : // constant
889 761652 : case 0:
890 :
891 : // linear
892 : case 1:
893 : case 2:
894 : case 3:
895 761652 : return 0.;
896 :
897 : // quadratic
898 133548 : case 4:
899 133548 : return 0.;
900 :
901 521357 : case 5:
902 521357 : return 1./distxy;
903 :
904 534192 : case 6:
905 : case 7:
906 : case 8:
907 : case 9:
908 534192 : return 0.;
909 :
910 : // cubic
911 38583 : case 10:
912 38583 : return 0.;
913 :
914 144390 : case 11:
915 144390 : return 2.*dx/distxy;
916 :
917 144390 : case 12:
918 144390 : return 2.*dy/distxy;
919 :
920 77166 : case 13:
921 : case 14:
922 77166 : return 0.;
923 :
924 144390 : case 15:
925 144390 : return dz/distxy;
926 :
927 154332 : case 16:
928 : case 17:
929 : case 18:
930 : case 19:
931 154332 : return 0.;
932 :
933 : // quartics
934 19635 : case 20:
935 19635 : return 0.;
936 :
937 73490 : case 21:
938 73490 : return 3.*dx*dx/distxy;
939 :
940 73490 : case 22:
941 73490 : return 4.*dx*dy/distxy;
942 :
943 73490 : case 23:
944 73490 : return 3.*dy*dy/distxy;
945 :
946 39270 : case 24:
947 : case 25:
948 39270 : return 0.;
949 :
950 73490 : case 26:
951 73490 : return 2.*dx*dz/distxy;
952 :
953 73490 : case 27:
954 73490 : return 2.*dy*dz/distxy;
955 :
956 39270 : case 28:
957 : case 29:
958 39270 : return 0.;
959 :
960 73490 : case 30:
961 73490 : return dz*dz/distxy;
962 :
963 78540 : case 31:
964 : case 32:
965 : case 33:
966 : case 34:
967 78540 : return 0.;
968 :
969 0 : default:
970 0 : unsigned int o = 0;
971 0 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
972 0 : unsigned int i2 = i - (o*(o+1)*(o+2)/6);
973 0 : unsigned int block=o, nz = 0;
974 0 : for (; block < i2; block += (o-nz+1)) { nz++; }
975 0 : const unsigned int nx = block - i2;
976 0 : const unsigned int ny = o - nx - nz;
977 0 : Real val = nx * ny;
978 0 : for (unsigned int index=1; index < nx; index++)
979 0 : val *= dx;
980 0 : for (unsigned int index=1; index < ny; index++)
981 0 : val *= dy;
982 0 : for (unsigned int index=0; index != nz; index++)
983 0 : val *= dz;
984 0 : return val/distxy;
985 : }
986 : }
987 :
988 :
989 : // d^2()/dy^2
990 8663584 : case 2:
991 : {
992 8663584 : switch (i)
993 : {
994 : // constant
995 761652 : case 0:
996 :
997 : // linear
998 : case 1:
999 : case 2:
1000 : case 3:
1001 761652 : return 0.;
1002 :
1003 : // quadratic
1004 267096 : case 4:
1005 : case 5:
1006 267096 : return 0.;
1007 :
1008 521357 : case 6:
1009 521357 : return 2./dist2y;
1010 :
1011 400644 : case 7:
1012 : case 8:
1013 : case 9:
1014 400644 : return 0.;
1015 :
1016 : // cubic
1017 77166 : case 10:
1018 : case 11:
1019 77166 : return 0.;
1020 :
1021 144390 : case 12:
1022 144390 : return 2.*dx/dist2y;
1023 144390 : case 13:
1024 144390 : return 6.*dy/dist2y;
1025 :
1026 77166 : case 14:
1027 : case 15:
1028 77166 : return 0.;
1029 :
1030 144390 : case 16:
1031 144390 : return 2.*dz/dist2y;
1032 :
1033 115749 : case 17:
1034 : case 18:
1035 : case 19:
1036 115749 : return 0.;
1037 :
1038 : // quartics
1039 39270 : case 20:
1040 : case 21:
1041 39270 : return 0.;
1042 :
1043 73490 : case 22:
1044 73490 : return 2.*dx*dx/dist2y;
1045 :
1046 73490 : case 23:
1047 73490 : return 6.*dx*dy/dist2y;
1048 :
1049 73490 : case 24:
1050 73490 : return 12.*dy*dy/dist2y;
1051 :
1052 39270 : case 25:
1053 : case 26:
1054 39270 : return 0.;
1055 :
1056 73490 : case 27:
1057 73490 : return 2.*dx*dz/dist2y;
1058 :
1059 73490 : case 28:
1060 73490 : return 6.*dy*dz/dist2y;
1061 :
1062 39270 : case 29:
1063 : case 30:
1064 39270 : return 0.;
1065 :
1066 73490 : case 31:
1067 73490 : return 2.*dz*dz/dist2y;
1068 :
1069 58905 : case 32:
1070 : case 33:
1071 : case 34:
1072 58905 : return 0.;
1073 :
1074 0 : default:
1075 0 : unsigned int o = 0;
1076 0 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1077 0 : unsigned int i2 = i - (o*(o+1)*(o+2)/6);
1078 0 : unsigned int block=o, nz = 0;
1079 0 : for (; block < i2; block += (o-nz+1)) { nz++; }
1080 0 : const unsigned int nx = block - i2;
1081 0 : const unsigned int ny = o - nx - nz;
1082 0 : Real val = ny * (ny - 1);
1083 0 : for (unsigned int index=0; index != nx; index++)
1084 0 : val *= dx;
1085 0 : for (unsigned int index=2; index < ny; index++)
1086 0 : val *= dy;
1087 0 : for (unsigned int index=0; index != nz; index++)
1088 0 : val *= dz;
1089 0 : return val/dist2y;
1090 : }
1091 : }
1092 :
1093 :
1094 : // d^2()/dxdz
1095 8663584 : case 3:
1096 : {
1097 8663584 : switch (i)
1098 : {
1099 : // constant
1100 761652 : case 0:
1101 :
1102 : // linear
1103 : case 1:
1104 : case 2:
1105 : case 3:
1106 761652 : return 0.;
1107 :
1108 : // quadratic
1109 400644 : case 4:
1110 : case 5:
1111 : case 6:
1112 400644 : return 0.;
1113 :
1114 521357 : case 7:
1115 521357 : return 1./distxz;
1116 :
1117 267096 : case 8:
1118 : case 9:
1119 267096 : return 0.;
1120 :
1121 : // cubic
1122 154332 : case 10:
1123 : case 11:
1124 : case 12:
1125 : case 13:
1126 154332 : return 0.;
1127 :
1128 144390 : case 14:
1129 144390 : return 2.*dx/distxz;
1130 :
1131 144390 : case 15:
1132 144390 : return dy/distxz;
1133 :
1134 38583 : case 16:
1135 38583 : return 0.;
1136 :
1137 144390 : case 17:
1138 144390 : return 2.*dz/distxz;
1139 :
1140 77166 : case 18:
1141 : case 19:
1142 77166 : return 0.;
1143 :
1144 : // quartics
1145 98175 : case 20:
1146 : case 21:
1147 : case 22:
1148 : case 23:
1149 : case 24:
1150 98175 : return 0.;
1151 :
1152 73490 : case 25:
1153 73490 : return 3.*dx*dx/distxz;
1154 :
1155 73490 : case 26:
1156 73490 : return 2.*dx*dy/distxz;
1157 :
1158 73490 : case 27:
1159 73490 : return dy*dy/distxz;
1160 :
1161 19635 : case 28:
1162 19635 : return 0.;
1163 :
1164 73490 : case 29:
1165 73490 : return 4.*dx*dz/distxz;
1166 :
1167 73490 : case 30:
1168 73490 : return 2.*dy*dz/distxz;
1169 :
1170 19635 : case 31:
1171 19635 : return 0.;
1172 :
1173 73490 : case 32:
1174 73490 : return 3.*dz*dz/distxz;
1175 :
1176 39270 : case 33:
1177 : case 34:
1178 39270 : return 0.;
1179 :
1180 0 : default:
1181 0 : unsigned int o = 0;
1182 0 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1183 0 : unsigned int i2 = i - (o*(o+1)*(o+2)/6);
1184 0 : unsigned int block=o, nz = 0;
1185 0 : for (; block < i2; block += (o-nz+1)) { nz++; }
1186 0 : const unsigned int nx = block - i2;
1187 0 : const unsigned int ny = o - nx - nz;
1188 0 : Real val = nx * nz;
1189 0 : for (unsigned int index=1; index < nx; index++)
1190 0 : val *= dx;
1191 0 : for (unsigned int index=0; index != ny; index++)
1192 0 : val *= dy;
1193 0 : for (unsigned int index=1; index < nz; index++)
1194 0 : val *= dz;
1195 0 : return val/distxz;
1196 : }
1197 : }
1198 :
1199 : // d^2()/dydz
1200 8663584 : case 4:
1201 : {
1202 8663584 : switch (i)
1203 : {
1204 : // constant
1205 761652 : case 0:
1206 :
1207 : // linear
1208 : case 1:
1209 : case 2:
1210 : case 3:
1211 761652 : return 0.;
1212 :
1213 : // quadratic
1214 534192 : case 4:
1215 : case 5:
1216 : case 6:
1217 : case 7:
1218 534192 : return 0.;
1219 :
1220 521357 : case 8:
1221 521357 : return 1./distyz;
1222 :
1223 133548 : case 9:
1224 133548 : return 0.;
1225 :
1226 : // cubic
1227 192915 : case 10:
1228 : case 11:
1229 : case 12:
1230 : case 13:
1231 : case 14:
1232 192915 : return 0.;
1233 :
1234 144390 : case 15:
1235 144390 : return dx/distyz;
1236 :
1237 144390 : case 16:
1238 144390 : return 2.*dy/distyz;
1239 :
1240 38583 : case 17:
1241 38583 : return 0.;
1242 :
1243 144390 : case 18:
1244 144390 : return 2.*dz/distyz;
1245 :
1246 38583 : case 19:
1247 38583 : return 0.;
1248 :
1249 : // quartics
1250 117810 : case 20:
1251 : case 21:
1252 : case 22:
1253 : case 23:
1254 : case 24:
1255 : case 25:
1256 117810 : return 0.;
1257 :
1258 73490 : case 26:
1259 73490 : return dx*dx/distyz;
1260 :
1261 73490 : case 27:
1262 73490 : return 2.*dx*dy/distyz;
1263 :
1264 73490 : case 28:
1265 73490 : return 3.*dy*dy/distyz;
1266 :
1267 19635 : case 29:
1268 19635 : return 0.;
1269 :
1270 73490 : case 30:
1271 73490 : return 2.*dx*dz/distyz;
1272 :
1273 73490 : case 31:
1274 73490 : return 4.*dy*dz/distyz;
1275 :
1276 19635 : case 32:
1277 19635 : return 0.;
1278 :
1279 73490 : case 33:
1280 73490 : return 3.*dz*dz/distyz;
1281 :
1282 19635 : case 34:
1283 19635 : return 0.;
1284 :
1285 0 : default:
1286 0 : unsigned int o = 0;
1287 0 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1288 0 : unsigned int i2 = i - (o*(o+1)*(o+2)/6);
1289 0 : unsigned int block=o, nz = 0;
1290 0 : for (; block < i2; block += (o-nz+1)) { nz++; }
1291 0 : const unsigned int nx = block - i2;
1292 0 : const unsigned int ny = o - nx - nz;
1293 0 : Real val = ny * nz;
1294 0 : for (unsigned int index=0; index != nx; index++)
1295 0 : val *= dx;
1296 0 : for (unsigned int index=1; index < ny; index++)
1297 0 : val *= dy;
1298 0 : for (unsigned int index=1; index < nz; index++)
1299 0 : val *= dz;
1300 0 : return val/distyz;
1301 : }
1302 : }
1303 :
1304 :
1305 : // d^2()/dz^2
1306 8663584 : case 5:
1307 : {
1308 8663584 : switch (i)
1309 : {
1310 : // constant
1311 761652 : case 0:
1312 :
1313 : // linear
1314 : case 1:
1315 : case 2:
1316 : case 3:
1317 761652 : return 0.;
1318 :
1319 : // quadratic
1320 667740 : case 4:
1321 : case 5:
1322 : case 6:
1323 : case 7:
1324 : case 8:
1325 667740 : return 0.;
1326 :
1327 521357 : case 9:
1328 521357 : return 2./dist2z;
1329 :
1330 : // cubic
1331 270081 : case 10:
1332 : case 11:
1333 : case 12:
1334 : case 13:
1335 : case 14:
1336 : case 15:
1337 : case 16:
1338 270081 : return 0.;
1339 :
1340 144390 : case 17:
1341 144390 : return 2.*dx/dist2z;
1342 :
1343 144390 : case 18:
1344 144390 : return 2.*dy/dist2z;
1345 :
1346 144390 : case 19:
1347 144390 : return 6.*dz/dist2z;
1348 :
1349 : // quartics
1350 176715 : case 20:
1351 : case 21:
1352 : case 22:
1353 : case 23:
1354 : case 24:
1355 : case 25:
1356 : case 26:
1357 : case 27:
1358 : case 28:
1359 176715 : return 0.;
1360 :
1361 73490 : case 29:
1362 73490 : return 2.*dx*dx/dist2z;
1363 :
1364 73490 : case 30:
1365 73490 : return 2.*dx*dy/dist2z;
1366 :
1367 73490 : case 31:
1368 73490 : return 2.*dy*dy/dist2z;
1369 :
1370 73490 : case 32:
1371 73490 : return 6.*dx*dz/dist2z;
1372 :
1373 73490 : case 33:
1374 73490 : return 6.*dy*dz/dist2z;
1375 :
1376 73490 : case 34:
1377 73490 : return 12.*dz*dz/dist2z;
1378 :
1379 0 : default:
1380 0 : unsigned int o = 0;
1381 0 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1382 0 : unsigned int i2 = i - (o*(o+1)*(o+2)/6);
1383 0 : unsigned int block=o, nz = 0;
1384 0 : for (; block < i2; block += (o-nz+1)) { nz++; }
1385 0 : const unsigned int nx = block - i2;
1386 0 : const unsigned int ny = o - nx - nz;
1387 0 : Real val = nz * (nz - 1);
1388 0 : for (unsigned int index=0; index != nx; index++)
1389 0 : val *= dx;
1390 0 : for (unsigned int index=0; index != ny; index++)
1391 0 : val *= dy;
1392 0 : for (unsigned int index=2; index < nz; index++)
1393 0 : val *= dz;
1394 0 : return val/dist2z;
1395 : }
1396 : }
1397 :
1398 :
1399 0 : default:
1400 0 : libmesh_error_msg("Invalid j = " << j);
1401 : }
1402 :
1403 : #else // LIBMESH_DIM != 3
1404 : libmesh_assert(true || order || add_p_level);
1405 : libmesh_ignore(elem, i, j, point_in);
1406 : libmesh_not_implemented();
1407 : #endif
1408 : }
1409 :
1410 :
1411 : template <>
1412 0 : Real FE<3,XYZ>::shape_second_deriv(const ElemType,
1413 : const Order,
1414 : const unsigned int,
1415 : const unsigned int,
1416 : const Point &)
1417 : {
1418 0 : libmesh_error_msg("XYZ polynomials require the element.");
1419 : return 0.;
1420 : }
1421 :
1422 :
1423 : template <>
1424 0 : Real FE<3,XYZ>::shape_second_deriv(const FEType fet,
1425 : const Elem * elem,
1426 : const unsigned int i,
1427 : const unsigned int j,
1428 : const Point & p,
1429 : const bool add_p_level)
1430 : {
1431 0 : return FE<3,XYZ>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
1432 : }
1433 :
1434 :
1435 :
1436 : #endif
1437 :
1438 : } // namespace libMesh
|
