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 : #include "libmesh/elem.h"
24 :
25 :
26 : namespace libMesh
27 : {
28 :
29 :
30 4169651 : LIBMESH_DEFAULT_VECTORIZED_FE(2,MONOMIAL)
31 :
32 :
33 : template <>
34 5353931 : Real FE<2,MONOMIAL>::shape(const ElemType,
35 : const Order libmesh_dbg_var(order),
36 : const unsigned int i,
37 : const Point & p)
38 : {
39 : #if LIBMESH_DIM > 1
40 :
41 462656 : libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)*
42 : (static_cast<unsigned int>(order)+2)/2);
43 :
44 5353931 : const Real xi = p(0);
45 5353931 : const Real eta = p(1);
46 :
47 5353931 : switch (i)
48 : {
49 : // constant
50 154520 : case 0:
51 154520 : return 1.;
52 :
53 : // linear
54 959318 : case 1:
55 959318 : return xi;
56 :
57 959318 : case 2:
58 959318 : return eta;
59 :
60 : // quadratics
61 144416 : case 3:
62 144416 : return xi*xi;
63 :
64 144416 : case 4:
65 144416 : return xi*eta;
66 :
67 144416 : case 5:
68 144416 : return eta*eta;
69 :
70 : // cubics
71 118292 : case 6:
72 118292 : return xi*xi*xi;
73 :
74 118292 : case 7:
75 118292 : return xi*xi*eta;
76 :
77 118292 : case 8:
78 118292 : return xi*eta*eta;
79 :
80 118292 : case 9:
81 118292 : return eta*eta*eta;
82 :
83 : // quartics
84 86944 : case 10:
85 86944 : return xi*xi*xi*xi;
86 :
87 86944 : case 11:
88 86944 : return xi*xi*xi*eta;
89 :
90 86944 : case 12:
91 86944 : return xi*xi*eta*eta;
92 :
93 86944 : case 13:
94 86944 : return xi*eta*eta*eta;
95 :
96 86944 : case 14:
97 86944 : return eta*eta*eta*eta;
98 :
99 24876 : default:
100 24876 : unsigned int o = 0;
101 1747152 : for (; i >= (o+1)*(o+2)/2; o++) { }
102 291192 : const int ny = i - (o*(o+1)/2);
103 291192 : const int nx = o - ny;
104 24876 : Real val = 1.;
105 1019172 : for (int index=0; index != nx; index++)
106 777732 : val *= xi;
107 1019172 : for (int index=0; index != ny; index++)
108 777732 : val *= eta;
109 24876 : return val;
110 : }
111 :
112 : #else // LIBMESH_DIM == 1
113 : libmesh_ignore(i, p);
114 : libmesh_assert(order);
115 : libmesh_not_implemented();
116 : #endif
117 : }
118 :
119 :
120 :
121 : template <>
122 4241501 : Real FE<2,MONOMIAL>::shape(const Elem * elem,
123 : const Order order,
124 : const unsigned int i,
125 : const Point & p,
126 : const bool add_p_level)
127 : {
128 361526 : libmesh_assert(elem);
129 :
130 : // by default call the orientation-independent shape functions
131 4964553 : return FE<2,MONOMIAL>::shape(elem->type(), order + add_p_level*elem->p_level(), i, p);
132 : }
133 :
134 :
135 : template <>
136 0 : Real FE<2,MONOMIAL>::shape(const FEType fet,
137 : const Elem * elem,
138 : const unsigned int i,
139 : const Point & p,
140 : const bool add_p_level)
141 : {
142 0 : libmesh_assert(elem);
143 : // by default call the orientation-independent shape functions
144 0 : return FE<2,MONOMIAL>::shape(elem->type(), fet.order + add_p_level*elem->p_level(), i, p);
145 : }
146 :
147 :
148 :
149 : template <>
150 4450494 : Real FE<2,MONOMIAL>::shape_deriv(const ElemType,
151 : const Order libmesh_dbg_var(order),
152 : const unsigned int i,
153 : const unsigned int j,
154 : const Point & p)
155 : {
156 : #if LIBMESH_DIM > 1
157 :
158 :
159 391704 : libmesh_assert_less (j, 2);
160 :
161 391704 : libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)*
162 : (static_cast<unsigned int>(order)+2)/2);
163 :
164 4450494 : const Real xi = p(0);
165 4450494 : const Real eta = p(1);
166 :
167 : // monomials. since they are hierarchic we only need one case block.
168 :
169 4450494 : switch (j)
170 : {
171 : // d()/dxi
172 2225247 : case 0:
173 : {
174 2225247 : switch (i)
175 : {
176 : // constants
177 41206 : case 0:
178 41206 : return 0.;
179 :
180 : // linears
181 71648 : case 1:
182 71648 : return 1.;
183 :
184 35824 : case 2:
185 35824 : return 0.;
186 :
187 : // quadratics
188 91408 : case 3:
189 91408 : return 2.*xi;
190 :
191 15864 : case 4:
192 15864 : return eta;
193 :
194 7932 : case 5:
195 7932 : return 0.;
196 :
197 : // cubics
198 70776 : case 6:
199 70776 : return 3.*xi*xi;
200 :
201 70776 : case 7:
202 70776 : return 2.*xi*eta;
203 :
204 70776 : case 8:
205 70776 : return eta*eta;
206 :
207 6144 : case 9:
208 6144 : return 0.;
209 :
210 : // quartics
211 48980 : case 10:
212 48980 : return 4.*xi*xi*xi;
213 :
214 48980 : case 11:
215 48980 : return 3.*xi*xi*eta;
216 :
217 48980 : case 12:
218 48980 : return 2.*xi*eta*eta;
219 :
220 48980 : case 13:
221 48980 : return eta*eta*eta;
222 :
223 4254 : case 14:
224 4254 : return 0.;
225 :
226 13356 : default:
227 13356 : unsigned int o = 0;
228 922032 : for (; i >= (o+1)*(o+2)/2; o++) { }
229 153672 : const int ny = i - (o*(o+1)/2);
230 153672 : const int nx = o - ny;
231 153672 : Real val = nx;
232 409792 : for (int index=1; index < nx; index++)
233 282832 : val *= xi;
234 537852 : for (int index=0; index != ny; index++)
235 410892 : val *= eta;
236 13356 : return val;
237 : }
238 : }
239 :
240 :
241 : // d()/deta
242 2225247 : case 1:
243 : {
244 2225247 : switch (i)
245 : {
246 : // constants
247 41206 : case 0:
248 41206 : return 0.;
249 :
250 : // linears
251 35824 : case 1:
252 35824 : return 0.;
253 :
254 71648 : case 2:
255 71648 : return 1.;
256 :
257 : // quadratics
258 7932 : case 3:
259 7932 : return 0.;
260 :
261 15864 : case 4:
262 15864 : return xi;
263 :
264 91408 : case 5:
265 91408 : return 2.*eta;
266 :
267 : // cubics
268 6144 : case 6:
269 6144 : return 0.;
270 :
271 70776 : case 7:
272 70776 : return xi*xi;
273 :
274 70776 : case 8:
275 70776 : return 2.*xi*eta;
276 :
277 70776 : case 9:
278 70776 : return 3.*eta*eta;
279 :
280 : // quartics
281 4254 : case 10:
282 4254 : return 0.;
283 :
284 48980 : case 11:
285 48980 : return xi*xi*xi;
286 :
287 48980 : case 12:
288 48980 : return 2.*xi*xi*eta;
289 :
290 48980 : case 13:
291 48980 : return 3.*xi*eta*eta;
292 :
293 48980 : case 14:
294 48980 : return 4.*eta*eta*eta;
295 :
296 13356 : default:
297 13356 : unsigned int o = 0;
298 922032 : for (; i >= (o+1)*(o+2)/2; o++) { }
299 153672 : const int ny = i - (o*(o+1)/2);
300 153672 : const int nx = o - ny;
301 153672 : Real val = ny;
302 537852 : for (int index=0; index != nx; index++)
303 410892 : val *= xi;
304 409792 : for (int index=1; index < ny; index++)
305 282832 : val *= eta;
306 13356 : return val;
307 : }
308 : }
309 :
310 0 : default:
311 0 : libmesh_error_msg("Invalid shape function derivative j = " << j);
312 : }
313 :
314 : #else // LIBMESH_DIM == 1
315 : libmesh_ignore(i, j, p);
316 : libmesh_assert(order);
317 : libmesh_not_implemented();
318 : #endif
319 : }
320 :
321 :
322 :
323 : template <>
324 2574574 : Real FE<2,MONOMIAL>::shape_deriv(const Elem * elem,
325 : const Order order,
326 : const unsigned int i,
327 : const unsigned int j,
328 : const Point & p,
329 : const bool add_p_level)
330 : {
331 222248 : libmesh_assert(elem);
332 :
333 : // by default call the orientation-independent shape functions
334 3019070 : return FE<2,MONOMIAL>::shape_deriv(elem->type(), order + add_p_level*elem->p_level(), i, j, p);
335 : }
336 :
337 :
338 :
339 : template <>
340 0 : Real FE<2,MONOMIAL>::shape_deriv(const FEType fet,
341 : const Elem * elem,
342 : const unsigned int i,
343 : const unsigned int j,
344 : const Point & p,
345 : const bool add_p_level)
346 : {
347 0 : libmesh_assert(elem);
348 : // by default call the orientation-independent shape functions
349 0 : return FE<2,MONOMIAL>::shape_deriv(elem->type(), fet.order + add_p_level*elem->p_level(), i, j, p);
350 : }
351 :
352 :
353 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
354 :
355 : template <>
356 3861648 : Real FE<2,MONOMIAL>::shape_second_deriv(const ElemType,
357 : const Order libmesh_dbg_var(order),
358 : const unsigned int i,
359 : const unsigned int j,
360 : const Point & p)
361 : {
362 : #if LIBMESH_DIM > 1
363 :
364 :
365 333366 : libmesh_assert_less_equal (j, 2);
366 :
367 333366 : libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)*
368 : (static_cast<unsigned int>(order)+2)/2);
369 :
370 3861648 : const Real xi = p(0);
371 3861648 : const Real eta = p(1);
372 :
373 : // monomials. since they are hierarchic we only need one case block.
374 :
375 3861648 : switch (j)
376 : {
377 : // d^2()/dxi^2
378 1287216 : case 0:
379 : {
380 1287216 : switch (i)
381 : {
382 : // constants
383 28124 : case 0:
384 : // linears
385 : case 1:
386 : case 2:
387 28124 : return 0.;
388 :
389 : // quadratics
390 15864 : case 3:
391 15864 : return 2.;
392 :
393 15864 : case 4:
394 : case 5:
395 15864 : return 0.;
396 :
397 : // cubics
398 70776 : case 6:
399 70776 : return 6.*xi;
400 :
401 70776 : case 7:
402 70776 : return 2.*eta;
403 :
404 12288 : case 8:
405 : case 9:
406 12288 : return 0.;
407 :
408 : // quartics
409 48980 : case 10:
410 48980 : return 12.*xi*xi;
411 :
412 48980 : case 11:
413 48980 : return 6.*xi*eta;
414 :
415 48980 : case 12:
416 48980 : return 2.*eta*eta;
417 :
418 8508 : case 13:
419 : case 14:
420 8508 : return 0.;
421 :
422 13356 : default:
423 13356 : unsigned int o = 0;
424 922032 : for (; i >= (o+1)*(o+2)/2; o++) { }
425 153672 : const int ny = i - (o*(o+1)/2);
426 153672 : const int nx = o - ny;
427 153672 : Real val = nx * (nx - 1);
428 307344 : for (int index=2; index < nx; index++)
429 180384 : val *= xi;
430 537852 : for (int index=0; index != ny; index++)
431 410892 : val *= eta;
432 13356 : return val;
433 : }
434 : }
435 :
436 : // d^2()/dxideta
437 1287216 : case 1:
438 : {
439 1287216 : switch (i)
440 : {
441 : // constants
442 28124 : case 0:
443 :
444 : // linears
445 : case 1:
446 : case 2:
447 28124 : return 0.;
448 :
449 : // quadratics
450 7932 : case 3:
451 7932 : return 0.;
452 :
453 15864 : case 4:
454 15864 : return 1.;
455 :
456 7932 : case 5:
457 7932 : return 0.;
458 :
459 : // cubics
460 6144 : case 6:
461 6144 : return 0.;
462 70776 : case 7:
463 70776 : return 2.*xi;
464 :
465 70776 : case 8:
466 70776 : return 2.*eta;
467 :
468 6144 : case 9:
469 6144 : return 0.;
470 :
471 : // quartics
472 4254 : case 10:
473 4254 : return 0.;
474 :
475 48980 : case 11:
476 48980 : return 3.*xi*xi;
477 :
478 48980 : case 12:
479 48980 : return 4.*xi*eta;
480 :
481 48980 : case 13:
482 48980 : return 3.*eta*eta;
483 :
484 4254 : case 14:
485 4254 : return 0.;
486 :
487 13356 : default:
488 13356 : unsigned int o = 0;
489 922032 : for (; i >= (o+1)*(o+2)/2; o++) { }
490 153672 : const int ny = i - (o*(o+1)/2);
491 153672 : const int nx = o - ny;
492 153672 : Real val = nx * ny;
493 409792 : for (int index=1; index < nx; index++)
494 282832 : val *= xi;
495 409792 : for (int index=1; index < ny; index++)
496 282832 : val *= eta;
497 13356 : return val;
498 : }
499 : }
500 :
501 : // d^2()/deta^2
502 1287216 : case 2:
503 : {
504 1287216 : switch (i)
505 : {
506 : // constants
507 28124 : case 0:
508 :
509 : // linears
510 : case 1:
511 : case 2:
512 28124 : return 0.;
513 :
514 : // quadratics
515 15864 : case 3:
516 : case 4:
517 15864 : return 0.;
518 :
519 15864 : case 5:
520 15864 : return 2.;
521 :
522 : // cubics
523 6144 : case 6:
524 6144 : return 0.;
525 :
526 6144 : case 7:
527 6144 : return 0.;
528 :
529 70776 : case 8:
530 70776 : return 2.*xi;
531 :
532 70776 : case 9:
533 70776 : return 6.*eta;
534 :
535 : // quartics
536 8508 : case 10:
537 : case 11:
538 8508 : return 0.;
539 :
540 48980 : case 12:
541 48980 : return 2.*xi*xi;
542 :
543 48980 : case 13:
544 48980 : return 6.*xi*eta;
545 :
546 48980 : case 14:
547 48980 : return 12.*eta*eta;
548 :
549 13356 : default:
550 13356 : unsigned int o = 0;
551 922032 : for (; i >= (o+1)*(o+2)/2; o++) { }
552 153672 : const int ny = i - (o*(o+1)/2);
553 153672 : const int nx = o - ny;
554 153672 : Real val = ny * (ny - 1);
555 537852 : for (int index=0; index != nx; index++)
556 410892 : val *= xi;
557 307344 : for (int index=2; index < ny; index++)
558 180384 : val *= eta;
559 13356 : return val;
560 : }
561 : }
562 :
563 0 : default:
564 0 : libmesh_error_msg("Invalid shape function derivative j = " << j);
565 : }
566 :
567 : #else // LIBMESH_DIM == 1
568 : libmesh_assert(order);
569 : libmesh_ignore(i, j, p);
570 : libmesh_not_implemented();
571 : #endif
572 : }
573 :
574 :
575 :
576 : template <>
577 3861648 : Real FE<2,MONOMIAL>::shape_second_deriv(const Elem * elem,
578 : const Order order,
579 : const unsigned int i,
580 : const unsigned int j,
581 : const Point & p,
582 : const bool add_p_level)
583 : {
584 333366 : libmesh_assert(elem);
585 :
586 : // by default call the orientation-independent shape functions
587 4528380 : return FE<2,MONOMIAL>::shape_second_deriv(elem->type(), order + add_p_level*elem->p_level(), i, j, p);
588 : }
589 :
590 :
591 : template <>
592 0 : Real FE<2,MONOMIAL>::shape_second_deriv(const FEType fet,
593 : const Elem * elem,
594 : const unsigned int i,
595 : const unsigned int j,
596 : const Point & p,
597 : const bool add_p_level)
598 : {
599 0 : libmesh_assert(elem);
600 : // by default call the orientation-independent shape functions
601 0 : return FE<2,MONOMIAL>::shape_second_deriv(elem->type(), fet.order + add_p_level*elem->p_level(), i, j, p);
602 : }
603 :
604 : #endif
605 :
606 : } // namespace libMesh
|
