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