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 2580 : LIBMESH_DEFAULT_VECTORIZED_FE(2,XYZ)
33 :
34 :
35 : template <>
36 4747966 : Real FE<2,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 > 1
43 :
44 398976 : libmesh_assert(elem);
45 :
46 4747966 : Point avg = elem->vertex_average();
47 398976 : Point max_distance = Point(0.,0.,0.);
48 34834938 : for (const Point & p : elem->node_ref_range())
49 90260916 : for (unsigned int d = 0; d < 2; d++)
50 : {
51 65230466 : const Real distance = std::abs(avg(d) - p(d));
52 60173944 : max_distance(d) = std::max(distance, max_distance(d));
53 : }
54 :
55 4747966 : const Real x = point_in(0);
56 4747966 : const Real y = point_in(1);
57 4747966 : const Real xc = avg(0);
58 4747966 : const Real yc = avg(1);
59 4747966 : const Real distx = max_distance(0);
60 4747966 : const Real disty = max_distance(1);
61 4747966 : const Real dx = (x - xc)/distx;
62 4747966 : const Real dy = (y - yc)/disty;
63 :
64 : #ifndef NDEBUG
65 : // totalorder is only used in the assertion below, so
66 : // we avoid declaring it when asserts are not active.
67 398976 : const unsigned int totalorder = order + add_p_level * elem->p_level();
68 : #endif
69 398976 : libmesh_assert_less (i, (totalorder+1)*(totalorder+2)/2);
70 :
71 :
72 : // monomials. since they are hierarchic we only need one case block.
73 4747966 : switch (i)
74 : {
75 : // constant
76 84033 : case 0:
77 84033 : return 1.;
78 :
79 : // linear
80 1002910 : case 1:
81 1002910 : return dx;
82 :
83 1002910 : case 2:
84 1002910 : return dy;
85 :
86 : // quadratics
87 391360 : case 3:
88 391360 : return dx*dx;
89 :
90 391360 : case 4:
91 391360 : return dx*dy;
92 :
93 391360 : case 5:
94 391360 : return dy*dy;
95 :
96 : // cubics
97 83416 : case 6:
98 83416 : return dx*dx*dx;
99 :
100 83416 : case 7:
101 83416 : return dx*dx*dy;
102 :
103 83416 : case 8:
104 83416 : return dx*dy*dy;
105 :
106 83416 : case 9:
107 83416 : return dy*dy*dy;
108 :
109 : // quartics
110 46068 : case 10:
111 46068 : return dx*dx*dx*dx;
112 :
113 46068 : case 11:
114 46068 : return dx*dx*dx*dy;
115 :
116 46068 : case 12:
117 46068 : return dx*dx*dy*dy;
118 :
119 46068 : case 13:
120 46068 : return dx*dy*dy*dy;
121 :
122 46068 : case 14:
123 46068 : return dy*dy*dy*dy;
124 :
125 0 : default:
126 0 : unsigned int o = 0;
127 0 : for (; i >= (o+1)*(o+2)/2; o++) { }
128 0 : unsigned int i2 = i - (o*(o+1)/2);
129 0 : Real val = 1.;
130 0 : for (unsigned int index=i2; index != o; index++)
131 0 : val *= dx;
132 0 : for (unsigned int index=0; index != i2; index++)
133 0 : val *= dy;
134 0 : return val;
135 : }
136 :
137 : #else // LIBMESH_DIM <= 1
138 : libmesh_assert(true || order || add_p_level);
139 : libmesh_ignore(elem, i, point_in);
140 : libmesh_not_implemented();
141 : #endif
142 : }
143 :
144 :
145 :
146 :
147 : template <>
148 0 : Real FE<2,XYZ>::shape(const ElemType,
149 : const Order,
150 : const unsigned int,
151 : const Point &)
152 : {
153 0 : libmesh_error_msg("XYZ polynomials require the element.");
154 : return 0.;
155 : }
156 :
157 :
158 :
159 : template <>
160 0 : Real FE<2,XYZ>::shape(const FEType fet,
161 : const Elem * elem,
162 : const unsigned int i,
163 : const Point & p,
164 : const bool add_p_level)
165 : {
166 0 : return FE<2,XYZ>::shape(elem, fet.order, i, p, add_p_level);
167 : }
168 :
169 :
170 :
171 :
172 :
173 : template <>
174 7926884 : Real FE<2,XYZ>::shape_deriv(const Elem * elem,
175 : const Order libmesh_dbg_var(order),
176 : const unsigned int i,
177 : const unsigned int j,
178 : const Point & point_in,
179 : const bool libmesh_dbg_var(add_p_level))
180 : {
181 : #if LIBMESH_DIM > 1
182 :
183 :
184 666952 : libmesh_assert_less (j, 2);
185 666952 : libmesh_assert(elem);
186 :
187 7926884 : Point avg = elem->vertex_average();
188 666952 : Point max_distance = Point(0.,0.,0.);
189 57575520 : for (const Point & p : elem->node_ref_range())
190 148945908 : for (unsigned int d = 0; d < 2; d++)
191 : {
192 107652648 : const Real distance = std::abs(avg(d) - p(d));
193 99297272 : max_distance(d) = std::max(distance, max_distance(d));
194 : }
195 :
196 7926884 : const Real x = point_in(0);
197 7926884 : const Real y = point_in(1);
198 7926884 : const Real xc = avg(0);
199 7926884 : const Real yc = avg(1);
200 7926884 : const Real distx = max_distance(0);
201 7926884 : const Real disty = max_distance(1);
202 7926884 : const Real dx = (x - xc)/distx;
203 7926884 : const Real dy = (y - yc)/disty;
204 :
205 : #ifndef NDEBUG
206 : // totalorder is only used in the assertion below, so
207 : // we avoid declaring it when asserts are not active.
208 666952 : const unsigned int totalorder = order + add_p_level * elem->p_level();
209 : #endif
210 666952 : libmesh_assert_less (i, (totalorder+1)*(totalorder+2)/2);
211 :
212 : // monomials. since they are hierarchic we only need one case block.
213 :
214 7926884 : switch (j)
215 : {
216 : // d()/dx
217 3963442 : case 0:
218 : {
219 3963442 : switch (i)
220 : {
221 : // constants
222 74704 : case 0:
223 74704 : return 0.;
224 :
225 : // linears
226 891938 : case 1:
227 891938 : return 1./distx;
228 :
229 74704 : case 2:
230 74704 : return 0.;
231 :
232 : // quadratics
233 316524 : case 3:
234 316524 : return 2.*dx/distx;
235 :
236 316524 : case 4:
237 316524 : return dy/distx;
238 :
239 26672 : case 5:
240 26672 : return 0.;
241 :
242 : // cubics
243 51404 : case 6:
244 51404 : return 3.*dx*dx/distx;
245 :
246 51404 : case 7:
247 51404 : return 2.*dx*dy/distx;
248 :
249 51404 : case 8:
250 51404 : return dy*dy/distx;
251 :
252 4462 : case 9:
253 4462 : return 0.;
254 :
255 : // quartics
256 26488 : case 10:
257 26488 : return 4.*dx*dx*dx/distx;
258 :
259 26488 : case 11:
260 26488 : return 3.*dx*dx*dy/distx;
261 :
262 26488 : case 12:
263 26488 : return 2.*dx*dy*dy/distx;
264 :
265 26488 : case 13:
266 26488 : return dy*dy*dy/distx;
267 :
268 2300 : case 14:
269 2300 : return 0.;
270 :
271 0 : default:
272 0 : unsigned int o = 0;
273 0 : for (; i >= (o+1)*(o+2)/2; o++) { }
274 0 : unsigned int i2 = i - (o*(o+1)/2);
275 0 : Real val = o - i2;
276 0 : for (unsigned int index=i2+1; index < o; index++)
277 0 : val *= dx;
278 0 : for (unsigned int index=0; index != i2; index++)
279 0 : val *= dy;
280 0 : return val/distx;
281 : }
282 : }
283 :
284 :
285 : // d()/dy
286 3963442 : case 1:
287 : {
288 3963442 : switch (i)
289 : {
290 : // constants
291 74704 : case 0:
292 74704 : return 0.;
293 :
294 : // linears
295 74704 : case 1:
296 74704 : return 0.;
297 :
298 891938 : case 2:
299 891938 : return 1./disty;
300 :
301 : // quadratics
302 26672 : case 3:
303 26672 : return 0.;
304 :
305 316524 : case 4:
306 316524 : return dx/disty;
307 :
308 316524 : case 5:
309 316524 : return 2.*dy/disty;
310 :
311 : // cubics
312 4462 : case 6:
313 4462 : return 0.;
314 :
315 51404 : case 7:
316 51404 : return dx*dx/disty;
317 :
318 51404 : case 8:
319 51404 : return 2.*dx*dy/disty;
320 :
321 51404 : case 9:
322 51404 : return 3.*dy*dy/disty;
323 :
324 : // quartics
325 2300 : case 10:
326 2300 : return 0.;
327 :
328 26488 : case 11:
329 26488 : return dx*dx*dx/disty;
330 :
331 26488 : case 12:
332 26488 : return 2.*dx*dx*dy/disty;
333 :
334 26488 : case 13:
335 26488 : return 3.*dx*dy*dy/disty;
336 :
337 26488 : case 14:
338 26488 : return 4.*dy*dy*dy/disty;
339 :
340 0 : default:
341 0 : unsigned int o = 0;
342 0 : for (; i >= (o+1)*(o+2)/2; o++) { }
343 0 : unsigned int i2 = i - (o*(o+1)/2);
344 0 : Real val = i2;
345 0 : for (unsigned int index=i2; index != o; index++)
346 0 : val *= dx;
347 0 : for (unsigned int index=1; index <= i2; index++)
348 0 : val *= dy;
349 0 : return val/disty;
350 : }
351 : }
352 :
353 :
354 0 : default:
355 0 : libmesh_error_msg("Invalid j = " << j);
356 : }
357 :
358 : #else // LIBMESH_DIM <= 1
359 : libmesh_assert(true || order || add_p_level);
360 : libmesh_ignore(elem, i, j, point_in);
361 : libmesh_not_implemented();
362 : #endif
363 : }
364 :
365 :
366 : template <>
367 0 : Real FE<2,XYZ>::shape_deriv(const ElemType,
368 : const Order,
369 : const unsigned int,
370 : const unsigned int,
371 : const Point &)
372 : {
373 0 : libmesh_error_msg("XYZ polynomials require the element.");
374 : return 0.;
375 : }
376 :
377 :
378 : template <>
379 0 : Real FE<2,XYZ>::shape_deriv(const FEType fet,
380 : const Elem * elem,
381 : const unsigned int i,
382 : const unsigned int j,
383 : const Point & p,
384 : const bool add_p_level)
385 : {
386 0 : return FE<2,XYZ>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
387 : }
388 :
389 :
390 :
391 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
392 :
393 : template <>
394 5496294 : Real FE<2,XYZ>::shape_second_deriv(const Elem * elem,
395 : const Order libmesh_dbg_var(order),
396 : const unsigned int i,
397 : const unsigned int j,
398 : const Point & point_in,
399 : const bool libmesh_dbg_var(add_p_level))
400 : {
401 : #if LIBMESH_DIM > 1
402 :
403 467376 : libmesh_assert_less_equal (j, 2);
404 467376 : libmesh_assert(elem);
405 :
406 5496294 : Point avg = elem->vertex_average();
407 467376 : Point max_distance = Point(0.,0.,0.);
408 41605056 : for (auto p : make_range(elem->n_nodes()))
409 108326286 : for (unsigned int d = 0; d < 2; d++)
410 : {
411 84490404 : const Real distance = std::abs(avg(d) - elem->point(p)(d));
412 72217524 : max_distance(d) = std::max(distance, max_distance(d));
413 : }
414 :
415 5496294 : const Real x = point_in(0);
416 5496294 : const Real y = point_in(1);
417 5496294 : const Real xc = avg(0);
418 5496294 : const Real yc = avg(1);
419 5496294 : const Real distx = max_distance(0);
420 5496294 : const Real disty = max_distance(1);
421 5496294 : const Real dx = (x - xc)/distx;
422 5496294 : const Real dy = (y - yc)/disty;
423 5496294 : const Real dist2x = pow(distx,2.);
424 5496294 : const Real dist2y = pow(disty,2.);
425 5496294 : const Real distxy = distx * disty;
426 :
427 : #ifndef NDEBUG
428 : // totalorder is only used in the assertion below, so
429 : // we avoid declaring it when asserts are not active.
430 467376 : const unsigned int totalorder = order + add_p_level * elem->p_level();
431 : #endif
432 467376 : libmesh_assert_less (i, (totalorder+1)*(totalorder+2)/2);
433 :
434 : // monomials. since they are hierarchic we only need one case block.
435 :
436 5496294 : switch (j)
437 : {
438 : // d^2()/dx^2
439 1832098 : case 0:
440 : {
441 1832098 : switch (i)
442 : {
443 : // constants
444 73428 : case 0:
445 : // linears
446 : case 1:
447 : case 2:
448 73428 : return 0.;
449 :
450 : // quadratics
451 208524 : case 3:
452 208524 : return 2./dist2x;
453 :
454 35344 : case 4:
455 : case 5:
456 35344 : return 0.;
457 :
458 : // cubics
459 51404 : case 6:
460 51404 : return 6.*dx/dist2x;
461 :
462 51404 : case 7:
463 51404 : return 2.*dy/dist2x;
464 :
465 8924 : case 8:
466 : case 9:
467 8924 : return 0.;
468 :
469 : // quartics
470 26488 : case 10:
471 26488 : return 12.*dx*dx/dist2x;
472 :
473 26488 : case 11:
474 26488 : return 6.*dx*dy/dist2x;
475 :
476 26488 : case 12:
477 26488 : return 2.*dy*dy/dist2x;
478 :
479 4600 : case 13:
480 : case 14:
481 4600 : return 0.;
482 :
483 0 : default:
484 0 : unsigned int o = 0;
485 0 : for (; i >= (o+1)*(o+2)/2; o++) { }
486 0 : unsigned int i2 = i - (o*(o+1)/2);
487 0 : Real val = (o - i2) * (o - i2 - 1);
488 0 : for (unsigned int index=i2+2; index < o; index++)
489 0 : val *= dx;
490 0 : for (unsigned int index=0; index != i2; index++)
491 0 : val *= dy;
492 0 : return val/dist2x;
493 : }
494 : }
495 :
496 : // d^2()/dxdy
497 1832098 : case 1:
498 : {
499 1832098 : switch (i)
500 : {
501 : // constants
502 73428 : case 0:
503 :
504 : // linears
505 : case 1:
506 : case 2:
507 73428 : return 0.;
508 :
509 : // quadratics
510 17672 : case 3:
511 17672 : return 0.;
512 :
513 208524 : case 4:
514 208524 : return 1./distxy;
515 :
516 17672 : case 5:
517 17672 : return 0.;
518 :
519 : // cubics
520 4462 : case 6:
521 4462 : return 0.;
522 51404 : case 7:
523 51404 : return 2.*dx/distxy;
524 :
525 51404 : case 8:
526 51404 : return 2.*dy/distxy;
527 :
528 4462 : case 9:
529 4462 : return 0.;
530 :
531 : // quartics
532 2300 : case 10:
533 2300 : return 0.;
534 :
535 26488 : case 11:
536 26488 : return 3.*dx*dx/distxy;
537 :
538 26488 : case 12:
539 26488 : return 4.*dx*dy/distxy;
540 :
541 26488 : case 13:
542 26488 : return 3.*dy*dy/distxy;
543 :
544 2300 : case 14:
545 2300 : return 0.;
546 :
547 0 : default:
548 0 : unsigned int o = 0;
549 0 : for (; i >= (o+1)*(o+2)/2; o++) { }
550 0 : unsigned int i2 = i - (o*(o+1)/2);
551 0 : Real val = (o - i2) * i2;
552 0 : for (unsigned int index=i2+1; index < o; index++)
553 0 : val *= dx;
554 0 : for (unsigned int index=1; index < i2; index++)
555 0 : val *= dy;
556 0 : return val/distxy;
557 : }
558 : }
559 :
560 : // d^2()/dy^2
561 1832098 : case 2:
562 : {
563 1832098 : switch (i)
564 : {
565 : // constants
566 73428 : case 0:
567 :
568 : // linears
569 : case 1:
570 : case 2:
571 73428 : return 0.;
572 :
573 : // quadratics
574 35344 : case 3:
575 : case 4:
576 35344 : return 0.;
577 :
578 208524 : case 5:
579 208524 : return 2./dist2y;
580 :
581 : // cubics
582 4462 : case 6:
583 4462 : return 0.;
584 :
585 4462 : case 7:
586 4462 : return 0.;
587 :
588 51404 : case 8:
589 51404 : return 2.*dx/dist2y;
590 :
591 51404 : case 9:
592 51404 : return 6.*dy/dist2y;
593 :
594 : // quartics
595 4600 : case 10:
596 : case 11:
597 4600 : return 0.;
598 :
599 26488 : case 12:
600 26488 : return 2.*dx*dx/dist2y;
601 :
602 26488 : case 13:
603 26488 : return 6.*dx*dy/dist2y;
604 :
605 26488 : case 14:
606 26488 : return 12.*dy*dy/dist2y;
607 :
608 0 : default:
609 0 : unsigned int o = 0;
610 0 : for (; i >= (o+1)*(o+2)/2; o++) { }
611 0 : unsigned int i2 = i - (o*(o+1)/2);
612 0 : Real val = i2 * (i2 - 1);
613 0 : for (unsigned int index=i2; index != o; index++)
614 0 : val *= dx;
615 0 : for (unsigned int index=2; index < i2; index++)
616 0 : val *= dy;
617 0 : return val/dist2y;
618 : }
619 : }
620 :
621 0 : default:
622 0 : libmesh_error_msg("Invalid shape function derivative j = " << j);
623 : }
624 :
625 : #else // LIBMESH_DIM <= 1
626 : libmesh_assert(true || order || add_p_level);
627 : libmesh_ignore(elem, i, j, point_in);
628 : libmesh_not_implemented();
629 : #endif
630 : }
631 :
632 :
633 : template <>
634 0 : Real FE<2,XYZ>::shape_second_deriv(const ElemType,
635 : const Order,
636 : const unsigned int,
637 : const unsigned int,
638 : const Point &)
639 : {
640 0 : libmesh_error_msg("XYZ polynomials require the element.");
641 : return 0.;
642 : }
643 :
644 :
645 :
646 : template <>
647 0 : Real FE<2,XYZ>::shape_second_deriv(const FEType fet,
648 : const Elem * elem,
649 : const unsigned int i,
650 : const unsigned int j,
651 : const Point & p,
652 : const bool add_p_level)
653 : {
654 0 : return FE<2,XYZ>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
655 : }
656 :
657 :
658 : #endif
659 :
660 : } // namespace libMesh
|
