libMesh
fe_hermite_shape_3D.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 
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  libmesh_error_msg("We'll never get here!");
444  return 0.;
445 }
446 
447 
448 
449 template <>
451  const Order,
452  const unsigned int,
453  const unsigned int,
454  const Point &)
455 {
456  libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
457  return 0.;
458 }
459 
460 
461 
462 template <>
464  const Order order,
465  const unsigned int i,
466  const unsigned int j,
467  const Point & p)
468 {
469  libmesh_assert(elem);
470  libmesh_assert (j == 0 || j == 1 || j == 2);
471 
472  std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
473 
474 #ifdef DEBUG
475  std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
476  std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
477 #endif //DEBUG
478 
479  hermite_compute_coefs(elem, dxdxi
480 #ifdef DEBUG
481  , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
482 #endif
483  );
484 
485  const ElemType type = elem->type();
486 
487  const Order totalorder = static_cast<Order>(order + elem->p_level());
488 
489  switch (totalorder)
490  {
491  // 3rd-order tricubic Hermite functions
492  case THIRD:
493  {
494  switch (type)
495  {
496  case HEX8:
497  case HEX20:
498  case HEX27:
499  {
500  libmesh_assert_less (i, 64);
501 
502  std::vector<unsigned int> bases1D;
503 
504  Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
505 
506  switch (j) // Derivative type
507  {
508  case 0:
509  return coef *
510  FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
511  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
512  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
513  break;
514  case 1:
515  return coef *
516  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
517  FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) *
518  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
519  break;
520  case 2:
521  return coef *
522  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
523  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
524  FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
525  break;
526  default:
527  libmesh_error_msg("Invalid shape function derivative j = " << j);
528  }
529 
530  }
531  default:
532  libmesh_error_msg("ERROR: Unsupported element type " << type);
533  }
534  }
535  // by default throw an error
536  default:
537  libmesh_error_msg("ERROR: Unsupported polynomial order " << totalorder);
538  }
539 
540  libmesh_error_msg("We'll never get here!");
541  return 0.;
542 }
543 
544 
545 
546 template <>
548  const Order order,
549  const unsigned int i,
550  const unsigned int j,
551  const Point & p)
552 {
553  libmesh_assert(elem);
554 
555  std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
556 
557 #ifdef DEBUG
558  std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
559  std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
560 #endif //DEBUG
561 
562  hermite_compute_coefs(elem, dxdxi
563 #ifdef DEBUG
564  , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
565 #endif
566  );
567 
568  const ElemType type = elem->type();
569 
570  const Order totalorder = static_cast<Order>(order + elem->p_level());
571 
572  switch (totalorder)
573  {
574  // 3rd-order tricubic Hermite functions
575  case THIRD:
576  {
577  switch (type)
578  {
579  case HEX8:
580  case HEX20:
581  case HEX27:
582  {
583  libmesh_assert_less (i, 64);
584 
585  std::vector<unsigned int> bases1D;
586 
587  Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
588 
589  switch (j) // Derivative type
590  {
591  case 0:
592  return coef *
594  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
595  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
596  break;
597  case 1:
598  return coef *
599  FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
600  FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) *
601  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
602  break;
603  case 2:
604  return coef *
605  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
607  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
608  break;
609  case 3:
610  return coef *
611  FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
612  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
613  FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
614  break;
615  case 4:
616  return coef *
617  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
618  FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) *
619  FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
620  break;
621  case 5:
622  return coef *
623  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
624  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
626  break;
627  default:
628  libmesh_error_msg("Invalid shape function derivative j = " << j);
629  }
630 
631  }
632  default:
633  libmesh_error_msg("ERROR: Unsupported element type " << type);
634  }
635  }
636  // by default throw an error
637  default:
638  libmesh_error_msg("ERROR: Unsupported polynomial order " << totalorder);
639  }
640 
641  libmesh_error_msg("We'll never get here!");
642  return 0.;
643 }
644 
645 } // 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
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)
static const Real TOLERANCE
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 cube_number_row[]
const unsigned char square_number_column[]
const unsigned char cube_number_page[]
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)