Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : // 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 12622 : LIBMESH_DEFAULT_VECTORIZED_FE(3,XYZ)
33 :
34 :
35 : template <>
36 77456206 : 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 6396911 : libmesh_assert(elem);
44 :
45 77456206 : Point avg = elem->vertex_average();
46 6396911 : Point max_distance = Point(0.,0.,0.);
47 1249458294 : for (auto p : make_range(elem->n_nodes()))
48 4688008352 : for (unsigned int d = 0; d < 3; d++)
49 : {
50 4099045506 : const Real distance = std::abs(avg(d) - elem->point(p)(d));
51 3516006264 : max_distance(d) = std::max(distance, max_distance(d));
52 : }
53 :
54 77456206 : const Real x = point_in(0);
55 77456206 : const Real y = point_in(1);
56 77456206 : const Real z = point_in(2);
57 77456206 : const Real xc = avg(0);
58 77456206 : const Real yc = avg(1);
59 77456206 : const Real zc = avg(2);
60 77456206 : const Real distx = max_distance(0);
61 77456206 : const Real disty = max_distance(1);
62 77456206 : const Real distz = max_distance(2);
63 77456206 : const Real dx = (x - xc)/distx;
64 77456206 : const Real dy = (y - yc)/disty;
65 77456206 : 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 6396911 : const unsigned int totalorder = order + add_p_level * elem->p_level();
71 : #endif
72 6396911 : 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 77456206 : switch (i)
77 : {
78 : // constant
79 820821 : case 0:
80 820821 : return 1.;
81 :
82 : // linears
83 9882689 : case 1:
84 9882689 : return dx;
85 :
86 9882689 : case 2:
87 9882689 : return dy;
88 :
89 9882689 : case 3:
90 9882689 : return dz;
91 :
92 : // quadratics
93 3247160 : case 4:
94 3247160 : return dx*dx;
95 :
96 3247160 : case 5:
97 3247160 : return dx*dy;
98 :
99 3247160 : case 6:
100 3247160 : return dy*dy;
101 :
102 3247160 : case 7:
103 3247160 : return dx*dz;
104 :
105 3247160 : case 8:
106 3247160 : return dz*dy;
107 :
108 3247160 : case 9:
109 3247160 : return dz*dz;
110 :
111 : // cubics
112 1056374 : case 10:
113 1056374 : return dx*dx*dx;
114 :
115 1056374 : case 11:
116 1056374 : return dx*dx*dy;
117 :
118 1056374 : case 12:
119 1056374 : return dx*dy*dy;
120 :
121 1056374 : case 13:
122 1056374 : return dy*dy*dy;
123 :
124 1056374 : case 14:
125 1056374 : return dx*dx*dz;
126 :
127 1056374 : case 15:
128 1056374 : return dx*dy*dz;
129 :
130 1056374 : case 16:
131 1056374 : return dy*dy*dz;
132 :
133 1056374 : case 17:
134 1056374 : return dx*dz*dz;
135 :
136 1056374 : case 18:
137 1056374 : return dy*dz*dz;
138 :
139 1056374 : case 19:
140 1056374 : return dz*dz*dz;
141 :
142 : // quartics
143 525250 : case 20:
144 525250 : return dx*dx*dx*dx;
145 :
146 525250 : case 21:
147 525250 : return dx*dx*dx*dy;
148 :
149 525250 : case 22:
150 525250 : return dx*dx*dy*dy;
151 :
152 525250 : case 23:
153 525250 : return dx*dy*dy*dy;
154 :
155 525250 : case 24:
156 525250 : return dy*dy*dy*dy;
157 :
158 525250 : case 25:
159 525250 : return dx*dx*dx*dz;
160 :
161 525250 : case 26:
162 525250 : return dx*dx*dy*dz;
163 :
164 525250 : case 27:
165 525250 : return dx*dy*dy*dz;
166 :
167 525250 : case 28:
168 525250 : return dy*dy*dy*dz;
169 :
170 525250 : case 29:
171 525250 : return dx*dx*dz*dz;
172 :
173 525250 : case 30:
174 525250 : return dx*dy*dz*dz;
175 :
176 525250 : case 31:
177 525250 : return dy*dy*dz*dz;
178 :
179 525250 : case 32:
180 525250 : return dx*dz*dz*dz;
181 :
182 525250 : case 33:
183 525250 : return dy*dz*dz*dz;
184 :
185 525250 : case 34:
186 525250 : 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 172272204 : 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 14224545 : libmesh_assert(elem);
252 14224545 : libmesh_assert_less (j, 3);
253 :
254 172272204 : Point avg = elem->vertex_average();
255 14224545 : Point max_distance = Point(0.,0.,0.);
256 2805709596 : for (auto p : make_range(elem->n_nodes()))
257 10533749568 : for (unsigned int d = 0; d < 3; d++)
258 : {
259 9211216770 : const Real distance = std::abs(avg(d) - elem->point(p)(d));
260 7900312176 : max_distance(d) = std::max(distance, max_distance(d));
261 : }
262 :
263 172272204 : const Real x = point_in(0);
264 172272204 : const Real y = point_in(1);
265 172272204 : const Real z = point_in(2);
266 172272204 : const Real xc = avg(0);
267 172272204 : const Real yc = avg(1);
268 172272204 : const Real zc = avg(2);
269 172272204 : const Real distx = max_distance(0);
270 172272204 : const Real disty = max_distance(1);
271 172272204 : const Real distz = max_distance(2);
272 172272204 : const Real dx = (x - xc)/distx;
273 172272204 : const Real dy = (y - yc)/disty;
274 172272204 : 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 14224545 : const unsigned int totalorder = order + add_p_level*elem->p_level();
280 : #endif
281 14224545 : libmesh_assert_less (i, (totalorder+1) * (totalorder+2) *
282 : (totalorder+3)/6);
283 :
284 172272204 : switch (j)
285 : {
286 : // d()/dx
287 57424068 : case 0:
288 : {
289 57424068 : switch (i)
290 : {
291 : // constant
292 723861 : case 0:
293 723861 : return 0.;
294 :
295 : // linear
296 8712261 : case 1:
297 8712261 : return 1./distx;
298 :
299 723861 : case 2:
300 723861 : return 0.;
301 :
302 723861 : case 3:
303 723861 : return 0.;
304 :
305 : // quadratic
306 2353934 : case 4:
307 2353934 : return 2.*dx/distx;
308 :
309 2353934 : case 5:
310 2353934 : return dy/distx;
311 :
312 194286 : case 6:
313 194286 : return 0.;
314 :
315 2353934 : case 7:
316 2353934 : return dz/distx;
317 :
318 194286 : case 8:
319 194286 : return 0.;
320 :
321 194286 : case 9:
322 194286 : return 0.;
323 :
324 : // cubic
325 483312 : case 10:
326 483312 : return 3.*dx*dx/distx;
327 :
328 483312 : case 11:
329 483312 : return 2.*dx*dy/distx;
330 :
331 483312 : case 12:
332 483312 : return dy*dy/distx;
333 :
334 38583 : case 13:
335 38583 : return 0.;
336 :
337 483312 : case 14:
338 483312 : return 2.*dx*dz/distx;
339 :
340 483312 : case 15:
341 483312 : return dy*dz/distx;
342 :
343 38583 : case 16:
344 38583 : return 0.;
345 :
346 483312 : case 17:
347 483312 : 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 241220 : case 20:
357 241220 : return 4.*dx*dx*dx/distx;
358 :
359 241220 : case 21:
360 241220 : return 3.*dx*dx*dy/distx;
361 :
362 241220 : case 22:
363 241220 : return 2.*dx*dy*dy/distx;
364 :
365 241220 : case 23:
366 241220 : return dy*dy*dy/distx;
367 :
368 19635 : case 24:
369 19635 : return 0.;
370 :
371 241220 : case 25:
372 241220 : return 3.*dx*dx*dz/distx;
373 :
374 241220 : case 26:
375 241220 : return 2.*dx*dy*dz/distx;
376 :
377 241220 : case 27:
378 241220 : return dy*dy*dz/distx;
379 :
380 19635 : case 28:
381 19635 : return 0.;
382 :
383 241220 : case 29:
384 241220 : return 2.*dx*dz*dz/distx;
385 :
386 241220 : case 30:
387 241220 : return dy*dz*dz/distx;
388 :
389 19635 : case 31:
390 19635 : return 0.;
391 :
392 241220 : case 32:
393 241220 : 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 57424068 : case 1:
423 : {
424 57424068 : switch (i)
425 : {
426 : // constant
427 723861 : case 0:
428 723861 : return 0.;
429 :
430 : // linear
431 723861 : case 1:
432 723861 : return 0.;
433 :
434 8712261 : case 2:
435 8712261 : return 1./disty;
436 :
437 723861 : case 3:
438 723861 : return 0.;
439 :
440 : // quadratic
441 194286 : case 4:
442 194286 : return 0.;
443 :
444 2353934 : case 5:
445 2353934 : return dx/disty;
446 :
447 2353934 : case 6:
448 2353934 : return 2.*dy/disty;
449 :
450 194286 : case 7:
451 194286 : return 0.;
452 :
453 2353934 : case 8:
454 2353934 : return dz/disty;
455 :
456 194286 : case 9:
457 194286 : return 0.;
458 :
459 : // cubic
460 38583 : case 10:
461 38583 : return 0.;
462 :
463 483312 : case 11:
464 483312 : return dx*dx/disty;
465 :
466 483312 : case 12:
467 483312 : return 2.*dx*dy/disty;
468 :
469 483312 : case 13:
470 483312 : return 3.*dy*dy/disty;
471 :
472 38583 : case 14:
473 38583 : return 0.;
474 :
475 483312 : case 15:
476 483312 : return dx*dz/disty;
477 :
478 483312 : case 16:
479 483312 : return 2.*dy*dz/disty;
480 :
481 38583 : case 17:
482 38583 : return 0.;
483 :
484 483312 : case 18:
485 483312 : 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 241220 : case 21:
495 241220 : return dx*dx*dx/disty;
496 :
497 241220 : case 22:
498 241220 : return 2.*dx*dx*dy/disty;
499 :
500 241220 : case 23:
501 241220 : return 3.*dx*dy*dy/disty;
502 :
503 241220 : case 24:
504 241220 : return 4.*dy*dy*dy/disty;
505 :
506 19635 : case 25:
507 19635 : return 0.;
508 :
509 241220 : case 26:
510 241220 : return dx*dx*dz/disty;
511 :
512 241220 : case 27:
513 241220 : return 2.*dx*dy*dz/disty;
514 :
515 241220 : case 28:
516 241220 : return 3.*dy*dy*dz/disty;
517 :
518 19635 : case 29:
519 19635 : return 0.;
520 :
521 241220 : case 30:
522 241220 : return dx*dz*dz/disty;
523 :
524 241220 : case 31:
525 241220 : return 2.*dy*dz*dz/disty;
526 :
527 19635 : case 32:
528 19635 : return 0.;
529 :
530 241220 : case 33:
531 241220 : 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 57424068 : case 2:
558 : {
559 57424068 : switch (i)
560 : {
561 : // constant
562 723861 : case 0:
563 723861 : return 0.;
564 :
565 : // linear
566 723861 : case 1:
567 723861 : return 0.;
568 :
569 723861 : case 2:
570 723861 : return 0.;
571 :
572 8712261 : case 3:
573 8712261 : return 1./distz;
574 :
575 : // quadratic
576 194286 : case 4:
577 194286 : return 0.;
578 :
579 194286 : case 5:
580 194286 : return 0.;
581 :
582 194286 : case 6:
583 194286 : return 0.;
584 :
585 2353934 : case 7:
586 2353934 : return dx/distz;
587 :
588 2353934 : case 8:
589 2353934 : return dy/distz;
590 :
591 2353934 : case 9:
592 2353934 : 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 483312 : case 14:
608 483312 : return dx*dx/distz;
609 :
610 483312 : case 15:
611 483312 : return dx*dy/distz;
612 :
613 483312 : case 16:
614 483312 : return dy*dy/distz;
615 :
616 483312 : case 17:
617 483312 : return 2.*dx*dz/distz;
618 :
619 483312 : case 18:
620 483312 : return 2.*dy*dz/distz;
621 :
622 483312 : case 19:
623 483312 : 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 241220 : case 25:
642 241220 : return dx*dx*dx/distz;
643 :
644 241220 : case 26:
645 241220 : return dx*dx*dy/distz;
646 :
647 241220 : case 27:
648 241220 : return dx*dy*dy/distz;
649 :
650 241220 : case 28:
651 241220 : return dy*dy*dy/distz;
652 :
653 241220 : case 29:
654 241220 : return 2.*dx*dx*dz/distz;
655 :
656 241220 : case 30:
657 241220 : return 2.*dx*dy*dz/distz;
658 :
659 241220 : case 31:
660 241220 : return 2.*dy*dy*dz/distz;
661 :
662 241220 : case 32:
663 241220 : return 3.*dx*dz*dz/distz;
664 :
665 241220 : case 33:
666 241220 : return 3.*dy*dz*dz/distz;
667 :
668 241220 : case 34:
669 241220 : 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 160961688 : 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 13150530 : libmesh_assert(elem);
742 13150530 : libmesh_assert_less (j, 6);
743 :
744 160961688 : Point avg = elem->vertex_average();
745 13150530 : Point max_distance = Point(0.,0.,0.);
746 2857678392 : for (const Point & p : elem->node_ref_range())
747 10786866816 : for (unsigned int d = 0; d < 3; d++)
748 : {
749 8758515186 : const Real distance = std::abs(avg(d) - p(d));
750 8090150112 : max_distance(d) = std::max(distance, max_distance(d));
751 : }
752 :
753 160961688 : const Real x = point_in(0);
754 160961688 : const Real y = point_in(1);
755 160961688 : const Real z = point_in(2);
756 160961688 : const Real xc = avg(0);
757 160961688 : const Real yc = avg(1);
758 160961688 : const Real zc = avg(2);
759 160961688 : const Real distx = max_distance(0);
760 160961688 : const Real disty = max_distance(1);
761 160961688 : const Real distz = max_distance(2);
762 160961688 : const Real dx = (x - xc)/distx;
763 160961688 : const Real dy = (y - yc)/disty;
764 160961688 : const Real dz = (z - zc)/distz;
765 160961688 : const Real dist2x = pow(distx,2.);
766 160961688 : const Real dist2y = pow(disty,2.);
767 160961688 : const Real dist2z = pow(distz,2.);
768 160961688 : const Real distxy = distx * disty;
769 160961688 : const Real distxz = distx * distz;
770 160961688 : 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 13150530 : const unsigned int totalorder = order + add_p_level*elem->p_level();
776 : #endif
777 13150530 : 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 160961688 : switch (j)
782 : {
783 : // d^2()/dx^2
784 26826948 : case 0:
785 : {
786 26826948 : switch (i)
787 : {
788 : // constant
789 732756 : case 0:
790 :
791 : // linear
792 : case 1:
793 : case 2:
794 : case 3:
795 732756 : return 0.;
796 :
797 : // quadratic
798 1579790 : case 4:
799 1579790 : return 2./dist2x;
800 :
801 648870 : case 5:
802 : case 6:
803 : case 7:
804 : case 8:
805 : case 9:
806 648870 : return 0.;
807 :
808 : // cubic
809 483312 : case 10:
810 483312 : return 6.*dx/dist2x;
811 :
812 483312 : case 11:
813 483312 : return 2.*dy/dist2x;
814 :
815 77166 : case 12:
816 : case 13:
817 77166 : return 0.;
818 :
819 483312 : case 14:
820 483312 : 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 241220 : case 20:
831 241220 : return 12.*dx*dx/dist2x;
832 :
833 241220 : case 21:
834 241220 : return 6.*dx*dy/dist2x;
835 :
836 241220 : case 22:
837 241220 : return 2.*dy*dy/dist2x;
838 :
839 39270 : case 23:
840 : case 24:
841 39270 : return 0.;
842 :
843 241220 : case 25:
844 241220 : return 6.*dx*dz/dist2x;
845 :
846 241220 : case 26:
847 241220 : return 2.*dy*dz/dist2x;
848 :
849 39270 : case 27:
850 : case 28:
851 39270 : return 0.;
852 :
853 241220 : case 29:
854 241220 : 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 26826948 : case 1:
885 : {
886 26826948 : switch (i)
887 : {
888 : // constant
889 732756 : case 0:
890 :
891 : // linear
892 : case 1:
893 : case 2:
894 : case 3:
895 732756 : return 0.;
896 :
897 : // quadratic
898 129774 : case 4:
899 129774 : return 0.;
900 :
901 1579790 : case 5:
902 1579790 : return 1./distxy;
903 :
904 519096 : case 6:
905 : case 7:
906 : case 8:
907 : case 9:
908 519096 : return 0.;
909 :
910 : // cubic
911 38583 : case 10:
912 38583 : return 0.;
913 :
914 483312 : case 11:
915 483312 : return 2.*dx/distxy;
916 :
917 483312 : case 12:
918 483312 : return 2.*dy/distxy;
919 :
920 77166 : case 13:
921 : case 14:
922 77166 : return 0.;
923 :
924 483312 : case 15:
925 483312 : 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 241220 : case 21:
938 241220 : return 3.*dx*dx/distxy;
939 :
940 241220 : case 22:
941 241220 : return 4.*dx*dy/distxy;
942 :
943 241220 : case 23:
944 241220 : return 3.*dy*dy/distxy;
945 :
946 39270 : case 24:
947 : case 25:
948 39270 : return 0.;
949 :
950 241220 : case 26:
951 241220 : return 2.*dx*dz/distxy;
952 :
953 241220 : case 27:
954 241220 : return 2.*dy*dz/distxy;
955 :
956 39270 : case 28:
957 : case 29:
958 39270 : return 0.;
959 :
960 241220 : case 30:
961 241220 : 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 26826948 : case 2:
991 : {
992 26826948 : switch (i)
993 : {
994 : // constant
995 732756 : case 0:
996 :
997 : // linear
998 : case 1:
999 : case 2:
1000 : case 3:
1001 732756 : return 0.;
1002 :
1003 : // quadratic
1004 259548 : case 4:
1005 : case 5:
1006 259548 : return 0.;
1007 :
1008 1579790 : case 6:
1009 1579790 : return 2./dist2y;
1010 :
1011 389322 : case 7:
1012 : case 8:
1013 : case 9:
1014 389322 : return 0.;
1015 :
1016 : // cubic
1017 77166 : case 10:
1018 : case 11:
1019 77166 : return 0.;
1020 :
1021 483312 : case 12:
1022 483312 : return 2.*dx/dist2y;
1023 483312 : case 13:
1024 483312 : return 6.*dy/dist2y;
1025 :
1026 77166 : case 14:
1027 : case 15:
1028 77166 : return 0.;
1029 :
1030 483312 : case 16:
1031 483312 : 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 241220 : case 22:
1044 241220 : return 2.*dx*dx/dist2y;
1045 :
1046 241220 : case 23:
1047 241220 : return 6.*dx*dy/dist2y;
1048 :
1049 241220 : case 24:
1050 241220 : return 12.*dy*dy/dist2y;
1051 :
1052 39270 : case 25:
1053 : case 26:
1054 39270 : return 0.;
1055 :
1056 241220 : case 27:
1057 241220 : return 2.*dx*dz/dist2y;
1058 :
1059 241220 : case 28:
1060 241220 : return 6.*dy*dz/dist2y;
1061 :
1062 39270 : case 29:
1063 : case 30:
1064 39270 : return 0.;
1065 :
1066 241220 : case 31:
1067 241220 : 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 26826948 : case 3:
1096 : {
1097 26826948 : switch (i)
1098 : {
1099 : // constant
1100 732756 : case 0:
1101 :
1102 : // linear
1103 : case 1:
1104 : case 2:
1105 : case 3:
1106 732756 : return 0.;
1107 :
1108 : // quadratic
1109 389322 : case 4:
1110 : case 5:
1111 : case 6:
1112 389322 : return 0.;
1113 :
1114 1579790 : case 7:
1115 1579790 : return 1./distxz;
1116 :
1117 259548 : case 8:
1118 : case 9:
1119 259548 : return 0.;
1120 :
1121 : // cubic
1122 154332 : case 10:
1123 : case 11:
1124 : case 12:
1125 : case 13:
1126 154332 : return 0.;
1127 :
1128 483312 : case 14:
1129 483312 : return 2.*dx/distxz;
1130 :
1131 483312 : case 15:
1132 483312 : return dy/distxz;
1133 :
1134 38583 : case 16:
1135 38583 : return 0.;
1136 :
1137 483312 : case 17:
1138 483312 : 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 241220 : case 25:
1153 241220 : return 3.*dx*dx/distxz;
1154 :
1155 241220 : case 26:
1156 241220 : return 2.*dx*dy/distxz;
1157 :
1158 241220 : case 27:
1159 241220 : return dy*dy/distxz;
1160 :
1161 19635 : case 28:
1162 19635 : return 0.;
1163 :
1164 241220 : case 29:
1165 241220 : return 4.*dx*dz/distxz;
1166 :
1167 241220 : case 30:
1168 241220 : return 2.*dy*dz/distxz;
1169 :
1170 19635 : case 31:
1171 19635 : return 0.;
1172 :
1173 241220 : case 32:
1174 241220 : 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 26826948 : case 4:
1201 : {
1202 26826948 : switch (i)
1203 : {
1204 : // constant
1205 732756 : case 0:
1206 :
1207 : // linear
1208 : case 1:
1209 : case 2:
1210 : case 3:
1211 732756 : return 0.;
1212 :
1213 : // quadratic
1214 519096 : case 4:
1215 : case 5:
1216 : case 6:
1217 : case 7:
1218 519096 : return 0.;
1219 :
1220 1579790 : case 8:
1221 1579790 : return 1./distyz;
1222 :
1223 129774 : case 9:
1224 129774 : 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 483312 : case 15:
1235 483312 : return dx/distyz;
1236 :
1237 483312 : case 16:
1238 483312 : return 2.*dy/distyz;
1239 :
1240 38583 : case 17:
1241 38583 : return 0.;
1242 :
1243 483312 : case 18:
1244 483312 : 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 241220 : case 26:
1259 241220 : return dx*dx/distyz;
1260 :
1261 241220 : case 27:
1262 241220 : return 2.*dx*dy/distyz;
1263 :
1264 241220 : case 28:
1265 241220 : return 3.*dy*dy/distyz;
1266 :
1267 19635 : case 29:
1268 19635 : return 0.;
1269 :
1270 241220 : case 30:
1271 241220 : return 2.*dx*dz/distyz;
1272 :
1273 241220 : case 31:
1274 241220 : return 4.*dy*dz/distyz;
1275 :
1276 19635 : case 32:
1277 19635 : return 0.;
1278 :
1279 241220 : case 33:
1280 241220 : 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 26826948 : case 5:
1307 : {
1308 26826948 : switch (i)
1309 : {
1310 : // constant
1311 732756 : case 0:
1312 :
1313 : // linear
1314 : case 1:
1315 : case 2:
1316 : case 3:
1317 732756 : return 0.;
1318 :
1319 : // quadratic
1320 648870 : case 4:
1321 : case 5:
1322 : case 6:
1323 : case 7:
1324 : case 8:
1325 648870 : return 0.;
1326 :
1327 1579790 : case 9:
1328 1579790 : 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 483312 : case 17:
1341 483312 : return 2.*dx/dist2z;
1342 :
1343 483312 : case 18:
1344 483312 : return 2.*dy/dist2z;
1345 :
1346 483312 : case 19:
1347 483312 : 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 241220 : case 29:
1362 241220 : return 2.*dx*dx/dist2z;
1363 :
1364 241220 : case 30:
1365 241220 : return 2.*dx*dy/dist2z;
1366 :
1367 241220 : case 31:
1368 241220 : return 2.*dy*dy/dist2z;
1369 :
1370 241220 : case 32:
1371 241220 : return 6.*dx*dz/dist2z;
1372 :
1373 241220 : case 33:
1374 241220 : return 6.*dy*dz/dist2z;
1375 :
1376 241220 : case 34:
1377 241220 : 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
|
