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 4719119 : LIBMESH_DEFAULT_VECTORIZED_FE(2,MONOMIAL)
31 :
32 :
33 : template <>
34 5786704 : 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 499304 : libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)*
42 : (static_cast<unsigned int>(order)+2)/2);
43 :
44 5786704 : const Real xi = p(0);
45 5786704 : const Real eta = p(1);
46 :
47 5786704 : switch (i)
48 : {
49 : // constant
50 157890 : case 0:
51 157890 : return 1.;
52 :
53 : // linear
54 991070 : case 1:
55 991070 : return xi;
56 :
57 991070 : case 2:
58 991070 : return eta;
59 :
60 : // quadratics
61 172880 : case 3:
62 172880 : return xi*xi;
63 :
64 172880 : case 4:
65 172880 : return xi*eta;
66 :
67 172880 : case 5:
68 172880 : return eta*eta;
69 :
70 : // cubics
71 141980 : case 6:
72 141980 : return xi*xi*xi;
73 :
74 141980 : case 7:
75 141980 : return xi*xi*eta;
76 :
77 141980 : case 8:
78 141980 : return xi*eta*eta;
79 :
80 141980 : case 9:
81 141980 : return eta*eta*eta;
82 :
83 : // quartics
84 104632 : case 10:
85 104632 : return xi*xi*xi*xi;
86 :
87 104632 : case 11:
88 104632 : return xi*xi*xi*eta;
89 :
90 104632 : case 12:
91 104632 : return xi*xi*eta*eta;
92 :
93 104632 : case 13:
94 104632 : return xi*eta*eta*eta;
95 :
96 104632 : case 14:
97 104632 : return eta*eta*eta*eta;
98 :
99 29964 : default:
100 29964 : unsigned int o = 0;
101 2108304 : for (; i >= (o+1)*(o+2)/2; o++) { }
102 351384 : const int ny = i - (o*(o+1)/2);
103 351384 : const int nx = o - ny;
104 29964 : Real val = 1.;
105 1229844 : for (int index=0; index != nx; index++)
106 938388 : val *= xi;
107 1229844 : for (int index=0; index != ny; index++)
108 938388 : val *= eta;
109 29964 : 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 4674274 : 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 398174 : libmesh_assert(elem);
129 :
130 : // by default call the orientation-independent shape functions
131 5470622 : 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 4784470 : 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 420816 : libmesh_assert_less (j, 2);
160 :
161 420816 : libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)*
162 : (static_cast<unsigned int>(order)+2)/2);
163 :
164 4784470 : const Real xi = p(0);
165 4784470 : const Real eta = p(1);
166 :
167 : // monomials. since they are hierarchic we only need one case block.
168 :
169 4784470 : switch (j)
170 : {
171 : // d()/dxi
172 2392235 : case 0:
173 : {
174 2392235 : switch (i)
175 : {
176 : // constants
177 42434 : case 0:
178 42434 : return 0.;
179 :
180 : // linears
181 74096 : case 1:
182 74096 : return 1.;
183 :
184 37048 : case 2:
185 37048 : return 0.;
186 :
187 : // quadratics
188 103888 : case 3:
189 103888 : return 2.*xi;
190 :
191 18040 : case 4:
192 18040 : return eta;
193 :
194 9020 : case 5:
195 9020 : return 0.;
196 :
197 : // cubics
198 80136 : case 6:
199 80136 : return 3.*xi*xi;
200 :
201 80136 : case 7:
202 80136 : return 2.*xi*eta;
203 :
204 80136 : case 8:
205 80136 : return eta*eta;
206 :
207 6960 : case 9:
208 6960 : return 0.;
209 :
210 : // quartics
211 55220 : case 10:
212 55220 : return 4.*xi*xi*xi;
213 :
214 55220 : case 11:
215 55220 : return 3.*xi*xi*eta;
216 :
217 55220 : case 12:
218 55220 : return 2.*xi*eta*eta;
219 :
220 55220 : case 13:
221 55220 : return eta*eta*eta;
222 :
223 4798 : case 14:
224 4798 : return 0.;
225 :
226 14988 : default:
227 14988 : unsigned int o = 0;
228 1034352 : for (; i >= (o+1)*(o+2)/2; o++) { }
229 172392 : const int ny = i - (o*(o+1)/2);
230 172392 : const int nx = o - ny;
231 172392 : Real val = nx;
232 459712 : for (int index=1; index < nx; index++)
233 317296 : val *= xi;
234 603372 : for (int index=0; index != ny; index++)
235 460956 : val *= eta;
236 14988 : return val;
237 : }
238 : }
239 :
240 :
241 : // d()/deta
242 2392235 : case 1:
243 : {
244 2392235 : switch (i)
245 : {
246 : // constants
247 42434 : case 0:
248 42434 : return 0.;
249 :
250 : // linears
251 37048 : case 1:
252 37048 : return 0.;
253 :
254 74096 : case 2:
255 74096 : return 1.;
256 :
257 : // quadratics
258 9020 : case 3:
259 9020 : return 0.;
260 :
261 18040 : case 4:
262 18040 : return xi;
263 :
264 103888 : case 5:
265 103888 : return 2.*eta;
266 :
267 : // cubics
268 6960 : case 6:
269 6960 : return 0.;
270 :
271 80136 : case 7:
272 80136 : return xi*xi;
273 :
274 80136 : case 8:
275 80136 : return 2.*xi*eta;
276 :
277 80136 : case 9:
278 80136 : return 3.*eta*eta;
279 :
280 : // quartics
281 4798 : case 10:
282 4798 : return 0.;
283 :
284 55220 : case 11:
285 55220 : return xi*xi*xi;
286 :
287 55220 : case 12:
288 55220 : return 2.*xi*xi*eta;
289 :
290 55220 : case 13:
291 55220 : return 3.*xi*eta*eta;
292 :
293 55220 : case 14:
294 55220 : return 4.*eta*eta*eta;
295 :
296 14988 : default:
297 14988 : unsigned int o = 0;
298 1034352 : for (; i >= (o+1)*(o+2)/2; o++) { }
299 172392 : const int ny = i - (o*(o+1)/2);
300 172392 : const int nx = o - ny;
301 172392 : Real val = ny;
302 603372 : for (int index=0; index != nx; index++)
303 460956 : val *= xi;
304 459712 : for (int index=1; index < ny; index++)
305 317296 : val *= eta;
306 14988 : 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 2908550 : 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 251360 : libmesh_assert(elem);
332 :
333 : // by default call the orientation-independent shape functions
334 3411270 : 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 4362513 : 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 377028 : libmesh_assert_less_equal (j, 2);
366 :
367 377028 : libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)*
368 : (static_cast<unsigned int>(order)+2)/2);
369 :
370 4362513 : const Real xi = p(0);
371 4362513 : const Real eta = p(1);
372 :
373 : // monomials. since they are hierarchic we only need one case block.
374 :
375 4362513 : switch (j)
376 : {
377 : // d^2()/dxi^2
378 1454171 : case 0:
379 : {
380 1454171 : switch (i)
381 : {
382 : // constants
383 31798 : case 0:
384 : // linears
385 : case 1:
386 : case 2:
387 31798 : return 0.;
388 :
389 : // quadratics
390 18040 : case 3:
391 18040 : return 2.;
392 :
393 18040 : case 4:
394 : case 5:
395 18040 : return 0.;
396 :
397 : // cubics
398 80136 : case 6:
399 80136 : return 6.*xi;
400 :
401 80136 : case 7:
402 80136 : return 2.*eta;
403 :
404 13920 : case 8:
405 : case 9:
406 13920 : return 0.;
407 :
408 : // quartics
409 55220 : case 10:
410 55220 : return 12.*xi*xi;
411 :
412 55220 : case 11:
413 55220 : return 6.*xi*eta;
414 :
415 55220 : case 12:
416 55220 : return 2.*eta*eta;
417 :
418 9596 : case 13:
419 : case 14:
420 9596 : return 0.;
421 :
422 14988 : default:
423 14988 : unsigned int o = 0;
424 1034352 : for (; i >= (o+1)*(o+2)/2; o++) { }
425 172392 : const int ny = i - (o*(o+1)/2);
426 172392 : const int nx = o - ny;
427 172392 : Real val = nx * (nx - 1);
428 344784 : for (int index=2; index < nx; index++)
429 202368 : val *= xi;
430 603372 : for (int index=0; index != ny; index++)
431 460956 : val *= eta;
432 14988 : return val;
433 : }
434 : }
435 :
436 : // d^2()/dxideta
437 1454171 : case 1:
438 : {
439 1454171 : switch (i)
440 : {
441 : // constants
442 31798 : case 0:
443 :
444 : // linears
445 : case 1:
446 : case 2:
447 31798 : return 0.;
448 :
449 : // quadratics
450 9020 : case 3:
451 9020 : return 0.;
452 :
453 18040 : case 4:
454 18040 : return 1.;
455 :
456 9020 : case 5:
457 9020 : return 0.;
458 :
459 : // cubics
460 6960 : case 6:
461 6960 : return 0.;
462 80136 : case 7:
463 80136 : return 2.*xi;
464 :
465 80136 : case 8:
466 80136 : return 2.*eta;
467 :
468 6960 : case 9:
469 6960 : return 0.;
470 :
471 : // quartics
472 4798 : case 10:
473 4798 : return 0.;
474 :
475 55220 : case 11:
476 55220 : return 3.*xi*xi;
477 :
478 55220 : case 12:
479 55220 : return 4.*xi*eta;
480 :
481 55220 : case 13:
482 55220 : return 3.*eta*eta;
483 :
484 4798 : case 14:
485 4798 : return 0.;
486 :
487 14988 : default:
488 14988 : unsigned int o = 0;
489 1034352 : for (; i >= (o+1)*(o+2)/2; o++) { }
490 172392 : const int ny = i - (o*(o+1)/2);
491 172392 : const int nx = o - ny;
492 172392 : Real val = nx * ny;
493 459712 : for (int index=1; index < nx; index++)
494 317296 : val *= xi;
495 459712 : for (int index=1; index < ny; index++)
496 317296 : val *= eta;
497 14988 : return val;
498 : }
499 : }
500 :
501 : // d^2()/deta^2
502 1454171 : case 2:
503 : {
504 1454171 : switch (i)
505 : {
506 : // constants
507 31798 : case 0:
508 :
509 : // linears
510 : case 1:
511 : case 2:
512 31798 : return 0.;
513 :
514 : // quadratics
515 18040 : case 3:
516 : case 4:
517 18040 : return 0.;
518 :
519 18040 : case 5:
520 18040 : return 2.;
521 :
522 : // cubics
523 6960 : case 6:
524 6960 : return 0.;
525 :
526 6960 : case 7:
527 6960 : return 0.;
528 :
529 80136 : case 8:
530 80136 : return 2.*xi;
531 :
532 80136 : case 9:
533 80136 : return 6.*eta;
534 :
535 : // quartics
536 9596 : case 10:
537 : case 11:
538 9596 : return 0.;
539 :
540 55220 : case 12:
541 55220 : return 2.*xi*xi;
542 :
543 55220 : case 13:
544 55220 : return 6.*xi*eta;
545 :
546 55220 : case 14:
547 55220 : return 12.*eta*eta;
548 :
549 14988 : default:
550 14988 : unsigned int o = 0;
551 1034352 : for (; i >= (o+1)*(o+2)/2; o++) { }
552 172392 : const int ny = i - (o*(o+1)/2);
553 172392 : const int nx = o - ny;
554 172392 : Real val = ny * (ny - 1);
555 603372 : for (int index=0; index != nx; index++)
556 460956 : val *= xi;
557 344784 : for (int index=2; index < ny; index++)
558 202368 : val *= eta;
559 14988 : 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 4362513 : 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 377028 : libmesh_assert(elem);
585 :
586 : // by default call the orientation-independent shape functions
587 5116569 : 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
|
