1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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 LIBMESH_DEFAULT_VECTORIZED_FE(1,BERNSTEIN)
38
39
40 template <>
shape(const ElemType,const Order order,const unsigned int i,const Point & p)41 Real FE<1,BERNSTEIN>::shape(const ElemType,
42 const Order order,
43 const unsigned int i,
44 const Point & p)
45 {
46 const Real xi = p(0);
47 using Utility::pow;
48
49 switch (order)
50 {
51 case FIRST:
52
53 switch(i)
54 {
55 case 0:
56 return (1.-xi)/2.;
57 case 1:
58 return (1.+xi)/2.;
59 default:
60 libmesh_error_msg("Invalid shape function index i = " << i);
61 }
62
63 case SECOND:
64
65 switch(i)
66 {
67 case 0:
68 return (1./4.)*pow<2>(1.-xi);
69 case 1:
70 return (1./4.)*pow<2>(1.+xi);
71 case 2:
72 return (1./2.)*(1.-xi)*(1.+xi);
73 default:
74 libmesh_error_msg("Invalid shape function index i = " << i);
75 }
76
77 case THIRD:
78
79 switch(i)
80 {
81 case 0:
82 return (1./8.)*pow<3>(1.-xi);
83 case 1:
84 return (1./8.)*pow<3>(1.+xi);
85 case 2:
86 return (3./8.)*(1.+xi)*pow<2>(1.-xi);
87 case 3:
88 return (3./8.)*pow<2>(1.+xi)*(1.-xi);
89 default:
90 libmesh_error_msg("Invalid shape function index i = " << i);
91 }
92
93 case FOURTH:
94
95 switch(i)
96 {
97 case 0:
98 return (1./16.)*pow<4>(1.-xi);
99 case 1:
100 return (1./16.)*pow<4>(1.+xi);
101 case 2:
102 return (1./ 4.)*(1.+xi)*pow<3>(1.-xi);
103 case 3:
104 return (3./ 8.)*pow<2>(1.+xi)*pow<2>(1.-xi);
105 case 4:
106 return (1./ 4.)*pow<3>(1.+xi)*(1.-xi);
107 default:
108 libmesh_error_msg("Invalid shape function index i = " << i);
109 }
110
111
112 case FIFTH:
113
114 switch(i)
115 {
116 case 0:
117 return (1./32.)*pow<5>(1.-xi);
118 case 1:
119 return (1./32.)*pow<5>(1.+xi);
120 case 2:
121 return (5./32.)*(1.+xi)*pow<4>(1.-xi);
122 case 3:
123 return (5./16.)*pow<2>(1.+xi)*pow<3>(1.-xi);
124 case 4:
125 return (5./16.)*pow<3>(1.+xi)*pow<2>(1.-xi);
126 case 5:
127 return (5./32.)*pow<4>(1.+xi)*(1.-xi);
128 default:
129 libmesh_error_msg("Invalid shape function index i = " << i);
130 }
131
132
133 case SIXTH:
134
135 switch (i)
136 {
137 case 0:
138 return ( 1./64.)*pow<6>(1.-xi);
139 case 1:
140 return ( 1./64.)*pow<6>(1.+xi);
141 case 2:
142 return ( 3./32.)*(1.+xi)*pow<5>(1.-xi);
143 case 3:
144 return (15./64.)*pow<2>(1.+xi)*pow<4>(1.-xi);
145 case 4:
146 return ( 5./16.)*pow<3>(1.+xi)*pow<3>(1.-xi);
147 case 5:
148 return (15./64.)*pow<4>(1.+xi)*pow<2>(1.-xi);
149 case 6:
150 return ( 3./32.)*pow<5>(1.+xi)*(1.-xi);
151 default:
152 libmesh_error_msg("Invalid shape function index i = " << i);
153 }
154
155 default:
156 {
157 libmesh_assert (order>6);
158
159 // Use this for arbitrary orders.
160 // Note that this implementation is less efficient.
161 const int p_order = static_cast<int>(order);
162 const int m = p_order-i+1;
163 const int n = (i-1);
164
165 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 if (i>1)
171 binomial_p_i = Real(Utility::binomial(static_cast<unsigned long>(p_order),
172 static_cast<unsigned long>(n)));
173
174 switch(i)
175 {
176 case 0:
177 return binomial_p_i * std::pow((1-xi)/2, p_order);
178 case 1:
179 return binomial_p_i * std::pow((1+xi)/2, p_order);
180 default:
181 {
182 return binomial_p_i * std::pow((1+xi)/2,n)
183 * std::pow((1-xi)/2,m);
184 }
185 }
186 }
187 }
188 }
189
190
191
192 template <>
shape(const Elem * elem,const Order order,const unsigned int i,const Point & p,const bool add_p_level)193 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 libmesh_assert(elem);
200
201 return FE<1,BERNSTEIN>::shape
202 (elem->type(),
203 static_cast<Order>(order + add_p_level*elem->p_level()), i, p);
204 }
205
206
207 template <>
shape(const FEType fet,const Elem * elem,const unsigned int i,const Point & p,const bool add_p_level)208 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 libmesh_assert(elem);
215 return FE<1,BERNSTEIN>::shape
216 (elem->type(),
217 static_cast<Order>(fet.order + add_p_level * elem->p_level()), i, p);
218 }
219
220
221 template <>
shape_deriv(const ElemType,const Order order,const unsigned int i,const unsigned int libmesh_dbg_var (j),const Point & p)222 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 libmesh_assert_equal_to (j, 0);
231
232 const Real xi = p(0);
233
234 using Utility::pow;
235
236 switch (order)
237 {
238 case FIRST:
239
240 switch(i)
241 {
242 case 0:
243 return -.5;
244 case 1:
245 return .5;
246 default:
247 libmesh_error_msg("Invalid shape function index i = " << i);
248 }
249
250 case SECOND:
251
252 switch(i)
253 {
254 case 0:
255 return (xi-1.)*.5;
256 case 1:
257 return (xi+1.)*.5;
258 case 2:
259 return -xi;
260 default:
261 libmesh_error_msg("Invalid shape function index i = " << i);
262 }
263
264 case THIRD:
265
266 switch(i)
267 {
268 case 0:
269 return -0.375*pow<2>(1.-xi);
270 case 1:
271 return 0.375*pow<2>(1.+xi);
272 case 2:
273 return -0.375 -.75*xi +1.125*pow<2>(xi);
274 case 3:
275 return 0.375 -.75*xi -1.125*pow<2>(xi);
276 default:
277 libmesh_error_msg("Invalid shape function index i = " << i);
278 }
279
280 case FOURTH:
281
282 switch(i)
283 {
284 case 0:
285 return -0.25*pow<3>(1.-xi);
286 case 1:
287 return 0.25*pow<3>(1.+xi);
288 case 2:
289 return -0.5 +1.5*pow<2>(xi)-pow<3>(xi);
290 case 3:
291 return 1.5*(pow<3>(xi)-xi);
292 case 4:
293 return 0.5 -1.5*pow<2>(xi)-pow<3>(xi);
294 default:
295 libmesh_error_msg("Invalid shape function index i = " << i);
296 }
297
298 case FIFTH:
299
300 switch(i)
301 {
302 case 0:
303 return -(5./32.)*pow<4>(xi-1.);
304 case 1:
305 return (5./32.)*pow<4>(xi+1.);
306 case 2:
307 return (5./32.)*pow<4>(1.-xi) -(5./8.)*(1.+xi)*pow<3>(1.-xi);
308 case 3:
309 return (5./ 8.)*(1.+xi)*pow<3>(1.-xi) -(15./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
310 case 4:
311 return -(5./ 8.)*pow<3>(1.+xi)*(1.-xi) +(15./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
312 case 5:
313 return (5./ 8.)*pow<3>(1.+xi)*(1.-xi) -(5./32.)*pow<4>(1.+xi);
314 default:
315 libmesh_error_msg("Invalid shape function index i = " << i);
316 }
317
318 case SIXTH:
319
320 switch(i)
321 {
322 case 0:
323 return -( 3./32.)*pow<5>(1.-xi);
324 case 1:
325 return ( 3./32.)*pow<5>(1.+xi);
326 case 2:
327 return ( 3./32.)*pow<5>(1.-xi)-(15./32.)*(1.+xi)*pow<4>(1.-xi);
328 case 3:
329 return (15./32.)*(1.+xi)*pow<4>(1.-xi)-(15./16.)*pow<2>(1.+xi)*pow<3>(1.-xi);
330 case 4:
331 return -(15./ 8.)*xi +(15./4.)*pow<3>(xi)-(15./8.)*pow<5>(xi);
332 case 5:
333 return -(15./32.)*(1.-xi)*pow<4>(1.+xi)+(15./16.)*pow<2>(1.-xi)*pow<3>(1.+xi);
334 case 6:
335 return (15./32.)*pow<4>(1.+xi)*(1.-xi)-(3./32.)*pow<5>(1.+xi);
336 default:
337 libmesh_error_msg("Invalid shape function index i = " << i);
338 }
339
340
341 default:
342 {
343 libmesh_assert (order>6);
344
345 // Use this for arbitrary orders
346 const int p_order = static_cast<int>(order);
347 const int m = p_order-(i-1);
348 const int n = (i-1);
349
350 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 if (i>1)
356 binomial_p_i = Real(Utility::binomial(static_cast<unsigned long>(p_order),
357 static_cast<unsigned long>(n)));
358
359 switch(i)
360 {
361 case 0:
362 return binomial_p_i * (-1./2.) * p_order * std::pow((1-xi)/2, p_order-1);
363 case 1:
364 return binomial_p_i * ( 1./2.) * p_order * std::pow((1+xi)/2, p_order-1);
365
366 default:
367 {
368 return binomial_p_i * (1./2. * n * std::pow((1+xi)/2,n-1) * std::pow((1-xi)/2,m)
369 - 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 <>
shape_deriv(const Elem * elem,const Order order,const unsigned int i,const unsigned int j,const Point & p,const bool add_p_level)380 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 libmesh_assert(elem);
388
389 return FE<1,BERNSTEIN>::shape_deriv
390 (elem->type(),
391 static_cast<Order>(order + add_p_level*elem->p_level()), i, j, p);
392 }
393
394 template <>
shape_deriv(const FEType fet,const Elem * elem,const unsigned int i,const unsigned int j,const Point & p,const bool add_p_level)395 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 libmesh_assert(elem);
403 return FE<1,BERNSTEIN>::shape_deriv
404 (elem->type(),
405 static_cast<Order>(fet.order + add_p_level * elem->p_level()), i, j, p);
406 }
407
408
409 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
410
411 template <>
shape_second_deriv(const ElemType,const Order order,const unsigned int i,const unsigned int libmesh_dbg_var (j),const Point & p)412 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 libmesh_assert_equal_to (j, 0);
421
422 const Real xi = p(0);
423
424 using Utility::pow;
425
426 switch (order)
427 {
428 case FIRST:
429
430 switch(i)
431 {
432 case 0:
433 case 1:
434 return 0;
435 default:
436 libmesh_error_msg("Invalid shape function index i = " << i);
437 }
438
439 case SECOND:
440
441 switch(i)
442 {
443 case 0:
444 case 1:
445 return .5;
446 case 2:
447 return -1;
448 default:
449 libmesh_error_msg("Invalid shape function index i = " << i);
450 }
451
452 case THIRD:
453
454 switch(i)
455 {
456 case 0:
457 return 0.75*(1.-xi);
458 case 1:
459 return 0.75*(1.+xi);
460 case 2:
461 return -.75 + 2.25*xi;
462 case 3:
463 return -.75 - 2.25*xi;
464 default:
465 libmesh_error_msg("Invalid shape function index i = " << i);
466 }
467
468 case FOURTH:
469
470 switch(i)
471 {
472 case 0:
473 return 0.75*pow<2>(1.-xi);
474 case 1:
475 return 0.75*pow<2>(1.+xi);
476 case 2:
477 return 3*(xi - pow<2>(xi));
478 case 3:
479 return 1.5*(3*pow<2>(xi)-1);
480 case 4:
481 return -3*xi-3*pow<2>(xi);
482 default:
483 libmesh_error_msg("Invalid shape function index i = " << i);
484 }
485
486 case FIFTH:
487
488 switch(i)
489 {
490 case 0:
491 return -(5./8.)*pow<3>(xi-1.);
492 case 1:
493 return (5./8.)*pow<3>(xi+1.);
494 case 2:
495 return -(5./4.)*pow<3>(1.-xi) + (15./8.)*(1.+xi)*pow<2>(1.-xi);
496 case 3:
497 return -(15./ 4.)*(1.+xi)*pow<2>(1.-xi) + (5./ 8.)*pow<3>(1.-xi)
498 + (15./8.)*pow<2>(1.+xi)*(1.-xi);
499 case 4:
500 return (5./ 8.)*pow<3>(1.+xi) - (15./ 4.)*pow<2>(1.+xi)*(1.-xi)
501 +(15./8.)*(1.+xi)*pow<2>(1.-xi);
502 case 5:
503 return -(5./ 8.)*pow<3>(1.+xi) + (15./ 8.)*pow<2>(1.+xi)*(1.-xi)
504 -(5./8.)*pow<3>(1.+xi);
505 default:
506 libmesh_error_msg("Invalid shape function index i = " << i);
507 }
508
509 case SIXTH:
510
511 switch(i)
512 {
513 case 0:
514 return ( 15./32.)*pow<4>(1.-xi);
515 case 1:
516 return ( 15./32.)*pow<4>(1.+xi);
517 case 2:
518 return -( 15./8.)*pow<4>(1.-xi) +
519 ( 15./8.)*(1.+xi)*pow<3>(1.-xi);
520 case 3:
521 return -(15./4.)*(1.+xi)*pow<3>(1.-xi)
522 + (15./32.)*pow<4>(1.-xi)
523 + (45./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
524 case 4:
525 return -(15./ 8.) +(45./4.)*pow<2>(xi) - (75./8.)*pow<4>(xi);
526 case 5:
527 return -(15./4.)*(1.-xi)*pow<3>(1.+xi)
528 + (15./32.)*pow<4>(1.+xi)
529 + (45./16.)*pow<2>(1.-xi)*pow<2>(1.+xi);
530 case 6:
531 return -(15./16.)*pow<4>(1.+xi)
532 + (15./8.)*pow<3>(1.+xi)*(1.-xi);
533 default:
534 libmesh_error_msg("Invalid shape function index i = " << i);
535 }
536
537
538 default:
539 {
540 libmesh_assert (order>6);
541
542 // Use this for arbitrary orders
543 const int p_order = static_cast<int>(order);
544 const int m = p_order-(i-1);
545 const int n = (i-1);
546
547 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 if (i>1)
553 binomial_p_i = Real(Utility::binomial(static_cast<unsigned long>(p_order),
554 static_cast<unsigned long>(n)));
555
556 switch(i)
557 {
558 case 0:
559 return binomial_p_i * (1./4.) * p_order * (p_order-1) * std::pow((1-xi)/2, p_order-2);
560 case 1:
561 return binomial_p_i * (1./4.) * p_order * (p_order-1) * std::pow((1+xi)/2, p_order-2);
562
563 default:
564 {
565 Real val = 0;
566
567 if (n == 1)
568 val +=
569 binomial_p_i * (-1./4. * m * std::pow((1-xi)/2,m-1));
570 else
571 val +=
572 binomial_p_i * (-1./4. * n * m * std::pow((1+xi)/2,n-1) * std::pow((1-xi)/2,m-1) +
573 1./4. * n * (n-1) * std::pow((1+xi)/2,n-2) * std::pow((1-xi)/2,m));
574
575 if (m == 1)
576 val += binomial_p_i * (-1./4. * n * std::pow((1+xi)/2,n-1));
577 else
578 val +=
579 binomial_p_i * (1./4. * m * (m-1) * std::pow((1+xi)/2,n) * std::pow((1-xi)/2,m-2)
580 - 1./4. * m * n * std::pow((1+xi)/2,n-1) * std::pow((1-xi)/2,m-1));
581
582 return val;
583 }
584 }
585 }
586
587 }
588 }
589
590
591
592 template <>
shape_second_deriv(const Elem * elem,const Order order,const unsigned int i,const unsigned int j,const Point & p,const bool add_p_level)593 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 libmesh_assert(elem);
601
602 return FE<1,BERNSTEIN>::shape_second_deriv
603 (elem->type(),
604 static_cast<Order>(order + add_p_level*elem->p_level()), i, j, p);
605 }
606
607 template <>
shape_second_deriv(const FEType fet,const Elem * elem,const unsigned int i,const unsigned int j,const Point & p,const bool add_p_level)608 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 libmesh_assert(elem);
616 return FE<1,BERNSTEIN>::shape_second_deriv
617 (elem->type(),
618 static_cast<Order>(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
627