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 : #include "libmesh/elem.h"
24 :
25 :
26 : namespace libMesh
27 : {
28 :
29 :
30 1360557 : LIBMESH_DEFAULT_VECTORIZED_FE(2,MONOMIAL)
31 :
32 :
33 : template <>
34 1676230 : 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 462919 : libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)*
42 : (static_cast<unsigned int>(order)+2)/2);
43 :
44 1676230 : const Real xi = p(0);
45 1676230 : const Real eta = p(1);
46 :
47 1676230 : switch (i)
48 : {
49 : // constant
50 154783 : case 0:
51 154783 : return 1.;
52 :
53 : // linear
54 287631 : case 1:
55 287631 : return xi;
56 :
57 287631 : case 2:
58 287631 : return eta;
59 :
60 : // quadratics
61 47464 : case 3:
62 47464 : return xi*xi;
63 :
64 47464 : case 4:
65 47464 : return xi*eta;
66 :
67 47464 : case 5:
68 47464 : return eta*eta;
69 :
70 : // cubics
71 38898 : case 6:
72 38898 : return xi*xi*xi;
73 :
74 38898 : case 7:
75 38898 : return xi*xi*eta;
76 :
77 38898 : case 8:
78 38898 : return xi*eta*eta;
79 :
80 38898 : case 9:
81 38898 : return eta*eta*eta;
82 :
83 : // quartics
84 28604 : case 10:
85 28604 : return xi*xi*xi*xi;
86 :
87 28604 : case 11:
88 28604 : return xi*xi*xi*eta;
89 :
90 28604 : case 12:
91 28604 : return xi*xi*eta*eta;
92 :
93 28604 : case 13:
94 28604 : return xi*eta*eta*eta;
95 :
96 28604 : case 14:
97 28604 : return eta*eta*eta*eta;
98 :
99 24876 : default:
100 24876 : unsigned int o = 0;
101 575064 : for (; i >= (o+1)*(o+2)/2; o++) { }
102 95844 : const int ny = i - (o*(o+1)/2);
103 95844 : const int nx = o - ny;
104 24876 : Real val = 1.;
105 335454 : for (int index=0; index != nx; index++)
106 289362 : val *= xi;
107 335454 : for (int index=0; index != ny; index++)
108 289362 : 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 1372840 : 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 361789 : libmesh_assert(elem);
129 :
130 : // by default call the orientation-independent shape functions
131 2096418 : 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 1344420 : 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 391712 : libmesh_assert_less (j, 2);
160 :
161 391712 : libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)*
162 : (static_cast<unsigned int>(order)+2)/2);
163 :
164 1344420 : const Real xi = p(0);
165 1344420 : const Real eta = p(1);
166 :
167 : // monomials. since they are hierarchic we only need one case block.
168 :
169 1344420 : switch (j)
170 : {
171 : // d()/dxi
172 672210 : case 0:
173 : {
174 672210 : switch (i)
175 : {
176 : // constants
177 41210 : case 0:
178 41210 : 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 29840 : case 3:
189 29840 : 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 23100 : case 6:
199 23100 : return 3.*xi*xi;
200 :
201 23100 : case 7:
202 23100 : return 2.*xi*eta;
203 :
204 23100 : case 8:
205 23100 : return eta*eta;
206 :
207 6144 : case 9:
208 6144 : return 0.;
209 :
210 : // quartics
211 15982 : case 10:
212 15982 : return 4.*xi*xi*xi;
213 :
214 15982 : case 11:
215 15982 : return 3.*xi*xi*eta;
216 :
217 15982 : case 12:
218 15982 : return 2.*xi*eta*eta;
219 :
220 15982 : case 13:
221 15982 : return eta*eta*eta;
222 :
223 4254 : case 14:
224 4254 : return 0.;
225 :
226 13356 : default:
227 13356 : unsigned int o = 0;
228 300744 : for (; i >= (o+1)*(o+2)/2; o++) { }
229 50124 : const int ny = i - (o*(o+1)/2);
230 50124 : const int nx = o - ny;
231 50124 : Real val = nx;
232 133664 : for (int index=1; index < nx; index++)
233 110252 : val *= xi;
234 175434 : for (int index=0; index != ny; index++)
235 152022 : val *= eta;
236 13356 : return val;
237 : }
238 : }
239 :
240 :
241 : // d()/deta
242 672210 : case 1:
243 : {
244 672210 : switch (i)
245 : {
246 : // constants
247 41210 : case 0:
248 41210 : 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 29840 : case 5:
265 29840 : return 2.*eta;
266 :
267 : // cubics
268 6144 : case 6:
269 6144 : return 0.;
270 :
271 23100 : case 7:
272 23100 : return xi*xi;
273 :
274 23100 : case 8:
275 23100 : return 2.*xi*eta;
276 :
277 23100 : case 9:
278 23100 : return 3.*eta*eta;
279 :
280 : // quartics
281 4254 : case 10:
282 4254 : return 0.;
283 :
284 15982 : case 11:
285 15982 : return xi*xi*xi;
286 :
287 15982 : case 12:
288 15982 : return 2.*xi*xi*eta;
289 :
290 15982 : case 13:
291 15982 : return 3.*xi*eta*eta;
292 :
293 15982 : case 14:
294 15982 : return 4.*eta*eta*eta;
295 :
296 13356 : default:
297 13356 : unsigned int o = 0;
298 300744 : for (; i >= (o+1)*(o+2)/2; o++) { }
299 50124 : const int ny = i - (o*(o+1)/2);
300 50124 : const int nx = o - ny;
301 50124 : Real val = ny;
302 175434 : for (int index=0; index != nx; index++)
303 152022 : val *= xi;
304 133664 : for (int index=1; index < ny; index++)
305 110252 : 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 836052 : 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 222256 : libmesh_assert(elem);
332 :
333 : // by default call the orientation-independent shape functions
334 1280564 : 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 1254036 : 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 333372 : libmesh_assert_less_equal (j, 2);
366 :
367 333372 : libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)*
368 : (static_cast<unsigned int>(order)+2)/2);
369 :
370 1254036 : const Real xi = p(0);
371 1254036 : const Real eta = p(1);
372 :
373 : // monomials. since they are hierarchic we only need one case block.
374 :
375 1254036 : switch (j)
376 : {
377 : // d^2()/dxi^2
378 418012 : case 0:
379 : {
380 418012 : switch (i)
381 : {
382 : // constants
383 28126 : case 0:
384 : // linears
385 : case 1:
386 : case 2:
387 28126 : 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 23100 : case 6:
399 23100 : return 6.*xi;
400 :
401 23100 : case 7:
402 23100 : return 2.*eta;
403 :
404 12288 : case 8:
405 : case 9:
406 12288 : return 0.;
407 :
408 : // quartics
409 15982 : case 10:
410 15982 : return 12.*xi*xi;
411 :
412 15982 : case 11:
413 15982 : return 6.*xi*eta;
414 :
415 15982 : case 12:
416 15982 : 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 300744 : for (; i >= (o+1)*(o+2)/2; o++) { }
425 50124 : const int ny = i - (o*(o+1)/2);
426 50124 : const int nx = o - ny;
427 50124 : Real val = nx * (nx - 1);
428 100248 : for (int index=2; index < nx; index++)
429 76836 : val *= xi;
430 175434 : for (int index=0; index != ny; index++)
431 152022 : val *= eta;
432 13356 : return val;
433 : }
434 : }
435 :
436 : // d^2()/dxideta
437 418012 : case 1:
438 : {
439 418012 : switch (i)
440 : {
441 : // constants
442 28126 : case 0:
443 :
444 : // linears
445 : case 1:
446 : case 2:
447 28126 : 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 23100 : case 7:
463 23100 : return 2.*xi;
464 :
465 23100 : case 8:
466 23100 : 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 15982 : case 11:
476 15982 : return 3.*xi*xi;
477 :
478 15982 : case 12:
479 15982 : return 4.*xi*eta;
480 :
481 15982 : case 13:
482 15982 : return 3.*eta*eta;
483 :
484 4254 : case 14:
485 4254 : return 0.;
486 :
487 13356 : default:
488 13356 : unsigned int o = 0;
489 300744 : for (; i >= (o+1)*(o+2)/2; o++) { }
490 50124 : const int ny = i - (o*(o+1)/2);
491 50124 : const int nx = o - ny;
492 50124 : Real val = nx * ny;
493 133664 : for (int index=1; index < nx; index++)
494 110252 : val *= xi;
495 133664 : for (int index=1; index < ny; index++)
496 110252 : val *= eta;
497 13356 : return val;
498 : }
499 : }
500 :
501 : // d^2()/deta^2
502 418012 : case 2:
503 : {
504 418012 : switch (i)
505 : {
506 : // constants
507 28126 : case 0:
508 :
509 : // linears
510 : case 1:
511 : case 2:
512 28126 : 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 23100 : case 8:
530 23100 : return 2.*xi;
531 :
532 23100 : case 9:
533 23100 : return 6.*eta;
534 :
535 : // quartics
536 8508 : case 10:
537 : case 11:
538 8508 : return 0.;
539 :
540 15982 : case 12:
541 15982 : return 2.*xi*xi;
542 :
543 15982 : case 13:
544 15982 : return 6.*xi*eta;
545 :
546 15982 : case 14:
547 15982 : return 12.*eta*eta;
548 :
549 13356 : default:
550 13356 : unsigned int o = 0;
551 300744 : for (; i >= (o+1)*(o+2)/2; o++) { }
552 50124 : const int ny = i - (o*(o+1)/2);
553 50124 : const int nx = o - ny;
554 50124 : Real val = ny * (ny - 1);
555 175434 : for (int index=0; index != nx; index++)
556 152022 : val *= xi;
557 100248 : for (int index=2; index < ny; index++)
558 76836 : 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 1254036 : 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 333372 : libmesh_assert(elem);
585 :
586 : // by default call the orientation-independent shape functions
587 1920780 : 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
|
