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 28447587 : 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 7206067 : libmesh_assert(elem);
44 :
45 28447587 : Point avg = elem->vertex_average();
46 7206067 : Point max_distance = Point(0.,0.,0.);
47 459094569 : for (auto p : make_range(elem->n_nodes()))
48 1722587928 : for (unsigned int d = 0; d < 3; d++)
49 : {
50 1947536796 : const Real distance = std::abs(avg(d) - elem->point(p)(d));
51 1291940946 : max_distance(d) = std::max(distance, max_distance(d));
52 : }
53 :
54 28447587 : const Real x = point_in(0);
55 28447587 : const Real y = point_in(1);
56 28447587 : const Real z = point_in(2);
57 28447587 : const Real xc = avg(0);
58 28447587 : const Real yc = avg(1);
59 28447587 : const Real zc = avg(2);
60 28447587 : const Real distx = max_distance(0);
61 28447587 : const Real disty = max_distance(1);
62 28447587 : const Real distz = max_distance(2);
63 28447587 : const Real dx = (x - xc)/distx;
64 28447587 : const Real dy = (y - yc)/disty;
65 28447587 : 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 7206067 : const unsigned int totalorder = order + add_p_level * elem->p_level();
71 : #endif
72 7206067 : 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 28447587 : switch (i)
77 : {
78 : // constant
79 861929 : case 0:
80 861929 : return 1.;
81 :
82 : // linears
83 3430280 : case 1:
84 3430280 : return dx;
85 :
86 3430280 : case 2:
87 3430280 : return dy;
88 :
89 3430280 : case 3:
90 3430280 : return dz;
91 :
92 : // quadratics
93 1197592 : case 4:
94 1197592 : return dx*dx;
95 :
96 1197592 : case 5:
97 1197592 : return dx*dy;
98 :
99 1197592 : case 6:
100 1197592 : return dy*dy;
101 :
102 1197592 : case 7:
103 1197592 : return dx*dz;
104 :
105 1197592 : case 8:
106 1197592 : return dz*dy;
107 :
108 1197592 : case 9:
109 1197592 : return dz*dz;
110 :
111 : // cubics
112 431209 : case 10:
113 431209 : return dx*dx*dx;
114 :
115 431209 : case 11:
116 431209 : return dx*dx*dy;
117 :
118 431209 : case 12:
119 431209 : return dx*dy*dy;
120 :
121 431209 : case 13:
122 431209 : return dy*dy*dy;
123 :
124 431209 : case 14:
125 431209 : return dx*dx*dz;
126 :
127 431209 : case 15:
128 431209 : return dx*dy*dz;
129 :
130 431209 : case 16:
131 431209 : return dy*dy*dz;
132 :
133 431209 : case 17:
134 431209 : return dx*dz*dz;
135 :
136 431209 : case 18:
137 431209 : return dy*dz*dz;
138 :
139 431209 : case 19:
140 431209 : return dz*dz*dz;
141 :
142 : // quartics
143 215255 : case 20:
144 215255 : return dx*dx*dx*dx;
145 :
146 215255 : case 21:
147 215255 : return dx*dx*dx*dy;
148 :
149 215255 : case 22:
150 215255 : return dx*dx*dy*dy;
151 :
152 215255 : case 23:
153 215255 : return dx*dy*dy*dy;
154 :
155 215255 : case 24:
156 215255 : return dy*dy*dy*dy;
157 :
158 215255 : case 25:
159 215255 : return dx*dx*dx*dz;
160 :
161 215255 : case 26:
162 215255 : return dx*dx*dy*dz;
163 :
164 215255 : case 27:
165 215255 : return dx*dy*dy*dz;
166 :
167 215255 : case 28:
168 215255 : return dy*dy*dy*dz;
169 :
170 215255 : case 29:
171 215255 : return dx*dx*dz*dz;
172 :
173 215255 : case 30:
174 215255 : return dx*dy*dz*dz;
175 :
176 215255 : case 31:
177 215255 : return dy*dy*dz*dz;
178 :
179 215255 : case 32:
180 215255 : return dx*dz*dz*dz;
181 :
182 215255 : case 33:
183 215255 : return dy*dz*dz*dz;
184 :
185 215255 : case 34:
186 215255 : 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 58242072 : 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 58242072 : Point avg = elem->vertex_average();
255 14829741 : Point max_distance = Point(0.,0.,0.);
256 953340366 : for (auto p : make_range(elem->n_nodes()))
257 3580393176 : for (unsigned int d = 0; d < 3; d++)
258 : {
259 4055486940 : const Real distance = std::abs(avg(d) - elem->point(p)(d));
260 2685294882 : max_distance(d) = std::max(distance, max_distance(d));
261 : }
262 :
263 58242072 : const Real x = point_in(0);
264 58242072 : const Real y = point_in(1);
265 58242072 : const Real z = point_in(2);
266 58242072 : const Real xc = avg(0);
267 58242072 : const Real yc = avg(1);
268 58242072 : const Real zc = avg(2);
269 58242072 : const Real distx = max_distance(0);
270 58242072 : const Real disty = max_distance(1);
271 58242072 : const Real distz = max_distance(2);
272 58242072 : const Real dx = (x - xc)/distx;
273 58242072 : const Real dy = (y - yc)/disty;
274 58242072 : 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 58242072 : switch (j)
285 : {
286 : // d()/dx
287 19414024 : case 0:
288 : {
289 19414024 : switch (i)
290 : {
291 : // constant
292 738621 : case 0:
293 738621 : return 0.;
294 :
295 : // linear
296 2937706 : case 1:
297 2937706 : 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 803405 : case 4:
307 803405 : return 2.*dx/distx;
308 :
309 803405 : case 5:
310 803405 : return dy/distx;
311 :
312 204588 : case 6:
313 204588 : return 0.;
314 :
315 803405 : case 7:
316 803405 : 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 161622 : case 10:
326 161622 : return 3.*dx*dx/distx;
327 :
328 161622 : case 11:
329 161622 : return 2.*dx*dy/distx;
330 :
331 161622 : case 12:
332 161622 : return dy*dy/distx;
333 :
334 43287 : case 13:
335 43287 : return 0.;
336 :
337 161622 : case 14:
338 161622 : return 2.*dx*dz/distx;
339 :
340 161622 : case 15:
341 161622 : return dy*dz/distx;
342 :
343 43287 : case 16:
344 43287 : return 0.;
345 :
346 161622 : case 17:
347 161622 : 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 81770 : case 20:
357 81770 : return 4.*dx*dx*dx/distx;
358 :
359 81770 : case 21:
360 81770 : return 3.*dx*dx*dy/distx;
361 :
362 81770 : case 22:
363 81770 : return 2.*dx*dy*dy/distx;
364 :
365 81770 : case 23:
366 81770 : return dy*dy*dy/distx;
367 :
368 21891 : case 24:
369 21891 : return 0.;
370 :
371 81770 : case 25:
372 81770 : return 3.*dx*dx*dz/distx;
373 :
374 81770 : case 26:
375 81770 : return 2.*dx*dy*dz/distx;
376 :
377 81770 : case 27:
378 81770 : return dy*dy*dz/distx;
379 :
380 21891 : case 28:
381 21891 : return 0.;
382 :
383 81770 : case 29:
384 81770 : return 2.*dx*dz*dz/distx;
385 :
386 81770 : case 30:
387 81770 : return dy*dz*dz/distx;
388 :
389 21891 : case 31:
390 21891 : return 0.;
391 :
392 81770 : case 32:
393 81770 : 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 19414024 : case 1:
423 : {
424 19414024 : 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 2937706 : case 2:
435 2937706 : 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 803405 : case 5:
445 803405 : return dx/disty;
446 :
447 803405 : case 6:
448 803405 : return 2.*dy/disty;
449 :
450 204588 : case 7:
451 204588 : return 0.;
452 :
453 803405 : case 8:
454 803405 : 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 161622 : case 11:
464 161622 : return dx*dx/disty;
465 :
466 161622 : case 12:
467 161622 : return 2.*dx*dy/disty;
468 :
469 161622 : case 13:
470 161622 : return 3.*dy*dy/disty;
471 :
472 43287 : case 14:
473 43287 : return 0.;
474 :
475 161622 : case 15:
476 161622 : return dx*dz/disty;
477 :
478 161622 : case 16:
479 161622 : return 2.*dy*dz/disty;
480 :
481 43287 : case 17:
482 43287 : return 0.;
483 :
484 161622 : case 18:
485 161622 : 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 81770 : case 21:
495 81770 : return dx*dx*dx/disty;
496 :
497 81770 : case 22:
498 81770 : return 2.*dx*dx*dy/disty;
499 :
500 81770 : case 23:
501 81770 : return 3.*dx*dy*dy/disty;
502 :
503 81770 : case 24:
504 81770 : return 4.*dy*dy*dy/disty;
505 :
506 21891 : case 25:
507 21891 : return 0.;
508 :
509 81770 : case 26:
510 81770 : return dx*dx*dz/disty;
511 :
512 81770 : case 27:
513 81770 : return 2.*dx*dy*dz/disty;
514 :
515 81770 : case 28:
516 81770 : return 3.*dy*dy*dz/disty;
517 :
518 21891 : case 29:
519 21891 : return 0.;
520 :
521 81770 : case 30:
522 81770 : return dx*dz*dz/disty;
523 :
524 81770 : case 31:
525 81770 : return 2.*dy*dz*dz/disty;
526 :
527 21891 : case 32:
528 21891 : return 0.;
529 :
530 81770 : case 33:
531 81770 : 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 19414024 : case 2:
558 : {
559 19414024 : 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 2937706 : case 3:
573 2937706 : 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 803405 : case 7:
586 803405 : return dx/distz;
587 :
588 803405 : case 8:
589 803405 : return dy/distz;
590 :
591 803405 : case 9:
592 803405 : 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 161622 : case 14:
608 161622 : return dx*dx/distz;
609 :
610 161622 : case 15:
611 161622 : return dx*dy/distz;
612 :
613 161622 : case 16:
614 161622 : return dy*dy/distz;
615 :
616 161622 : case 17:
617 161622 : return 2.*dx*dz/distz;
618 :
619 161622 : case 18:
620 161622 : return 2.*dy*dz/distz;
621 :
622 161622 : case 19:
623 161622 : 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 81770 : case 25:
642 81770 : return dx*dx*dx/distz;
643 :
644 81770 : case 26:
645 81770 : return dx*dx*dy/distz;
646 :
647 81770 : case 27:
648 81770 : return dx*dy*dy/distz;
649 :
650 81770 : case 28:
651 81770 : return dy*dy*dy/distz;
652 :
653 81770 : case 29:
654 81770 : return 2.*dx*dx*dz/distz;
655 :
656 81770 : case 30:
657 81770 : return 2.*dx*dy*dz/distz;
658 :
659 81770 : case 31:
660 81770 : return 2.*dy*dy*dz/distz;
661 :
662 81770 : case 32:
663 81770 : return 3.*dx*dz*dz/distz;
664 :
665 81770 : case 33:
666 81770 : return 3.*dy*dz*dz/distz;
667 :
668 81770 : case 34:
669 81770 : 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 55289904 : 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 55289904 : Point avg = elem->vertex_average();
745 14360922 : Point max_distance = Point(0.,0.,0.);
746 988767132 : for (const Point & p : elem->node_ref_range())
747 3733908912 : for (unsigned int d = 0; d < 3; d++)
748 : {
749 3528084222 : const Real distance = std::abs(avg(d) - p(d));
750 2800431684 : max_distance(d) = std::max(distance, max_distance(d));
751 : }
752 :
753 55289904 : const Real x = point_in(0);
754 55289904 : const Real y = point_in(1);
755 55289904 : const Real z = point_in(2);
756 55289904 : const Real xc = avg(0);
757 55289904 : const Real yc = avg(1);
758 55289904 : const Real zc = avg(2);
759 55289904 : const Real distx = max_distance(0);
760 55289904 : const Real disty = max_distance(1);
761 55289904 : const Real distz = max_distance(2);
762 55289904 : const Real dx = (x - xc)/distx;
763 55289904 : const Real dy = (y - yc)/disty;
764 55289904 : const Real dz = (z - zc)/distz;
765 55289904 : const Real dist2x = pow(distx,2.);
766 55289904 : const Real dist2y = pow(disty,2.);
767 55289904 : const Real dist2z = pow(distz,2.);
768 55289904 : const Real distxy = distx * disty;
769 55289904 : const Real distxz = distx * distz;
770 55289904 : 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 55289904 : switch (j)
782 : {
783 : // d^2()/dx^2
784 9214984 : case 0:
785 : {
786 9214984 : 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 545357 : case 4:
799 545357 : 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 161622 : case 10:
810 161622 : return 6.*dx/dist2x;
811 :
812 161622 : case 11:
813 161622 : return 2.*dy/dist2x;
814 :
815 86574 : case 12:
816 : case 13:
817 86574 : return 0.;
818 :
819 161622 : case 14:
820 161622 : 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 81770 : case 20:
831 81770 : return 12.*dx*dx/dist2x;
832 :
833 81770 : case 21:
834 81770 : return 6.*dx*dy/dist2x;
835 :
836 81770 : case 22:
837 81770 : return 2.*dy*dy/dist2x;
838 :
839 43782 : case 23:
840 : case 24:
841 43782 : return 0.;
842 :
843 81770 : case 25:
844 81770 : return 6.*dx*dz/dist2x;
845 :
846 81770 : case 26:
847 81770 : return 2.*dy*dz/dist2x;
848 :
849 43782 : case 27:
850 : case 28:
851 43782 : return 0.;
852 :
853 81770 : case 29:
854 81770 : 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 9214984 : case 1:
885 : {
886 9214984 : 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 545357 : case 5:
902 545357 : 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 161622 : case 11:
915 161622 : return 2.*dx/distxy;
916 :
917 161622 : case 12:
918 161622 : return 2.*dy/distxy;
919 :
920 86574 : case 13:
921 : case 14:
922 86574 : return 0.;
923 :
924 161622 : case 15:
925 161622 : 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 81770 : case 21:
938 81770 : return 3.*dx*dx/distxy;
939 :
940 81770 : case 22:
941 81770 : return 4.*dx*dy/distxy;
942 :
943 81770 : case 23:
944 81770 : return 3.*dy*dy/distxy;
945 :
946 43782 : case 24:
947 : case 25:
948 43782 : return 0.;
949 :
950 81770 : case 26:
951 81770 : return 2.*dx*dz/distxy;
952 :
953 81770 : case 27:
954 81770 : return 2.*dy*dz/distxy;
955 :
956 43782 : case 28:
957 : case 29:
958 43782 : return 0.;
959 :
960 81770 : case 30:
961 81770 : 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 9214984 : case 2:
991 : {
992 9214984 : 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 545357 : case 6:
1009 545357 : 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 161622 : case 12:
1022 161622 : return 2.*dx/dist2y;
1023 161622 : case 13:
1024 161622 : return 6.*dy/dist2y;
1025 :
1026 86574 : case 14:
1027 : case 15:
1028 86574 : return 0.;
1029 :
1030 161622 : case 16:
1031 161622 : 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 81770 : case 22:
1044 81770 : return 2.*dx*dx/dist2y;
1045 :
1046 81770 : case 23:
1047 81770 : return 6.*dx*dy/dist2y;
1048 :
1049 81770 : case 24:
1050 81770 : return 12.*dy*dy/dist2y;
1051 :
1052 43782 : case 25:
1053 : case 26:
1054 43782 : return 0.;
1055 :
1056 81770 : case 27:
1057 81770 : return 2.*dx*dz/dist2y;
1058 :
1059 81770 : case 28:
1060 81770 : return 6.*dy*dz/dist2y;
1061 :
1062 43782 : case 29:
1063 : case 30:
1064 43782 : return 0.;
1065 :
1066 81770 : case 31:
1067 81770 : 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 9214984 : case 3:
1096 : {
1097 9214984 : 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 545357 : case 7:
1115 545357 : 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 161622 : case 14:
1129 161622 : return 2.*dx/distxz;
1130 :
1131 161622 : case 15:
1132 161622 : return dy/distxz;
1133 :
1134 43287 : case 16:
1135 43287 : return 0.;
1136 :
1137 161622 : case 17:
1138 161622 : 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 81770 : case 25:
1153 81770 : return 3.*dx*dx/distxz;
1154 :
1155 81770 : case 26:
1156 81770 : return 2.*dx*dy/distxz;
1157 :
1158 81770 : case 27:
1159 81770 : return dy*dy/distxz;
1160 :
1161 21891 : case 28:
1162 21891 : return 0.;
1163 :
1164 81770 : case 29:
1165 81770 : return 4.*dx*dz/distxz;
1166 :
1167 81770 : case 30:
1168 81770 : return 2.*dy*dz/distxz;
1169 :
1170 21891 : case 31:
1171 21891 : return 0.;
1172 :
1173 81770 : case 32:
1174 81770 : 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 9214984 : case 4:
1201 : {
1202 9214984 : 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 545357 : case 8:
1221 545357 : 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 161622 : case 15:
1235 161622 : return dx/distyz;
1236 :
1237 161622 : case 16:
1238 161622 : return 2.*dy/distyz;
1239 :
1240 43287 : case 17:
1241 43287 : return 0.;
1242 :
1243 161622 : case 18:
1244 161622 : 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 81770 : case 26:
1259 81770 : return dx*dx/distyz;
1260 :
1261 81770 : case 27:
1262 81770 : return 2.*dx*dy/distyz;
1263 :
1264 81770 : case 28:
1265 81770 : return 3.*dy*dy/distyz;
1266 :
1267 21891 : case 29:
1268 21891 : return 0.;
1269 :
1270 81770 : case 30:
1271 81770 : return 2.*dx*dz/distyz;
1272 :
1273 81770 : case 31:
1274 81770 : return 4.*dy*dz/distyz;
1275 :
1276 21891 : case 32:
1277 21891 : return 0.;
1278 :
1279 81770 : case 33:
1280 81770 : 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 9214984 : case 5:
1307 : {
1308 9214984 : 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 545357 : case 9:
1328 545357 : 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 161622 : case 17:
1341 161622 : return 2.*dx/dist2z;
1342 :
1343 161622 : case 18:
1344 161622 : return 2.*dy/dist2z;
1345 :
1346 161622 : case 19:
1347 161622 : 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 81770 : case 29:
1362 81770 : return 2.*dx*dx/dist2z;
1363 :
1364 81770 : case 30:
1365 81770 : return 2.*dx*dy/dist2z;
1366 :
1367 81770 : case 31:
1368 81770 : return 2.*dy*dy/dist2z;
1369 :
1370 81770 : case 32:
1371 81770 : return 6.*dx*dz/dist2z;
1372 :
1373 81770 : case 33:
1374 81770 : return 6.*dy*dz/dist2z;
1375 :
1376 81770 : case 34:
1377 81770 : 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
|
