fe_hierarchic_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++ 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::err << "Hierarchic polynomials require the element type\n"
651  << "because edge and face orientation is needed."
652  << std::endl;
653 
654  libmesh_error();
655  return 0.;
656 }
657 
658 
659 
660 template <>
662  const Order order,
663  const unsigned int i,
664  const Point& p)
665 {
666 #if LIBMESH_DIM == 3
667 
668  libmesh_assert(elem);
669  const ElemType type = elem->type();
670 
671  const Order totalorder = static_cast<Order>(order+elem->p_level());
672 
673  switch (type)
674  {
675  case HEX8:
676  case HEX20:
677  libmesh_assert_less (totalorder, 2);
678  case HEX27:
679  {
680  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+1u));
681 
682  // Compute hex shape functions as a tensor-product
683  Real xi = p(0);
684  Real eta = p(1);
685  Real zeta = p(2);
686 
687  unsigned int i0, i1, i2;
688 
689  cube_indices(elem, totalorder, i, xi, eta, zeta, i0, i1, i2);
690 
691  return (FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i0, xi)*
692  FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i1, eta)*
693  FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i2, zeta));
694  }
695 
696  default:
697  libmesh_error();
698  }
699 
700 #endif
701 
702  libmesh_error();
703  return 0.;
704 }
705 
706 
707 
708 
709 template <>
711  const Order,
712  const unsigned int,
713  const unsigned int,
714  const Point& )
715 {
716  libMesh::err << "Hierarchic polynomials require the element type\n"
717  << "because edge and face orientation is needed."
718  << std::endl;
719  libmesh_error();
720 
721  return 0.;
722 }
723 
724 
725 
726 template <>
728  const Order order,
729  const unsigned int i,
730  const unsigned int j,
731  const Point& p)
732 {
733 #if LIBMESH_DIM == 3
734  libmesh_assert(elem);
735 
736  libmesh_assert_less (j, 3);
737 
738  // cheat by using finite difference approximations:
739  const Real eps = 1.e-6;
740  Point pp, pm;
741 
742  switch (j)
743  {
744  // d()/dxi
745  case 0:
746  {
747  pp = Point(p(0)+eps, p(1), p(2));
748  pm = Point(p(0)-eps, p(1), p(2));
749  break;
750  }
751 
752  // d()/deta
753  case 1:
754  {
755  pp = Point(p(0), p(1)+eps, p(2));
756  pm = Point(p(0), p(1)-eps, p(2));
757  break;
758  }
759 
760  // d()/dzeta
761  case 2:
762  {
763  pp = Point(p(0), p(1), p(2)+eps);
764  pm = Point(p(0), p(1), p(2)-eps);
765  break;
766  }
767 
768  default:
769  libmesh_error();
770  }
771 
772  return (FE<3,HIERARCHIC>::shape(elem, order, i, pp) -
773  FE<3,HIERARCHIC>::shape(elem, order, i, pm))/2./eps;
774 #endif
775 
776  libmesh_error();
777  return 0.;
778 }
779 
780 
781 
782 template <>
784  const Order,
785  const unsigned int,
786  const unsigned int,
787  const Point& )
788 {
789  libMesh::err << "Hierarchic polynomials require the element type\n"
790  << "because edge and face orientation is needed."
791  << std::endl;
792  libmesh_error();
793 
794  return 0.;
795 }
796 
797 
798 
799 template <>
801  const Order order,
802  const unsigned int i,
803  const unsigned int j,
804  const Point& p)
805 {
806  libmesh_assert(elem);
807 
808  const Real eps = 1.e-6;
809  Point pp, pm;
810  unsigned int prevj = libMesh::invalid_uint;
811 
812  switch (j)
813  {
814  // d^2()/dxi^2
815  case 0:
816  {
817  pp = Point(p(0)+eps, p(1), p(2));
818  pm = Point(p(0)-eps, p(1), p(2));
819  prevj = 0;
820  break;
821  }
822 
823  // d^2()/dxideta
824  case 1:
825  {
826  pp = Point(p(0), p(1)+eps, p(2));
827  pm = Point(p(0), p(1)-eps, p(2));
828  prevj = 0;
829  break;
830  }
831 
832  // d^2()/deta^2
833  case 2:
834  {
835  pp = Point(p(0), p(1)+eps, p(2));
836  pm = Point(p(0), p(1)-eps, p(2));
837  prevj = 1;
838  break;
839  }
840 
841  // d^2()/dxidzeta
842  case 3:
843  {
844  pp = Point(p(0), p(1), p(2)+eps);
845  pm = Point(p(0), p(1), p(2)-eps);
846  prevj = 0;
847  break;
848  }
849 
850  // d^2()/detadzeta
851  case 4:
852  {
853  pp = Point(p(0), p(1), p(2)+eps);
854  pm = Point(p(0), p(1), p(2)-eps);
855  prevj = 1;
856  break;
857  }
858 
859  // d^2()/dzeta^2
860  case 5:
861  {
862  pp = Point(p(0), p(1), p(2)+eps);
863  pm = Point(p(0), p(1), p(2)-eps);
864  prevj = 2;
865  break;
866  }
867  default:
868  libmesh_error();
869  }
870  return (FE<3,HIERARCHIC>::shape_deriv(elem, order, i, prevj, pp) -
871  FE<3,HIERARCHIC>::shape_deriv(elem, order, i, prevj, pm))
872  / 2. / eps;
873 }
874 
875 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo