fe_hermite_shape_3D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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++ inlcludes
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
105  (std::vector<unsigned int> &bases1D,
106  const std::vector<std::vector<Real> > &dxdxi,
107  const Order &o,
108  unsigned int i)
109 {
110  bases1D.clear();
111  bases1D.resize(3,0);
112  Real coef = 1.0;
113 
114  unsigned int e = o-2;
115 
116  // Nodes
117  if (i < 64)
118  {
119  switch (i / 8)
120  {
121  case 0:
122  break;
123  case 1:
124  bases1D[0] = 1;
125  break;
126  case 2:
127  bases1D[0] = 1;
128  bases1D[1] = 1;
129  break;
130  case 3:
131  bases1D[1] = 1;
132  break;
133  case 4:
134  bases1D[2] = 1;
135  break;
136  case 5:
137  bases1D[0] = 1;
138  bases1D[2] = 1;
139  break;
140  case 6:
141  bases1D[0] = 1;
142  bases1D[1] = 1;
143  bases1D[2] = 1;
144  break;
145  case 7:
146  bases1D[1] = 1;
147  bases1D[2] = 1;
148  break;
149  }
150 
151  unsigned int basisnum = i%8;
152  switch (basisnum) // DoF type
153  {
154  case 0: // DoF = value at node
155  coef = 1.0;
156  break;
157  case 1: // DoF = x derivative at node
158  coef = dxdxi[0][bases1D[0]];
159  bases1D[0] += 2; break;
160  case 2: // DoF = y derivative at node
161  coef = dxdxi[1][bases1D[1]];
162  bases1D[1] += 2; break;
163  case 3: // DoF = xy derivative at node
164  coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]];
165  bases1D[0] += 2; bases1D[1] += 2; break;
166  case 4: // DoF = z derivative at node
167  coef = dxdxi[2][bases1D[2]];
168  bases1D[2] += 2; break;
169  case 5: // DoF = xz derivative at node
170  coef = dxdxi[0][bases1D[0]] * dxdxi[2][bases1D[2]];
171  bases1D[0] += 2; bases1D[2] += 2; break;
172  case 6: // DoF = yz derivative at node
173  coef = dxdxi[1][bases1D[1]] * dxdxi[2][bases1D[2]];
174  bases1D[1] += 2; bases1D[2] += 2; break;
175  case 7: // DoF = xyz derivative at node
176  coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]] * dxdxi[2][bases1D[2]];
177  bases1D[0] += 2; bases1D[1] += 2; bases1D[2] += 2; break;
178  }
179  }
180  // Edges
181  else if (i < 64 + 12*4*e)
182  {
183  unsigned int basisnum = (i - 64) % (4*e);
184  switch ((i - 64) / (4*e))
185  {
186  case 0:
187  bases1D[0] = basisnum / 4 + 4;
188  bases1D[1] = basisnum % 4 / 2 * 2;
189  bases1D[2] = basisnum % 2 * 2;
190  if (basisnum % 4 / 2)
191  coef *= dxdxi[1][0];
192  if (basisnum % 2)
193  coef *= dxdxi[2][0];
194  break;
195  case 1:
196  bases1D[0] = basisnum % 4 / 2 * 2 + 1;
197  bases1D[1] = basisnum / 4 + 4;
198  bases1D[2] = basisnum % 2 * 2;
199  if (basisnum % 4 / 2)
200  coef *= dxdxi[0][1];
201  if (basisnum % 2)
202  coef *= dxdxi[2][0];
203  break;
204  case 2:
205  bases1D[0] = basisnum / 4 + 4;
206  bases1D[1] = basisnum % 4 / 2 * 2 + 1;
207  bases1D[2] = basisnum % 2 * 2;
208  if (basisnum % 4 / 2)
209  coef *= dxdxi[1][1];
210  if (basisnum % 2)
211  coef *= dxdxi[2][0];
212  break;
213  case 3:
214  bases1D[0] = basisnum % 4 / 2 * 2;
215  bases1D[1] = basisnum / 4 + 4;
216  bases1D[2] = basisnum % 2 * 2;
217  if (basisnum % 4 / 2)
218  coef *= dxdxi[0][0];
219  if (basisnum % 2)
220  coef *= dxdxi[2][0];
221  break;
222  case 4:
223  bases1D[0] = basisnum % 4 / 2 * 2;
224  bases1D[1] = basisnum % 2 * 2;
225  bases1D[2] = basisnum / 4 + 4;
226  if (basisnum % 4 / 2)
227  coef *= dxdxi[0][0];
228  if (basisnum % 2)
229  coef *= dxdxi[1][0];
230  break;
231  case 5:
232  bases1D[0] = basisnum % 4 / 2 * 2 + 1;
233  bases1D[1] = basisnum % 2 * 2;
234  bases1D[2] = basisnum / 4 + 4;
235  if (basisnum % 4 / 2)
236  coef *= dxdxi[0][1];
237  if (basisnum % 2)
238  coef *= dxdxi[1][0];
239  break;
240  case 6:
241  bases1D[0] = basisnum % 4 / 2 * 2 + 1;
242  bases1D[1] = basisnum % 2 * 2 + 1;
243  bases1D[2] = basisnum / 4 + 4;
244  if (basisnum % 4 / 2)
245  coef *= dxdxi[0][1];
246  if (basisnum % 2)
247  coef *= dxdxi[1][1];
248  break;
249  case 7:
250  bases1D[0] = basisnum % 4 / 2 * 2;
251  bases1D[1] = basisnum % 2 * 2 + 1;
252  bases1D[2] = basisnum / 4 + 4;
253  if (basisnum % 4 / 2)
254  coef *= dxdxi[0][0];
255  if (basisnum % 2)
256  coef *= dxdxi[1][1];
257  break;
258  case 8:
259  bases1D[0] = basisnum / 4 + 4;
260  bases1D[1] = basisnum % 4 / 2 * 2;
261  bases1D[2] = basisnum % 2 * 2 + 1;
262  if (basisnum % 4 / 2)
263  coef *= dxdxi[1][0];
264  if (basisnum % 2)
265  coef *= dxdxi[2][1];
266  break;
267  case 9:
268  bases1D[0] = basisnum % 4 / 2 * 2 + 1;
269  bases1D[1] = basisnum / 4 + 4;
270  bases1D[2] = basisnum % 2 * 2;
271  if (basisnum % 4 / 2)
272  coef *= dxdxi[0][1];
273  if (basisnum % 2)
274  coef *= dxdxi[2][1];
275  break;
276  case 10:
277  bases1D[0] = basisnum / 4 + 4;
278  bases1D[1] = basisnum % 4 / 2 * 2 + 1;
279  bases1D[2] = basisnum % 2 * 2 + 1;
280  if (basisnum % 4 / 2)
281  coef *= dxdxi[1][1];
282  if (basisnum % 2)
283  coef *= dxdxi[2][1];
284  break;
285  case 11:
286  bases1D[0] = basisnum % 4 / 2 * 2;
287  bases1D[1] = basisnum / 4 + 4;
288  bases1D[2] = basisnum % 2 * 2 + 1;
289  if (basisnum % 4 / 2)
290  coef *= dxdxi[0][0];
291  if (basisnum % 2)
292  coef *= dxdxi[2][1];
293  break;
294  }
295  }
296  // Faces
297  else if (i < 64 + 12*4*e + 6*2*e*e)
298  {
299  unsigned int basisnum = (i - 64 - 12*4*e) % (2*e*e);
300  switch ((i - 64 - 12*4*e) / (2*e*e))
301  {
302  case 0:
303  bases1D[0] = square_number_column[basisnum / 2];
304  bases1D[1] = square_number_row[basisnum / 2];
305  bases1D[2] = basisnum % 2 * 2;
306  if (basisnum % 2)
307  coef *= dxdxi[2][0];
308  break;
309  case 1:
310  bases1D[0] = square_number_column[basisnum / 2];
311  bases1D[1] = basisnum % 2 * 2;
312  bases1D[2] = square_number_row[basisnum / 2];
313  if (basisnum % 2)
314  coef *= dxdxi[1][0];
315  break;
316  case 2:
317  bases1D[0] = basisnum % 2 * 2 + 1;
318  bases1D[1] = square_number_column[basisnum / 2];
319  bases1D[2] = square_number_row[basisnum / 2];
320  if (basisnum % 2)
321  coef *= dxdxi[0][1];
322  break;
323  case 3:
324  bases1D[0] = square_number_column[basisnum / 2];
325  bases1D[1] = basisnum % 2 * 2 + 1;
326  bases1D[2] = square_number_row[basisnum / 2];
327  if (basisnum % 2)
328  coef *= dxdxi[1][1];
329  break;
330  case 4:
331  bases1D[0] = basisnum % 2 * 2;
332  bases1D[1] = square_number_column[basisnum / 2];
333  bases1D[2] = square_number_row[basisnum / 2];
334  if (basisnum % 2)
335  coef *= dxdxi[0][0];
336  break;
337  case 5:
338  bases1D[0] = square_number_column[basisnum / 2];
339  bases1D[1] = square_number_row[basisnum / 2];
340  bases1D[2] = basisnum % 2 * 2 + 1;
341  if (basisnum % 2)
342  coef *= dxdxi[2][1];
343  break;
344  }
345  }
346  // Interior
347  else
348  {
349  unsigned int basisnum = i - 64 - 12*4*e;
350  bases1D[0] = cube_number_column[basisnum] + 4;
351  bases1D[1] = cube_number_row[basisnum] + 4;
352  bases1D[2] = cube_number_page[basisnum] + 4;
353  }
354 
355  // No singular elements
356  libmesh_assert(coef);
357  return coef;
358 }
359 
360 
361 } // end anonymous namespace
362 
363 
364 namespace libMesh
365 {
366 
367 
368 template <>
370  const Order,
371  const unsigned int,
372  const Point&)
373 {
374  libMesh::err << "Hermite elements require the real element\n"
375  << "to construct gradient-based degrees of freedom."
376  << std::endl;
377 
378  libmesh_error();
379  return 0.;
380 }
381 
382 
383 
384 template <>
386  const Order order,
387  const unsigned int i,
388  const Point& p)
389 {
390  libmesh_assert(elem);
391 
392  std::vector<std::vector<Real> > dxdxi(3, std::vector<Real>(2, 0));
393 
394 #ifdef DEBUG
395  std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
396  std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
397 #endif //DEBUG
398 
399  hermite_compute_coefs(elem, dxdxi
400 #ifdef DEBUG
401  , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
402 #endif
403  );
404 
405  const ElemType type = elem->type();
406 
407  const Order totalorder = static_cast<Order>(order + elem->p_level());
408 
409  switch (totalorder)
410  {
411  // 3rd-order tricubic Hermite functions
412  case THIRD:
413  {
414  switch (type)
415  {
416  case HEX8:
417  case HEX20:
418  case HEX27:
419  {
420  libmesh_assert_less (i, 64);
421 
422  std::vector<unsigned int> bases1D;
423 
424  Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
425 
426  return coef *
427  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
428  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
429  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
430  }
431  default:
432  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
433  libmesh_error();
434  }
435  }
436  // by default throw an error
437  default:
438  libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl;
439  libmesh_error();
440  }
441 
442  libmesh_error();
443  return 0.;
444 }
445 
446 
447 
448 template <>
450  const Order,
451  const unsigned int,
452  const unsigned int,
453  const Point&)
454 {
455  libMesh::err << "Hermite elements require the real element\n"
456  << "to construct gradient-based degrees of freedom."
457  << std::endl;
458 
459  libmesh_error();
460  return 0.;
461 }
462 
463 
464 
465 template <>
467  const Order order,
468  const unsigned int i,
469  const unsigned int j,
470  const Point& p)
471 {
472  libmesh_assert(elem);
473  libmesh_assert (j == 0 || j == 1 || j == 2);
474 
475  std::vector<std::vector<Real> > dxdxi(3, std::vector<Real>(2, 0));
476 
477 #ifdef DEBUG
478  std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
479  std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
480 #endif //DEBUG
481 
482  hermite_compute_coefs(elem, dxdxi
483 #ifdef DEBUG
484  , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
485 #endif
486  );
487 
488  const ElemType type = elem->type();
489 
490  const Order totalorder = static_cast<Order>(order + elem->p_level());
491 
492  switch (totalorder)
493  {
494  // 3rd-order tricubic Hermite functions
495  case THIRD:
496  {
497  switch (type)
498  {
499  case HEX8:
500  case HEX20:
501  case HEX27:
502  {
503  libmesh_assert_less (i, 64);
504 
505  std::vector<unsigned int> bases1D;
506 
507  Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
508 
509  switch (j) // Derivative type
510  {
511  case 0:
512  return coef *
513  FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
514  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
515  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
516  break;
517  case 1:
518  return coef *
519  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
520  FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) *
521  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
522  break;
523  case 2:
524  return coef *
525  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
526  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
527  FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
528  break;
529  }
530 
531  }
532  default:
533  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
534  libmesh_error();
535  }
536  }
537  // by default throw an error
538  default:
539  libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl;
540  libmesh_error();
541  }
542 
543  libmesh_error();
544  return 0.;
545 }
546 
547 
548 
549 template <>
551  const Order order,
552  const unsigned int i,
553  const unsigned int j,
554  const Point& p)
555 {
556  libmesh_assert(elem);
557 
558  std::vector<std::vector<Real> > dxdxi(3, std::vector<Real>(2, 0));
559 
560 #ifdef DEBUG
561  std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
562  std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
563 #endif //DEBUG
564 
565  hermite_compute_coefs(elem, dxdxi
566 #ifdef DEBUG
567  , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
568 #endif
569  );
570 
571  const ElemType type = elem->type();
572 
573  const Order totalorder = static_cast<Order>(order + elem->p_level());
574 
575  switch (totalorder)
576  {
577  // 3rd-order tricubic Hermite functions
578  case THIRD:
579  {
580  switch (type)
581  {
582  case HEX8:
583  case HEX20:
584  case HEX27:
585  {
586  libmesh_assert_less (i, 64);
587 
588  std::vector<unsigned int> bases1D;
589 
590  Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
591 
592  switch (j) // Derivative type
593  {
594  case 0:
595  return coef *
597  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
598  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
599  break;
600  case 1:
601  return coef *
602  FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
603  FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) *
604  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
605  break;
606  case 2:
607  return coef *
608  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
610  FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
611  break;
612  case 3:
613  return coef *
614  FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
615  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
616  FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
617  break;
618  case 4:
619  return coef *
620  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
621  FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) *
622  FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
623  break;
624  case 5:
625  return coef *
626  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
627  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
629  break;
630  }
631 
632  }
633  default:
634  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
635  libmesh_error();
636  }
637  }
638  // by default throw an error
639  default:
640  libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl;
641  libmesh_error();
642  }
643 
644  libmesh_error();
645  return 0.;
646 }
647 
648 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo