libMesh
fe_hermite_shape_2D.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 // C++ includes
20 
21 // Local includes
22 #include "libmesh/fe.h"
23 #include "libmesh/elem.h"
24 #include "libmesh/fe_interface.h"
25 #include "libmesh/number_lookups.h"
26 
27 namespace
28 {
29 using namespace libMesh;
30 
31 // Compute the static coefficients for an element
32 void hermite_compute_coefs(const Elem * elem, std::vector<std::vector<Real>> & dxdxi
33 #ifdef DEBUG
34  , std::vector<Real> & dxdeta, std::vector<Real> & dydxi
35 #endif
36  )
37 {
38  const FEFamily mapping_family = FEMap::map_fe_type(*elem);
39  const Order mapping_order (elem->default_order());
40  const ElemType mapping_elem_type (elem->type());
41 
42  const FEType map_fe_type(mapping_order, mapping_family);
43 
44  const int n_mapping_shape_functions =
45  FEInterface::n_shape_functions (2, map_fe_type,
46  mapping_elem_type);
47 
48  std::vector<Point> dofpt;
49  dofpt.push_back(Point(-1,-1));
50  dofpt.push_back(Point(1,1));
51 
52  FEInterface::shape_deriv_ptr shape_deriv_ptr =
53  FEInterface::shape_deriv_function(2, map_fe_type);
54 
55  for (int p = 0; p != 2; ++p)
56  {
57  dxdxi[0][p] = 0;
58  dxdxi[1][p] = 0;
59 #ifdef DEBUG
60  dxdeta[p] = 0;
61  dydxi[p] = 0;
62 #endif
63  for (int i = 0; i != n_mapping_shape_functions; ++i)
64  {
65  const Real ddxi = shape_deriv_ptr
66  (elem, mapping_order, i, 0, dofpt[p], false);
67  const Real ddeta = shape_deriv_ptr
68  (elem, mapping_order, i, 1, dofpt[p], false);
69 
70  dxdxi[0][p] += elem->point(i)(0) * ddxi;
71  dxdxi[1][p] += elem->point(i)(1) * ddeta;
72  // dxdeta and dydxi should be 0!
73 #ifdef DEBUG
74  dxdeta[p] += elem->point(i)(0) * ddeta;
75  dydxi[p] += elem->point(i)(1) * ddxi;
76 #endif
77  }
78  // No singular elements!
79  libmesh_assert(dxdxi[0][p]);
80  libmesh_assert(dxdxi[1][p]);
81  // No non-rectilinear or non-axis-aligned elements!
82 #ifdef DEBUG
83  libmesh_assert_less (std::abs(dxdeta[p]), 1e-9);
84  libmesh_assert_less (std::abs(dydxi[p]), 1e-9);
85 #endif
86  }
87 }
88 
89 
90 
91 Real hermite_bases_2D (std::vector<unsigned int> & bases1D,
92  const std::vector<std::vector<Real>> & dxdxi,
93  const Order & o,
94  unsigned int i)
95 {
96  bases1D.clear();
97  bases1D.resize(2,0);
98  Real coef = 1.0;
99 
100  unsigned int e = o-3;
101 
102  // Nodes
103  if (i < 16)
104  {
105  switch (i / 4)
106  {
107  case 0:
108  break;
109  case 1:
110  bases1D[0] = 1;
111  break;
112  case 2:
113  bases1D[0] = 1;
114  bases1D[1] = 1;
115  break;
116  case 3:
117  bases1D[1] = 1;
118  break;
119  default:
120  libmesh_error_msg("Invalid basis node = " << i/4);
121  }
122 
123  unsigned int basisnum = i%4;
124  switch (basisnum)
125  {
126  case 0: // DoF = value at node
127  coef = 1.0;
128  break;
129  case 1: // DoF = x derivative at node
130  coef = dxdxi[0][bases1D[0]];
131  bases1D[0] += 2; break;
132  case 2: // DoF = y derivative at node
133  coef = dxdxi[1][bases1D[1]];
134  bases1D[1] += 2; break;
135  case 3: // DoF = xy derivative at node
136  coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]];
137  bases1D[0] += 2; bases1D[1] += 2; break;
138  default:
139  libmesh_error_msg("Invalid basisnum = " << basisnum);
140  }
141  }
142  // Edges
143  else if (i < 16 + 4*2*e)
144  {
145  unsigned int basisnum = (i - 16) % (2*e);
146  switch ((i - 16) / (2*e))
147  {
148  case 0:
149  bases1D[0] = basisnum/2 + 4;
150  bases1D[1] = basisnum%2*2;
151  if (basisnum%2)
152  coef = dxdxi[1][0];
153  break;
154  case 1:
155  bases1D[0] = basisnum%2*2 + 1;
156  bases1D[1] = basisnum/2 + 4;
157  if (basisnum%2)
158  coef = dxdxi[0][1];
159  break;
160  case 2:
161  bases1D[0] = basisnum/2 + 4;
162  bases1D[1] = basisnum%2*2 + 1;
163  if (basisnum%2)
164  coef = dxdxi[1][1];
165  break;
166  case 3:
167  bases1D[0] = basisnum%2*2;
168  bases1D[1] = basisnum/2 + 4;
169  if (basisnum%2)
170  coef = dxdxi[0][0];
171  break;
172  default:
173  libmesh_error_msg("Invalid basisnum = " << basisnum);
174  }
175  }
176  // Interior
177  else
178  {
179  unsigned int basisnum = i - 16 - 8*e;
180  bases1D[0] = square_number_row[basisnum]+4;
181  bases1D[1] = square_number_column[basisnum]+4;
182  }
183 
184  // No singular elements
185  libmesh_assert(coef);
186  return coef;
187 }
188 
189 } // end anonymous namespace
190 
191 
192 namespace libMesh
193 {
194 
195 
196 template <>
197 Real FE<2,HERMITE>::shape(const ElemType,
198  const Order,
199  const unsigned int,
200  const Point &)
201 {
202  libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
203  return 0.;
204 }
205 
206 
207 
208 template <>
209 Real FE<2,HERMITE>::shape(const Elem * elem,
210  const Order order,
211  const unsigned int i,
212  const Point & p,
213  const bool add_p_level)
214 {
215  libmesh_assert(elem);
216 
217  std::vector<std::vector<Real>> dxdxi(2, std::vector<Real>(2, 0));
218 
219 #ifdef DEBUG
220  std::vector<Real> dxdeta(2), dydxi(2);
221 #endif
222 
223  hermite_compute_coefs(elem,dxdxi
224 #ifdef DEBUG
225  ,dxdeta,dydxi
226 #endif
227  );
228 
229  const ElemType type = elem->type();
230 
231  const Order totalorder =
232  static_cast<Order>(order + add_p_level * elem->p_level());
233 
234  switch (type)
235  {
236  case QUAD4:
237  case QUADSHELL4:
238  libmesh_assert_less (totalorder, 4);
239  libmesh_fallthrough();
240  case QUAD8:
241  case QUADSHELL8:
242  case QUAD9:
243  {
244  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
245 
246  std::vector<unsigned int> bases1D;
247 
248  Real coef = hermite_bases_2D(bases1D, dxdxi, totalorder, i);
249 
250  return coef * FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
251  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1));
252  }
253  default:
254  libmesh_error_msg("ERROR: Unsupported element type = " << type);
255  }
256 }
257 
258 
259 
260 template <>
261 Real FE<2,HERMITE>::shape_deriv(const ElemType,
262  const Order,
263  const unsigned int,
264  const unsigned int,
265  const Point &)
266 {
267  libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
268  return 0.;
269 }
270 
271 
272 
273 template <>
274 Real FE<2,HERMITE>::shape_deriv(const Elem * elem,
275  const Order order,
276  const unsigned int i,
277  const unsigned int j,
278  const Point & p,
279  const bool add_p_level)
280 {
281  libmesh_assert(elem);
282  libmesh_assert (j == 0 || j == 1);
283 
284  std::vector<std::vector<Real>> dxdxi(2, std::vector<Real>(2, 0));
285 
286 #ifdef DEBUG
287  std::vector<Real> dxdeta(2), dydxi(2);
288 #endif
289 
290  hermite_compute_coefs(elem,dxdxi
291 #ifdef DEBUG
292  ,dxdeta,dydxi
293 #endif
294  );
295 
296  const ElemType type = elem->type();
297 
298  const Order totalorder =
299  static_cast<Order>(order + add_p_level * elem->p_level());
300 
301  switch (type)
302  {
303  case QUAD4:
304  case QUADSHELL4:
305  libmesh_assert_less (totalorder, 4);
306  libmesh_fallthrough();
307  case QUAD8:
308  case QUADSHELL8:
309  case QUAD9:
310  {
311  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
312 
313  std::vector<unsigned int> bases1D;
314 
315  Real coef = hermite_bases_2D(bases1D, dxdxi, totalorder, i);
316 
317  switch (j)
318  {
319  case 0:
320  return coef *
321  FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
322  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1));
323  case 1:
324  return coef *
325  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
326  FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1));
327  default:
328  libmesh_error_msg("Invalid derivative index j = " << j);
329  }
330  }
331  default:
332  libmesh_error_msg("ERROR: Unsupported element type = " << type);
333  }
334 }
335 
336 
337 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
338 
339 template <>
340 Real FE<2,HERMITE>::shape_second_deriv(const ElemType,
341  const Order,
342  const unsigned int,
343  const unsigned int,
344  const Point &)
345 {
346  libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
347  return 0.;
348 }
349 
350 
351 
352 template <>
353 Real FE<2,HERMITE>::shape_second_deriv(const Elem * elem,
354  const Order order,
355  const unsigned int i,
356  const unsigned int j,
357  const Point & p,
358  const bool add_p_level)
359 {
360  libmesh_assert(elem);
361  libmesh_assert (j == 0 || j == 1 || j == 2);
362 
363  std::vector<std::vector<Real>> dxdxi(2, std::vector<Real>(2, 0));
364 
365 #ifdef DEBUG
366  std::vector<Real> dxdeta(2), dydxi(2);
367 #endif
368 
369  hermite_compute_coefs(elem,dxdxi
370 #ifdef DEBUG
371  ,dxdeta,dydxi
372 #endif
373  );
374 
375  const ElemType type = elem->type();
376 
377  const Order totalorder =
378  static_cast<Order>(order + add_p_level * elem->p_level());
379 
380  switch (type)
381  {
382  case QUAD4:
383  case QUADSHELL4:
384  libmesh_assert_less (totalorder, 4);
385  libmesh_fallthrough();
386  case QUAD8:
387  case QUADSHELL8:
388  case QUAD9:
389  {
390  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
391 
392  std::vector<unsigned int> bases1D;
393 
394  Real coef = hermite_bases_2D(bases1D, dxdxi, totalorder, i);
395 
396  switch (j)
397  {
398  case 0:
399  return coef *
400  FEHermite<1>::hermite_raw_shape_second_deriv(bases1D[0],p(0)) *
401  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1));
402  case 1:
403  return coef *
404  FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
405  FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1));
406  case 2:
407  return coef *
408  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
409  FEHermite<1>::hermite_raw_shape_second_deriv(bases1D[1],p(1));
410  default:
411  libmesh_error_msg("Invalid derivative index j = " << j);
412  }
413  }
414  default:
415  libmesh_error_msg("ERROR: Unsupported element type = " << type);
416  }
417 }
418 
419 #endif
420 
421 } // namespace libMesh
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::FEInterface::n_shape_functions
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:446
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::FEInterface::shape_deriv_function
static shape_deriv_ptr shape_deriv_function(const unsigned int dim, const FEType &fe_t)
Definition: fe_interface.C:916
libMesh::FEHermite::hermite_raw_shape_second_deriv
static Real hermite_raw_shape_second_deriv(const unsigned int basis_num, const Real xi)
1D hermite functions on unit interval
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2512
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:1955
libMesh::square_number_column
const unsigned char square_number_column[]
Definition: number_lookups.C:56
libMesh::libmesh_assert
libmesh_assert(ctx)
std::abs
MetaPhysicL::DualNumber< T, D > abs(const MetaPhysicL::DualNumber< T, D > &in)
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::FEHermite::hermite_raw_shape
static Real hermite_raw_shape(const unsigned int basis_num, const Real xi)
libMesh::QUAD4
Definition: enum_elem_type.h:41
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
libMesh::Elem::default_order
virtual Order default_order() const =0
libMesh::FEMap::map_fe_type
static FEFamily map_fe_type(const Elem &elem)
Definition: fe_map.C:47
libMesh::FEInterface::shape_deriv_ptr
Real(* shape_deriv_ptr)(const Elem *elem, const Order o, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level)
Definition: fe_interface.h:345
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
libMesh::FEFamily
FEFamily
Definition: enum_fe_family.h:34
libMesh::QUADSHELL8
Definition: enum_elem_type.h:73
libMesh::square_number_row
const unsigned char square_number_row[]
Definition: number_lookups.C:69
libMesh::QUADSHELL4
Definition: enum_elem_type.h:72
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:100
libMesh::QUAD9
Definition: enum_elem_type.h:43
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::FEHermite::hermite_raw_shape_deriv
static Real hermite_raw_shape_deriv(const unsigned int basis_num, const Real xi)
libMesh::QUAD8
Definition: enum_elem_type.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