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 12814 : LIBMESH_DEFAULT_VECTORIZED_FE(3,XYZ)
33 :
34 :
35 : template <>
36 85694926 : 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 7059011 : libmesh_assert(elem);
44 :
45 85694926 : Point avg = elem->vertex_average();
46 7059011 : Point max_distance = Point(0.,0.,0.);
47 1380079794 : for (auto p : make_range(elem->n_nodes()))
48 5177539472 : for (unsigned int d = 0; d < 3; d++)
49 : {
50 4525811124 : const Real distance = std::abs(avg(d) - elem->point(p)(d));
51 3883154604 : max_distance(d) = std::max(distance, max_distance(d));
52 : }
53 :
54 85694926 : const Real x = point_in(0);
55 85694926 : const Real y = point_in(1);
56 85694926 : const Real z = point_in(2);
57 85694926 : const Real xc = avg(0);
58 85694926 : const Real yc = avg(1);
59 85694926 : const Real zc = avg(2);
60 85694926 : const Real distx = max_distance(0);
61 85694926 : const Real disty = max_distance(1);
62 85694926 : const Real distz = max_distance(2);
63 85694926 : const Real dx = (x - xc)/distx;
64 85694926 : const Real dy = (y - yc)/disty;
65 85694926 : 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 7059011 : const unsigned int totalorder = order + add_p_level * elem->p_level();
71 : #endif
72 7059011 : 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 85694926 : switch (i)
77 : {
78 : // constant
79 855912 : case 0:
80 855912 : return 1.;
81 :
82 : // linears
83 10320821 : case 1:
84 10320821 : return dx;
85 :
86 10320821 : case 2:
87 10320821 : return dy;
88 :
89 10320821 : case 3:
90 10320821 : return dz;
91 :
92 : // quadratics
93 3612752 : case 4:
94 3612752 : return dx*dx;
95 :
96 3612752 : case 5:
97 3612752 : return dx*dy;
98 :
99 3612752 : case 6:
100 3612752 : return dy*dy;
101 :
102 3612752 : case 7:
103 3612752 : return dx*dz;
104 :
105 3612752 : case 8:
106 3612752 : return dz*dy;
107 :
108 3612752 : case 9:
109 3612752 : return dz*dz;
110 :
111 : // cubics
112 1305638 : case 10:
113 1305638 : return dx*dx*dx;
114 :
115 1305638 : case 11:
116 1305638 : return dx*dx*dy;
117 :
118 1305638 : case 12:
119 1305638 : return dx*dy*dy;
120 :
121 1305638 : case 13:
122 1305638 : return dy*dy*dy;
123 :
124 1305638 : case 14:
125 1305638 : return dx*dx*dz;
126 :
127 1305638 : case 15:
128 1305638 : return dx*dy*dz;
129 :
130 1305638 : case 16:
131 1305638 : return dy*dy*dz;
132 :
133 1305638 : case 17:
134 1305638 : return dx*dz*dz;
135 :
136 1305638 : case 18:
137 1305638 : return dy*dz*dz;
138 :
139 1305638 : case 19:
140 1305638 : return dz*dz*dz;
141 :
142 : // quartics
143 645250 : case 20:
144 645250 : return dx*dx*dx*dx;
145 :
146 645250 : case 21:
147 645250 : return dx*dx*dx*dy;
148 :
149 645250 : case 22:
150 645250 : return dx*dx*dy*dy;
151 :
152 645250 : case 23:
153 645250 : return dx*dy*dy*dy;
154 :
155 645250 : case 24:
156 645250 : return dy*dy*dy*dy;
157 :
158 645250 : case 25:
159 645250 : return dx*dx*dx*dz;
160 :
161 645250 : case 26:
162 645250 : return dx*dx*dy*dz;
163 :
164 645250 : case 27:
165 645250 : return dx*dy*dy*dz;
166 :
167 645250 : case 28:
168 645250 : return dy*dy*dy*dz;
169 :
170 645250 : case 29:
171 645250 : return dx*dx*dz*dz;
172 :
173 645250 : case 30:
174 645250 : return dx*dy*dz*dz;
175 :
176 645250 : case 31:
177 645250 : return dy*dy*dz*dz;
178 :
179 645250 : case 32:
180 645250 : return dx*dz*dz*dz;
181 :
182 645250 : case 33:
183 645250 : return dy*dz*dz*dz;
184 :
185 645250 : case 34:
186 645250 : 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 180415116 : 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 14829741 : libmesh_assert(elem);
252 14829741 : libmesh_assert_less (j, 3);
253 :
254 180415116 : Point avg = elem->vertex_average();
255 14829741 : Point max_distance = Point(0.,0.,0.);
256 2941872108 : for (auto p : make_range(elem->n_nodes()))
257 11045827968 : for (unsigned int d = 0; d < 3; d++)
258 : {
259 9654563034 : const Real distance = std::abs(avg(d) - elem->point(p)(d));
260 8284370976 : max_distance(d) = std::max(distance, max_distance(d));
261 : }
262 :
263 180415116 : const Real x = point_in(0);
264 180415116 : const Real y = point_in(1);
265 180415116 : const Real z = point_in(2);
266 180415116 : const Real xc = avg(0);
267 180415116 : const Real yc = avg(1);
268 180415116 : const Real zc = avg(2);
269 180415116 : const Real distx = max_distance(0);
270 180415116 : const Real disty = max_distance(1);
271 180415116 : const Real distz = max_distance(2);
272 180415116 : const Real dx = (x - xc)/distx;
273 180415116 : const Real dy = (y - yc)/disty;
274 180415116 : 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 14829741 : const unsigned int totalorder = order + add_p_level*elem->p_level();
280 : #endif
281 14829741 : libmesh_assert_less (i, (totalorder+1) * (totalorder+2) *
282 : (totalorder+3)/6);
283 :
284 180415116 : switch (j)
285 : {
286 : // d()/dx
287 60138372 : case 0:
288 : {
289 60138372 : switch (i)
290 : {
291 : // constant
292 738621 : case 0:
293 738621 : return 0.;
294 :
295 : // linear
296 8906421 : case 1:
297 8906421 : return 1./distx;
298 :
299 738621 : case 2:
300 738621 : return 0.;
301 :
302 738621 : case 3:
303 738621 : return 0.;
304 :
305 : // quadratic
306 2490998 : case 4:
307 2490998 : return 2.*dx/distx;
308 :
309 2490998 : case 5:
310 2490998 : return dy/distx;
311 :
312 204588 : case 6:
313 204588 : return 0.;
314 :
315 2490998 : case 7:
316 2490998 : return dz/distx;
317 :
318 204588 : case 8:
319 204588 : return 0.;
320 :
321 204588 : case 9:
322 204588 : return 0.;
323 :
324 : // cubic
325 549840 : case 10:
326 549840 : return 3.*dx*dx/distx;
327 :
328 549840 : case 11:
329 549840 : return 2.*dx*dy/distx;
330 :
331 549840 : case 12:
332 549840 : return dy*dy/distx;
333 :
334 43287 : case 13:
335 43287 : return 0.;
336 :
337 549840 : case 14:
338 549840 : return 2.*dx*dz/distx;
339 :
340 549840 : case 15:
341 549840 : return dy*dz/distx;
342 :
343 43287 : case 16:
344 43287 : return 0.;
345 :
346 549840 : case 17:
347 549840 : return dz*dz/distx;
348 :
349 43287 : case 18:
350 43287 : return 0.;
351 :
352 43287 : case 19:
353 43287 : return 0.;
354 :
355 : // quartics
356 271220 : case 20:
357 271220 : return 4.*dx*dx*dx/distx;
358 :
359 271220 : case 21:
360 271220 : return 3.*dx*dx*dy/distx;
361 :
362 271220 : case 22:
363 271220 : return 2.*dx*dy*dy/distx;
364 :
365 271220 : case 23:
366 271220 : return dy*dy*dy/distx;
367 :
368 21891 : case 24:
369 21891 : return 0.;
370 :
371 271220 : case 25:
372 271220 : return 3.*dx*dx*dz/distx;
373 :
374 271220 : case 26:
375 271220 : return 2.*dx*dy*dz/distx;
376 :
377 271220 : case 27:
378 271220 : return dy*dy*dz/distx;
379 :
380 21891 : case 28:
381 21891 : return 0.;
382 :
383 271220 : case 29:
384 271220 : return 2.*dx*dz*dz/distx;
385 :
386 271220 : case 30:
387 271220 : return dy*dz*dz/distx;
388 :
389 21891 : case 31:
390 21891 : return 0.;
391 :
392 271220 : case 32:
393 271220 : return dz*dz*dz/distx;
394 :
395 21891 : case 33:
396 21891 : return 0.;
397 :
398 21891 : case 34:
399 21891 : 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 60138372 : case 1:
423 : {
424 60138372 : switch (i)
425 : {
426 : // constant
427 738621 : case 0:
428 738621 : return 0.;
429 :
430 : // linear
431 738621 : case 1:
432 738621 : return 0.;
433 :
434 8906421 : case 2:
435 8906421 : return 1./disty;
436 :
437 738621 : case 3:
438 738621 : return 0.;
439 :
440 : // quadratic
441 204588 : case 4:
442 204588 : return 0.;
443 :
444 2490998 : case 5:
445 2490998 : return dx/disty;
446 :
447 2490998 : case 6:
448 2490998 : return 2.*dy/disty;
449 :
450 204588 : case 7:
451 204588 : return 0.;
452 :
453 2490998 : case 8:
454 2490998 : return dz/disty;
455 :
456 204588 : case 9:
457 204588 : return 0.;
458 :
459 : // cubic
460 43287 : case 10:
461 43287 : return 0.;
462 :
463 549840 : case 11:
464 549840 : return dx*dx/disty;
465 :
466 549840 : case 12:
467 549840 : return 2.*dx*dy/disty;
468 :
469 549840 : case 13:
470 549840 : return 3.*dy*dy/disty;
471 :
472 43287 : case 14:
473 43287 : return 0.;
474 :
475 549840 : case 15:
476 549840 : return dx*dz/disty;
477 :
478 549840 : case 16:
479 549840 : return 2.*dy*dz/disty;
480 :
481 43287 : case 17:
482 43287 : return 0.;
483 :
484 549840 : case 18:
485 549840 : return dz*dz/disty;
486 :
487 43287 : case 19:
488 43287 : return 0.;
489 :
490 : // quartics
491 21891 : case 20:
492 21891 : return 0.;
493 :
494 271220 : case 21:
495 271220 : return dx*dx*dx/disty;
496 :
497 271220 : case 22:
498 271220 : return 2.*dx*dx*dy/disty;
499 :
500 271220 : case 23:
501 271220 : return 3.*dx*dy*dy/disty;
502 :
503 271220 : case 24:
504 271220 : return 4.*dy*dy*dy/disty;
505 :
506 21891 : case 25:
507 21891 : return 0.;
508 :
509 271220 : case 26:
510 271220 : return dx*dx*dz/disty;
511 :
512 271220 : case 27:
513 271220 : return 2.*dx*dy*dz/disty;
514 :
515 271220 : case 28:
516 271220 : return 3.*dy*dy*dz/disty;
517 :
518 21891 : case 29:
519 21891 : return 0.;
520 :
521 271220 : case 30:
522 271220 : return dx*dz*dz/disty;
523 :
524 271220 : case 31:
525 271220 : return 2.*dy*dz*dz/disty;
526 :
527 21891 : case 32:
528 21891 : return 0.;
529 :
530 271220 : case 33:
531 271220 : return dz*dz*dz/disty;
532 :
533 21891 : case 34:
534 21891 : 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 60138372 : case 2:
558 : {
559 60138372 : switch (i)
560 : {
561 : // constant
562 738621 : case 0:
563 738621 : return 0.;
564 :
565 : // linear
566 738621 : case 1:
567 738621 : return 0.;
568 :
569 738621 : case 2:
570 738621 : return 0.;
571 :
572 8906421 : case 3:
573 8906421 : return 1./distz;
574 :
575 : // quadratic
576 204588 : case 4:
577 204588 : return 0.;
578 :
579 204588 : case 5:
580 204588 : return 0.;
581 :
582 204588 : case 6:
583 204588 : return 0.;
584 :
585 2490998 : case 7:
586 2490998 : return dx/distz;
587 :
588 2490998 : case 8:
589 2490998 : return dy/distz;
590 :
591 2490998 : case 9:
592 2490998 : return 2.*dz/distz;
593 :
594 : // cubic
595 43287 : case 10:
596 43287 : return 0.;
597 :
598 43287 : case 11:
599 43287 : return 0.;
600 :
601 43287 : case 12:
602 43287 : return 0.;
603 :
604 43287 : case 13:
605 43287 : return 0.;
606 :
607 549840 : case 14:
608 549840 : return dx*dx/distz;
609 :
610 549840 : case 15:
611 549840 : return dx*dy/distz;
612 :
613 549840 : case 16:
614 549840 : return dy*dy/distz;
615 :
616 549840 : case 17:
617 549840 : return 2.*dx*dz/distz;
618 :
619 549840 : case 18:
620 549840 : return 2.*dy*dz/distz;
621 :
622 549840 : case 19:
623 549840 : return 3.*dz*dz/distz;
624 :
625 : // quartics
626 21891 : case 20:
627 21891 : return 0.;
628 :
629 21891 : case 21:
630 21891 : return 0.;
631 :
632 21891 : case 22:
633 21891 : return 0.;
634 :
635 21891 : case 23:
636 21891 : return 0.;
637 :
638 21891 : case 24:
639 21891 : return 0.;
640 :
641 271220 : case 25:
642 271220 : return dx*dx*dx/distz;
643 :
644 271220 : case 26:
645 271220 : return dx*dx*dy/distz;
646 :
647 271220 : case 27:
648 271220 : return dx*dy*dy/distz;
649 :
650 271220 : case 28:
651 271220 : return dy*dy*dy/distz;
652 :
653 271220 : case 29:
654 271220 : return 2.*dx*dx*dz/distz;
655 :
656 271220 : case 30:
657 271220 : return 2.*dx*dy*dz/distz;
658 :
659 271220 : case 31:
660 271220 : return 2.*dy*dy*dz/distz;
661 :
662 271220 : case 32:
663 271220 : return 3.*dx*dz*dz/distz;
664 :
665 271220 : case 33:
666 271220 : return 3.*dy*dz*dz/distz;
667 :
668 271220 : case 34:
669 271220 : 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 177247512 : 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 14360922 : libmesh_assert(elem);
742 14360922 : libmesh_assert_less (j, 6);
743 :
744 177247512 : Point avg = elem->vertex_average();
745 14360922 : Point max_distance = Point(0.,0.,0.);
746 3130003416 : for (const Point & p : elem->node_ref_range())
747 11811023616 : for (unsigned int d = 0; d < 3; d++)
748 : {
749 9585920250 : const Real distance = std::abs(avg(d) - p(d));
750 8858267712 : max_distance(d) = std::max(distance, max_distance(d));
751 : }
752 :
753 177247512 : const Real x = point_in(0);
754 177247512 : const Real y = point_in(1);
755 177247512 : const Real z = point_in(2);
756 177247512 : const Real xc = avg(0);
757 177247512 : const Real yc = avg(1);
758 177247512 : const Real zc = avg(2);
759 177247512 : const Real distx = max_distance(0);
760 177247512 : const Real disty = max_distance(1);
761 177247512 : const Real distz = max_distance(2);
762 177247512 : const Real dx = (x - xc)/distx;
763 177247512 : const Real dy = (y - yc)/disty;
764 177247512 : const Real dz = (z - zc)/distz;
765 177247512 : const Real dist2x = pow(distx,2.);
766 177247512 : const Real dist2y = pow(disty,2.);
767 177247512 : const Real dist2z = pow(distz,2.);
768 177247512 : const Real distxy = distx * disty;
769 177247512 : const Real distxz = distx * distz;
770 177247512 : 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 14360922 : const unsigned int totalorder = order + add_p_level*elem->p_level();
776 : #endif
777 14360922 : 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 177247512 : switch (j)
782 : {
783 : // d^2()/dx^2
784 29541252 : case 0:
785 : {
786 29541252 : switch (i)
787 : {
788 : // constant
789 791796 : case 0:
790 :
791 : // linear
792 : case 1:
793 : case 2:
794 : case 3:
795 791796 : return 0.;
796 :
797 : // quadratic
798 1716854 : case 4:
799 1716854 : return 2./dist2x;
800 :
801 700380 : case 5:
802 : case 6:
803 : case 7:
804 : case 8:
805 : case 9:
806 700380 : return 0.;
807 :
808 : // cubic
809 549840 : case 10:
810 549840 : return 6.*dx/dist2x;
811 :
812 549840 : case 11:
813 549840 : return 2.*dy/dist2x;
814 :
815 86574 : case 12:
816 : case 13:
817 86574 : return 0.;
818 :
819 549840 : case 14:
820 549840 : return 2.*dz/dist2x;
821 :
822 216435 : case 15:
823 : case 16:
824 : case 17:
825 : case 18:
826 : case 19:
827 216435 : return 0.;
828 :
829 : // quartics
830 271220 : case 20:
831 271220 : return 12.*dx*dx/dist2x;
832 :
833 271220 : case 21:
834 271220 : return 6.*dx*dy/dist2x;
835 :
836 271220 : case 22:
837 271220 : return 2.*dy*dy/dist2x;
838 :
839 43782 : case 23:
840 : case 24:
841 43782 : return 0.;
842 :
843 271220 : case 25:
844 271220 : return 6.*dx*dz/dist2x;
845 :
846 271220 : case 26:
847 271220 : return 2.*dy*dz/dist2x;
848 :
849 43782 : case 27:
850 : case 28:
851 43782 : return 0.;
852 :
853 271220 : case 29:
854 271220 : return 2.*dz*dz/dist2x;
855 :
856 109455 : case 30:
857 : case 31:
858 : case 32:
859 : case 33:
860 : case 34:
861 109455 : 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 29541252 : case 1:
885 : {
886 29541252 : switch (i)
887 : {
888 : // constant
889 791796 : case 0:
890 :
891 : // linear
892 : case 1:
893 : case 2:
894 : case 3:
895 791796 : return 0.;
896 :
897 : // quadratic
898 140076 : case 4:
899 140076 : return 0.;
900 :
901 1716854 : case 5:
902 1716854 : return 1./distxy;
903 :
904 560304 : case 6:
905 : case 7:
906 : case 8:
907 : case 9:
908 560304 : return 0.;
909 :
910 : // cubic
911 43287 : case 10:
912 43287 : return 0.;
913 :
914 549840 : case 11:
915 549840 : return 2.*dx/distxy;
916 :
917 549840 : case 12:
918 549840 : return 2.*dy/distxy;
919 :
920 86574 : case 13:
921 : case 14:
922 86574 : return 0.;
923 :
924 549840 : case 15:
925 549840 : return dz/distxy;
926 :
927 173148 : case 16:
928 : case 17:
929 : case 18:
930 : case 19:
931 173148 : return 0.;
932 :
933 : // quartics
934 21891 : case 20:
935 21891 : return 0.;
936 :
937 271220 : case 21:
938 271220 : return 3.*dx*dx/distxy;
939 :
940 271220 : case 22:
941 271220 : return 4.*dx*dy/distxy;
942 :
943 271220 : case 23:
944 271220 : return 3.*dy*dy/distxy;
945 :
946 43782 : case 24:
947 : case 25:
948 43782 : return 0.;
949 :
950 271220 : case 26:
951 271220 : return 2.*dx*dz/distxy;
952 :
953 271220 : case 27:
954 271220 : return 2.*dy*dz/distxy;
955 :
956 43782 : case 28:
957 : case 29:
958 43782 : return 0.;
959 :
960 271220 : case 30:
961 271220 : return dz*dz/distxy;
962 :
963 87564 : case 31:
964 : case 32:
965 : case 33:
966 : case 34:
967 87564 : 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 29541252 : case 2:
991 : {
992 29541252 : switch (i)
993 : {
994 : // constant
995 791796 : case 0:
996 :
997 : // linear
998 : case 1:
999 : case 2:
1000 : case 3:
1001 791796 : return 0.;
1002 :
1003 : // quadratic
1004 280152 : case 4:
1005 : case 5:
1006 280152 : return 0.;
1007 :
1008 1716854 : case 6:
1009 1716854 : return 2./dist2y;
1010 :
1011 420228 : case 7:
1012 : case 8:
1013 : case 9:
1014 420228 : return 0.;
1015 :
1016 : // cubic
1017 86574 : case 10:
1018 : case 11:
1019 86574 : return 0.;
1020 :
1021 549840 : case 12:
1022 549840 : return 2.*dx/dist2y;
1023 549840 : case 13:
1024 549840 : return 6.*dy/dist2y;
1025 :
1026 86574 : case 14:
1027 : case 15:
1028 86574 : return 0.;
1029 :
1030 549840 : case 16:
1031 549840 : return 2.*dz/dist2y;
1032 :
1033 129861 : case 17:
1034 : case 18:
1035 : case 19:
1036 129861 : return 0.;
1037 :
1038 : // quartics
1039 43782 : case 20:
1040 : case 21:
1041 43782 : return 0.;
1042 :
1043 271220 : case 22:
1044 271220 : return 2.*dx*dx/dist2y;
1045 :
1046 271220 : case 23:
1047 271220 : return 6.*dx*dy/dist2y;
1048 :
1049 271220 : case 24:
1050 271220 : return 12.*dy*dy/dist2y;
1051 :
1052 43782 : case 25:
1053 : case 26:
1054 43782 : return 0.;
1055 :
1056 271220 : case 27:
1057 271220 : return 2.*dx*dz/dist2y;
1058 :
1059 271220 : case 28:
1060 271220 : return 6.*dy*dz/dist2y;
1061 :
1062 43782 : case 29:
1063 : case 30:
1064 43782 : return 0.;
1065 :
1066 271220 : case 31:
1067 271220 : return 2.*dz*dz/dist2y;
1068 :
1069 65673 : case 32:
1070 : case 33:
1071 : case 34:
1072 65673 : 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 29541252 : case 3:
1096 : {
1097 29541252 : switch (i)
1098 : {
1099 : // constant
1100 791796 : case 0:
1101 :
1102 : // linear
1103 : case 1:
1104 : case 2:
1105 : case 3:
1106 791796 : return 0.;
1107 :
1108 : // quadratic
1109 420228 : case 4:
1110 : case 5:
1111 : case 6:
1112 420228 : return 0.;
1113 :
1114 1716854 : case 7:
1115 1716854 : return 1./distxz;
1116 :
1117 280152 : case 8:
1118 : case 9:
1119 280152 : return 0.;
1120 :
1121 : // cubic
1122 173148 : case 10:
1123 : case 11:
1124 : case 12:
1125 : case 13:
1126 173148 : return 0.;
1127 :
1128 549840 : case 14:
1129 549840 : return 2.*dx/distxz;
1130 :
1131 549840 : case 15:
1132 549840 : return dy/distxz;
1133 :
1134 43287 : case 16:
1135 43287 : return 0.;
1136 :
1137 549840 : case 17:
1138 549840 : return 2.*dz/distxz;
1139 :
1140 86574 : case 18:
1141 : case 19:
1142 86574 : return 0.;
1143 :
1144 : // quartics
1145 109455 : case 20:
1146 : case 21:
1147 : case 22:
1148 : case 23:
1149 : case 24:
1150 109455 : return 0.;
1151 :
1152 271220 : case 25:
1153 271220 : return 3.*dx*dx/distxz;
1154 :
1155 271220 : case 26:
1156 271220 : return 2.*dx*dy/distxz;
1157 :
1158 271220 : case 27:
1159 271220 : return dy*dy/distxz;
1160 :
1161 21891 : case 28:
1162 21891 : return 0.;
1163 :
1164 271220 : case 29:
1165 271220 : return 4.*dx*dz/distxz;
1166 :
1167 271220 : case 30:
1168 271220 : return 2.*dy*dz/distxz;
1169 :
1170 21891 : case 31:
1171 21891 : return 0.;
1172 :
1173 271220 : case 32:
1174 271220 : return 3.*dz*dz/distxz;
1175 :
1176 43782 : case 33:
1177 : case 34:
1178 43782 : 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 29541252 : case 4:
1201 : {
1202 29541252 : switch (i)
1203 : {
1204 : // constant
1205 791796 : case 0:
1206 :
1207 : // linear
1208 : case 1:
1209 : case 2:
1210 : case 3:
1211 791796 : return 0.;
1212 :
1213 : // quadratic
1214 560304 : case 4:
1215 : case 5:
1216 : case 6:
1217 : case 7:
1218 560304 : return 0.;
1219 :
1220 1716854 : case 8:
1221 1716854 : return 1./distyz;
1222 :
1223 140076 : case 9:
1224 140076 : return 0.;
1225 :
1226 : // cubic
1227 216435 : case 10:
1228 : case 11:
1229 : case 12:
1230 : case 13:
1231 : case 14:
1232 216435 : return 0.;
1233 :
1234 549840 : case 15:
1235 549840 : return dx/distyz;
1236 :
1237 549840 : case 16:
1238 549840 : return 2.*dy/distyz;
1239 :
1240 43287 : case 17:
1241 43287 : return 0.;
1242 :
1243 549840 : case 18:
1244 549840 : return 2.*dz/distyz;
1245 :
1246 43287 : case 19:
1247 43287 : return 0.;
1248 :
1249 : // quartics
1250 131346 : case 20:
1251 : case 21:
1252 : case 22:
1253 : case 23:
1254 : case 24:
1255 : case 25:
1256 131346 : return 0.;
1257 :
1258 271220 : case 26:
1259 271220 : return dx*dx/distyz;
1260 :
1261 271220 : case 27:
1262 271220 : return 2.*dx*dy/distyz;
1263 :
1264 271220 : case 28:
1265 271220 : return 3.*dy*dy/distyz;
1266 :
1267 21891 : case 29:
1268 21891 : return 0.;
1269 :
1270 271220 : case 30:
1271 271220 : return 2.*dx*dz/distyz;
1272 :
1273 271220 : case 31:
1274 271220 : return 4.*dy*dz/distyz;
1275 :
1276 21891 : case 32:
1277 21891 : return 0.;
1278 :
1279 271220 : case 33:
1280 271220 : return 3.*dz*dz/distyz;
1281 :
1282 21891 : case 34:
1283 21891 : 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 29541252 : case 5:
1307 : {
1308 29541252 : switch (i)
1309 : {
1310 : // constant
1311 791796 : case 0:
1312 :
1313 : // linear
1314 : case 1:
1315 : case 2:
1316 : case 3:
1317 791796 : return 0.;
1318 :
1319 : // quadratic
1320 700380 : case 4:
1321 : case 5:
1322 : case 6:
1323 : case 7:
1324 : case 8:
1325 700380 : return 0.;
1326 :
1327 1716854 : case 9:
1328 1716854 : return 2./dist2z;
1329 :
1330 : // cubic
1331 303009 : case 10:
1332 : case 11:
1333 : case 12:
1334 : case 13:
1335 : case 14:
1336 : case 15:
1337 : case 16:
1338 303009 : return 0.;
1339 :
1340 549840 : case 17:
1341 549840 : return 2.*dx/dist2z;
1342 :
1343 549840 : case 18:
1344 549840 : return 2.*dy/dist2z;
1345 :
1346 549840 : case 19:
1347 549840 : return 6.*dz/dist2z;
1348 :
1349 : // quartics
1350 197019 : 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 197019 : return 0.;
1360 :
1361 271220 : case 29:
1362 271220 : return 2.*dx*dx/dist2z;
1363 :
1364 271220 : case 30:
1365 271220 : return 2.*dx*dy/dist2z;
1366 :
1367 271220 : case 31:
1368 271220 : return 2.*dy*dy/dist2z;
1369 :
1370 271220 : case 32:
1371 271220 : return 6.*dx*dz/dist2z;
1372 :
1373 271220 : case 33:
1374 271220 : return 6.*dy*dz/dist2z;
1375 :
1376 271220 : case 34:
1377 271220 : 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
|
