libMesh
fe_hierarchic_shape_3D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2019 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 // Local includes
20 #include "libmesh/fe.h"
21 #include "libmesh/elem.h"
22 #include "libmesh/number_lookups.h"
23 
24 // Anonymous namespace for functions shared by HIERARCHIC and
25 // L2_HIERARCHIC implementations. Implementations appear at the bottom
26 // of this file.
27 namespace
28 {
29 using namespace libMesh;
30 
31 Real fe_hierarchic_3D_shape(const Elem * elem,
32  const Order order,
33  const unsigned int i,
34  const Point & p,
35  const bool add_p_level);
36 
37 Real fe_hierarchic_3D_shape_deriv(const Elem * elem,
38  const Order order,
39  const unsigned int i,
40  const unsigned int j,
41  const Point & p,
42  const bool add_p_level);
43 
44 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
45 
46 Real fe_hierarchic_3D_shape_second_deriv(const Elem * elem,
47  const Order order,
48  const unsigned int i,
49  const unsigned int j,
50  const Point & p,
51  const bool add_p_level);
52 
53 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
54 
55 #if LIBMESH_DIM > 2
56 Point get_min_point(const Elem * elem,
57  unsigned int a,
58  unsigned int b,
59  unsigned int c,
60  unsigned int d)
61 {
62  return std::min(std::min(elem->point(a),elem->point(b)),
63  std::min(elem->point(c),elem->point(d)));
64 }
65 
66 void cube_indices(const Elem * elem,
67  const unsigned int totalorder,
68  const unsigned int i,
69  Real & xi, Real & eta, Real & zeta,
70  unsigned int & i0,
71  unsigned int & i1,
72  unsigned int & i2)
73 {
74  // The only way to make any sense of this
75  // is to look at the mgflo/mg2/mgf documentation
76  // and make the cut-out cube!
77  // Example i0 and i1 values for totalorder = 3:
78  // FIXME - these examples are incorrect now that we've got truly
79  // hierarchic basis functions
80  // Nodes 0 1 2 3 4 5 6 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 20 20 21 21 21 21 22 22 22 22 23 23 23 23 24 24 24 24 25 25 25 25 26 26 26 26 26 26 26 26
81  // DOFS 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 18 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 60 62 63
82  // static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3};
83  // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3, 0, 0, 0, 0, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3};
84  // static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3};
85 
86  // the number of DoFs per edge appears everywhere:
87  const unsigned int e = totalorder - 1u;
88 
89  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+1u));
90 
91  Real xi_saved = xi, eta_saved = eta, zeta_saved = zeta;
92 
93  // Vertices:
94  if (i == 0)
95  {
96  i0 = 0;
97  i1 = 0;
98  i2 = 0;
99  }
100  else if (i == 1)
101  {
102  i0 = 1;
103  i1 = 0;
104  i2 = 0;
105  }
106  else if (i == 2)
107  {
108  i0 = 1;
109  i1 = 1;
110  i2 = 0;
111  }
112  else if (i == 3)
113  {
114  i0 = 0;
115  i1 = 1;
116  i2 = 0;
117  }
118  else if (i == 4)
119  {
120  i0 = 0;
121  i1 = 0;
122  i2 = 1;
123  }
124  else if (i == 5)
125  {
126  i0 = 1;
127  i1 = 0;
128  i2 = 1;
129  }
130  else if (i == 6)
131  {
132  i0 = 1;
133  i1 = 1;
134  i2 = 1;
135  }
136  else if (i == 7)
137  {
138  i0 = 0;
139  i1 = 1;
140  i2 = 1;
141  }
142  // Edge 0
143  else if (i < 8 + e)
144  {
145  i0 = i - 6;
146  i1 = 0;
147  i2 = 0;
148  if (elem->point(0) > elem->point(1))
149  xi = -xi_saved;
150  }
151  // Edge 1
152  else if (i < 8 + 2*e)
153  {
154  i0 = 1;
155  i1 = i - e - 6;
156  i2 = 0;
157  if (elem->point(1) > elem->point(2))
158  eta = -eta_saved;
159  }
160  // Edge 2
161  else if (i < 8 + 3*e)
162  {
163  i0 = i - 2*e - 6;
164  i1 = 1;
165  i2 = 0;
166  if (elem->point(3) > elem->point(2))
167  xi = -xi_saved;
168  }
169  // Edge 3
170  else if (i < 8 + 4*e)
171  {
172  i0 = 0;
173  i1 = i - 3*e - 6;
174  i2 = 0;
175  if (elem->point(0) > elem->point(3))
176  eta = -eta_saved;
177  }
178  // Edge 4
179  else if (i < 8 + 5*e)
180  {
181  i0 = 0;
182  i1 = 0;
183  i2 = i - 4*e - 6;
184  if (elem->point(0) > elem->point(4))
185  zeta = -zeta_saved;
186  }
187  // Edge 5
188  else if (i < 8 + 6*e)
189  {
190  i0 = 1;
191  i1 = 0;
192  i2 = i - 5*e - 6;
193  if (elem->point(1) > elem->point(5))
194  zeta = -zeta_saved;
195  }
196  // Edge 6
197  else if (i < 8 + 7*e)
198  {
199  i0 = 1;
200  i1 = 1;
201  i2 = i - 6*e - 6;
202  if (elem->point(2) > elem->point(6))
203  zeta = -zeta_saved;
204  }
205  // Edge 7
206  else if (i < 8 + 8*e)
207  {
208  i0 = 0;
209  i1 = 1;
210  i2 = i - 7*e - 6;
211  if (elem->point(3) > elem->point(7))
212  zeta = -zeta_saved;
213  }
214  // Edge 8
215  else if (i < 8 + 9*e)
216  {
217  i0 = i - 8*e - 6;
218  i1 = 0;
219  i2 = 1;
220  if (elem->point(4) > elem->point(5))
221  xi = -xi_saved;
222  }
223  // Edge 9
224  else if (i < 8 + 10*e)
225  {
226  i0 = 1;
227  i1 = i - 9*e - 6;
228  i2 = 1;
229  if (elem->point(5) > elem->point(6))
230  eta = -eta_saved;
231  }
232  // Edge 10
233  else if (i < 8 + 11*e)
234  {
235  i0 = i - 10*e - 6;
236  i1 = 1;
237  i2 = 1;
238  if (elem->point(7) > elem->point(6))
239  xi = -xi_saved;
240  }
241  // Edge 11
242  else if (i < 8 + 12*e)
243  {
244  i0 = 0;
245  i1 = i - 11*e - 6;
246  i2 = 1;
247  if (elem->point(4) > elem->point(7))
248  eta = -eta_saved;
249  }
250  // Face 0
251  else if (i < 8 + 12*e + e*e)
252  {
253  unsigned int basisnum = i - 8 - 12*e;
254  i0 = square_number_row[basisnum] + 2;
255  i1 = square_number_column[basisnum] + 2;
256  i2 = 0;
257  const Point min_point = get_min_point(elem, 1, 2, 0, 3);
258 
259  if (elem->point(0) == min_point)
260  if (elem->point(1) == std::min(elem->point(1), elem->point(3)))
261  {
262  // Case 1
263  xi = xi_saved;
264  eta = eta_saved;
265  }
266  else
267  {
268  // Case 2
269  xi = eta_saved;
270  eta = xi_saved;
271  }
272 
273  else if (elem->point(3) == min_point)
274  if (elem->point(0) == std::min(elem->point(0), elem->point(2)))
275  {
276  // Case 3
277  xi = -eta_saved;
278  eta = xi_saved;
279  }
280  else
281  {
282  // Case 4
283  xi = xi_saved;
284  eta = -eta_saved;
285  }
286 
287  else if (elem->point(2) == min_point)
288  if (elem->point(3) == std::min(elem->point(3), elem->point(1)))
289  {
290  // Case 5
291  xi = -xi_saved;
292  eta = -eta_saved;
293  }
294  else
295  {
296  // Case 6
297  xi = -eta_saved;
298  eta = -xi_saved;
299  }
300 
301  else if (elem->point(1) == min_point)
302  {
303  if (elem->point(2) == std::min(elem->point(2), elem->point(0)))
304  {
305  // Case 7
306  xi = eta_saved;
307  eta = -xi_saved;
308  }
309  else
310  {
311  // Case 8
312  xi = -xi_saved;
313  eta = eta_saved;
314  }
315  }
316  }
317  // Face 1
318  else if (i < 8 + 12*e + 2*e*e)
319  {
320  unsigned int basisnum = i - 8 - 12*e - e*e;
321  i0 = square_number_row[basisnum] + 2;
322  i1 = 0;
323  i2 = square_number_column[basisnum] + 2;
324  const Point min_point = get_min_point(elem, 0, 1, 5, 4);
325 
326  if (elem->point(0) == min_point)
327  if (elem->point(1) == std::min(elem->point(1), elem->point(4)))
328  {
329  // Case 1
330  xi = xi_saved;
331  zeta = zeta_saved;
332  }
333  else
334  {
335  // Case 2
336  xi = zeta_saved;
337  zeta = xi_saved;
338  }
339 
340  else if (elem->point(1) == min_point)
341  if (elem->point(5) == std::min(elem->point(5), elem->point(0)))
342  {
343  // Case 3
344  xi = zeta_saved;
345  zeta = -xi_saved;
346  }
347  else
348  {
349  // Case 4
350  xi = -xi_saved;
351  zeta = zeta_saved;
352  }
353 
354  else if (elem->point(5) == min_point)
355  if (elem->point(4) == std::min(elem->point(4), elem->point(1)))
356  {
357  // Case 5
358  xi = -xi_saved;
359  zeta = -zeta_saved;
360  }
361  else
362  {
363  // Case 6
364  xi = -zeta_saved;
365  zeta = -xi_saved;
366  }
367 
368  else if (elem->point(4) == min_point)
369  {
370  if (elem->point(0) == std::min(elem->point(0), elem->point(5)))
371  {
372  // Case 7
373  xi = -xi_saved;
374  zeta = zeta_saved;
375  }
376  else
377  {
378  // Case 8
379  xi = xi_saved;
380  zeta = -zeta_saved;
381  }
382  }
383  }
384  // Face 2
385  else if (i < 8 + 12*e + 3*e*e)
386  {
387  unsigned int basisnum = i - 8 - 12*e - 2*e*e;
388  i0 = 1;
389  i1 = square_number_row[basisnum] + 2;
390  i2 = square_number_column[basisnum] + 2;
391  const Point min_point = get_min_point(elem, 1, 2, 6, 5);
392 
393  if (elem->point(1) == min_point)
394  if (elem->point(2) == std::min(elem->point(2), elem->point(5)))
395  {
396  // Case 1
397  eta = eta_saved;
398  zeta = zeta_saved;
399  }
400  else
401  {
402  // Case 2
403  eta = zeta_saved;
404  zeta = eta_saved;
405  }
406 
407  else if (elem->point(2) == min_point)
408  if (elem->point(6) == std::min(elem->point(6), elem->point(1)))
409  {
410  // Case 3
411  eta = zeta_saved;
412  zeta = -eta_saved;
413  }
414  else
415  {
416  // Case 4
417  eta = -eta_saved;
418  zeta = zeta_saved;
419  }
420 
421  else if (elem->point(6) == min_point)
422  if (elem->point(5) == std::min(elem->point(5), elem->point(2)))
423  {
424  // Case 5
425  eta = -eta_saved;
426  zeta = -zeta_saved;
427  }
428  else
429  {
430  // Case 6
431  eta = -zeta_saved;
432  zeta = -eta_saved;
433  }
434 
435  else if (elem->point(5) == min_point)
436  {
437  if (elem->point(1) == std::min(elem->point(1), elem->point(6)))
438  {
439  // Case 7
440  eta = -zeta_saved;
441  zeta = eta_saved;
442  }
443  else
444  {
445  // Case 8
446  eta = eta_saved;
447  zeta = -zeta_saved;
448  }
449  }
450  }
451  // Face 3
452  else if (i < 8 + 12*e + 4*e*e)
453  {
454  unsigned int basisnum = i - 8 - 12*e - 3*e*e;
455  i0 = square_number_row[basisnum] + 2;
456  i1 = 1;
457  i2 = square_number_column[basisnum] + 2;
458  const Point min_point = get_min_point(elem, 2, 3, 7, 6);
459 
460  if (elem->point(3) == min_point)
461  if (elem->point(2) == std::min(elem->point(2), elem->point(7)))
462  {
463  // Case 1
464  xi = xi_saved;
465  zeta = zeta_saved;
466  }
467  else
468  {
469  // Case 2
470  xi = zeta_saved;
471  zeta = xi_saved;
472  }
473 
474  else if (elem->point(7) == min_point)
475  if (elem->point(3) == std::min(elem->point(3), elem->point(6)))
476  {
477  // Case 3
478  xi = -zeta_saved;
479  zeta = xi_saved;
480  }
481  else
482  {
483  // Case 4
484  xi = xi_saved;
485  zeta = -zeta_saved;
486  }
487 
488  else if (elem->point(6) == min_point)
489  if (elem->point(7) == std::min(elem->point(7), elem->point(2)))
490  {
491  // Case 5
492  xi = -xi_saved;
493  zeta = -zeta_saved;
494  }
495  else
496  {
497  // Case 6
498  xi = -zeta_saved;
499  zeta = -xi_saved;
500  }
501 
502  else if (elem->point(2) == min_point)
503  {
504  if (elem->point(6) == std::min(elem->point(3), elem->point(6)))
505  {
506  // Case 7
507  xi = zeta_saved;
508  zeta = -xi_saved;
509  }
510  else
511  {
512  // Case 8
513  xi = -xi_saved;
514  zeta = zeta_saved;
515  }
516  }
517  }
518  // Face 4
519  else if (i < 8 + 12*e + 5*e*e)
520  {
521  unsigned int basisnum = i - 8 - 12*e - 4*e*e;
522  i0 = 0;
523  i1 = square_number_row[basisnum] + 2;
524  i2 = square_number_column[basisnum] + 2;
525  const Point min_point = get_min_point(elem, 3, 0, 4, 7);
526 
527  if (elem->point(0) == min_point)
528  if (elem->point(3) == std::min(elem->point(3), elem->point(4)))
529  {
530  // Case 1
531  eta = eta_saved;
532  zeta = zeta_saved;
533  }
534  else
535  {
536  // Case 2
537  eta = zeta_saved;
538  zeta = eta_saved;
539  }
540 
541  else if (elem->point(4) == min_point)
542  if (elem->point(0) == std::min(elem->point(0), elem->point(7)))
543  {
544  // Case 3
545  eta = -zeta_saved;
546  zeta = eta_saved;
547  }
548  else
549  {
550  // Case 4
551  eta = eta_saved;
552  zeta = -zeta_saved;
553  }
554 
555  else if (elem->point(7) == min_point)
556  if (elem->point(4) == std::min(elem->point(4), elem->point(3)))
557  {
558  // Case 5
559  eta = -eta_saved;
560  zeta = -zeta_saved;
561  }
562  else
563  {
564  // Case 6
565  eta = -zeta_saved;
566  zeta = -eta_saved;
567  }
568 
569  else if (elem->point(3) == min_point)
570  {
571  if (elem->point(7) == std::min(elem->point(7), elem->point(0)))
572  {
573  // Case 7
574  eta = zeta_saved;
575  zeta = -eta_saved;
576  }
577  else
578  {
579  // Case 8
580  eta = -eta_saved;
581  zeta = zeta_saved;
582  }
583  }
584  }
585  // Face 5
586  else if (i < 8 + 12*e + 6*e*e)
587  {
588  unsigned int basisnum = i - 8 - 12*e - 5*e*e;
589  i0 = square_number_row[basisnum] + 2;
590  i1 = square_number_column[basisnum] + 2;
591  i2 = 1;
592  const Point min_point = get_min_point(elem, 4, 5, 6, 7);
593 
594  if (elem->point(4) == min_point)
595  if (elem->point(5) == std::min(elem->point(5), elem->point(7)))
596  {
597  // Case 1
598  xi = xi_saved;
599  eta = eta_saved;
600  }
601  else
602  {
603  // Case 2
604  xi = eta_saved;
605  eta = xi_saved;
606  }
607 
608  else if (elem->point(5) == min_point)
609  if (elem->point(6) == std::min(elem->point(6), elem->point(4)))
610  {
611  // Case 3
612  xi = eta_saved;
613  eta = -xi_saved;
614  }
615  else
616  {
617  // Case 4
618  xi = -xi_saved;
619  eta = eta_saved;
620  }
621 
622  else if (elem->point(6) == min_point)
623  if (elem->point(7) == std::min(elem->point(7), elem->point(5)))
624  {
625  // Case 5
626  xi = -xi_saved;
627  eta = -eta_saved;
628  }
629  else
630  {
631  // Case 6
632  xi = -eta_saved;
633  eta = -xi_saved;
634  }
635 
636  else if (elem->point(7) == min_point)
637  {
638  if (elem->point(4) == std::min(elem->point(4), elem->point(6)))
639  {
640  // Case 7
641  xi = -eta_saved;
642  eta = xi_saved;
643  }
644  else
645  {
646  // Case 8
647  xi = xi_saved;
648  eta = eta_saved;
649  }
650  }
651  }
652 
653  // Internal DoFs
654  else
655  {
656  unsigned int basisnum = i - 8 - 12*e - 6*e*e;
657  i0 = cube_number_column[basisnum] + 2;
658  i1 = cube_number_row[basisnum] + 2;
659  i2 = cube_number_page[basisnum] + 2;
660  }
661 }
662 #endif // LIBMESH_DIM > 2
663 
664 } // end anonymous namespace
665 
666 
667 
668 namespace libMesh
669 {
670 
671 template <>
672 Real FE<3,HIERARCHIC>::shape(const ElemType,
673  const Order,
674  const unsigned int,
675  const Point &)
676 {
677  libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
678  return 0.;
679 }
680 
681 
682 
683 template <>
684 Real FE<3,L2_HIERARCHIC>::shape(const ElemType,
685  const Order,
686  const unsigned int,
687  const Point &)
688 {
689  libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
690  return 0.;
691 }
692 
693 
694 
695 template <>
696 Real FE<3,HIERARCHIC>::shape(const Elem * elem,
697  const Order order,
698  const unsigned int i,
699  const Point & p,
700  const bool add_p_level)
701 {
702  return fe_hierarchic_3D_shape(elem, order, i, p, add_p_level);
703 }
704 
705 
706 
707 template <>
708 Real FE<3,L2_HIERARCHIC>::shape(const Elem * elem,
709  const Order order,
710  const unsigned int i,
711  const Point & p,
712  const bool add_p_level)
713 {
714  return fe_hierarchic_3D_shape(elem, order, i, p, add_p_level);
715 }
716 
717 
718 
719 template <>
720 Real FE<3,HIERARCHIC>::shape_deriv(const ElemType,
721  const Order,
722  const unsigned int,
723  const unsigned int,
724  const Point & )
725 {
726  libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
727  return 0.;
728 }
729 
730 
731 
732 template <>
733 Real FE<3,L2_HIERARCHIC>::shape_deriv(const ElemType,
734  const Order,
735  const unsigned int,
736  const unsigned int,
737  const Point & )
738 {
739  libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
740  return 0.;
741 }
742 
743 
744 
745 template <>
746 Real FE<3,HIERARCHIC>::shape_deriv(const Elem * elem,
747  const Order order,
748  const unsigned int i,
749  const unsigned int j,
750  const Point & p,
751  const bool add_p_level)
752 {
753  return fe_hierarchic_3D_shape_deriv(elem, order, i, j, p, add_p_level);
754 }
755 
756 
757 template <>
758 Real FE<3,L2_HIERARCHIC>::shape_deriv(const Elem * elem,
759  const Order order,
760  const unsigned int i,
761  const unsigned int j,
762  const Point & p,
763  const bool add_p_level)
764 {
765  return fe_hierarchic_3D_shape_deriv(elem, order, i, j, p, add_p_level);
766 }
767 
768 
769 
770 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
771 
772 template <>
773 Real FE<3,HIERARCHIC>::shape_second_deriv(const ElemType,
774  const Order,
775  const unsigned int,
776  const unsigned int,
777  const Point & )
778 {
779  libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
780  return 0.;
781 }
782 
783 
784 
785 template <>
786 Real FE<3,L2_HIERARCHIC>::shape_second_deriv(const ElemType,
787  const Order,
788  const unsigned int,
789  const unsigned int,
790  const Point & )
791 {
792  libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
793  return 0.;
794 }
795 
796 
797 
798 template <>
799 Real FE<3,HIERARCHIC>::shape_second_deriv(const Elem * elem,
800  const Order order,
801  const unsigned int i,
802  const unsigned int j,
803  const Point & p,
804  const bool add_p_level)
805 {
806  return fe_hierarchic_3D_shape_second_deriv(elem, order, i, j, p, add_p_level);
807 }
808 
809 
810 
811 template <>
812 Real FE<3,L2_HIERARCHIC>::shape_second_deriv(const Elem * elem,
813  const Order order,
814  const unsigned int i,
815  const unsigned int j,
816  const Point & p,
817  const bool add_p_level)
818 {
819  return fe_hierarchic_3D_shape_second_deriv(elem, order, i, j, p, add_p_level);
820 }
821 
822 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
823 
824 } // namespace libMesh
825 
826 
827 
828 namespace
829 {
830 using namespace libMesh;
831 
832 Real fe_hierarchic_3D_shape(const Elem * elem,
833  const Order order,
834  const unsigned int i,
835  const Point & p,
836  const bool add_p_level)
837 {
838 #if LIBMESH_DIM == 3
839 
840  libmesh_assert(elem);
841  const ElemType type = elem->type();
842 
843  const Order totalorder =
844  static_cast<Order>(order+add_p_level*elem->p_level());
845 
846  switch (type)
847  {
848  case HEX8:
849  case HEX20:
850  libmesh_assert_less (totalorder, 2);
851  libmesh_fallthrough();
852  case HEX27:
853  {
854  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+1u));
855 
856  // Compute hex shape functions as a tensor-product
857  Real xi = p(0);
858  Real eta = p(1);
859  Real zeta = p(2);
860 
861  unsigned int i0, i1, i2;
862 
863  cube_indices(elem, totalorder, i, xi, eta, zeta, i0, i1, i2);
864 
865  return (FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i0, xi)*
866  FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i1, eta)*
867  FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i2, zeta));
868  }
869 
870  default:
871  libmesh_error_msg("Invalid element type = " << type);
872  }
873 
874 #else // LIBMESH_DIM != 3
875  libmesh_ignore(elem, order, i, p, add_p_level);
876  libmesh_not_implemented();
877 #endif
878 }
879 
880 
881 
882 Real fe_hierarchic_3D_shape_deriv(const Elem * elem,
883  const Order order,
884  const unsigned int i,
885  const unsigned int j,
886  const Point & p,
887  const bool add_p_level)
888 {
889 #if LIBMESH_DIM == 3
890  libmesh_assert(elem);
891 
892  libmesh_assert_less (j, 3);
893 
894  // cheat by using finite difference approximations:
895  const Real eps = 1.e-6;
896  Point pp, pm;
897 
898  switch (j)
899  {
900  // d()/dxi
901  case 0:
902  {
903  pp = Point(p(0)+eps, p(1), p(2));
904  pm = Point(p(0)-eps, p(1), p(2));
905  break;
906  }
907 
908  // d()/deta
909  case 1:
910  {
911  pp = Point(p(0), p(1)+eps, p(2));
912  pm = Point(p(0), p(1)-eps, p(2));
913  break;
914  }
915 
916  // d()/dzeta
917  case 2:
918  {
919  pp = Point(p(0), p(1), p(2)+eps);
920  pm = Point(p(0), p(1), p(2)-eps);
921  break;
922  }
923 
924  default:
925  libmesh_error_msg("Invalid derivative index j = " << j);
926  }
927 
928  return (FE<3,HIERARCHIC>::shape(elem, order, i, pp, add_p_level) -
929  FE<3,HIERARCHIC>::shape(elem, order, i, pm, add_p_level))/2./eps;
930 #else // LIBMESH_DIM != 3
931  libmesh_ignore(elem, order, i, j, p, add_p_level);
932  libmesh_not_implemented();
933 #endif
934 }
935 
936 
937 
938 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
939 
940 Real fe_hierarchic_3D_shape_second_deriv(const Elem * elem,
941  const Order order,
942  const unsigned int i,
943  const unsigned int j,
944  const Point & p,
945  const bool add_p_level)
946 {
947  libmesh_assert(elem);
948 
949  const Real eps = 1.e-6;
950  Point pp, pm;
951  unsigned int prevj = libMesh::invalid_uint;
952 
953  switch (j)
954  {
955  // d^2()/dxi^2
956  case 0:
957  {
958  pp = Point(p(0)+eps, p(1), p(2));
959  pm = Point(p(0)-eps, p(1), p(2));
960  prevj = 0;
961  break;
962  }
963 
964  // d^2()/dxideta
965  case 1:
966  {
967  pp = Point(p(0), p(1)+eps, p(2));
968  pm = Point(p(0), p(1)-eps, p(2));
969  prevj = 0;
970  break;
971  }
972 
973  // d^2()/deta^2
974  case 2:
975  {
976  pp = Point(p(0), p(1)+eps, p(2));
977  pm = Point(p(0), p(1)-eps, p(2));
978  prevj = 1;
979  break;
980  }
981 
982  // d^2()/dxidzeta
983  case 3:
984  {
985  pp = Point(p(0), p(1), p(2)+eps);
986  pm = Point(p(0), p(1), p(2)-eps);
987  prevj = 0;
988  break;
989  }
990 
991  // d^2()/detadzeta
992  case 4:
993  {
994  pp = Point(p(0), p(1), p(2)+eps);
995  pm = Point(p(0), p(1), p(2)-eps);
996  prevj = 1;
997  break;
998  }
999 
1000  // d^2()/dzeta^2
1001  case 5:
1002  {
1003  pp = Point(p(0), p(1), p(2)+eps);
1004  pm = Point(p(0), p(1), p(2)-eps);
1005  prevj = 2;
1006  break;
1007  }
1008  default:
1009  libmesh_error_msg("Invalid derivative index j = " << j);
1010  }
1011 
1012  return (FE<3,HIERARCHIC>::shape_deriv(elem, order, i, prevj, pp, add_p_level) -
1013  FE<3,HIERARCHIC>::shape_deriv(elem, order, i, prevj, pm, add_p_level))
1014  / 2. / eps;
1015 }
1016 
1017 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1018 
1019 
1020 } // anonymous namespace
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)
libMesh::HEX20
Definition: enum_elem_type.h:48
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value.
Definition: libmesh.h:249
libMesh::HEX8
Definition: enum_elem_type.h:47
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
Definition: factoryfunction.C:55
libMesh::Order
Order
Definition: enum_order.h:40
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2512
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:1955
libMesh::square_number_column
const unsigned char square_number_column[]
Definition: number_lookups.C:56
libMesh::cube_number_page
const unsigned char cube_number_page[]
Definition: number_lookups.C:308
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::cube_number_column
const unsigned char cube_number_column[]
Definition: number_lookups.C:84
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::HEX27
Definition: enum_elem_type.h:49
libMesh::cube_number_row
const unsigned char cube_number_row[]
Definition: number_lookups.C:196
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::FE
A specific instantiation of the FEBase class.
Definition: fe.h:95
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:526
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
libMesh::EDGE3
Definition: enum_elem_type.h:36
libMesh::square_number_row
const unsigned char square_number_row[]
Definition: number_lookups.C:69
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:100
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
Generated on Sat Jan 25 2020 12:06:57 for libMesh by   doxygen 1.8.16