libMesh
fe_bernstein_shape_1D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 <>
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  // Using an unsigned long here will work for any of the orders we support.
163  unsigned long binomial_p_i = 1;
164 
165  // the binomial coefficient (p choose n)
166  if (i>1)
167  binomial_p_i = Utility::binomial(static_cast<unsigned long>(p_order),
168  static_cast<unsigned long>(n));
169 
170  switch(i)
171  {
172  case 0:
173  return binomial_p_i * std::pow((1-xi)/2, p_order);
174  case 1:
175  return binomial_p_i * std::pow((1+xi)/2, p_order);
176  default:
177  {
178  return binomial_p_i * std::pow((1+xi)/2,n)
179  * std::pow((1-xi)/2,m);
180  }
181  }
182  }
183  }
184 
185  libmesh_error_msg("We'll never get here!");
186  return 0.;
187 }
188 
189 
190 
191 template <>
193  const Order order,
194  const unsigned int i,
195  const Point & p)
196 {
197  libmesh_assert(elem);
198 
199  return FE<1,BERNSTEIN>::shape(elem->type(), static_cast<Order>(order + elem->p_level()), i, p);
200 }
201 
202 
203 
204 template <>
206  const Order order,
207  const unsigned int i,
208  const unsigned int libmesh_dbg_var(j),
209  const Point & p)
210 {
211  // only d()/dxi in 1D!
212 
213  libmesh_assert_equal_to (j, 0);
214 
215  const Real xi = p(0);
216 
217  using Utility::pow;
218 
219  switch (order)
220  {
221  case FIRST:
222 
223  switch(i)
224  {
225  case 0:
226  return -.5;
227  case 1:
228  return .5;
229  default:
230  libmesh_error_msg("Invalid shape function index i = " << i);
231  }
232 
233  case SECOND:
234 
235  switch(i)
236  {
237  case 0:
238  return (xi-1.)*.5;
239  case 1:
240  return (xi+1.)*.5;
241  case 2:
242  return -xi;
243  default:
244  libmesh_error_msg("Invalid shape function index i = " << i);
245  }
246 
247  case THIRD:
248 
249  switch(i)
250  {
251  case 0:
252  return -0.375*pow<2>(1.-xi);
253  case 1:
254  return 0.375*pow<2>(1.+xi);
255  case 2:
256  return -0.375 -.75*xi +1.125*pow<2>(xi);
257  case 3:
258  return 0.375 -.75*xi -1.125*pow<2>(xi);
259  default:
260  libmesh_error_msg("Invalid shape function index i = " << i);
261  }
262 
263  case FOURTH:
264 
265  switch(i)
266  {
267  case 0:
268  return -0.25*pow<3>(1.-xi);
269  case 1:
270  return 0.25*pow<3>(1.+xi);
271  case 2:
272  return -0.5 +1.5*pow<2>(xi)-pow<3>(xi);
273  case 3:
274  return 1.5*(pow<3>(xi)-xi);
275  case 4:
276  return 0.5 -1.5*pow<2>(xi)-pow<3>(xi);
277  default:
278  libmesh_error_msg("Invalid shape function index i = " << i);
279  }
280 
281  case FIFTH:
282 
283  switch(i)
284  {
285  case 0:
286  return -(5./32.)*pow<4>(xi-1.);
287  case 1:
288  return (5./32.)*pow<4>(xi+1.);
289  case 2:
290  return (5./32.)*pow<4>(1.-xi) -(5./8.)*(1.+xi)*pow<3>(1.-xi);
291  case 3:
292  return (5./ 8.)*(1.+xi)*pow<3>(1.-xi) -(15./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
293  case 4:
294  return -(5./ 8.)*pow<3>(1.+xi)*(1.-xi) +(15./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
295  case 5:
296  return (5./ 8.)*pow<3>(1.+xi)*(1.-xi) -(5./32.)*pow<4>(1.+xi);
297  default:
298  libmesh_error_msg("Invalid shape function index i = " << i);
299  }
300 
301  case SIXTH:
302 
303  switch(i)
304  {
305  case 0:
306  return -( 3./32.)*pow<5>(1.-xi);
307  case 1:
308  return ( 3./32.)*pow<5>(1.+xi);
309  case 2:
310  return ( 3./32.)*pow<5>(1.-xi)-(15./32.)*(1.+xi)*pow<4>(1.-xi);
311  case 3:
312  return (15./32.)*(1.+xi)*pow<4>(1.-xi)-(15./16.)*pow<2>(1.+xi)*pow<3>(1.-xi);
313  case 4:
314  return -(15./ 8.)*xi +(15./4.)*pow<3>(xi)-(15./8.)*pow<5>(xi);
315  case 5:
316  return -(15./32.)*(1.-xi)*pow<4>(1.+xi)+(15./16.)*pow<2>(1.-xi)*pow<3>(1.+xi);
317  case 6:
318  return (15./32.)*pow<4>(1.+xi)*(1.-xi)-(3./32.)*pow<5>(1.+xi);
319  default:
320  libmesh_error_msg("Invalid shape function index i = " << i);
321  }
322 
323 
324  default:
325  {
326  libmesh_assert (order>6);
327 
328  // Use this for arbitrary orders
329  const int p_order = static_cast<int>(order);
330  const int m = p_order-(i-1);
331  const int n = (i-1);
332 
333  // Using an unsigned long here will work for any of the orders we support.
334  unsigned long binomial_p_i = 1;
335 
336  // the binomial coefficient (p choose n)
337  if (i>1)
338  binomial_p_i = Utility::binomial(static_cast<unsigned long>(p_order),
339  static_cast<unsigned long>(n));
340 
341  switch(i)
342  {
343  case 0:
344  return binomial_p_i * (-1./2.) * p_order * std::pow((1-xi)/2, p_order-1);
345  case 1:
346  return binomial_p_i * ( 1./2.) * p_order * std::pow((1+xi)/2, p_order-1);
347 
348  default:
349  {
350  return binomial_p_i * (1./2. * n * std::pow((1+xi)/2,n-1) * std::pow((1-xi)/2,m)
351  - 1./2. * m * std::pow((1+xi)/2,n) * std::pow((1-xi)/2,m-1));
352  }
353  }
354  }
355 
356  }
357 
358  libmesh_error_msg("We'll never get here!");
359  return 0.;
360 }
361 
362 
363 
364 template <>
366  const Order order,
367  const unsigned int i,
368  const unsigned int j,
369  const Point & p)
370 {
371  libmesh_assert(elem);
372 
373  return FE<1,BERNSTEIN>::shape_deriv(elem->type(),
374  static_cast<Order>(order + elem->p_level()), i, j, p);
375 }
376 
377 
378 
379 template <>
381  const Order,
382  const unsigned int,
383  const unsigned int,
384  const Point &)
385 {
386  static bool warning_given = false;
387 
388  if (!warning_given)
389  libMesh::err << "Second derivatives for Bernstein elements "
390  << "are not yet implemented!"
391  << std::endl;
392 
393  warning_given = true;
394  return 0.;
395 }
396 
397 
398 
399 
400 template <>
402  const Order,
403  const unsigned int,
404  const unsigned int,
405  const Point &)
406 {
407  static bool warning_given = false;
408 
409  if (!warning_given)
410  libMesh::err << "Second derivatives for Bernstein elements "
411  << "are not yet implemented!"
412  << std::endl;
413 
414  warning_given = true;
415  return 0.;
416 }
417 
418 } // namespace libMesh
419 
420 
421 #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
OStreamProxy err
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
virtual ElemType type() const =0
unsigned int p_level() const
Definition: elem.h:2422
ElemType
Defines an enum for geometric element types.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
The libMesh namespace provides an interface to certain functionality in the library.
libmesh_assert(j)
T pow(const T &x)
Definition: utility.h:195
PetscErrorCode Vec Mat libmesh_dbg_var(j)
double pow(double a, int b)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
T binomial(T n, T k)
Definition: utility.h:221