libMesh
fe_hermite_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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  {
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 <>
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 <>
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  libmesh_error_msg("We'll never get here!");
248  return 0.;
249 }
250 
251 
252 
253 template <>
255  const Order,
256  const unsigned int,
257  const unsigned int,
258  const Point &)
259 {
260  libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
261  return 0.;
262 }
263 
264 
265 
266 template <>
268  const Order order,
269  const unsigned int i,
270  const unsigned int j,
271  const Point & p)
272 {
273  libmesh_assert(elem);
274  libmesh_assert (j == 0 || j == 1);
275 
276  std::vector<std::vector<Real>> dxdxi(2, std::vector<Real>(2, 0));
277 
278 #ifdef DEBUG
279  std::vector<Real> dxdeta(2), dydxi(2);
280 #endif
281 
282  hermite_compute_coefs(elem,dxdxi
283 #ifdef DEBUG
284  ,dxdeta,dydxi
285 #endif
286  );
287 
288  const ElemType type = elem->type();
289 
290  const Order totalorder = static_cast<Order>(order + elem->p_level());
291 
292  switch (type)
293  {
294  case QUAD4:
295  case QUADSHELL4:
296  libmesh_assert_less (totalorder, 4);
297  libmesh_fallthrough();
298  case QUAD8:
299  case QUADSHELL8:
300  case QUAD9:
301  {
302  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
303 
304  std::vector<unsigned int> bases1D;
305 
306  Real coef = hermite_bases_2D(bases1D, dxdxi, totalorder, i);
307 
308  switch (j)
309  {
310  case 0:
311  return coef *
312  FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
313  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1));
314  case 1:
315  return coef *
316  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
317  FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1));
318  default:
319  libmesh_error_msg("Invalid derivative index j = " << j);
320  }
321  }
322  default:
323  libmesh_error_msg("ERROR: Unsupported element type = " << type);
324  }
325 
326  libmesh_error_msg("We'll never get here!");
327  return 0.;
328 }
329 
330 
331 
332 template <>
334  const Order order,
335  const unsigned int i,
336  const unsigned int j,
337  const Point & p)
338 {
339  libmesh_assert(elem);
340  libmesh_assert (j == 0 || j == 1 || j == 2);
341 
342  std::vector<std::vector<Real>> dxdxi(2, std::vector<Real>(2, 0));
343 
344 #ifdef DEBUG
345  std::vector<Real> dxdeta(2), dydxi(2);
346 #endif
347 
348  hermite_compute_coefs(elem,dxdxi
349 #ifdef DEBUG
350  ,dxdeta,dydxi
351 #endif
352  );
353 
354  const ElemType type = elem->type();
355 
356  const Order totalorder = static_cast<Order>(order + elem->p_level());
357 
358  switch (type)
359  {
360  case QUAD4:
361  case QUADSHELL4:
362  libmesh_assert_less (totalorder, 4);
363  libmesh_fallthrough();
364  case QUAD8:
365  case QUADSHELL8:
366  case QUAD9:
367  {
368  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
369 
370  std::vector<unsigned int> bases1D;
371 
372  Real coef = hermite_bases_2D(bases1D, dxdxi, totalorder, i);
373 
374  switch (j)
375  {
376  case 0:
377  return coef *
379  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1));
380  case 1:
381  return coef *
382  FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
383  FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1));
384  case 2:
385  return coef *
386  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
388  default:
389  libmesh_error_msg("Invalid derivative index j = " << j);
390  }
391  }
392  default:
393  libmesh_error_msg("ERROR: Unsupported element type = " << type);
394  }
395 
396  libmesh_error_msg("We'll never get here!");
397  return 0.;
398 }
399 
400 } // namespace libMesh
double abs(double a)
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
virtual ElemType type() const =0
unsigned int p_level() const
Definition: elem.h:2422
ElemType
Defines an enum for geometric element types.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
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_assert(j)
const unsigned char square_number_column[]
virtual Order default_order() const =0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
virtual unsigned int n_shape_functions() const libmesh_override
Definition: fe.C:36
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
const unsigned char square_number_row[]
static Real hermite_raw_shape_deriv(const unsigned int basis_num, const Real xi)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)