fe_hierarchic_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 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 <>
645 Real FE<3,HIERARCHIC>::shape(const ElemType,
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 <>
657 Real FE<3,HIERARCHIC>::shape(const Elem * elem,
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 
700 
701 
702 
703 template <>
704 Real FE<3,HIERARCHIC>::shape_deriv(const ElemType,
705  const Order,
706  const unsigned int,
707  const unsigned int,
708  const Point & )
709 {
710  libmesh_error_msg("Hierarchic polynomials require the element type \nbecause edge and face orientation is needed.");
711  return 0.;
712 }
713 
714 
715 
716 template <>
717 Real FE<3,HIERARCHIC>::shape_deriv(const Elem * elem,
718  const Order order,
719  const unsigned int i,
720  const unsigned int j,
721  const Point & p)
722 {
723 #if LIBMESH_DIM == 3
724  libmesh_assert(elem);
725 
726  libmesh_assert_less (j, 3);
727 
728  // cheat by using finite difference approximations:
729  const Real eps = 1.e-6;
730  Point pp, pm;
731 
732  switch (j)
733  {
734  // d()/dxi
735  case 0:
736  {
737  pp = Point(p(0)+eps, p(1), p(2));
738  pm = Point(p(0)-eps, p(1), p(2));
739  break;
740  }
741 
742  // d()/deta
743  case 1:
744  {
745  pp = Point(p(0), p(1)+eps, p(2));
746  pm = Point(p(0), p(1)-eps, p(2));
747  break;
748  }
749 
750  // d()/dzeta
751  case 2:
752  {
753  pp = Point(p(0), p(1), p(2)+eps);
754  pm = Point(p(0), p(1), p(2)-eps);
755  break;
756  }
757 
758  default:
759  libmesh_error_msg("Invalid derivative index j = " << j);
760  }
761 
762  return (FE<3,HIERARCHIC>::shape(elem, order, i, pp) -
763  FE<3,HIERARCHIC>::shape(elem, order, i, pm))/2./eps;
764 #endif
765 }
766 
767 
768 
769 template <>
770 Real FE<3,HIERARCHIC>::shape_second_deriv(const ElemType,
771  const Order,
772  const unsigned int,
773  const unsigned int,
774  const Point & )
775 {
776  libmesh_error_msg("Hierarchic polynomials require the element type \nbecause edge and face orientation is needed.");
777  return 0.;
778 }
779 
780 
781 
782 template <>
783 Real FE<3,HIERARCHIC>::shape_second_deriv(const Elem * elem,
784  const Order order,
785  const unsigned int i,
786  const unsigned int j,
787  const Point & p)
788 {
789  libmesh_assert(elem);
790 
791  const Real eps = 1.e-6;
792  Point pp, pm;
793  unsigned int prevj = libMesh::invalid_uint;
794 
795  switch (j)
796  {
797  // d^2()/dxi^2
798  case 0:
799  {
800  pp = Point(p(0)+eps, p(1), p(2));
801  pm = Point(p(0)-eps, p(1), p(2));
802  prevj = 0;
803  break;
804  }
805 
806  // d^2()/dxideta
807  case 1:
808  {
809  pp = Point(p(0), p(1)+eps, p(2));
810  pm = Point(p(0), p(1)-eps, p(2));
811  prevj = 0;
812  break;
813  }
814 
815  // d^2()/deta^2
816  case 2:
817  {
818  pp = Point(p(0), p(1)+eps, p(2));
819  pm = Point(p(0), p(1)-eps, p(2));
820  prevj = 1;
821  break;
822  }
823 
824  // d^2()/dxidzeta
825  case 3:
826  {
827  pp = Point(p(0), p(1), p(2)+eps);
828  pm = Point(p(0), p(1), p(2)-eps);
829  prevj = 0;
830  break;
831  }
832 
833  // d^2()/detadzeta
834  case 4:
835  {
836  pp = Point(p(0), p(1), p(2)+eps);
837  pm = Point(p(0), p(1), p(2)-eps);
838  prevj = 1;
839  break;
840  }
841 
842  // d^2()/dzeta^2
843  case 5:
844  {
845  pp = Point(p(0), p(1), p(2)+eps);
846  pm = Point(p(0), p(1), p(2)-eps);
847  prevj = 2;
848  break;
849  }
850  default:
851  libmesh_error_msg("Invalid derivative index j = " << j);
852  }
853 
854  return (FE<3,HIERARCHIC>::shape_deriv(elem, order, i, prevj, pp) -
855  FE<3,HIERARCHIC>::shape_deriv(elem, order, i, prevj, pm))
856  / 2. / eps;
857 }
858 
859 } // namespace libMesh
libMesh::invalid_uint
const unsigned int invalid_uint
Definition: libmesh.h:245
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Elem
The base class for all geometric element types.
Definition: elem.h:100
libMesh::cube_number_column
const unsigned char cube_number_column[]
Definition: number_lookups.C:84
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::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2555
libMesh::cube_number_row
const unsigned char cube_number_row[]
Definition: number_lookups.C:196
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::FE
Template class which generates the different FE families and orders.
Definition: fe.h:89
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::square_number_row
const unsigned char square_number_row[]
Definition: number_lookups.C:69
number_lookups.h
libMesh::Elem::type
virtual ElemType type() const =0
std::min
long double min(long double a, double b)
Definition: libmesh_augment_std_namespace.h:29
libMesh::Point
A geometric point in (x,y,z) space.
Definition: point.h:38
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)
generated by doxygen