libMesh
fe_bernstein_shape_1D.C
Go to the documentation of this file.
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 LIBMESH_DEFAULT_VECTORIZED_FE(1,BERNSTEIN)
38 
39 
40 template <>
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 <>
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  order + add_p_level*elem->p_level(), i, p);
204 }
205 
206 
207 template <>
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  fet.order + add_p_level*elem->p_level(), i, p);
218 }
219 
220 
221 template <>
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 <>
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  order + add_p_level*elem->p_level(), i, j, p);
392 }
393 
394 template <>
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  fet.order + add_p_level*elem->p_level(), i, j, p);
406 }
407 
408 
409 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
410 
411 template <>
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 <>
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  order + add_p_level*elem->p_level(), i, j, p);
605 }
606 
607 template <>
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  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
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::Order
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
libMesh::FE::shape_deriv
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:3109
libMesh::FEType::order
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:215
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::LIBMESH_DEFAULT_VECTORIZED_FE
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
Definition: fe_bernstein_shape_0D.C:30
libMesh::Utility::pow
T pow(const T &x)
Definition: utility.h:328
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::Utility::binomial
T binomial(T n, T k)
Definition: utility.h:354
Generated on Thu Oct 30 2025 06:31:47 for libMesh by   doxygen 1.8.14