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 : // Local includes
20 : #include "libmesh/fe.h"
21 : #include "libmesh/elem.h"
22 : #include "libmesh/utility.h"
23 :
24 :
25 : // Anonymous namespace for functions shared by HIERARCHIC and
26 : // L2_HIERARCHIC implementations. Implementations appear at the bottom
27 : // of this file.
28 : namespace
29 : {
30 : using namespace libMesh;
31 :
32 : Real fe_hierarchic_1D_shape(const ElemType,
33 : const Order libmesh_dbg_var(order),
34 : const unsigned int i,
35 : const Point & p);
36 :
37 : Real fe_hierarchic_1D_shape_deriv(const ElemType,
38 : const Order libmesh_dbg_var(order),
39 : const unsigned int i,
40 : const unsigned int libmesh_dbg_var(j),
41 : const Point & p);
42 :
43 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
44 :
45 : Real fe_hierarchic_1D_shape_second_deriv(const ElemType,
46 : const Order libmesh_dbg_var(order),
47 : const unsigned int i,
48 : const unsigned int libmesh_dbg_var(j),
49 : const Point & p);
50 :
51 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
52 :
53 : } // anonymous namespace
54 :
55 :
56 :
57 : namespace libMesh
58 : {
59 :
60 :
61 126399 : LIBMESH_DEFAULT_VECTORIZED_FE(1,HIERARCHIC)
62 30624 : LIBMESH_DEFAULT_VECTORIZED_FE(1,L2_HIERARCHIC)
63 7200 : LIBMESH_DEFAULT_VECTORIZED_FE(1,SIDE_HIERARCHIC)
64 :
65 :
66 : template <>
67 2663749945 : Real FE<1,HIERARCHIC>::shape(const ElemType elem_type,
68 : const Order order,
69 : const unsigned int i,
70 : const Point & p)
71 : {
72 2663749945 : return fe_hierarchic_1D_shape(elem_type, order, i, p);
73 : }
74 :
75 :
76 :
77 : template <>
78 2712501930 : Real FE<1,L2_HIERARCHIC>::shape(const ElemType elem_type,
79 : const Order order,
80 : const unsigned int i,
81 : const Point & p)
82 : {
83 2712501930 : return fe_hierarchic_1D_shape(elem_type, order, i, p);
84 : }
85 :
86 :
87 :
88 : template <>
89 0 : Real FE<1,SIDE_HIERARCHIC>::shape(const ElemType,
90 : const Order,
91 : const unsigned int i,
92 : const Point & p)
93 : {
94 0 : unsigned int right_side = p(0) > 0; // 0 false, 1 true
95 0 : return (right_side == i);
96 : }
97 :
98 :
99 :
100 : template <>
101 126386 : Real FE<1,HIERARCHIC>::shape(const Elem * elem,
102 : const Order order,
103 : const unsigned int i,
104 : const Point & p,
105 : const bool add_p_level)
106 : {
107 9974 : libmesh_assert(elem);
108 :
109 146302 : return fe_hierarchic_1D_shape(elem->type(), order + add_p_level*elem->p_level(), i, p);
110 : }
111 :
112 :
113 :
114 : template <>
115 0 : Real FE<1,HIERARCHIC>::shape(const FEType fet,
116 : const Elem * elem,
117 : const unsigned int i,
118 : const Point & p,
119 : const bool add_p_level)
120 : {
121 0 : libmesh_assert(elem);
122 0 : return fe_hierarchic_1D_shape(elem->type(), fet.order + add_p_level*elem->p_level(), i, p);
123 : }
124 :
125 :
126 :
127 :
128 :
129 : template <>
130 20028 : Real FE<1,L2_HIERARCHIC>::shape(const Elem * elem,
131 : const Order order,
132 : const unsigned int i,
133 : const Point & p,
134 : const bool add_p_level)
135 : {
136 1669 : libmesh_assert(elem);
137 :
138 23366 : return fe_hierarchic_1D_shape(elem->type(), order + add_p_level*elem->p_level(), i, p);
139 : }
140 :
141 : template <>
142 0 : Real FE<1,L2_HIERARCHIC>::shape(const FEType fet,
143 : const Elem * elem,
144 : const unsigned int i,
145 : const Point & p,
146 : const bool add_p_level)
147 : {
148 0 : libmesh_assert(elem);
149 0 : return fe_hierarchic_1D_shape(elem->type(), fet.order + add_p_level*elem->p_level(), i, p);
150 : }
151 :
152 :
153 :
154 : template <>
155 2400 : Real FE<1,SIDE_HIERARCHIC>::shape(const Elem *,
156 : const Order,
157 : const unsigned int i,
158 : const Point & p,
159 : const bool)
160 : {
161 2400 : unsigned int right_side = p(0) > 0; // 0 false, 1 true
162 2400 : return (right_side == i);
163 : }
164 :
165 : template <>
166 0 : Real FE<1,SIDE_HIERARCHIC>::shape(const FEType,
167 : const Elem *,
168 : const unsigned int i,
169 : const Point & p,
170 : const bool)
171 : {
172 0 : unsigned int right_side = p(0) > 0; // 0 false, 1 true
173 0 : return (right_side == i);
174 : }
175 :
176 :
177 : template <>
178 12298886 : Real FE<1,HIERARCHIC>::shape_deriv(const ElemType elem_type,
179 : const Order order,
180 : const unsigned int i,
181 : const unsigned int j,
182 : const Point & p)
183 : {
184 12298886 : return fe_hierarchic_1D_shape_deriv(elem_type, order, i, j, p);
185 : }
186 :
187 :
188 :
189 : template <>
190 5062176 : Real FE<1,L2_HIERARCHIC>::shape_deriv(const ElemType elem_type,
191 : const Order order,
192 : const unsigned int i,
193 : const unsigned int j,
194 : const Point & p)
195 : {
196 5062176 : return fe_hierarchic_1D_shape_deriv(elem_type, order, i, j, p);
197 : }
198 :
199 :
200 :
201 : template <>
202 0 : Real FE<1,SIDE_HIERARCHIC>::shape_deriv(const ElemType,
203 : const Order,
204 : const unsigned int,
205 : const unsigned int,
206 : const Point &)
207 : {
208 0 : return 0;
209 : }
210 :
211 :
212 :
213 : template <>
214 59372 : Real FE<1,HIERARCHIC>::shape_deriv(const Elem * elem,
215 : const Order order,
216 : const unsigned int i,
217 : const unsigned int j,
218 : const Point & p,
219 : const bool add_p_level)
220 : {
221 4574 : libmesh_assert(elem);
222 :
223 63946 : return fe_hierarchic_1D_shape_deriv(elem->type(),
224 63946 : order + add_p_level*elem->p_level(), i, j, p);
225 : }
226 :
227 :
228 :
229 : template <>
230 0 : Real FE<1,HIERARCHIC>::shape_deriv(const FEType fet,
231 : const Elem * elem,
232 : const unsigned int i,
233 : const unsigned int j,
234 : const Point & p,
235 : const bool add_p_level)
236 : {
237 0 : libmesh_assert(elem);
238 0 : return fe_hierarchic_1D_shape_deriv(elem->type(), fet.order + add_p_level*elem->p_level(), i, j, p);
239 : }
240 :
241 :
242 :
243 :
244 : template <>
245 12468 : Real FE<1,L2_HIERARCHIC>::shape_deriv(const Elem * elem,
246 : const Order order,
247 : const unsigned int i,
248 : const unsigned int j,
249 : const Point & p,
250 : const bool add_p_level)
251 : {
252 1039 : libmesh_assert(elem);
253 :
254 13507 : return fe_hierarchic_1D_shape_deriv(elem->type(),
255 13507 : order + add_p_level*elem->p_level(), i, j, p);
256 : }
257 :
258 :
259 :
260 : template <>
261 0 : Real FE<1,L2_HIERARCHIC>::shape_deriv(const FEType fet,
262 : const Elem * elem,
263 : const unsigned int i,
264 : const unsigned int j,
265 : const Point & p,
266 : const bool add_p_level)
267 : {
268 0 : libmesh_assert(elem);
269 0 : return fe_hierarchic_1D_shape_deriv(elem->type(), fet.order + add_p_level*elem->p_level(), i, j, p);
270 : }
271 :
272 :
273 :
274 : template <>
275 2400 : Real FE<1,SIDE_HIERARCHIC>::shape_deriv(const Elem *,
276 : const Order,
277 : const unsigned int,
278 : const unsigned int,
279 : const Point &,
280 : const bool)
281 : {
282 2400 : return 0;
283 : }
284 :
285 :
286 :
287 : template <>
288 0 : Real FE<1,SIDE_HIERARCHIC>::shape_deriv(const FEType,
289 : const Elem *,
290 : const unsigned int,
291 : const unsigned int,
292 : const Point &,
293 : const bool)
294 : {
295 0 : return 0;
296 : }
297 :
298 :
299 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
300 :
301 : template <>
302 103680 : Real FE<1,HIERARCHIC>::shape_second_deriv(const ElemType elem_type,
303 : const Order order,
304 : const unsigned int i,
305 : const unsigned int j,
306 : const Point & p)
307 : {
308 103680 : return fe_hierarchic_1D_shape_second_deriv(elem_type, order, i, j, p);
309 : }
310 :
311 :
312 :
313 :
314 : template <>
315 0 : Real FE<1,L2_HIERARCHIC>::shape_second_deriv(const ElemType elem_type,
316 : const Order order,
317 : const unsigned int i,
318 : const unsigned int j,
319 : const Point & p)
320 : {
321 0 : return fe_hierarchic_1D_shape_second_deriv(elem_type, order, i, j, p);
322 : }
323 :
324 :
325 :
326 : template <>
327 0 : Real FE<1,SIDE_HIERARCHIC>::shape_second_deriv(const ElemType,
328 : const Order,
329 : const unsigned int,
330 : const unsigned int,
331 : const Point &)
332 : {
333 0 : return 0;
334 : }
335 :
336 :
337 :
338 : template <>
339 27130 : Real FE<1,HIERARCHIC>::shape_second_deriv(const Elem * elem,
340 : const Order order,
341 : const unsigned int i,
342 : const unsigned int j,
343 : const Point & p,
344 : const bool add_p_level)
345 : {
346 2101 : libmesh_assert(elem);
347 :
348 29231 : return fe_hierarchic_1D_shape_second_deriv(elem->type(),
349 29231 : order + add_p_level*elem->p_level(), i, j, p);
350 : }
351 :
352 :
353 :
354 : template <>
355 0 : Real FE<1,HIERARCHIC>::shape_second_deriv(const FEType fet,
356 : const Elem * elem,
357 : const unsigned int i,
358 : const unsigned int j,
359 : const Point & p,
360 : const bool add_p_level)
361 : {
362 0 : libmesh_assert(elem);
363 0 : return fe_hierarchic_1D_shape_second_deriv(elem->type(),
364 0 : fet.order + add_p_level*elem->p_level(), i, j, p);
365 : }
366 :
367 :
368 : template <>
369 12468 : Real FE<1,L2_HIERARCHIC>::shape_second_deriv(const Elem * elem,
370 : const Order order,
371 : const unsigned int i,
372 : const unsigned int j,
373 : const Point & p,
374 : const bool add_p_level)
375 : {
376 1039 : libmesh_assert(elem);
377 :
378 13507 : return fe_hierarchic_1D_shape_second_deriv(elem->type(),
379 13507 : order + add_p_level*elem->p_level(), i, j, p);
380 : }
381 :
382 :
383 : template <>
384 0 : Real FE<1,L2_HIERARCHIC>::shape_second_deriv(const FEType fet,
385 : const Elem * elem,
386 : const unsigned int i,
387 : const unsigned int j,
388 : const Point & p,
389 : const bool add_p_level)
390 : {
391 0 : libmesh_assert(elem);
392 0 : return fe_hierarchic_1D_shape_second_deriv(elem->type(),
393 0 : fet.order + add_p_level*elem->p_level(), i, j, p);
394 : }
395 :
396 :
397 : template <>
398 2400 : Real FE<1,SIDE_HIERARCHIC>::shape_second_deriv(const Elem *,
399 : const Order,
400 : const unsigned int,
401 : const unsigned int,
402 : const Point &,
403 : const bool)
404 : {
405 2400 : return 0.;
406 : }
407 :
408 :
409 : template <>
410 0 : Real FE<1,SIDE_HIERARCHIC>::shape_second_deriv(const FEType,
411 : const Elem *,
412 : const unsigned int,
413 : const unsigned int,
414 : const Point &,
415 : const bool)
416 : {
417 0 : return 0.;
418 : }
419 :
420 : #endif
421 :
422 : } // namespace libMesh
423 :
424 :
425 :
426 : namespace
427 : {
428 : using namespace libMesh;
429 :
430 4930423361 : Real fe_hierarchic_1D_shape(const ElemType,
431 : const Order order,
432 : const unsigned int i,
433 : const Point & p)
434 : {
435 445978112 : libmesh_assert_less (i, order+1u);
436 :
437 : // If we were to define p=0 here, it wouldn't be hierarchic
438 4930423361 : libmesh_error_msg_if (order <= 0,
439 : "HIERARCHIC FE families do not support p=0");
440 :
441 : // Declare that we are using our own special power function
442 : // from the Utility namespace. This saves typing later.
443 : using Utility::pow;
444 :
445 4930423361 : const Real xi = p(0);
446 :
447 445978112 : Real returnval = 1.;
448 :
449 4930423361 : switch (i)
450 : {
451 912485645 : case 0:
452 912485645 : returnval = .5*(1. - xi);
453 912485645 : break;
454 912485645 : case 1:
455 912485645 : returnval = .5*(1. + xi);
456 912485645 : break;
457 : // All even-terms have the same form.
458 : // (xi^p - 1.)/p!
459 1508419400 : case 2:
460 1508419400 : returnval = (xi*xi - 1.)/2.;
461 1508419400 : break;
462 51056723 : case 4:
463 553377424 : returnval = (pow<4>(xi) - 1.)/24.;
464 553377424 : break;
465 132 : case 6:
466 1452 : returnval = (pow<6>(xi) - 1.)/720.;
467 1452 : break;
468 :
469 : // All odd-terms have the same form.
470 : // (xi^p - xi)/p!
471 1042900686 : case 3:
472 1042900686 : returnval = (xi*xi*xi - xi)/6.;
473 1042900686 : break;
474 51453 : case 5:
475 750843 : returnval = (pow<5>(xi) - xi)/120.;
476 750843 : break;
477 103 : case 7:
478 1133 : returnval = (pow<7>(xi) - xi)/5040.;
479 1133 : break;
480 103 : default:
481 103 : Real denominator = 1.;
482 10560 : for (unsigned int n=1; n <= i; ++n)
483 : {
484 9427 : returnval *= xi;
485 9427 : denominator *= n;
486 103 : }
487 : // Odd:
488 1133 : if (i % 2)
489 363 : returnval = (returnval - xi)/denominator;
490 : // Even:
491 : else
492 770 : returnval = (returnval - 1.)/denominator;
493 103 : break;
494 : }
495 :
496 4930423361 : return returnval;
497 : }
498 :
499 :
500 :
501 16190599 : Real fe_hierarchic_1D_shape_deriv(const ElemType,
502 : const Order order,
503 : const unsigned int i,
504 : const unsigned int libmesh_dbg_var(j),
505 : const Point & p)
506 : {
507 : // only d()/dxi in 1D!
508 1242303 : libmesh_assert_equal_to (j, 0);
509 1242303 : libmesh_assert_less (i, order+1u);
510 :
511 : // If we were to define p=0 here, it wouldn't be hierarchic
512 16190599 : libmesh_error_msg_if (order <= 0,
513 : "HIERARCHIC FE families do not support p=0");
514 :
515 : // Declare that we are using our own special power function
516 : // from the Utility namespace. This saves typing later.
517 : using Utility::pow;
518 :
519 16190599 : const Real xi = p(0);
520 :
521 1242303 : Real returnval = 1.;
522 :
523 16190599 : switch (i)
524 : {
525 395107 : case 0:
526 395107 : returnval = -.5;
527 395107 : break;
528 5022783 : case 1:
529 395107 : returnval = .5;
530 5022783 : break;
531 : // All even-terms have the same form.
532 : // xi^(p-1)/(p-1)!
533 3655663 : case 2:
534 271923 : returnval = xi;
535 3655663 : break;
536 63844 : case 4:
537 881630 : returnval = pow<3>(xi)/6.;
538 881630 : break;
539 68 : case 6:
540 748 : returnval = pow<5>(xi)/120.;
541 748 : break;
542 : // All odd-terms have the same form.
543 : // (p*xi^(p-1) - 1.)/p!
544 1402332 : case 3:
545 1402332 : returnval = (3*xi*xi - 1.)/6.;
546 1402332 : break;
547 10472 : case 5:
548 203428 : returnval = (5.*pow<4>(xi) - 1.)/120.;
549 203428 : break;
550 54 : case 7:
551 594 : returnval = (7.*pow<6>(xi) - 1.)/5040.;
552 594 : break;
553 58 : default:
554 58 : Real denominator = 1.;
555 5324 : for (unsigned int n=1; n != i; ++n)
556 : {
557 4686 : returnval *= xi;
558 4686 : denominator *= n;
559 58 : }
560 : // Odd:
561 638 : if (i % 2)
562 220 : returnval = (i * returnval - 1.)/denominator/i;
563 : // Even:
564 : else
565 418 : returnval = returnval/denominator;
566 58 : break;
567 : }
568 :
569 16190599 : return returnval;
570 : }
571 :
572 :
573 :
574 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
575 :
576 131498 : Real fe_hierarchic_1D_shape_second_deriv(const ElemType,
577 : const Order order,
578 : const unsigned int i,
579 : const unsigned int libmesh_dbg_var(j),
580 : const Point & p)
581 : {
582 : // only d2()/d2xi in 1D!
583 11780 : libmesh_assert_equal_to (j, 0);
584 11780 : libmesh_assert_less (i, order+1u);
585 :
586 : // If we were to define p=0 here, it wouldn't be hierarchic
587 131498 : libmesh_error_msg_if (order <= 0,
588 : "HIERARCHIC FE families do not support p=0");
589 :
590 : // Declare that we are using our own special power function
591 : // from the Utility namespace. This saves typing later.
592 : using Utility::pow;
593 :
594 131498 : const Real xi = p(0);
595 :
596 11780 : Real returnval = 1.;
597 :
598 131498 : switch (i)
599 : {
600 6164 : case 0:
601 : case 1:
602 6164 : returnval = 0;
603 6164 : break;
604 : // All terms have the same form.
605 : // xi^(p-2)/(p-2)!
606 29708 : case 2:
607 2650 : returnval = 1;
608 29708 : break;
609 20123 : case 3:
610 1822 : returnval = xi;
611 20123 : break;
612 1002 : case 4:
613 11131 : returnval = pow<2>(xi)/2.;
614 11131 : break;
615 52 : case 5:
616 626 : returnval = pow<3>(xi)/6.;
617 626 : break;
618 34 : case 6:
619 374 : returnval = pow<4>(xi)/24.;
620 374 : break;
621 27 : case 7:
622 297 : returnval = pow<5>(xi)/120.;
623 297 : break;
624 :
625 29 : default:
626 29 : Real denominator = 1.;
627 2662 : for (unsigned int n=1; n != i; ++n)
628 : {
629 2343 : returnval *= xi;
630 2343 : denominator *= n;
631 29 : }
632 : // Odd:
633 319 : if (i % 2)
634 110 : returnval = (i * returnval - 1.)/denominator/i;
635 : // Even:
636 : else
637 209 : returnval = returnval/denominator;
638 29 : break;
639 : }
640 :
641 131498 : return returnval;
642 : }
643 :
644 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
645 :
646 : } // anonymous namespace
|
