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 2580 : LIBMESH_DEFAULT_VECTORIZED_FE(2,XYZ)
33 :
34 :
35 : template <>
36 4531510 : 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 380576 : libmesh_assert(elem);
45 :
46 4531510 : Point avg = elem->vertex_average();
47 380576 : Point max_distance = Point(0.,0.,0.);
48 33061278 : for (const Point & p : elem->node_ref_range())
49 85589304 : for (unsigned int d = 0; d < 2; d++)
50 : {
51 61851804 : const Real distance = std::abs(avg(d) - p(d));
52 57059536 : max_distance(d) = std::max(distance, max_distance(d));
53 : }
54 :
55 4531510 : const Real x = point_in(0);
56 4531510 : const Real y = point_in(1);
57 4531510 : const Real xc = avg(0);
58 4531510 : const Real yc = avg(1);
59 4531510 : const Real distx = max_distance(0);
60 4531510 : const Real disty = max_distance(1);
61 4531510 : const Real dx = (x - xc)/distx;
62 4531510 : 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 380576 : const unsigned int totalorder = order + add_p_level * elem->p_level();
68 : #endif
69 380576 : libmesh_assert_less (i, (totalorder+1)*(totalorder+2)/2);
70 :
71 :
72 : // monomials. since they are hierarchic we only need one case block.
73 4531510 : switch (i)
74 : {
75 : // constant
76 82122 : case 0:
77 82122 : return 1.;
78 :
79 : // linear
80 981298 : case 1:
81 981298 : return dx;
82 :
83 981298 : case 2:
84 981298 : return dy;
85 :
86 : // quadratics
87 372172 : case 3:
88 372172 : return dx*dx;
89 :
90 372172 : case 4:
91 372172 : return dx*dy;
92 :
93 372172 : case 5:
94 372172 : return dy*dy;
95 :
96 : // cubics
97 69760 : case 6:
98 69760 : return dx*dx*dx;
99 :
100 69760 : case 7:
101 69760 : return dx*dx*dy;
102 :
103 69760 : case 8:
104 69760 : return dx*dy*dy;
105 :
106 69760 : case 9:
107 69760 : return dy*dy*dy;
108 :
109 : // quartics
110 38412 : case 10:
111 38412 : return dx*dx*dx*dx;
112 :
113 38412 : case 11:
114 38412 : return dx*dx*dx*dy;
115 :
116 38412 : case 12:
117 38412 : return dx*dx*dy*dy;
118 :
119 38412 : case 13:
120 38412 : return dx*dy*dy*dy;
121 :
122 38412 : case 14:
123 38412 : 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 7724084 : 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 649272 : libmesh_assert_less (j, 2);
185 649272 : libmesh_assert(elem);
186 :
187 7724084 : Point avg = elem->vertex_average();
188 649272 : Point max_distance = Point(0.,0.,0.);
189 55892208 : for (const Point & p : elem->node_ref_range())
190 144504372 : for (unsigned int d = 0; d < 2; d++)
191 : {
192 104434984 : const Real distance = std::abs(avg(d) - p(d));
193 96336248 : max_distance(d) = std::max(distance, max_distance(d));
194 : }
195 :
196 7724084 : const Real x = point_in(0);
197 7724084 : const Real y = point_in(1);
198 7724084 : const Real xc = avg(0);
199 7724084 : const Real yc = avg(1);
200 7724084 : const Real distx = max_distance(0);
201 7724084 : const Real disty = max_distance(1);
202 7724084 : const Real dx = (x - xc)/distx;
203 7724084 : 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 649272 : const unsigned int totalorder = order + add_p_level * elem->p_level();
209 : #endif
210 649272 : libmesh_assert_less (i, (totalorder+1)*(totalorder+2)/2);
211 :
212 : // monomials. since they are hierarchic we only need one case block.
213 :
214 7724084 : switch (j)
215 : {
216 : // d()/dx
217 3862042 : case 0:
218 : {
219 3862042 : switch (i)
220 : {
221 : // constants
222 73752 : case 0:
223 73752 : return 0.;
224 :
225 : // linears
226 881018 : case 1:
227 881018 : return 1./distx;
228 :
229 73752 : case 2:
230 73752 : return 0.;
231 :
232 : // quadratics
233 307164 : case 3:
234 307164 : return 2.*dx/distx;
235 :
236 307164 : case 4:
237 307164 : return dy/distx;
238 :
239 25856 : case 5:
240 25856 : return 0.;
241 :
242 : // cubics
243 45164 : case 6:
244 45164 : return 3.*dx*dx/distx;
245 :
246 45164 : case 7:
247 45164 : return 2.*dx*dy/distx;
248 :
249 45164 : case 8:
250 45164 : return dy*dy/distx;
251 :
252 3918 : case 9:
253 3918 : return 0.;
254 :
255 : // quartics
256 23368 : case 10:
257 23368 : return 4.*dx*dx*dx/distx;
258 :
259 23368 : case 11:
260 23368 : return 3.*dx*dx*dy/distx;
261 :
262 23368 : case 12:
263 23368 : return 2.*dx*dy*dy/distx;
264 :
265 23368 : case 13:
266 23368 : return dy*dy*dy/distx;
267 :
268 2028 : case 14:
269 2028 : 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 3862042 : case 1:
287 : {
288 3862042 : switch (i)
289 : {
290 : // constants
291 73752 : case 0:
292 73752 : return 0.;
293 :
294 : // linears
295 73752 : case 1:
296 73752 : return 0.;
297 :
298 881018 : case 2:
299 881018 : return 1./disty;
300 :
301 : // quadratics
302 25856 : case 3:
303 25856 : return 0.;
304 :
305 307164 : case 4:
306 307164 : return dx/disty;
307 :
308 307164 : case 5:
309 307164 : return 2.*dy/disty;
310 :
311 : // cubics
312 3918 : case 6:
313 3918 : return 0.;
314 :
315 45164 : case 7:
316 45164 : return dx*dx/disty;
317 :
318 45164 : case 8:
319 45164 : return 2.*dx*dy/disty;
320 :
321 45164 : case 9:
322 45164 : return 3.*dy*dy/disty;
323 :
324 : // quartics
325 2028 : case 10:
326 2028 : return 0.;
327 :
328 23368 : case 11:
329 23368 : return dx*dx*dx/disty;
330 :
331 23368 : case 12:
332 23368 : return 2.*dx*dx*dy/disty;
333 :
334 23368 : case 13:
335 23368 : return 3.*dx*dy*dy/disty;
336 :
337 23368 : case 14:
338 23368 : 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 5192094 : 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 440856 : libmesh_assert_less_equal (j, 2);
404 440856 : libmesh_assert(elem);
405 :
406 5192094 : Point avg = elem->vertex_average();
407 440856 : Point max_distance = Point(0.,0.,0.);
408 39080088 : for (auto p : make_range(elem->n_nodes()))
409 101663982 : for (unsigned int d = 0; d < 2; d++)
410 : {
411 79278948 : const Real distance = std::abs(avg(d) - elem->point(p)(d));
412 67775988 : max_distance(d) = std::max(distance, max_distance(d));
413 : }
414 :
415 5192094 : const Real x = point_in(0);
416 5192094 : const Real y = point_in(1);
417 5192094 : const Real xc = avg(0);
418 5192094 : const Real yc = avg(1);
419 5192094 : const Real distx = max_distance(0);
420 5192094 : const Real disty = max_distance(1);
421 5192094 : const Real dx = (x - xc)/distx;
422 5192094 : const Real dy = (y - yc)/disty;
423 5192094 : const Real dist2x = pow(distx,2.);
424 5192094 : const Real dist2y = pow(disty,2.);
425 5192094 : 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 440856 : const unsigned int totalorder = order + add_p_level * elem->p_level();
431 : #endif
432 440856 : libmesh_assert_less (i, (totalorder+1)*(totalorder+2)/2);
433 :
434 : // monomials. since they are hierarchic we only need one case block.
435 :
436 5192094 : switch (j)
437 : {
438 : // d^2()/dx^2
439 1730698 : case 0:
440 : {
441 1730698 : switch (i)
442 : {
443 : // constants
444 70572 : case 0:
445 : // linears
446 : case 1:
447 : case 2:
448 70572 : return 0.;
449 :
450 : // quadratics
451 199164 : case 3:
452 199164 : return 2./dist2x;
453 :
454 33712 : case 4:
455 : case 5:
456 33712 : return 0.;
457 :
458 : // cubics
459 45164 : case 6:
460 45164 : return 6.*dx/dist2x;
461 :
462 45164 : case 7:
463 45164 : return 2.*dy/dist2x;
464 :
465 7836 : case 8:
466 : case 9:
467 7836 : return 0.;
468 :
469 : // quartics
470 23368 : case 10:
471 23368 : return 12.*dx*dx/dist2x;
472 :
473 23368 : case 11:
474 23368 : return 6.*dx*dy/dist2x;
475 :
476 23368 : case 12:
477 23368 : return 2.*dy*dy/dist2x;
478 :
479 4056 : case 13:
480 : case 14:
481 4056 : 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 1730698 : case 1:
498 : {
499 1730698 : switch (i)
500 : {
501 : // constants
502 70572 : case 0:
503 :
504 : // linears
505 : case 1:
506 : case 2:
507 70572 : return 0.;
508 :
509 : // quadratics
510 16856 : case 3:
511 16856 : return 0.;
512 :
513 199164 : case 4:
514 199164 : return 1./distxy;
515 :
516 16856 : case 5:
517 16856 : return 0.;
518 :
519 : // cubics
520 3918 : case 6:
521 3918 : return 0.;
522 45164 : case 7:
523 45164 : return 2.*dx/distxy;
524 :
525 45164 : case 8:
526 45164 : return 2.*dy/distxy;
527 :
528 3918 : case 9:
529 3918 : return 0.;
530 :
531 : // quartics
532 2028 : case 10:
533 2028 : return 0.;
534 :
535 23368 : case 11:
536 23368 : return 3.*dx*dx/distxy;
537 :
538 23368 : case 12:
539 23368 : return 4.*dx*dy/distxy;
540 :
541 23368 : case 13:
542 23368 : return 3.*dy*dy/distxy;
543 :
544 2028 : case 14:
545 2028 : 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 1730698 : case 2:
562 : {
563 1730698 : switch (i)
564 : {
565 : // constants
566 70572 : case 0:
567 :
568 : // linears
569 : case 1:
570 : case 2:
571 70572 : return 0.;
572 :
573 : // quadratics
574 33712 : case 3:
575 : case 4:
576 33712 : return 0.;
577 :
578 199164 : case 5:
579 199164 : return 2./dist2y;
580 :
581 : // cubics
582 3918 : case 6:
583 3918 : return 0.;
584 :
585 3918 : case 7:
586 3918 : return 0.;
587 :
588 45164 : case 8:
589 45164 : return 2.*dx/dist2y;
590 :
591 45164 : case 9:
592 45164 : return 6.*dy/dist2y;
593 :
594 : // quartics
595 4056 : case 10:
596 : case 11:
597 4056 : return 0.;
598 :
599 23368 : case 12:
600 23368 : return 2.*dx*dx/dist2y;
601 :
602 23368 : case 13:
603 23368 : return 6.*dx*dy/dist2y;
604 :
605 23368 : case 14:
606 23368 : 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
|
