libMesh
fe_hermite_shape_1D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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 // Local includes
19 #include "libmesh/fe.h"
20 #include "libmesh/elem.h"
21 #include "libmesh/fe_interface.h"
22 #include "libmesh/utility.h"
23 #include "libmesh/enum_to_string.h"
24 
25 namespace
26 {
27 using namespace libMesh;
28 
29 // Compute the static coefficients for an element
30 void hermite_compute_coefs(const Elem * elem, Real & d1xd1x, Real & d2xd2x)
31 {
32  const FEFamily mapping_family = FEMap::map_fe_type(*elem);
33  const FEType map_fe_type(elem->default_order(), mapping_family);
34 
35  // Note: we explicitly don't consider the elem->p_level() when
36  // computing the number of mapping shape functions.
37  const int n_mapping_shape_functions =
38  FEInterface::n_shape_functions (map_fe_type, /*extra_order=*/0, elem);
39 
40  // Degrees of freedom are at vertices and edge midpoints
41  static const Point dofpt[2] = {Point(-1), Point(1)};
42 
43  // Mapping functions - first derivatives at each dofpt
44  std::vector<Real> dxdxi(2);
45  std::vector<Real> dxidx(2);
46 
47  FEInterface::shape_deriv_ptr shape_deriv_ptr =
48  FEInterface::shape_deriv_function(map_fe_type, elem);
49 
50  for (int p = 0; p != 2; ++p)
51  {
52  dxdxi[p] = 0;
53  for (int i = 0; i != n_mapping_shape_functions; ++i)
54  {
55  const Real ddxi = shape_deriv_ptr
56  (map_fe_type, elem, i, 0, dofpt[p], /*add_p_level=*/false);
57  dxdxi[p] += elem->point(i)(0) * ddxi;
58  }
59  }
60 
61  // Calculate derivative scaling factors
62  d1xd1x = dxdxi[0];
63  d2xd2x = dxdxi[1];
64 }
65 
66 
67 } // end anonymous namespace
68 
69 
70 namespace libMesh
71 {
72 
73 
75 
76 
77 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
78 
79 template<>
80 Real FEHermite<1>::hermite_raw_shape_second_deriv (const unsigned int i, const Real xi)
81 {
82  using Utility::pow;
83 
84  switch (i)
85  {
86  case 0:
87  return 1.5 * xi;
88  case 1:
89  return -1.5 * xi;
90  case 2:
91  return 0.5 * (-1. + 3.*xi);
92  case 3:
93  return 0.5 * (1. + 3.*xi);
94  case 4:
95  return (8.*xi*xi + 4.*(xi*xi-1.))/24.;
96  case 5:
97  return (8.*xi*xi*xi + 12.*xi*(xi*xi-1.))/120.;
98  // case 6:
99  // return (8.*pow<4>(xi) + 20.*xi*xi*(xi*xi-1.) +
100  // 2.*(xi*xi-1)*(xi*xi-1))/720.;
101  default:
102  Real denominator = 720., xipower = 1.;
103  for (unsigned n=6; n != i; ++n)
104  {
105  xipower *= xi;
106  denominator *= (n+1);
107  }
108  return (8.*pow<4>(xi)*xipower +
109  (8.*(i-4)+4.)*xi*xi*xipower*(xi*xi-1.) +
110  (i-4)*(i-5)*xipower*(xi*xi-1.)*(xi*xi-1.))/denominator;
111  }
112 }
113 
114 #endif
115 
116 
117 template<>
118 Real FEHermite<1>::hermite_raw_shape_deriv(const unsigned int i, const Real xi)
119 {
120  switch (i)
121  {
122  case 0:
123  return 0.75 * (-1. + xi*xi);
124  case 1:
125  return 0.75 * (1. - xi*xi);
126  case 2:
127  return 0.25 * (-1. - 2.*xi + 3.*xi*xi);
128  case 3:
129  return 0.25 * (-1. + 2.*xi + 3.*xi*xi);
130  case 4:
131  return 4.*xi * (xi*xi-1.)/24.;
132  case 5:
133  return (4*xi*xi*(xi*xi-1.) + (xi*xi-1.)*(xi*xi-1.))/120.;
134  // case 6:
135  // return (4*xi*xi*xi*(xi*xi-1.) + 2*xi*(xi*xi-1.)*(xi*xi-1.))/720.;
136  default:
137  Real denominator = 720., xipower = 1.;
138  for (unsigned n=6; n != i; ++n)
139  {
140  xipower *= xi;
141  denominator *= (n+1);
142  }
143  return (4*xi*xi*xi*xipower*(xi*xi-1.) +
144  (i-4)*xi*xipower*(xi*xi-1.)*(xi*xi-1.))/denominator;
145  }
146 }
147 
148 template<>
149 Real FEHermite<1>::hermite_raw_shape(const unsigned int i, const Real xi)
150 {
151  switch (i)
152  {
153  case 0:
154  return 0.25 * (2. - 3.*xi + xi*xi*xi);
155  case 1:
156  return 0.25 * (2. + 3.*xi - xi*xi*xi);
157  case 2:
158  return 0.25 * (1. - xi - xi*xi + xi*xi*xi);
159  case 3:
160  return 0.25 * (-1. - xi + xi*xi + xi*xi*xi);
161  // All high order terms have the form x^(p-4)(x^2-1)^2/p!
162  case 4:
163  return (xi*xi-1.) * (xi*xi-1.)/24.;
164  case 5:
165  return xi * (xi*xi-1.) * (xi*xi-1.)/120.;
166  // case 6:
167  // return xi*xi * (xi*xi-1.) * (xi*xi-1.)/720.;
168  default:
169  Real denominator = 720., xipower = 1.;
170  for (unsigned n=6; n != i; ++n)
171  {
172  xipower *= xi;
173  denominator *= (n+1);
174  }
175  return (xi*xi*xipower*(xi*xi-1.)*(xi*xi-1.))/denominator;
176 
177  }
178 }
179 
180 
181 template <>
183  const Order libmesh_dbg_var(order),
184  const unsigned int i,
185  const Point & p,
186  const bool libmesh_dbg_var(add_p_level))
187 {
188  libmesh_assert(elem);
189 
190  // Coefficient naming: d(1)d(2n) is the coefficient of the
191  // global shape function corresponding to value 1 in terms of the
192  // local shape function corresponding to normal derivative 2
193  Real d1xd1x, d2xd2x;
194 
195  hermite_compute_coefs(elem, d1xd1x, d2xd2x);
196 
197  const ElemType type = elem->type();
198 
199 #ifndef NDEBUG
200  const unsigned int totalorder =
201  order + add_p_level * elem->p_level();
202 #endif
203 
204  switch (type)
205  {
206  // C1 functions on the C1 cubic edge
207  case EDGE2:
208  case EDGE3:
209  {
210  libmesh_assert_less (i, totalorder+1);
211 
212  switch (i)
213  {
214  case 0:
215  return FEHermite<1>::hermite_raw_shape(0, p(0));
216  case 1:
217  return d1xd1x * FEHermite<1>::hermite_raw_shape(2, p(0));
218  case 2:
219  return FEHermite<1>::hermite_raw_shape(1, p(0));
220  case 3:
221  return d2xd2x * FEHermite<1>::hermite_raw_shape(3, p(0));
222  default:
223  return FEHermite<1>::hermite_raw_shape(i, p(0));
224  }
225  }
226  default:
227  libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
228  }
229 }
230 
231 
232 
233 template <>
235  const Order,
236  const unsigned int,
237  const Point &)
238 {
239  libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
240  return 0.;
241 }
242 
243 
244 template <>
246  const Elem * elem,
247  const unsigned int i,
248  const Point & p,
249  const bool add_p_level)
250 {
251  return FE<1,HERMITE>::shape(elem, fet.order, i, p, add_p_level);
252 }
253 
254 
255 template <>
257  const Order libmesh_dbg_var(order),
258  const unsigned int i,
259  const unsigned int,
260  const Point & p,
261  const bool libmesh_dbg_var(add_p_level))
262 {
263  libmesh_assert(elem);
264 
265  // Coefficient naming: d(1)d(2n) is the coefficient of the
266  // global shape function corresponding to value 1 in terms of the
267  // local shape function corresponding to normal derivative 2
268  Real d1xd1x, d2xd2x;
269 
270  hermite_compute_coefs(elem, d1xd1x, d2xd2x);
271 
272  const ElemType type = elem->type();
273 
274 #ifndef NDEBUG
275  const unsigned int totalorder =
276  order + add_p_level * elem->p_level();
277 #endif
278 
279  switch (type)
280  {
281  // C1 functions on the C1 cubic edge
282  case EDGE2:
283  case EDGE3:
284  {
285  libmesh_assert_less (i, totalorder+1);
286 
287  switch (i)
288  {
289  case 0:
291  case 1:
292  return d1xd1x * FEHermite<1>::hermite_raw_shape_deriv(2, p(0));
293  case 2:
295  case 3:
296  return d2xd2x * FEHermite<1>::hermite_raw_shape_deriv(3, p(0));
297  default:
299  }
300  }
301  default:
302  libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
303  }
304 }
305 
306 
307 template <>
309  const Order,
310  const unsigned int,
311  const unsigned int,
312  const Point &)
313 {
314  libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
315  return 0.;
316 }
317 
318 template <>
320  const Elem * elem,
321  const unsigned int i,
322  const unsigned int j,
323  const Point & p,
324  const bool add_p_level)
325 {
326  return FE<1,HERMITE>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
327 }
328 
329 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
330 
331 
332 template <>
334  const Order libmesh_dbg_var(order),
335  const unsigned int i,
336  const unsigned int,
337  const Point & p,
338  const bool libmesh_dbg_var(add_p_level))
339 {
340  libmesh_assert(elem);
341 
342  // Coefficient naming: d(1)d(2n) is the coefficient of the
343  // global shape function corresponding to value 1 in terms of the
344  // local shape function corresponding to normal derivative 2
345  Real d1xd1x, d2xd2x;
346 
347  hermite_compute_coefs(elem, d1xd1x, d2xd2x);
348 
349  const ElemType type = elem->type();
350 
351 #ifndef NDEBUG
352  const unsigned int totalorder =
353  order + add_p_level * elem->p_level();
354 #endif
355 
356  switch (type)
357  {
358  // C1 functions on the C1 cubic edge
359  case EDGE2:
360  case EDGE3:
361  {
362  libmesh_assert_less (i, totalorder+1);
363 
364  switch (i)
365  {
366  case 0:
368  case 1:
369  return d1xd1x * FEHermite<1>::hermite_raw_shape_second_deriv(2, p(0));
370  case 2:
372  case 3:
373  return d2xd2x * FEHermite<1>::hermite_raw_shape_second_deriv(3, p(0));
374  default:
376  }
377  }
378  default:
379  libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
380  }
381 }
382 
383 
384 template <>
386  const Order,
387  const unsigned int,
388  const unsigned int,
389  const Point &)
390 {
391  libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
392  return 0.;
393 }
394 
395 template <>
397  const Elem * elem,
398  const unsigned int i,
399  const unsigned int j,
400  const Point & p,
401  const bool add_p_level)
402 {
403  return FE<1,HERMITE>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
404 }
405 
406 
407 #endif
408 
409 } // namespace libMesh
Real(* shape_deriv_ptr)(const FEType fet, const Elem *elem, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level)
Typedef for pointer to a function that returns FE shape function derivative values.
Definition: fe_interface.h:623
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
ElemType
Defines an enum for geometric element types.
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
unsigned int p_level() const
Definition: elem.h:2945
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:201
The libMesh namespace provides an interface to certain functionality in the library.
static Real hermite_raw_shape_second_deriv(const unsigned int basis_num, const Real xi)
1D hermite functions on unit interval
static Real hermite_raw_shape(const unsigned int basis_num, const Real xi)
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
T pow(const T &x)
Definition: utility.h:328
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:515
libmesh_assert(ctx)
static shape_deriv_ptr shape_deriv_function(const unsigned int dim, const FEType &fe_t, const ElemType t)
std::string enum_to_string(const T e)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static Real hermite_raw_shape_deriv(const unsigned int basis_num, const Real xi)
FEFamily
defines an enum for finite element families.
virtual Order default_order() const =0
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
const Point & point(const unsigned int i) const
Definition: elem.h:2277
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
static FEFamily map_fe_type(const Elem &elem)
Definition: fe_map.C:45