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