fe_hermite_shape_3D.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 
33 #ifdef DEBUG
34  , std::vector<Real> & dydxi, std::vector<Real> & dzdeta, std::vector<Real> & dxdzeta,
35  std::vector<Real> & dzdxi, std::vector<Real> & dxdeta, std::vector<Real> & dydzeta
36 #endif
37  )
38 {
39 
40  const Order mapping_order (elem->default_order());
41  const ElemType mapping_elem_type (elem->type());
42  const int n_mapping_shape_functions =
43  FE<3,LAGRANGE>::n_shape_functions(mapping_elem_type,
44  mapping_order);
45 
46  static const Point dofpt[2] = {Point(-1,-1,-1), Point(1,1,1)};
47 
48  for (int p = 0; p != 2; ++p)
49  {
50  dxdxi[0][p] = 0;
51  dxdxi[1][p] = 0;
52  dxdxi[2][p] = 0;
53 #ifdef DEBUG
54  dydxi[p] = 0;
55  dzdeta[p] = 0;
56  dxdzeta[p] = 0;
57  dzdxi[p] = 0;
58  dxdeta[p] = 0;
59  dydzeta[p] = 0;
60 #endif
61  for (int i = 0; i != n_mapping_shape_functions; ++i)
62  {
64  (mapping_elem_type, mapping_order, i, 0, dofpt[p]);
65  const Real ddeta = FE<3,LAGRANGE>::shape_deriv
66  (mapping_elem_type, mapping_order, i, 1, dofpt[p]);
67  const Real ddzeta = FE<3,LAGRANGE>::shape_deriv
68  (mapping_elem_type, mapping_order, i, 2, dofpt[p]);
69 
70  // dxdeta, dxdzeta, dydxi, dydzeta, dzdxi, dzdeta should all
71  // be 0!
72  const Point & point_i = elem->point(i);
73  dxdxi[0][p] += point_i(0) * ddxi;
74  dxdxi[1][p] += point_i(1) * ddeta;
75  dxdxi[2][p] += point_i(2) * ddzeta;
76 #ifdef DEBUG
77  dydxi[p] += point_i(1) * ddxi;
78  dzdeta[p] += point_i(2) * ddeta;
79  dxdzeta[p] += point_i(0) * ddzeta;
80  dzdxi[p] += point_i(2) * ddxi;
81  dxdeta[p] += point_i(0) * ddeta;
82  dydzeta[p] += point_i(1) * ddzeta;
83 #endif
84  }
85 
86  // No singular elements!
87  libmesh_assert(dxdxi[0][p]);
88  libmesh_assert(dxdxi[1][p]);
89  libmesh_assert(dxdxi[2][p]);
90  // No non-rectilinear or non-axis-aligned elements!
91 #ifdef DEBUG
92  libmesh_assert_less (std::abs(dydxi[p]), TOLERANCE);
93  libmesh_assert_less (std::abs(dzdeta[p]), TOLERANCE);
94  libmesh_assert_less (std::abs(dxdzeta[p]), TOLERANCE);
95  libmesh_assert_less (std::abs(dzdxi[p]), TOLERANCE);
96  libmesh_assert_less (std::abs(dxdeta[p]), TOLERANCE);
97  libmesh_assert_less (std::abs(dydzeta[p]), TOLERANCE);
98 #endif
99  }
100 }
101 
102 
103 
104 Real hermite_bases_3D (std::vector<unsigned int> & bases1D,
105  const std::vector<std::vector<Real>> & dxdxi,
106  const Order & o,
107  unsigned int i)
108 {
109  bases1D.clear();
110  bases1D.resize(3,0);
111  Real coef = 1.0;
112 
113  unsigned int e = o-2;
114 
115  // Nodes
116  if (i < 64)
117  {
118  switch (i / 8)
119  {
120  case 0:
121  break;
122  case 1:
123  bases1D[0] = 1;
124  break;
125  case 2:
126  bases1D[0] = 1;
127  bases1D[1] = 1;
128  break;
129  case 3:
130  bases1D[1] = 1;
131  break;
132  case 4:
133  bases1D[2] = 1;
134  break;
135  case 5:
136  bases1D[0] = 1;
137  bases1D[2] = 1;
138  break;
139  case 6:
140  bases1D[0] = 1;
141  bases1D[1] = 1;
142  bases1D[2] = 1;
143  break;
144  case 7:
145  bases1D[1] = 1;
146  bases1D[2] = 1;
147  break;
148  default:
149  libmesh_error_msg("Invalid basis node = " << i/8);
150  }
151 
152  unsigned int basisnum = i%8;
153  switch (basisnum) // DoF type
154  {
155  case 0: // DoF = value at node
156  coef = 1.0;
157  break;
158  case 1: // DoF = x derivative at node
159  coef = dxdxi[0][bases1D[0]];
160  bases1D[0] += 2; break;
161  case 2: // DoF = y derivative at node
162  coef = dxdxi[1][bases1D[1]];
163  bases1D[1] += 2; break;
164  case 3: // DoF = xy derivative at node
165  coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]];
166  bases1D[0] += 2; bases1D[1] += 2; break;
167  case 4: // DoF = z derivative at node
168  coef = dxdxi[2][bases1D[2]];
169  bases1D[2] += 2; break;
170  case 5: // DoF = xz derivative at node
171  coef = dxdxi[0][bases1D[0]] * dxdxi[2][bases1D[2]];
172  bases1D[0] += 2; bases1D[2] += 2; break;
173  case 6: // DoF = yz derivative at node
174  coef = dxdxi[1][bases1D[1]] * dxdxi[2][bases1D[2]];
175  bases1D[1] += 2; bases1D[2] += 2; break;
176  case 7: // DoF = xyz derivative at node
177  coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]] * dxdxi[2][bases1D[2]];
178  bases1D[0] += 2; bases1D[1] += 2; bases1D[2] += 2; break;
179  default:
180  libmesh_error_msg("Invalid basisnum = " << basisnum);
181  }
182  }
183  // Edges
184  else if (i < 64 + 12*4*e)
185  {
186  unsigned int basisnum = (i - 64) % (4*e);
187  switch ((i - 64) / (4*e))
188  {
189  case 0:
190  bases1D[0] = basisnum / 4 + 4;
191  bases1D[1] = basisnum % 4 / 2 * 2;
192  bases1D[2] = basisnum % 2 * 2;
193  if (basisnum % 4 / 2)
194  coef *= dxdxi[1][0];
195  if (basisnum % 2)
196  coef *= dxdxi[2][0];
197  break;
198  case 1:
199  bases1D[0] = basisnum % 4 / 2 * 2 + 1;
200  bases1D[1] = basisnum / 4 + 4;
201  bases1D[2] = basisnum % 2 * 2;
202  if (basisnum % 4 / 2)
203  coef *= dxdxi[0][1];
204  if (basisnum % 2)
205  coef *= dxdxi[2][0];
206  break;
207  case 2:
208  bases1D[0] = basisnum / 4 + 4;
209  bases1D[1] = basisnum % 4 / 2 * 2 + 1;
210  bases1D[2] = basisnum % 2 * 2;
211  if (basisnum % 4 / 2)
212  coef *= dxdxi[1][1];
213  if (basisnum % 2)
214  coef *= dxdxi[2][0];
215  break;
216  case 3:
217  bases1D[0] = basisnum % 4 / 2 * 2;
218  bases1D[1] = basisnum / 4 + 4;
219  bases1D[2] = basisnum % 2 * 2;
220  if (basisnum % 4 / 2)
221  coef *= dxdxi[0][0];
222  if (basisnum % 2)
223  coef *= dxdxi[2][0];
224  break;
225  case 4:
226  bases1D[0] = basisnum % 4 / 2 * 2;
227  bases1D[1] = basisnum % 2 * 2;
228  bases1D[2] = basisnum / 4 + 4;
229  if (basisnum % 4 / 2)
230  coef *= dxdxi[0][0];
231  if (basisnum % 2)
232  coef *= dxdxi[1][0];
233  break;
234  case 5:
235  bases1D[0] = basisnum % 4 / 2 * 2 + 1;
236  bases1D[1] = basisnum % 2 * 2;
237  bases1D[2] = basisnum / 4 + 4;
238  if (basisnum % 4 / 2)
239  coef *= dxdxi[0][1];
240  if (basisnum % 2)
241  coef *= dxdxi[1][0];
242  break;
243  case 6:
244  bases1D[0] = basisnum % 4 / 2 * 2 + 1;
245  bases1D[1] = basisnum % 2 * 2 + 1;
246  bases1D[2] = basisnum / 4 + 4;
247  if (basisnum % 4 / 2)
248  coef *= dxdxi[0][1];
249  if (basisnum % 2)
250  coef *= dxdxi[1][1];
251  break;
252  case 7:
253  bases1D[0] = basisnum % 4 / 2 * 2;
254  bases1D[1] = basisnum % 2 * 2 + 1;
255  bases1D[2] = basisnum / 4 + 4;
256  if (basisnum % 4 / 2)
257  coef *= dxdxi[0][0];
258  if (basisnum % 2)
259  coef *= dxdxi[1][1];
260  break;
261  case 8:
262  bases1D[0] = basisnum / 4 + 4;
263  bases1D[1] = basisnum % 4 / 2 * 2;
264  bases1D[2] = basisnum % 2 * 2 + 1;
265  if (basisnum % 4 / 2)
266  coef *= dxdxi[1][0];
267  if (basisnum % 2)
268  coef *= dxdxi[2][1];
269  break;
270  case 9:
271  bases1D[0] = basisnum % 4 / 2 * 2 + 1;
272  bases1D[1] = basisnum / 4 + 4;
273  bases1D[2] = basisnum % 2 * 2;
274  if (basisnum % 4 / 2)
275  coef *= dxdxi[0][1];
276  if (basisnum % 2)
277  coef *= dxdxi[2][1];
278  break;
279  case 10:
280  bases1D[0] = basisnum / 4 + 4;
281  bases1D[1] = basisnum % 4 / 2 * 2 + 1;
282  bases1D[2] = basisnum % 2 * 2 + 1;
283  if (basisnum % 4 / 2)
284  coef *= dxdxi[1][1];
285  if (basisnum % 2)
286  coef *= dxdxi[2][1];
287  break;
288  case 11:
289  bases1D[0] = basisnum % 4 / 2 * 2;
290  bases1D[1] = basisnum / 4 + 4;
291  bases1D[2] = basisnum % 2 * 2 + 1;
292  if (basisnum % 4 / 2)
293  coef *= dxdxi[0][0];
294  if (basisnum % 2)
295  coef *= dxdxi[2][1];
296  break;
297  default:
298  libmesh_error_msg("Invalid basis node = " << (i - 64) / (4*e));
299  }
300  }
301  // Faces
302  else if (i < 64 + 12*4*e + 6*2*e*e)
303  {
304  unsigned int basisnum = (i - 64 - 12*4*e) % (2*e*e);
305  switch ((i - 64 - 12*4*e) / (2*e*e))
306  {
307  case 0:
308  bases1D[0] = square_number_column[basisnum / 2];
309  bases1D[1] = square_number_row[basisnum / 2];
310  bases1D[2] = basisnum % 2 * 2;
311  if (basisnum % 2)
312  coef *= dxdxi[2][0];
313  break;
314  case 1:
315  bases1D[0] = square_number_column[basisnum / 2];
316  bases1D[1] = basisnum % 2 * 2;
317  bases1D[2] = square_number_row[basisnum / 2];
318  if (basisnum % 2)
319  coef *= dxdxi[1][0];
320  break;
321  case 2:
322  bases1D[0] = basisnum % 2 * 2 + 1;
323  bases1D[1] = square_number_column[basisnum / 2];
324  bases1D[2] = square_number_row[basisnum / 2];
325  if (basisnum % 2)
326  coef *= dxdxi[0][1];
327  break;
328  case 3:
329  bases1D[0] = square_number_column[basisnum / 2];
330  bases1D[1] = basisnum % 2 * 2 + 1;
331  bases1D[2] = square_number_row[basisnum / 2];
332  if (basisnum % 2)
333  coef *= dxdxi[1][1];
334  break;
335  case 4:
336  bases1D[0] = basisnum % 2 * 2;
337  bases1D[1] = square_number_column[basisnum / 2];
338  bases1D[2] = square_number_row[basisnum / 2];
339  if (basisnum % 2)
340  coef *= dxdxi[0][0];
341  break;
342  case 5:
343  bases1D[0] = square_number_column[basisnum / 2];
344  bases1D[1] = square_number_row[basisnum / 2];
345  bases1D[2] = basisnum % 2 * 2 + 1;
346  if (basisnum % 2)
347  coef *= dxdxi[2][1];
348  break;
349  default:
350  libmesh_error_msg("Invalid basis node = " << (i - 64 - 12*4*e) / (2*e*e));
351  }
352  }
353  // Interior
354  else
355  {
356  unsigned int basisnum = i - 64 - 12*4*e;
357  bases1D[0] = cube_number_column[basisnum] + 4;
358  bases1D[1] = cube_number_row[basisnum] + 4;
359  bases1D[2] = cube_number_page[basisnum] + 4;
360  }
361 
362  // No singular elements
363  libmesh_assert(coef);
364  return coef;
365 }
366 
367 
368 } // end anonymous namespace
369 
370 
371 namespace libMesh
372 {
373 
374 
375 template <>
377  const Order,
378  const unsigned int,
379  const Point &)
380 {
381  libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
382  return 0.;
383 }
384 
385 
386 
387 template <>
389  const Order order,
390  const unsigned int i,
391  const Point & p)
392 {
393  libmesh_assert(elem);
394 
395  std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
396 
397 #ifdef DEBUG
398  std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
399  std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
400 #endif //DEBUG
401 
402  hermite_compute_coefs(elem, dxdxi
403 #ifdef DEBUG
404  , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
405 #endif
406  );
407 
408  const ElemType type = elem->type();
409 
410  const Order totalorder = static_cast<Order>(order + elem->p_level());
411 
412  switch (totalorder)
413  {
414  // 3rd-order tricubic Hermite functions
415  case THIRD:
416  {
417  switch (type)
418  {
419  case HEX8:
420  case HEX20:
421  case HEX27:
422  {
423  libmesh_assert_less (i, 64);
424 
425  std::vector<unsigned int> bases1D;
426 
427  Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
428 
429  return coef *
430  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
431  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
432  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
433  }
434  default:
435  libmesh_error_msg("ERROR: Unsupported element type " << type);
436  }
437  }
438  // by default throw an error
439  default:
440  libmesh_error_msg("ERROR: Unsupported polynomial order " << totalorder);
441  }
442 }
443 
444 
445 
446 template <>
448  const Order,
449  const unsigned int,
450  const unsigned int,
451  const Point &)
452 {
453  libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
454  return 0.;
455 }
456 
457 
458 
459 template <>
461  const Order order,
462  const unsigned int i,
463  const unsigned int j,
464  const Point & p)
465 {
466  libmesh_assert(elem);
467  libmesh_assert (j == 0 || j == 1 || j == 2);
468 
469  std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
470 
471 #ifdef DEBUG
472  std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
473  std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
474 #endif //DEBUG
475 
476  hermite_compute_coefs(elem, dxdxi
477 #ifdef DEBUG
478  , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
479 #endif
480  );
481 
482  const ElemType type = elem->type();
483 
484  const Order totalorder = static_cast<Order>(order + elem->p_level());
485 
486  switch (totalorder)
487  {
488  // 3rd-order tricubic Hermite functions
489  case THIRD:
490  {
491  switch (type)
492  {
493  case HEX8:
494  case HEX20:
495  case HEX27:
496  {
497  libmesh_assert_less (i, 64);
498 
499  std::vector<unsigned int> bases1D;
500 
501  Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
502 
503  switch (j) // Derivative type
504  {
505  case 0:
506  return coef *
507  FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
508  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
509  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
510  break;
511  case 1:
512  return coef *
513  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
514  FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) *
515  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
516  break;
517  case 2:
518  return coef *
519  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
520  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
521  FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
522  break;
523  default:
524  libmesh_error_msg("Invalid shape function derivative j = " << j);
525  }
526 
527  }
528  default:
529  libmesh_error_msg("ERROR: Unsupported element type " << type);
530  }
531  }
532  // by default throw an error
533  default:
534  libmesh_error_msg("ERROR: Unsupported polynomial order " << totalorder);
535  }
536 }
537 
538 
539 
540 template <>
542  const Order order,
543  const unsigned int i,
544  const unsigned int j,
545  const Point & p)
546 {
547  libmesh_assert(elem);
548 
549  std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
550 
551 #ifdef DEBUG
552  std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
553  std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
554 #endif //DEBUG
555 
556  hermite_compute_coefs(elem, dxdxi
557 #ifdef DEBUG
558  , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
559 #endif
560  );
561 
562  const ElemType type = elem->type();
563 
564  const Order totalorder = static_cast<Order>(order + elem->p_level());
565 
566  switch (totalorder)
567  {
568  // 3rd-order tricubic Hermite functions
569  case THIRD:
570  {
571  switch (type)
572  {
573  case HEX8:
574  case HEX20:
575  case HEX27:
576  {
577  libmesh_assert_less (i, 64);
578 
579  std::vector<unsigned int> bases1D;
580 
581  Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
582 
583  switch (j) // Derivative type
584  {
585  case 0:
586  return coef *
588  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
589  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
590  break;
591  case 1:
592  return coef *
593  FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
594  FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) *
595  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
596  break;
597  case 2:
598  return coef *
599  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
601  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
602  break;
603  case 3:
604  return coef *
605  FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
606  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
607  FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
608  break;
609  case 4:
610  return coef *
611  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
612  FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) *
613  FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
614  break;
615  case 5:
616  return coef *
617  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
618  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
620  break;
621  default:
622  libmesh_error_msg("Invalid shape function derivative j = " << j);
623  }
624 
625  }
626  default:
627  libmesh_error_msg("ERROR: Unsupported element type " << type);
628  }
629  }
630  // by default throw an error
631  default:
632  libmesh_error_msg("ERROR: Unsupported polynomial order " << totalorder);
633  }
634 }
635 
636 } // namespace libMesh
double abs(double a)
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
The base class for all geometric element types.
Definition: elem.h:100
const unsigned char cube_number_column[]
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:2555
static const Real TOLERANCE
static Real hermite_raw_shape_second_deriv(const unsigned int basis_num, const Real xi)
static Real hermite_raw_shape(const unsigned int basis_num, const Real xi)
const unsigned char cube_number_row[]
const unsigned char square_number_column[]
const unsigned char cube_number_page[]
virtual unsigned int n_shape_functions() const override
Definition: fe.C:50
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const unsigned char square_number_row[]
static Real hermite_raw_shape_deriv(const unsigned int basis_num, const Real xi)
virtual Order default_order() const =0
virtual ElemType type() const =0
A geometric point in (x,y,z) space.
Definition: point.h:38
const Point & point(const unsigned int i) const
Definition: elem.h:1892
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)