fe_bernstein_shape_1D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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::err << "Invalid shape function index!" << std::endl;
58  libmesh_error();
59  }
60 
61  case SECOND:
62 
63  switch(i)
64  {
65  case 0:
66  return (1./4.)*pow<2>(1.-xi);
67  case 1:
68  return (1./4.)*pow<2>(1.+xi);
69  case 2:
70  return (1./2.)*(1.-xi)*(1.+xi);
71  default:
72  libMesh::err << "Invalid shape function index!" << std::endl;
73  libmesh_error();
74  }
75 
76  case THIRD:
77 
78  switch(i)
79  {
80  case 0:
81  return (1./8.)*pow<3>(1.-xi);
82  case 1:
83  return (1./8.)*pow<3>(1.+xi);
84  case 2:
85  return (3./8.)*(1.+xi)*pow<2>(1.-xi);
86  case 3:
87  return (3./8.)*pow<2>(1.+xi)*(1.-xi);
88  default:
89  libMesh::err << "Invalid shape function index!" << std::endl;
90  libmesh_error();
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::err << "Invalid shape function index!" << std::endl;
109  libmesh_error();
110  }
111 
112 
113  case FIFTH:
114 
115  switch(i)
116  {
117  case 0:
118  return (1./32.)*pow<5>(1.-xi);
119  case 1:
120  return (1./32.)*pow<5>(1.+xi);
121  case 2:
122  return (5./32.)*(1.+xi)*pow<4>(1.-xi);
123  case 3:
124  return (5./16.)*pow<2>(1.+xi)*pow<3>(1.-xi);
125  case 4:
126  return (5./16.)*pow<3>(1.+xi)*pow<2>(1.-xi);
127  case 5:
128  return (5./32.)*pow<4>(1.+xi)*(1.-xi);
129  default:
130  libMesh::err << "Invalid shape function index!" << std::endl;
131  libmesh_error();
132  }
133 
134 
135  case SIXTH:
136 
137  switch (i)
138  {
139  case 0:
140  return ( 1./64.)*pow<6>(1.-xi);
141  case 1:
142  return ( 1./64.)*pow<6>(1.+xi);
143  case 2:
144  return ( 3./32.)*(1.+xi)*pow<5>(1.-xi);
145  case 3:
146  return (15./64.)*pow<2>(1.+xi)*pow<4>(1.-xi);
147  case 4:
148  return ( 5./16.)*pow<3>(1.+xi)*pow<3>(1.-xi);
149  case 5:
150  return (15./64.)*pow<4>(1.+xi)*pow<2>(1.-xi);
151  case 6:
152  return ( 3./32.)*pow<5>(1.+xi)*(1.-xi);
153  default:
154  libMesh::err << "Invalid shape function index!" << std::endl;
155  libmesh_error();
156  }
157 
158  default:
159  {
160  libmesh_assert (order>6);
161 
162  // Use this for arbitrary orders.
163  // Note that this implementation is less efficient.
164  const int p_order = static_cast<unsigned int>(order);
165  const int m = p_order-i+1;
166  const int n = (i-1);
167 
168  unsigned int binomial_p_i = 1;
169 
170  // the binomial coefficient (p choose n)
171  if (i>1)
172  binomial_p_i = Utility::factorial(p_order)
173  / (Utility::factorial(n)*Utility::factorial(p_order-n));
174 
175 
176  switch(i)
177  {
178  case 0:
179  return binomial_p_i * std::pow((1-xi)/2,static_cast<int>(p_order));
180  case 1:
181  return binomial_p_i * std::pow((1+xi)/2,static_cast<int>(p_order));
182  default:
183  {
184  return binomial_p_i * std::pow((1+xi)/2,n)
185  * std::pow((1-xi)/2,m);
186  }
187  }
188 
189  // we should never get here
190  libmesh_error();
191  }
192  }
193 
194  libmesh_error();
195  return 0.;
196 }
197 
198 
199 
200 template <>
202  const Order order,
203  const unsigned int i,
204  const Point& p)
205 {
206  libmesh_assert(elem);
207 
208  return FE<1,BERNSTEIN>::shape(elem->type(), static_cast<Order>(order + elem->p_level()), i, p);
209 }
210 
211 
212 
213 template <>
215  const Order order,
216  const unsigned int i,
217  const unsigned int libmesh_dbg_var(j),
218  const Point& p)
219 {
220  // only d()/dxi in 1D!
221 
222  libmesh_assert_equal_to (j, 0);
223 
224  const Real xi = p(0);
225 
226  using Utility::pow;
227 
228  switch (order)
229  {
230  case FIRST:
231 
232  switch(i)
233  {
234  case 0:
235  return -.5;
236  case 1:
237  return .5;
238  default:
239  libMesh::err << "Invalid shape function index " << i << std::endl;
240  libmesh_error();
241  }
242 
243  case SECOND:
244 
245  switch(i)
246  {
247  case 0:
248  return (xi-1.)*.5;
249  case 1:
250  return (xi+1.)*.5;
251  case 2:
252  return -xi;
253  default:
254  libMesh::err << "Invalid shape function index!" << std::endl;
255  libmesh_error();
256  }
257 
258  case THIRD:
259 
260  switch(i)
261  {
262  case 0:
263  return -0.375*pow<2>(1.-xi);
264  case 1:
265  return 0.375*pow<2>(1.+xi);
266  case 2:
267  return -0.375 -.75*xi +1.125*pow<2>(xi);
268  case 3:
269  return 0.375 -.75*xi -1.125*pow<2>(xi);
270  default:
271  libMesh::err << "Invalid shape function index!" << std::endl;
272  libmesh_error();
273  }
274 
275  case FOURTH:
276 
277  switch(i)
278  {
279  case 0:
280  return -0.25*pow<3>(1.-xi);
281  case 1:
282  return 0.25*pow<3>(1.+xi);
283  case 2:
284  return -0.5 +1.5*pow<2>(xi)-pow<3>(xi);
285  case 3:
286  return 1.5*(pow<3>(xi)-xi);
287  case 4:
288  return 0.5 -1.5*pow<2>(xi)-pow<3>(xi);
289  default:
290  libMesh::err << "Invalid shape function index!" << std::endl;
291  libmesh_error();
292  }
293 
294  case FIFTH:
295 
296  switch(i)
297  {
298  case 0:
299  return -(5./32.)*pow<4>(xi-1.);
300  case 1:
301  return (5./32.)*pow<4>(xi+1.);
302  case 2:
303  return (5./32.)*pow<4>(1.-xi) -(5./8.)*(1.+xi)*pow<3>(1.-xi);
304  case 3:
305  return (5./ 8.)*(1.+xi)*pow<3>(1.-xi) -(15./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
306  case 4:
307  return -(5./ 8.)*pow<3>(1.+xi)*(1.-xi) +(15./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
308  case 5:
309  return (5./ 8.)*pow<3>(1.+xi)*(1.-xi) -(5./32.)*pow<4>(1.+xi);
310  default:
311  libMesh::err << "Invalid shape function index!" << std::endl;
312  libmesh_error();
313  }
314 
315  case SIXTH:
316 
317  switch(i)
318  {
319  case 0:
320  return -( 3./32.)*pow<5>(1.-xi);
321  case 1:
322  return ( 3./32.)*pow<5>(1.+xi);
323  case 2:
324  return ( 3./32.)*pow<5>(1.-xi)-(15./32.)*(1.+xi)*pow<4>(1.-xi);
325  case 3:
326  return (15./32.)*(1.+xi)*pow<4>(1.-xi)-(15./16.)*pow<2>(1.+xi)*pow<3>(1.-xi);
327  case 4:
328  return -(15./ 8.)*xi +(15./4.)*pow<3>(xi)-(15./8.)*pow<5>(xi);
329  case 5:
330  return -(15./32.)*(1.-xi)*pow<4>(1.+xi)+(15./16.)*pow<2>(1.-xi)*pow<3>(1.+xi);
331  case 6:
332  return (15./32.)*pow<4>(1.+xi)*(1.-xi)-(3./32.)*pow<5>(1.+xi);
333  default:
334  libMesh::err << "Invalid shape function index!" << std::endl;
335  libmesh_error();
336  }
337 
338 
339  default:
340  {
341  libmesh_assert (order>6);
342 
343  // Use this for arbitrary orders
344  const int p_order = static_cast<unsigned int>(order);
345  const int m = p_order-(i-1);
346  const int n = (i-1);
347 
348  unsigned int binomial_p_i = 1;
349 
350  // the binomial coefficient (p choose n)
351  if (i>1)
352  binomial_p_i = Utility::factorial(p_order)
353  / (Utility::factorial(n)*Utility::factorial(p_order-n));
354 
355 
356 
357  switch(i)
358  {
359  case 0:
360  return binomial_p_i * (-1./2.) * p_order * std::pow((1-xi)/2,static_cast<int>(p_order-1));
361  case 1:
362  return binomial_p_i * ( 1./2.) * p_order * std::pow((1+xi)/2,static_cast<int>(p_order-1));
363 
364  default:
365  {
366  return binomial_p_i * (1./2. * n * std::pow((1+xi)/2,n-1) * std::pow((1-xi)/2,m)
367  - 1./2. * m * std::pow((1+xi)/2,n) * std::pow((1-xi)/2,m-1));
368  }
369  }
370 
371  // we should never get here
372  libmesh_error();
373  }
374 
375  }
376 
377  libmesh_error();
378  return 0.;
379 }
380 
381 
382 
383 template <>
385  const Order order,
386  const unsigned int i,
387  const unsigned int j,
388  const Point& p)
389 {
390  libmesh_assert(elem);
391 
392  return FE<1,BERNSTEIN>::shape_deriv(elem->type(),
393  static_cast<Order>(order + elem->p_level()), i, j, p);
394 }
395 
396 
397 
398 template <>
400  const Order,
401  const unsigned int,
402  const unsigned int,
403  const Point&)
404 {
405  static bool warning_given = false;
406 
407  if (!warning_given)
408  libMesh::err << "Second derivatives for Bernstein elements "
409  << "are not yet implemented!"
410  << std::endl;
411 
412  warning_given = true;
413  return 0.;
414 }
415 
416 
417 
418 
419 template <>
421  const Order,
422  const unsigned int,
423  const unsigned int,
424  const Point&)
425 {
426  static bool warning_given = false;
427 
428  if (!warning_given)
429  libMesh::err << "Second derivatives for Bernstein elements "
430  << "are not yet implemented!"
431  << std::endl;
432 
433  warning_given = true;
434  return 0.;
435 }
436 
437 } // namespace libMesh
438 
439 
440 #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo