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