libMesh
fe_bernstein_shape_1D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2019 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 template <>
38 Real FE<1,BERNSTEIN>::shape(const ElemType,
39  const Order order,
40  const unsigned int i,
41  const Point & p)
42 {
43  const Real xi = p(0);
44  using Utility::pow;
45 
46  switch (order)
47  {
48  case FIRST:
49 
50  switch(i)
51  {
52  case 0:
53  return (1.-xi)/2.;
54  case 1:
55  return (1.+xi)/2.;
56  default:
57  libmesh_error_msg("Invalid shape function index i = " << i);
58  }
59 
60  case SECOND:
61 
62  switch(i)
63  {
64  case 0:
65  return (1./4.)*pow<2>(1.-xi);
66  case 1:
67  return (1./4.)*pow<2>(1.+xi);
68  case 2:
69  return (1./2.)*(1.-xi)*(1.+xi);
70  default:
71  libmesh_error_msg("Invalid shape function index i = " << i);
72  }
73 
74  case THIRD:
75 
76  switch(i)
77  {
78  case 0:
79  return (1./8.)*pow<3>(1.-xi);
80  case 1:
81  return (1./8.)*pow<3>(1.+xi);
82  case 2:
83  return (3./8.)*(1.+xi)*pow<2>(1.-xi);
84  case 3:
85  return (3./8.)*pow<2>(1.+xi)*(1.-xi);
86  default:
87  libmesh_error_msg("Invalid shape function index i = " << i);
88  }
89 
90  case FOURTH:
91 
92  switch(i)
93  {
94  case 0:
95  return (1./16.)*pow<4>(1.-xi);
96  case 1:
97  return (1./16.)*pow<4>(1.+xi);
98  case 2:
99  return (1./ 4.)*(1.+xi)*pow<3>(1.-xi);
100  case 3:
101  return (3./ 8.)*pow<2>(1.+xi)*pow<2>(1.-xi);
102  case 4:
103  return (1./ 4.)*pow<3>(1.+xi)*(1.-xi);
104  default:
105  libmesh_error_msg("Invalid shape function index i = " << i);
106  }
107 
108 
109  case FIFTH:
110 
111  switch(i)
112  {
113  case 0:
114  return (1./32.)*pow<5>(1.-xi);
115  case 1:
116  return (1./32.)*pow<5>(1.+xi);
117  case 2:
118  return (5./32.)*(1.+xi)*pow<4>(1.-xi);
119  case 3:
120  return (5./16.)*pow<2>(1.+xi)*pow<3>(1.-xi);
121  case 4:
122  return (5./16.)*pow<3>(1.+xi)*pow<2>(1.-xi);
123  case 5:
124  return (5./32.)*pow<4>(1.+xi)*(1.-xi);
125  default:
126  libmesh_error_msg("Invalid shape function index i = " << i);
127  }
128 
129 
130  case SIXTH:
131 
132  switch (i)
133  {
134  case 0:
135  return ( 1./64.)*pow<6>(1.-xi);
136  case 1:
137  return ( 1./64.)*pow<6>(1.+xi);
138  case 2:
139  return ( 3./32.)*(1.+xi)*pow<5>(1.-xi);
140  case 3:
141  return (15./64.)*pow<2>(1.+xi)*pow<4>(1.-xi);
142  case 4:
143  return ( 5./16.)*pow<3>(1.+xi)*pow<3>(1.-xi);
144  case 5:
145  return (15./64.)*pow<4>(1.+xi)*pow<2>(1.-xi);
146  case 6:
147  return ( 3./32.)*pow<5>(1.+xi)*(1.-xi);
148  default:
149  libmesh_error_msg("Invalid shape function index i = " << i);
150  }
151 
152  default:
153  {
154  libmesh_assert (order>6);
155 
156  // Use this for arbitrary orders.
157  // Note that this implementation is less efficient.
158  const int p_order = static_cast<int>(order);
159  const int m = p_order-i+1;
160  const int n = (i-1);
161 
162  Real binomial_p_i = 1;
163 
164  // the binomial coefficient (p choose n)
165  // Using an unsigned long here will work for any of the orders we support.
166  // Explicitly construct a Real to prevent conversion warnings
167  if (i>1)
168  binomial_p_i = Real(Utility::binomial(static_cast<unsigned long>(p_order),
169  static_cast<unsigned long>(n)));
170 
171  switch(i)
172  {
173  case 0:
174  return binomial_p_i * std::pow((1-xi)/2, p_order);
175  case 1:
176  return binomial_p_i * std::pow((1+xi)/2, p_order);
177  default:
178  {
179  return binomial_p_i * std::pow((1+xi)/2,n)
180  * std::pow((1-xi)/2,m);
181  }
182  }
183  }
184  }
185 }
186 
187 
188 
189 template <>
190 Real FE<1,BERNSTEIN>::shape(const Elem * elem,
191  const Order order,
192  const unsigned int i,
193  const Point & p,
194  const bool add_p_level)
195 {
196  libmesh_assert(elem);
197 
198  return FE<1,BERNSTEIN>::shape
199  (elem->type(),
200  static_cast<Order>(order + add_p_level*elem->p_level()), i, p);
201 }
202 
203 
204 
205 template <>
206 Real FE<1,BERNSTEIN>::shape_deriv(const ElemType,
207  const Order order,
208  const unsigned int i,
209  const unsigned int libmesh_dbg_var(j),
210  const Point & p)
211 {
212  // only d()/dxi in 1D!
213 
214  libmesh_assert_equal_to (j, 0);
215 
216  const Real xi = p(0);
217 
218  using Utility::pow;
219 
220  switch (order)
221  {
222  case FIRST:
223 
224  switch(i)
225  {
226  case 0:
227  return -.5;
228  case 1:
229  return .5;
230  default:
231  libmesh_error_msg("Invalid shape function index i = " << i);
232  }
233 
234  case SECOND:
235 
236  switch(i)
237  {
238  case 0:
239  return (xi-1.)*.5;
240  case 1:
241  return (xi+1.)*.5;
242  case 2:
243  return -xi;
244  default:
245  libmesh_error_msg("Invalid shape function index i = " << i);
246  }
247 
248  case THIRD:
249 
250  switch(i)
251  {
252  case 0:
253  return -0.375*pow<2>(1.-xi);
254  case 1:
255  return 0.375*pow<2>(1.+xi);
256  case 2:
257  return -0.375 -.75*xi +1.125*pow<2>(xi);
258  case 3:
259  return 0.375 -.75*xi -1.125*pow<2>(xi);
260  default:
261  libmesh_error_msg("Invalid shape function index i = " << i);
262  }
263 
264  case FOURTH:
265 
266  switch(i)
267  {
268  case 0:
269  return -0.25*pow<3>(1.-xi);
270  case 1:
271  return 0.25*pow<3>(1.+xi);
272  case 2:
273  return -0.5 +1.5*pow<2>(xi)-pow<3>(xi);
274  case 3:
275  return 1.5*(pow<3>(xi)-xi);
276  case 4:
277  return 0.5 -1.5*pow<2>(xi)-pow<3>(xi);
278  default:
279  libmesh_error_msg("Invalid shape function index i = " << i);
280  }
281 
282  case FIFTH:
283 
284  switch(i)
285  {
286  case 0:
287  return -(5./32.)*pow<4>(xi-1.);
288  case 1:
289  return (5./32.)*pow<4>(xi+1.);
290  case 2:
291  return (5./32.)*pow<4>(1.-xi) -(5./8.)*(1.+xi)*pow<3>(1.-xi);
292  case 3:
293  return (5./ 8.)*(1.+xi)*pow<3>(1.-xi) -(15./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
294  case 4:
295  return -(5./ 8.)*pow<3>(1.+xi)*(1.-xi) +(15./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
296  case 5:
297  return (5./ 8.)*pow<3>(1.+xi)*(1.-xi) -(5./32.)*pow<4>(1.+xi);
298  default:
299  libmesh_error_msg("Invalid shape function index i = " << i);
300  }
301 
302  case SIXTH:
303 
304  switch(i)
305  {
306  case 0:
307  return -( 3./32.)*pow<5>(1.-xi);
308  case 1:
309  return ( 3./32.)*pow<5>(1.+xi);
310  case 2:
311  return ( 3./32.)*pow<5>(1.-xi)-(15./32.)*(1.+xi)*pow<4>(1.-xi);
312  case 3:
313  return (15./32.)*(1.+xi)*pow<4>(1.-xi)-(15./16.)*pow<2>(1.+xi)*pow<3>(1.-xi);
314  case 4:
315  return -(15./ 8.)*xi +(15./4.)*pow<3>(xi)-(15./8.)*pow<5>(xi);
316  case 5:
317  return -(15./32.)*(1.-xi)*pow<4>(1.+xi)+(15./16.)*pow<2>(1.-xi)*pow<3>(1.+xi);
318  case 6:
319  return (15./32.)*pow<4>(1.+xi)*(1.-xi)-(3./32.)*pow<5>(1.+xi);
320  default:
321  libmesh_error_msg("Invalid shape function index i = " << i);
322  }
323 
324 
325  default:
326  {
327  libmesh_assert (order>6);
328 
329  // Use this for arbitrary orders
330  const int p_order = static_cast<int>(order);
331  const int m = p_order-(i-1);
332  const int n = (i-1);
333 
334  Real binomial_p_i = 1;
335 
336  // the binomial coefficient (p choose n)
337  // Using an unsigned long here will work for any of the orders we support.
338  // Explicitly construct a Real to prevent conversion warnings
339  if (i>1)
340  binomial_p_i = Real(Utility::binomial(static_cast<unsigned long>(p_order),
341  static_cast<unsigned long>(n)));
342 
343  switch(i)
344  {
345  case 0:
346  return binomial_p_i * (-1./2.) * p_order * std::pow((1-xi)/2, p_order-1);
347  case 1:
348  return binomial_p_i * ( 1./2.) * p_order * std::pow((1+xi)/2, p_order-1);
349 
350  default:
351  {
352  return binomial_p_i * (1./2. * n * std::pow((1+xi)/2,n-1) * std::pow((1-xi)/2,m)
353  - 1./2. * m * std::pow((1+xi)/2,n) * std::pow((1-xi)/2,m-1));
354  }
355  }
356  }
357 
358  }
359 }
360 
361 
362 
363 template <>
364 Real FE<1,BERNSTEIN>::shape_deriv(const Elem * elem,
365  const Order order,
366  const unsigned int i,
367  const unsigned int j,
368  const Point & p,
369  const bool add_p_level)
370 {
371  libmesh_assert(elem);
372 
373  return FE<1,BERNSTEIN>::shape_deriv
374  (elem->type(),
375  static_cast<Order>(order + add_p_level*elem->p_level()), i, j, p);
376 }
377 
378 
379 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
380 
381 template <>
382 Real FE<1,BERNSTEIN>::shape_second_deriv(const ElemType,
383  const Order,
384  const unsigned int,
385  const unsigned int,
386  const Point &)
387 {
388  static bool warning_given = false;
389 
390  if (!warning_given)
391  libMesh::err << "Second derivatives for Bernstein elements "
392  << "are not yet implemented!"
393  << std::endl;
394 
395  warning_given = true;
396  return 0.;
397 }
398 
399 
400 
401 
402 template <>
403 Real FE<1,BERNSTEIN>::shape_second_deriv(const Elem *,
404  const Order,
405  const unsigned int,
406  const unsigned int,
407  const Point &,
408  const bool)
409 {
410  static bool warning_given = false;
411 
412  if (!warning_given)
413  libMesh::err << "Second derivatives for Bernstein elements "
414  << "are not yet implemented!"
415  << std::endl;
416 
417  warning_given = true;
418  return 0.;
419 }
420 
421 #endif
422 
423 } // namespace libMesh
424 
425 
426 #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
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:272
libMesh::SIXTH
Definition: enum_order.h:47
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
Definition: factoryfunction.C:55
libMesh::Order
Order
Definition: enum_order.h:40
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2512
libMesh::FIFTH
Definition: enum_order.h:46
libMesh::SECOND
Definition: enum_order.h:43
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
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::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
std::pow
double pow(double a, int b)
Definition: libmesh_augment_std_namespace.h:58
libMesh::FOURTH
Definition: enum_order.h:45
libMesh::Utility::pow
T pow(const T &x)
Definition: utility.h:246
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:100
libMesh::THIRD
Definition: enum_order.h:44
libMesh::err
OStreamProxy err
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::FIRST
Definition: enum_order.h:42
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
Generated on Sat Jan 25 2020 12:06:57 for libMesh by   doxygen 1.8.16