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 :
20 : // Local includes
21 : #include "libmesh/libmesh_config.h"
22 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
23 :
24 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
25 : #include <cmath>
26 :
27 : #include "libmesh/libmesh_common.h"
28 : #include "libmesh/fe.h"
29 : #include "libmesh/elem.h"
30 : #include "libmesh/utility.h"
31 :
32 :
33 : namespace libMesh
34 : {
35 :
36 :
37 35928 : LIBMESH_DEFAULT_VECTORIZED_FE(1,BERNSTEIN)
38 :
39 :
40 : template <>
41 12428941152 : Real FE<1,BERNSTEIN>::shape(const ElemType,
42 : const Order order,
43 : const unsigned int i,
44 : const Point & p)
45 : {
46 12428941152 : const Real xi = p(0);
47 : using Utility::pow;
48 :
49 12428941152 : switch (order)
50 : {
51 249189792 : case FIRST:
52 :
53 249189792 : switch(i)
54 : {
55 124594896 : case 0:
56 124594896 : return (1.-xi)/2.;
57 124594896 : case 1:
58 124594896 : return (1.+xi)/2.;
59 0 : default:
60 0 : libmesh_error_msg("Invalid shape function index i = " << i);
61 : }
62 :
63 11403626280 : case SECOND:
64 :
65 11403626280 : switch(i)
66 : {
67 2935612152 : case 0:
68 3175801422 : return (1./4.)*pow<2>(1.-xi);
69 2935612152 : case 1:
70 3175801422 : return (1./4.)*pow<2>(1.+xi);
71 5532401976 : case 2:
72 5532401976 : return (1./2.)*(1.-xi)*(1.+xi);
73 0 : default:
74 0 : libmesh_error_msg("Invalid shape function index i = " << i);
75 : }
76 :
77 252379680 : case THIRD:
78 :
79 252379680 : switch(i)
80 : {
81 63094920 : case 0:
82 68352830 : return (1./8.)*pow<3>(1.-xi);
83 63094920 : case 1:
84 68352830 : return (1./8.)*pow<3>(1.+xi);
85 63094920 : case 2:
86 68352830 : return (3./8.)*(1.+xi)*pow<2>(1.-xi);
87 63094920 : case 3:
88 68352830 : return (3./8.)*pow<2>(1.+xi)*(1.-xi);
89 0 : default:
90 0 : libmesh_error_msg("Invalid shape function index i = " << i);
91 : }
92 :
93 523745400 : case FOURTH:
94 :
95 523745400 : switch(i)
96 : {
97 104749080 : case 0:
98 113478170 : return (1./16.)*pow<4>(1.-xi);
99 104749080 : case 1:
100 113478170 : return (1./16.)*pow<4>(1.+xi);
101 104749080 : case 2:
102 113478170 : return (1./ 4.)*(1.+xi)*pow<3>(1.-xi);
103 104749080 : case 3:
104 122207260 : return (3./ 8.)*pow<2>(1.+xi)*pow<2>(1.-xi);
105 104749080 : case 4:
106 113478170 : return (1./ 4.)*pow<3>(1.+xi)*(1.-xi);
107 0 : default:
108 0 : libmesh_error_msg("Invalid shape function index i = " << i);
109 : }
110 :
111 :
112 0 : case FIFTH:
113 :
114 0 : switch(i)
115 : {
116 0 : case 0:
117 0 : return (1./32.)*pow<5>(1.-xi);
118 0 : case 1:
119 0 : return (1./32.)*pow<5>(1.+xi);
120 0 : case 2:
121 0 : return (5./32.)*(1.+xi)*pow<4>(1.-xi);
122 0 : case 3:
123 0 : return (5./16.)*pow<2>(1.+xi)*pow<3>(1.-xi);
124 0 : case 4:
125 0 : return (5./16.)*pow<3>(1.+xi)*pow<2>(1.-xi);
126 0 : case 5:
127 0 : return (5./32.)*pow<4>(1.+xi)*(1.-xi);
128 0 : default:
129 0 : libmesh_error_msg("Invalid shape function index i = " << i);
130 : }
131 :
132 :
133 0 : case SIXTH:
134 :
135 0 : switch (i)
136 : {
137 0 : case 0:
138 0 : return ( 1./64.)*pow<6>(1.-xi);
139 0 : case 1:
140 0 : return ( 1./64.)*pow<6>(1.+xi);
141 0 : case 2:
142 0 : return ( 3./32.)*(1.+xi)*pow<5>(1.-xi);
143 0 : case 3:
144 0 : return (15./64.)*pow<2>(1.+xi)*pow<4>(1.-xi);
145 0 : case 4:
146 0 : return ( 5./16.)*pow<3>(1.+xi)*pow<3>(1.-xi);
147 0 : case 5:
148 0 : return (15./64.)*pow<4>(1.+xi)*pow<2>(1.-xi);
149 0 : case 6:
150 0 : return ( 3./32.)*pow<5>(1.+xi)*(1.-xi);
151 0 : default:
152 0 : libmesh_error_msg("Invalid shape function index i = " << i);
153 : }
154 :
155 0 : default:
156 : {
157 0 : libmesh_assert (order>6);
158 :
159 : // Use this for arbitrary orders.
160 : // Note that this implementation is less efficient.
161 0 : const int p_order = static_cast<int>(order);
162 0 : const int m = p_order-i+1;
163 0 : const int n = (i-1);
164 :
165 0 : Real binomial_p_i = 1;
166 :
167 : // the binomial coefficient (p choose n)
168 : // Using an unsigned long here will work for any of the orders we support.
169 : // Explicitly construct a Real to prevent conversion warnings
170 0 : if (i>1)
171 0 : binomial_p_i = Real(Utility::binomial(static_cast<unsigned long>(p_order),
172 : static_cast<unsigned long>(n)));
173 :
174 0 : switch(i)
175 : {
176 0 : case 0:
177 0 : return binomial_p_i * std::pow((1-xi)/2, p_order);
178 0 : case 1:
179 0 : return binomial_p_i * std::pow((1+xi)/2, p_order);
180 0 : default:
181 : {
182 0 : return binomial_p_i * std::pow((1+xi)/2,n)
183 0 : * std::pow((1-xi)/2,m);
184 : }
185 : }
186 : }
187 : }
188 : }
189 :
190 :
191 :
192 : template <>
193 2752740 : Real FE<1,BERNSTEIN>::shape(const Elem * elem,
194 : const Order order,
195 : const unsigned int i,
196 : const Point & p,
197 : const bool add_p_level)
198 : {
199 250032 : libmesh_assert(elem);
200 :
201 : return FE<1,BERNSTEIN>::shape
202 3002772 : (elem->type(),
203 3002772 : order + add_p_level*elem->p_level(), i, p);
204 : }
205 :
206 :
207 : template <>
208 0 : Real FE<1,BERNSTEIN>::shape(const FEType fet,
209 : const Elem * elem,
210 : const unsigned int i,
211 : const Point & p,
212 : const bool add_p_level)
213 : {
214 0 : libmesh_assert(elem);
215 : return FE<1,BERNSTEIN>::shape
216 0 : (elem->type(),
217 0 : fet.order + add_p_level*elem->p_level(), i, p);
218 : }
219 :
220 :
221 : template <>
222 2624108190 : Real FE<1,BERNSTEIN>::shape_deriv(const ElemType,
223 : const Order order,
224 : const unsigned int i,
225 : const unsigned int libmesh_dbg_var(j),
226 : const Point & p)
227 : {
228 : // only d()/dxi in 1D!
229 :
230 217043843 : libmesh_assert_equal_to (j, 0);
231 :
232 2624108190 : const Real xi = p(0);
233 :
234 : using Utility::pow;
235 :
236 2624108190 : switch (order)
237 : {
238 87942720 : case FIRST:
239 :
240 87942720 : switch(i)
241 : {
242 3664280 : case 0:
243 3664280 : return -.5;
244 43971360 : case 1:
245 43971360 : return .5;
246 0 : default:
247 0 : libmesh_error_msg("Invalid shape function index i = " << i);
248 : }
249 :
250 2535224850 : case SECOND:
251 :
252 2535224850 : switch(i)
253 : {
254 700712454 : case 0:
255 700712454 : return (xi-1.)*.5;
256 700712454 : case 1:
257 700712454 : return (xi+1.)*.5;
258 1133799942 : case 2:
259 1133799942 : return -xi;
260 0 : default:
261 0 : libmesh_error_msg("Invalid shape function index i = " << i);
262 : }
263 :
264 354720 : case THIRD:
265 :
266 354720 : switch(i)
267 : {
268 88680 : case 0:
269 96070 : return -0.375*pow<2>(1.-xi);
270 88680 : case 1:
271 96070 : return 0.375*pow<2>(1.+xi);
272 88680 : case 2:
273 96070 : return -0.375 -.75*xi +1.125*pow<2>(xi);
274 88680 : case 3:
275 96070 : return 0.375 -.75*xi -1.125*pow<2>(xi);
276 0 : default:
277 0 : libmesh_error_msg("Invalid shape function index i = " << i);
278 : }
279 :
280 585900 : case FOURTH:
281 :
282 585900 : switch(i)
283 : {
284 117180 : case 0:
285 126945 : return -0.25*pow<3>(1.-xi);
286 117180 : case 1:
287 126945 : return 0.25*pow<3>(1.+xi);
288 19530 : case 2:
289 126945 : return -0.5 +1.5*pow<2>(xi)-pow<3>(xi);
290 19530 : case 3:
291 117180 : return 1.5*(pow<3>(xi)-xi);
292 19530 : case 4:
293 234360 : return 0.5 -1.5*pow<2>(xi)-pow<3>(xi);
294 0 : default:
295 0 : libmesh_error_msg("Invalid shape function index i = " << i);
296 : }
297 :
298 0 : case FIFTH:
299 :
300 0 : switch(i)
301 : {
302 0 : case 0:
303 0 : return -(5./32.)*pow<4>(xi-1.);
304 0 : case 1:
305 0 : return (5./32.)*pow<4>(xi+1.);
306 0 : case 2:
307 0 : return (5./32.)*pow<4>(1.-xi) -(5./8.)*(1.+xi)*pow<3>(1.-xi);
308 0 : case 3:
309 0 : return (5./ 8.)*(1.+xi)*pow<3>(1.-xi) -(15./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
310 0 : case 4:
311 0 : return -(5./ 8.)*pow<3>(1.+xi)*(1.-xi) +(15./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
312 0 : case 5:
313 0 : return (5./ 8.)*pow<3>(1.+xi)*(1.-xi) -(5./32.)*pow<4>(1.+xi);
314 0 : default:
315 0 : libmesh_error_msg("Invalid shape function index i = " << i);
316 : }
317 :
318 0 : case SIXTH:
319 :
320 0 : switch(i)
321 : {
322 0 : case 0:
323 0 : return -( 3./32.)*pow<5>(1.-xi);
324 0 : case 1:
325 0 : return ( 3./32.)*pow<5>(1.+xi);
326 0 : case 2:
327 0 : return ( 3./32.)*pow<5>(1.-xi)-(15./32.)*(1.+xi)*pow<4>(1.-xi);
328 0 : case 3:
329 0 : return (15./32.)*(1.+xi)*pow<4>(1.-xi)-(15./16.)*pow<2>(1.+xi)*pow<3>(1.-xi);
330 0 : case 4:
331 0 : return -(15./ 8.)*xi +(15./4.)*pow<3>(xi)-(15./8.)*pow<5>(xi);
332 0 : case 5:
333 0 : return -(15./32.)*(1.-xi)*pow<4>(1.+xi)+(15./16.)*pow<2>(1.-xi)*pow<3>(1.+xi);
334 0 : case 6:
335 0 : return (15./32.)*pow<4>(1.+xi)*(1.-xi)-(3./32.)*pow<5>(1.+xi);
336 0 : default:
337 0 : libmesh_error_msg("Invalid shape function index i = " << i);
338 : }
339 :
340 :
341 0 : default:
342 : {
343 0 : libmesh_assert (order>6);
344 :
345 : // Use this for arbitrary orders
346 0 : const int p_order = static_cast<int>(order);
347 0 : const int m = p_order-(i-1);
348 0 : const int n = (i-1);
349 :
350 0 : Real binomial_p_i = 1;
351 :
352 : // the binomial coefficient (p choose n)
353 : // Using an unsigned long here will work for any of the orders we support.
354 : // Explicitly construct a Real to prevent conversion warnings
355 0 : if (i>1)
356 0 : binomial_p_i = Real(Utility::binomial(static_cast<unsigned long>(p_order),
357 : static_cast<unsigned long>(n)));
358 :
359 0 : switch(i)
360 : {
361 0 : case 0:
362 0 : return binomial_p_i * (-1./2.) * p_order * std::pow((1-xi)/2, p_order-1);
363 0 : case 1:
364 0 : return binomial_p_i * ( 1./2.) * p_order * std::pow((1+xi)/2, p_order-1);
365 :
366 0 : default:
367 : {
368 0 : return binomial_p_i * (1./2. * n * std::pow((1+xi)/2,n-1) * std::pow((1-xi)/2,m)
369 0 : - 1./2. * m * std::pow((1+xi)/2,n) * std::pow((1-xi)/2,m-1));
370 : }
371 : }
372 : }
373 :
374 : }
375 : }
376 :
377 :
378 :
379 : template <>
380 1381782 : Real FE<1,BERNSTEIN>::shape_deriv(const Elem * elem,
381 : const Order order,
382 : const unsigned int i,
383 : const unsigned int j,
384 : const Point & p,
385 : const bool add_p_level)
386 : {
387 125467 : libmesh_assert(elem);
388 :
389 : return FE<1,BERNSTEIN>::shape_deriv
390 1507249 : (elem->type(),
391 1507249 : order + add_p_level*elem->p_level(), i, j, p);
392 : }
393 :
394 : template <>
395 0 : Real FE<1,BERNSTEIN>::shape_deriv(const FEType fet,
396 : const Elem * elem,
397 : const unsigned int i,
398 : const unsigned int j,
399 : const Point & p,
400 : const bool add_p_level)
401 : {
402 0 : libmesh_assert(elem);
403 : return FE<1,BERNSTEIN>::shape_deriv
404 0 : (elem->type(),
405 0 : fet.order + add_p_level*elem->p_level(), i, j, p);
406 : }
407 :
408 :
409 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
410 :
411 : template <>
412 1842276 : Real FE<1,BERNSTEIN>::shape_second_deriv(const ElemType,
413 : const Order order,
414 : const unsigned int i,
415 : const unsigned int libmesh_dbg_var(j),
416 : const Point & p)
417 : {
418 : // only d^2()/dxi^2 in 1D!
419 :
420 153523 : libmesh_assert_equal_to (j, 0);
421 :
422 1842276 : const Real xi = p(0);
423 :
424 : using Utility::pow;
425 :
426 1842276 : switch (order)
427 : {
428 241152 : case FIRST:
429 :
430 241152 : switch(i)
431 : {
432 20096 : case 0:
433 : case 1:
434 20096 : return 0;
435 0 : default:
436 0 : libmesh_error_msg("Invalid shape function index i = " << i);
437 : }
438 :
439 1127376 : case SECOND:
440 :
441 1127376 : switch(i)
442 : {
443 62632 : case 0:
444 : case 1:
445 62632 : return .5;
446 375792 : case 2:
447 375792 : return -1;
448 0 : default:
449 0 : libmesh_error_msg("Invalid shape function index i = " << i);
450 : }
451 :
452 178848 : case THIRD:
453 :
454 178848 : switch(i)
455 : {
456 44712 : case 0:
457 44712 : return 0.75*(1.-xi);
458 44712 : case 1:
459 44712 : return 0.75*(1.+xi);
460 44712 : case 2:
461 44712 : return -.75 + 2.25*xi;
462 44712 : case 3:
463 44712 : return -.75 - 2.25*xi;
464 0 : default:
465 0 : libmesh_error_msg("Invalid shape function index i = " << i);
466 : }
467 :
468 294900 : case FOURTH:
469 :
470 294900 : switch(i)
471 : {
472 58980 : case 0:
473 63895 : return 0.75*pow<2>(1.-xi);
474 58980 : case 1:
475 63895 : return 0.75*pow<2>(1.+xi);
476 58980 : case 2:
477 58980 : return 3*(xi - pow<2>(xi));
478 9830 : case 3:
479 58980 : return 1.5*(3*pow<2>(xi)-1);
480 58980 : case 4:
481 63895 : return -3*xi-3*pow<2>(xi);
482 0 : default:
483 0 : libmesh_error_msg("Invalid shape function index i = " << i);
484 : }
485 :
486 0 : case FIFTH:
487 :
488 0 : switch(i)
489 : {
490 0 : case 0:
491 0 : return -(5./8.)*pow<3>(xi-1.);
492 0 : case 1:
493 0 : return (5./8.)*pow<3>(xi+1.);
494 0 : case 2:
495 0 : return -(5./4.)*pow<3>(1.-xi) + (15./8.)*(1.+xi)*pow<2>(1.-xi);
496 0 : case 3:
497 0 : return -(15./ 4.)*(1.+xi)*pow<2>(1.-xi) + (5./ 8.)*pow<3>(1.-xi)
498 0 : + (15./8.)*pow<2>(1.+xi)*(1.-xi);
499 0 : case 4:
500 0 : return (5./ 8.)*pow<3>(1.+xi) - (15./ 4.)*pow<2>(1.+xi)*(1.-xi)
501 0 : +(15./8.)*(1.+xi)*pow<2>(1.-xi);
502 0 : case 5:
503 0 : return -(5./ 8.)*pow<3>(1.+xi) + (15./ 8.)*pow<2>(1.+xi)*(1.-xi)
504 0 : -(5./8.)*pow<3>(1.+xi);
505 0 : default:
506 0 : libmesh_error_msg("Invalid shape function index i = " << i);
507 : }
508 :
509 0 : case SIXTH:
510 :
511 0 : switch(i)
512 : {
513 0 : case 0:
514 0 : return ( 15./32.)*pow<4>(1.-xi);
515 0 : case 1:
516 0 : return ( 15./32.)*pow<4>(1.+xi);
517 0 : case 2:
518 0 : return -( 15./8.)*pow<4>(1.-xi) +
519 0 : ( 15./8.)*(1.+xi)*pow<3>(1.-xi);
520 0 : case 3:
521 0 : return -(15./4.)*(1.+xi)*pow<3>(1.-xi)
522 0 : + (15./32.)*pow<4>(1.-xi)
523 0 : + (45./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
524 0 : case 4:
525 0 : return -(15./ 8.) +(45./4.)*pow<2>(xi) - (75./8.)*pow<4>(xi);
526 0 : case 5:
527 0 : return -(15./4.)*(1.-xi)*pow<3>(1.+xi)
528 0 : + (15./32.)*pow<4>(1.+xi)
529 0 : + (45./16.)*pow<2>(1.-xi)*pow<2>(1.+xi);
530 0 : case 6:
531 0 : return -(15./16.)*pow<4>(1.+xi)
532 0 : + (15./8.)*pow<3>(1.+xi)*(1.-xi);
533 0 : default:
534 0 : libmesh_error_msg("Invalid shape function index i = " << i);
535 : }
536 :
537 :
538 0 : default:
539 : {
540 0 : libmesh_assert (order>6);
541 :
542 : // Use this for arbitrary orders
543 0 : const int p_order = static_cast<int>(order);
544 0 : const int m = p_order-(i-1);
545 0 : const int n = (i-1);
546 :
547 0 : Real binomial_p_i = 1;
548 :
549 : // the binomial coefficient (p choose n)
550 : // Using an unsigned long here will work for any of the orders we support.
551 : // Explicitly construct a Real to prevent conversion warnings
552 0 : if (i>1)
553 0 : binomial_p_i = Real(Utility::binomial(static_cast<unsigned long>(p_order),
554 : static_cast<unsigned long>(n)));
555 :
556 0 : switch(i)
557 : {
558 0 : case 0:
559 0 : return binomial_p_i * (1./4.) * p_order * (p_order-1) * std::pow((1-xi)/2, p_order-2);
560 0 : case 1:
561 0 : return binomial_p_i * (1./4.) * p_order * (p_order-1) * std::pow((1+xi)/2, p_order-2);
562 :
563 0 : default:
564 : {
565 0 : Real val = 0;
566 :
567 0 : if (n == 1)
568 0 : val +=
569 0 : binomial_p_i * (-1./4. * m * std::pow((1-xi)/2,m-1));
570 : else
571 0 : val +=
572 0 : binomial_p_i * (-1./4. * n * m * std::pow((1+xi)/2,n-1) * std::pow((1-xi)/2,m-1) +
573 0 : 1./4. * n * (n-1) * std::pow((1+xi)/2,n-2) * std::pow((1-xi)/2,m));
574 :
575 0 : if (m == 1)
576 0 : val += binomial_p_i * (-1./4. * n * std::pow((1+xi)/2,n-1));
577 : else
578 0 : val +=
579 0 : binomial_p_i * (1./4. * m * (m-1) * std::pow((1+xi)/2,n) * std::pow((1-xi)/2,m-2)
580 0 : - 1./4. * m * n * std::pow((1+xi)/2,n-1) * std::pow((1-xi)/2,m-1));
581 :
582 0 : return val;
583 : }
584 : }
585 : }
586 :
587 : }
588 : }
589 :
590 :
591 :
592 : template <>
593 19404 : Real FE<1,BERNSTEIN>::shape_second_deriv(const Elem * elem,
594 : const Order order,
595 : const unsigned int i,
596 : const unsigned int j,
597 : const Point & p,
598 : const bool add_p_level)
599 : {
600 1617 : libmesh_assert(elem);
601 :
602 : return FE<1,BERNSTEIN>::shape_second_deriv
603 21021 : (elem->type(),
604 21021 : order + add_p_level*elem->p_level(), i, j, p);
605 : }
606 :
607 : template <>
608 0 : Real FE<1,BERNSTEIN>::shape_second_deriv(const FEType fet,
609 : const Elem * elem,
610 : const unsigned int i,
611 : const unsigned int j,
612 : const Point & p,
613 : const bool add_p_level)
614 : {
615 0 : libmesh_assert(elem);
616 : return FE<1,BERNSTEIN>::shape_second_deriv
617 0 : (elem->type(),
618 0 : fet.order + add_p_level*elem->p_level(), i, j, p);
619 : }
620 :
621 : #endif
622 :
623 : } // namespace libMesh
624 :
625 :
626 : #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
|
