fe_l2_lagrange_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
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
25 namespace libMesh
26 {
27
28
29
30
31 template <>
33  const Order order,
34  const unsigned int i,
35  const Point & p)
36 {
37 #if LIBMESH_DIM == 3
38
39
40  switch (order)
41  {
42  // linear Lagrange shape functions
43  case FIRST:
44  {
45  switch (type)
46  {
47  // trilinear hexahedral shape functions
48  case HEX8:
49  case HEX20:
50  case HEX27:
51  {
52  libmesh_assert_less (i, 8);
53
54  // Compute hex shape functions as a tensor-product
55  const Real xi = p(0);
56  const Real eta = p(1);
57  const Real zeta = p(2);
58
59  // 0 1 2 3 4 5 6 7
60  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
61  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
62  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
63
64  return (FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i0[i], xi)*
65  FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i1[i], eta)*
66  FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i2[i], zeta));
67  }
68
69  // linear tetrahedral shape functions
70  case TET4:
71  case TET10:
72  {
73  libmesh_assert_less (i, 4);
74
75  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
76  const Real zeta1 = p(0);
77  const Real zeta2 = p(1);
78  const Real zeta3 = p(2);
79  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
80
81  switch(i)
82  {
83  case 0:
84  return zeta0;
85
86  case 1:
87  return zeta1;
88
89  case 2:
90  return zeta2;
91
92  case 3:
93  return zeta3;
94
95  default:
96  libmesh_error_msg("Invalid i = " << i);
97  }
98  }
99
100  // linear prism shape functions
101  case PRISM6:
102  case PRISM15:
103  case PRISM18:
104  {
105  libmesh_assert_less (i, 6);
106
107  // Compute prism shape functions as a tensor-product
108  // of a triangle and an edge
109
110  Point p2d(p(0),p(1));
111  Point p1d(p(2));
112
113  // 0 1 2 3 4 5
114  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1};
115  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2};
116
117  return (FE<2,L2_LAGRANGE>::shape(TRI3, FIRST, i1[i], p2d)*
118  FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d));
119  }
120
121  // linear pyramid shape functions
122  case PYRAMID5:
123  {
124  libmesh_assert_less (i, 5);
125
126  const Real xi = p(0);
127  const Real eta = p(1);
128  const Real zeta = p(2);
129  const Real eps = 1.e-35;
130
131  switch(i)
132  {
133  case 0:
134  return .25*(zeta + xi - 1.)*(zeta + eta - 1.)/((1. - zeta) + eps);
135
136  case 1:
137  return .25*(zeta - xi - 1.)*(zeta + eta - 1.)/((1. - zeta) + eps);
138
139  case 2:
140  return .25*(zeta - xi - 1.)*(zeta - eta - 1.)/((1. - zeta) + eps);
141
142  case 3:
143  return .25*(zeta + xi - 1.)*(zeta - eta - 1.)/((1. - zeta) + eps);
144
145  case 4:
146  return zeta;
147
148  default:
149  libmesh_error_msg("Invalid i = " << i);
150  }
151  }
152
153
154  default:
155  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
156  }
157  }
158
159
160  // quadratic Lagrange shape functions
161  case SECOND:
162  {
163  switch (type)
164  {
165
166  // serendipity hexahedral quadratic shape functions
167  case HEX20:
168  {
169  libmesh_assert_less (i, 20);
170
171  const Real xi = p(0);
172  const Real eta = p(1);
173  const Real zeta = p(2);
174
175  // these functions are defined for (x,y,z) in [0,1]^3
176  // so transform the locations
177  const Real x = .5*(xi + 1.);
178  const Real y = .5*(eta + 1.);
179  const Real z = .5*(zeta + 1.);
180
181  switch (i)
182  {
183  case 0:
184  return (1. - x)*(1. - y)*(1. - z)*(1. - 2.*x - 2.*y - 2.*z);
185
186  case 1:
187  return x*(1. - y)*(1. - z)*(2.*x - 2.*y - 2.*z - 1.);
188
189  case 2:
190  return x*y*(1. - z)*(2.*x + 2.*y - 2.*z - 3.);
191
192  case 3:
193  return (1. - x)*y*(1. - z)*(2.*y - 2.*x - 2.*z - 1.);
194
195  case 4:
196  return (1. - x)*(1. - y)*z*(2.*z - 2.*x - 2.*y - 1.);
197
198  case 5:
199  return x*(1. - y)*z*(2.*x - 2.*y + 2.*z - 3.);
200
201  case 6:
202  return x*y*z*(2.*x + 2.*y + 2.*z - 5.);
203
204  case 7:
205  return (1. - x)*y*z*(2.*y - 2.*x + 2.*z - 3.);
206
207  case 8:
208  return 4.*x*(1. - x)*(1. - y)*(1. - z);
209
210  case 9:
211  return 4.*x*y*(1. - y)*(1. - z);
212
213  case 10:
214  return 4.*x*(1. - x)*y*(1. - z);
215
216  case 11:
217  return 4.*(1. - x)*y*(1. - y)*(1. - z);
218
219  case 12:
220  return 4.*(1. - x)*(1. - y)*z*(1. - z);
221
222  case 13:
223  return 4.*x*(1. - y)*z*(1. - z);
224
225  case 14:
226  return 4.*x*y*z*(1. - z);
227
228  case 15:
229  return 4.*(1. - x)*y*z*(1. - z);
230
231  case 16:
232  return 4.*x*(1. - x)*(1. - y)*z;
233
234  case 17:
235  return 4.*x*y*(1. - y)*z;
236
237  case 18:
238  return 4.*x*(1. - x)*y*z;
239
240  case 19:
241  return 4.*(1. - x)*y*(1. - y)*z;
242
243  default:
244  libmesh_error_msg("Invalid i = " << i);
245  }
246  }
247
248  // triquadratic hexahedral shape functions
249  case HEX27:
250  {
251  libmesh_assert_less (i, 27);
252
253  // Compute hex shape functions as a tensor-product
254  const Real xi = p(0);
255  const Real eta = p(1);
256  const Real zeta = p(2);
257
258  // The only way to make any sense of this
259  // is to look at the mgflo/mg2/mgf documentation
260  // and make the cut-out cube!
261  // 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
262  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
263  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
264  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2};
265
266  return (FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], xi)*
267  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i1[i], eta)*
268  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i2[i], zeta));
269  }
270
271  // quadratic tetrahedral shape functions
272  case TET10:
273  {
274  libmesh_assert_less (i, 10);
275
276  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
277  const Real zeta1 = p(0);
278  const Real zeta2 = p(1);
279  const Real zeta3 = p(2);
280  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
281
282  switch(i)
283  {
284  case 0:
285  return zeta0*(2.*zeta0 - 1.);
286
287  case 1:
288  return zeta1*(2.*zeta1 - 1.);
289
290  case 2:
291  return zeta2*(2.*zeta2 - 1.);
292
293  case 3:
294  return zeta3*(2.*zeta3 - 1.);
295
296  case 4:
297  return 4.*zeta0*zeta1;
298
299  case 5:
300  return 4.*zeta1*zeta2;
301
302  case 6:
303  return 4.*zeta2*zeta0;
304
305  case 7:
306  return 4.*zeta0*zeta3;
307
308  case 8:
309  return 4.*zeta1*zeta3;
310
311  case 9:
312  return 4.*zeta2*zeta3;
313
314  default:
315  libmesh_error_msg("Invalid i = " << i);
316  }
317  }
318
319  // quadratic prism shape functions
320  case PRISM18:
321  {
322  libmesh_assert_less (i, 18);
323
324  // Compute prism shape functions as a tensor-product
325  // of a triangle and an edge
326
327  Point p2d(p(0),p(1));
328  Point p1d(p(2));
329
330  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
331  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
332  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};
333
334  return (FE<2,L2_LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
335  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
336  }
337
338
339  default:
340  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
341  }
342  }
343
344
345  // unsupported order
346  default:
347  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order);
348  }
349
350 #endif
351 }
352
353
354
355 template <>
357  const Order order,
358  const unsigned int i,
359  const Point & p)
360 {
361  libmesh_assert(elem);
362
363  // call the orientation-independent shape functions
364  return FE<3,L2_LAGRANGE>::shape(elem->type(), static_cast<Order>(order + elem->p_level()), i, p);
365 }
366
367
368
369
370 template <>
372  const Order order,
373  const unsigned int i,
374  const unsigned int j,
375  const Point & p)
376 {
377 #if LIBMESH_DIM == 3
378
379  libmesh_assert_less (j, 3);
380
381  switch (order)
382  {
383  // linear Lagrange shape functions
384  case FIRST:
385  {
386  switch (type)
387  {
388  // trilinear hexahedral shape functions
389  case HEX8:
390  case HEX20:
391  case HEX27:
392  {
393  libmesh_assert_less (i, 8);
394
395  // Compute hex shape functions as a tensor-product
396  const Real xi = p(0);
397  const Real eta = p(1);
398  const Real zeta = p(2);
399
400  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
401  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
402  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
403
404  switch(j)
405  {
406  case 0:
407  return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, xi)*
408  FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i1[i], eta)*
409  FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i2[i], zeta));
410
411  case 1:
412  return (FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i0[i], xi)*
413  FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i1[i], 0, eta)*
414  FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i2[i], zeta));
415
416  case 2:
417  return (FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i0[i], xi)*
418  FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i1[i], eta)*
419  FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i2[i], 0, zeta));
420
421  default:
422  libmesh_error_msg("Invalid j = " << j);
423  }
424  }
425
426  // linear tetrahedral shape functions
427  case TET4:
428  case TET10:
429  {
430  libmesh_assert_less (i, 4);
431
432  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
433  const Real dzeta0dxi = -1.;
434  const Real dzeta1dxi = 1.;
435  const Real dzeta2dxi = 0.;
436  const Real dzeta3dxi = 0.;
437
438  const Real dzeta0deta = -1.;
439  const Real dzeta1deta = 0.;
440  const Real dzeta2deta = 1.;
441  const Real dzeta3deta = 0.;
442
443  const Real dzeta0dzeta = -1.;
444  const Real dzeta1dzeta = 0.;
445  const Real dzeta2dzeta = 0.;
446  const Real dzeta3dzeta = 1.;
447
448  switch (j)
449  {
450  // d()/dxi
451  case 0:
452  {
453  switch(i)
454  {
455  case 0:
456  return dzeta0dxi;
457
458  case 1:
459  return dzeta1dxi;
460
461  case 2:
462  return dzeta2dxi;
463
464  case 3:
465  return dzeta3dxi;
466
467  default:
468  libmesh_error_msg("Invalid i = " << i);
469  }
470  }
471
472  // d()/deta
473  case 1:
474  {
475  switch(i)
476  {
477  case 0:
478  return dzeta0deta;
479
480  case 1:
481  return dzeta1deta;
482
483  case 2:
484  return dzeta2deta;
485
486  case 3:
487  return dzeta3deta;
488
489  default:
490  libmesh_error_msg("Invalid i = " << i);
491  }
492  }
493
494  // d()/dzeta
495  case 2:
496  {
497  switch(i)
498  {
499  case 0:
500  return dzeta0dzeta;
501
502  case 1:
503  return dzeta1dzeta;
504
505  case 2:
506  return dzeta2dzeta;
507
508  case 3:
509  return dzeta3dzeta;
510
511  default:
512  libmesh_error_msg("Invalid i = " << i);
513  }
514  }
515
516  default:
517  libmesh_error_msg("Invalid shape function derivative j = " << j);
518  }
519  }
520
521  // linear prism shape functions
522  case PRISM6:
523  case PRISM15:
524  case PRISM18:
525  {
526  libmesh_assert_less (i, 6);
527
528  // Compute prism shape functions as a tensor-product
529  // of a triangle and an edge
530
531  Point p2d(p(0),p(1));
532  Point p1d(p(2));
533
534  // 0 1 2 3 4 5
535  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1};
536  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2};
537
538  switch (j)
539  {
540  // d()/dxi
541  case 0:
542  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 0, p2d)*
543  FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d));
544
545  // d()/deta
546  case 1:
547  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 1, p2d)*
548  FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d));
549
550  // d()/dzeta
551  case 2:
552  return (FE<2,L2_LAGRANGE>::shape(TRI3, FIRST, i1[i], p2d)*
553  FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, p1d));
554
555  default:
556  libmesh_error_msg("Invalid shape function derivative j = " << j);
557  }
558  }
559
560  // linear pyramid shape functions
561  case PYRAMID5:
562  {
563  libmesh_assert_less (i, 5);
564
565  const Real xi = p(0);
566  const Real eta = p(1);
567  const Real zeta = p(2);
568  const Real eps = 1.e-35;
569
570  switch (j)
571  {
572  // d/dxi
573  case 0:
574  switch(i)
575  {
576  case 0:
577  return .25*(zeta + eta - 1.)/((1. - zeta) + eps);
578
579  case 1:
580  return -.25*(zeta + eta - 1.)/((1. - zeta) + eps);
581
582  case 2:
583  return -.25*(zeta - eta - 1.)/((1. - zeta) + eps);
584
585  case 3:
586  return .25*(zeta - eta - 1.)/((1. - zeta) + eps);
587
588  case 4:
589  return 0;
590
591  default:
592  libmesh_error_msg("Invalid i = " << i);
593  }
594
595
596  // d/deta
597  case 1:
598  switch(i)
599  {
600  case 0:
601  return .25*(zeta + xi - 1.)/((1. - zeta) + eps);
602
603  case 1:
604  return .25*(zeta - xi - 1.)/((1. - zeta) + eps);
605
606  case 2:
607  return -.25*(zeta - xi - 1.)/((1. - zeta) + eps);
608
609  case 3:
610  return -.25*(zeta + xi - 1.)/((1. - zeta) + eps);
611
612  case 4:
613  return 0;
614
615  default:
616  libmesh_error_msg("Invalid i = " << i);
617  }
618
619
620  // d/dzeta
621  case 2:
622  switch(i)
623  {
624  case 0:
625  {
626  const Real a=1.;
627  const Real b=1.;
628
629  return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
630  (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
631  ((1. - zeta)*(1. - zeta) + eps));
632  }
633
634  case 1:
635  {
636  const Real a=-1.;
637  const Real b=1.;
638
639  return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
640  (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
641  ((1. - zeta)*(1. - zeta) + eps));
642  }
643
644  case 2:
645  {
646  const Real a=-1.;
647  const Real b=-1.;
648
649  return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
650  (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
651  ((1. - zeta)*(1. - zeta) + eps));
652  }
653
654  case 3:
655  {
656  const Real a=1.;
657  const Real b=-1.;
658
659  return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
660  (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
661  ((1. - zeta)*(1. - zeta) + eps));
662  }
663
664  case 4:
665  return 1.;
666
667  default:
668  libmesh_error_msg("Invalid i = " << i);
669  }
670
671
672  default:
673  libmesh_error_msg("Invalid j = " << j);
674  }
675  }
676
677
678  default:
679  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
680  }
681  }
682
683
684  // quadratic Lagrange shape functions
685  case SECOND:
686  {
687  switch (type)
688  {
689
690  // serendipity hexahedral quadratic shape functions
691  case HEX20:
692  {
693  libmesh_assert_less (i, 20);
694
695  const Real xi = p(0);
696  const Real eta = p(1);
697  const Real zeta = p(2);
698
699  // these functions are defined for (x,y,z) in [0,1]^3
700  // so transform the locations
701  const Real x = .5*(xi + 1.);
702  const Real y = .5*(eta + 1.);
703  const Real z = .5*(zeta + 1.);
704
705  // and don't forget the chain rule!
706
707  switch (j)
708  {
709
710  // d/dx*dx/dxi
711  case 0:
712  switch (i)
713  {
714  case 0:
715  return .5*(1. - y)*(1. - z)*((1. - x)*(-2.) +
716  (-1.)*(1. - 2.*x - 2.*y - 2.*z));
717
718  case 1:
719  return .5*(1. - y)*(1. - z)*(x*(2.) +
720  (1.)*(2.*x - 2.*y - 2.*z - 1.));
721
722  case 2:
723  return .5*y*(1. - z)*(x*(2.) +
724  (1.)*(2.*x + 2.*y - 2.*z - 3.));
725
726  case 3:
727  return .5*y*(1. - z)*((1. - x)*(-2.) +
728  (-1.)*(2.*y - 2.*x - 2.*z - 1.));
729
730  case 4:
731  return .5*(1. - y)*z*((1. - x)*(-2.) +
732  (-1.)*(2.*z - 2.*x - 2.*y - 1.));
733
734  case 5:
735  return .5*(1. - y)*z*(x*(2.) +
736  (1.)*(2.*x - 2.*y + 2.*z - 3.));
737
738  case 6:
739  return .5*y*z*(x*(2.) +
740  (1.)*(2.*x + 2.*y + 2.*z - 5.));
741
742  case 7:
743  return .5*y*z*((1. - x)*(-2.) +
744  (-1.)*(2.*y - 2.*x + 2.*z - 3.));
745
746  case 8:
747  return 2.*(1. - y)*(1. - z)*(1. - 2.*x);
748
749  case 9:
750  return 2.*y*(1. - y)*(1. - z);
751
752  case 10:
753  return 2.*y*(1. - z)*(1. - 2.*x);
754
755  case 11:
756  return 2.*y*(1. - y)*(1. - z)*(-1.);
757
758  case 12:
759  return 2.*(1. - y)*z*(1. - z)*(-1.);
760
761  case 13:
762  return 2.*(1. - y)*z*(1. - z);
763
764  case 14:
765  return 2.*y*z*(1. - z);
766
767  case 15:
768  return 2.*y*z*(1. - z)*(-1.);
769
770  case 16:
771  return 2.*(1. - y)*z*(1. - 2.*x);
772
773  case 17:
774  return 2.*y*(1. - y)*z;
775
776  case 18:
777  return 2.*y*z*(1. - 2.*x);
778
779  case 19:
780  return 2.*y*(1. - y)*z*(-1.);
781
782  default:
783  libmesh_error_msg("Invalid i = " << i);
784  }
785
786
787  // d/dy*dy/deta
788  case 1:
789  switch (i)
790  {
791  case 0:
792  return .5*(1. - x)*(1. - z)*((1. - y)*(-2.) +
793  (-1.)*(1. - 2.*x - 2.*y - 2.*z));
794
795  case 1:
796  return .5*x*(1. - z)*((1. - y)*(-2.) +
797  (-1.)*(2.*x - 2.*y - 2.*z - 1.));
798
799  case 2:
800  return .5*x*(1. - z)*(y*(2.) +
801  (1.)*(2.*x + 2.*y - 2.*z - 3.));
802
803  case 3:
804  return .5*(1. - x)*(1. - z)*(y*(2.) +
805  (1.)*(2.*y - 2.*x - 2.*z - 1.));
806
807  case 4:
808  return .5*(1. - x)*z*((1. - y)*(-2.) +
809  (-1.)*(2.*z - 2.*x - 2.*y - 1.));
810
811  case 5:
812  return .5*x*z*((1. - y)*(-2.) +
813  (-1.)*(2.*x - 2.*y + 2.*z - 3.));
814
815  case 6:
816  return .5*x*z*(y*(2.) +
817  (1.)*(2.*x + 2.*y + 2.*z - 5.));
818
819  case 7:
820  return .5*(1. - x)*z*(y*(2.) +
821  (1.)*(2.*y - 2.*x + 2.*z - 3.));
822
823  case 8:
824  return 2.*x*(1. - x)*(1. - z)*(-1.);
825
826  case 9:
827  return 2.*x*(1. - z)*(1. - 2.*y);
828
829  case 10:
830  return 2.*x*(1. - x)*(1. - z);
831
832  case 11:
833  return 2.*(1. - x)*(1. - z)*(1. - 2.*y);
834
835  case 12:
836  return 2.*(1. - x)*z*(1. - z)*(-1.);
837
838  case 13:
839  return 2.*x*z*(1. - z)*(-1.);
840
841  case 14:
842  return 2.*x*z*(1. - z);
843
844  case 15:
845  return 2.*(1. - x)*z*(1. - z);
846
847  case 16:
848  return 2.*x*(1. - x)*z*(-1.);
849
850  case 17:
851  return 2.*x*z*(1. - 2.*y);
852
853  case 18:
854  return 2.*x*(1. - x)*z;
855
856  case 19:
857  return 2.*(1. - x)*z*(1. - 2.*y);
858
859  default:
860  libmesh_error_msg("Invalid i = " << i);
861  }
862
863
864  // d/dz*dz/dzeta
865  case 2:
866  switch (i)
867  {
868  case 0:
869  return .5*(1. - x)*(1. - y)*((1. - z)*(-2.) +
870  (-1.)*(1. - 2.*x - 2.*y - 2.*z));
871
872  case 1:
873  return .5*x*(1. - y)*((1. - z)*(-2.) +
874  (-1.)*(2.*x - 2.*y - 2.*z - 1.));
875
876  case 2:
877  return .5*x*y*((1. - z)*(-2.) +
878  (-1.)*(2.*x + 2.*y - 2.*z - 3.));
879
880  case 3:
881  return .5*(1. - x)*y*((1. - z)*(-2.) +
882  (-1.)*(2.*y - 2.*x - 2.*z - 1.));
883
884  case 4:
885  return .5*(1. - x)*(1. - y)*(z*(2.) +
886  (1.)*(2.*z - 2.*x - 2.*y - 1.));
887
888  case 5:
889  return .5*x*(1. - y)*(z*(2.) +
890  (1.)*(2.*x - 2.*y + 2.*z - 3.));
891
892  case 6:
893  return .5*x*y*(z*(2.) +
894  (1.)*(2.*x + 2.*y + 2.*z - 5.));
895
896  case 7:
897  return .5*(1. - x)*y*(z*(2.) +
898  (1.)*(2.*y - 2.*x + 2.*z - 3.));
899
900  case 8:
901  return 2.*x*(1. - x)*(1. - y)*(-1.);
902
903  case 9:
904  return 2.*x*y*(1. - y)*(-1.);
905
906  case 10:
907  return 2.*x*(1. - x)*y*(-1.);
908
909  case 11:
910  return 2.*(1. - x)*y*(1. - y)*(-1.);
911
912  case 12:
913  return 2.*(1. - x)*(1. - y)*(1. - 2.*z);
914
915  case 13:
916  return 2.*x*(1. - y)*(1. - 2.*z);
917
918  case 14:
919  return 2.*x*y*(1. - 2.*z);
920
921  case 15:
922  return 2.*(1. - x)*y*(1. - 2.*z);
923
924  case 16:
925  return 2.*x*(1. - x)*(1. - y);
926
927  case 17:
928  return 2.*x*y*(1. - y);
929
930  case 18:
931  return 2.*x*(1. - x)*y;
932
933  case 19:
934  return 2.*(1. - x)*y*(1. - y);
935
936  default:
937  libmesh_error_msg("Invalid i = " << i);
938  }
939
940  default:
941  libmesh_error_msg("Invalid shape function derivative j = " << j);
942  }
943  }
944
945  // triquadratic hexahedral shape functions
946  case HEX27:
947  {
948  libmesh_assert_less (i, 27);
949
950  // Compute hex shape functions as a tensor-product
951  const Real xi = p(0);
952  const Real eta = p(1);
953  const Real zeta = p(2);
954
955  // The only way to make any sense of this
956  // is to look at the mgflo/mg2/mgf documentation
957  // and make the cut-out cube!
958  // 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
959  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
960  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
961  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2};
962
963  switch(j)
964  {
965  case 0:
966  return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)*
967  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)*
968  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta));
969
970  case 1:
971  return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)*
973  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta));
974
975  case 2:
976  return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)*
977  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)*
978  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta));
979
980  default:
981  libmesh_error_msg("Invalid j = " << j);
982  }
983  }
984
985  // quadratic tetrahedral shape functions
986  case TET10:
987  {
988  libmesh_assert_less (i, 10);
989
990  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
991  const Real zeta1 = p(0);
992  const Real zeta2 = p(1);
993  const Real zeta3 = p(2);
994  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
995
996  const Real dzeta0dxi = -1.;
997  const Real dzeta1dxi = 1.;
998  const Real dzeta2dxi = 0.;
999  const Real dzeta3dxi = 0.;
1000
1001  const Real dzeta0deta = -1.;
1002  const Real dzeta1deta = 0.;
1003  const Real dzeta2deta = 1.;
1004  const Real dzeta3deta = 0.;
1005
1006  const Real dzeta0dzeta = -1.;
1007  const Real dzeta1dzeta = 0.;
1008  const Real dzeta2dzeta = 0.;
1009  const Real dzeta3dzeta = 1.;
1010
1011  switch (j)
1012  {
1013  // d()/dxi
1014  case 0:
1015  {
1016  switch(i)
1017  {
1018  case 0:
1019  return (4.*zeta0 - 1.)*dzeta0dxi;
1020
1021  case 1:
1022  return (4.*zeta1 - 1.)*dzeta1dxi;
1023
1024  case 2:
1025  return (4.*zeta2 - 1.)*dzeta2dxi;
1026
1027  case 3:
1028  return (4.*zeta3 - 1.)*dzeta3dxi;
1029
1030  case 4:
1031  return 4.*(zeta0*dzeta1dxi + dzeta0dxi*zeta1);
1032
1033  case 5:
1034  return 4.*(zeta1*dzeta2dxi + dzeta1dxi*zeta2);
1035
1036  case 6:
1037  return 4.*(zeta0*dzeta2dxi + dzeta0dxi*zeta2);
1038
1039  case 7:
1040  return 4.*(zeta0*dzeta3dxi + dzeta0dxi*zeta3);
1041
1042  case 8:
1043  return 4.*(zeta1*dzeta3dxi + dzeta1dxi*zeta3);
1044
1045  case 9:
1046  return 4.*(zeta2*dzeta3dxi + dzeta2dxi*zeta3);
1047
1048  default:
1049  libmesh_error_msg("Invalid i = " << i);
1050  }
1051  }
1052
1053  // d()/deta
1054  case 1:
1055  {
1056  switch(i)
1057  {
1058  case 0:
1059  return (4.*zeta0 - 1.)*dzeta0deta;
1060
1061  case 1:
1062  return (4.*zeta1 - 1.)*dzeta1deta;
1063
1064  case 2:
1065  return (4.*zeta2 - 1.)*dzeta2deta;
1066
1067  case 3:
1068  return (4.*zeta3 - 1.)*dzeta3deta;
1069
1070  case 4:
1071  return 4.*(zeta0*dzeta1deta + dzeta0deta*zeta1);
1072
1073  case 5:
1074  return 4.*(zeta1*dzeta2deta + dzeta1deta*zeta2);
1075
1076  case 6:
1077  return 4.*(zeta0*dzeta2deta + dzeta0deta*zeta2);
1078
1079  case 7:
1080  return 4.*(zeta0*dzeta3deta + dzeta0deta*zeta3);
1081
1082  case 8:
1083  return 4.*(zeta1*dzeta3deta + dzeta1deta*zeta3);
1084
1085  case 9:
1086  return 4.*(zeta2*dzeta3deta + dzeta2deta*zeta3);
1087
1088  default:
1089  libmesh_error_msg("Invalid i = " << i);
1090  }
1091  }
1092
1093  // d()/dzeta
1094  case 2:
1095  {
1096  switch(i)
1097  {
1098  case 0:
1099  return (4.*zeta0 - 1.)*dzeta0dzeta;
1100
1101  case 1:
1102  return (4.*zeta1 - 1.)*dzeta1dzeta;
1103
1104  case 2:
1105  return (4.*zeta2 - 1.)*dzeta2dzeta;
1106
1107  case 3:
1108  return (4.*zeta3 - 1.)*dzeta3dzeta;
1109
1110  case 4:
1111  return 4.*(zeta0*dzeta1dzeta + dzeta0dzeta*zeta1);
1112
1113  case 5:
1114  return 4.*(zeta1*dzeta2dzeta + dzeta1dzeta*zeta2);
1115
1116  case 6:
1117  return 4.*(zeta0*dzeta2dzeta + dzeta0dzeta*zeta2);
1118
1119  case 7:
1120  return 4.*(zeta0*dzeta3dzeta + dzeta0dzeta*zeta3);
1121
1122  case 8:
1123  return 4.*(zeta1*dzeta3dzeta + dzeta1dzeta*zeta3);
1124
1125  case 9:
1126  return 4.*(zeta2*dzeta3dzeta + dzeta2dzeta*zeta3);
1127
1128  default:
1129  libmesh_error_msg("Invalid i = " << i);
1130  }
1131  }
1132
1133  default:
1134  libmesh_error_msg("Invalid j = " << j);
1135  }
1136  }
1137
1138
1139
1140  // quadratic prism shape functions
1141  case PRISM18:
1142  {
1143  libmesh_assert_less (i, 18);
1144
1145  // Compute prism shape functions as a tensor-product
1146  // of a triangle and an edge
1147
1148  Point p2d(p(0),p(1));
1149  Point p1d(p(2));
1150
1151  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
1152  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
1153  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};
1154
1155  switch (j)
1156  {
1157  // d()/dxi
1158  case 0:
1159  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)*
1160  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
1161
1162  // d()/deta
1163  case 1:
1164  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)*
1165  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
1166
1167  // d()/dzeta
1168  case 2:
1169  return (FE<2,L2_LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
1170  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d));
1171
1172  default:
1173  libmesh_error_msg("Invalid shape function derivative j = " << j);
1174  }
1175  }
1176
1177
1178  default:
1179  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
1180  }
1181  }
1182
1183
1184  // unsupported order
1185  default:
1186  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order);
1187  }
1188
1189 #endif
1190 }
1191
1192
1193
1194 template <>
1196  const Order order,
1197  const unsigned int i,
1198  const unsigned int j,
1199  const Point & p)
1200 {
1201  libmesh_assert(elem);
1202
1203  // call the orientation-independent shape function derivatives
1204  return FE<3,L2_LAGRANGE>::shape_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p);
1205 }
1206
1207
1208
1209 template <>
1211  const Order order,
1212  const unsigned int i,
1213  const unsigned int j,
1214  const Point & p)
1215 {
1216 #if LIBMESH_DIM == 3
1217
1218  libmesh_assert_less (j, 6);
1219
1220  switch (order)
1221  {
1222  // linear Lagrange shape functions
1223  case FIRST:
1224  {
1225  return 0.;
1226  }
1227
1228  // quadratic Lagrange shape functions
1229  case SECOND:
1230  {
1231  switch (type)
1232  {
1233
1234  // serendipity hexahedral quadratic shape functions
1235  case HEX20:
1236  {
1237  static bool warning_given_HEX20 = false;
1238
1239  if (!warning_given_HEX20)
1240  libMesh::err << "Second derivatives for 3D Lagrangian HEX20"
1241  << " elements are not yet implemented!"
1242  << std::endl;
1243  warning_given_HEX20 = true;
1244  }
1245  libmesh_fallthrough();
1246
1247  case HEX27:
1248  // triquadratic hexahedral shape functions
1249  {
1250  libmesh_assert_less (i, 27);
1251
1252  // Compute hex shape functions as a tensor-product
1253  const Real xi = p(0);
1254  const Real eta = p(1);
1255  const Real zeta = p(2);
1256
1257  // The only way to make any sense of this
1258  // is to look at the mgflo/mg2/mgf documentation
1259  // and make the cut-out cube!
1260  // 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
1261  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
1262  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
1263  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2};
1264
1265  switch(j)
1266  {
1267  // d^2()/dxi^2
1268  case 0:
1269  return (FE<1,L2_LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i0[i], 0, xi)*
1270  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)*
1271  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta));
1272
1273  // d^2()/dxideta
1274  case 1:
1275  return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)*
1276  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)*
1277  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta));
1278
1279  // d^2()/deta^2
1280  case 2:
1281  return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)*
1283  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta));
1284
1285  // d^2()/dxidzeta
1286  case 3:
1287  return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)*
1288  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)*
1289  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta));
1290
1292  case 4:
1293  return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)*
1294  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)*
1295  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta));
1296
1297  // d^2()/dzeta^2
1298  case 5:
1299  return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)*
1300  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)*
1302
1303  default:
1304  libmesh_error_msg("Invalid j = " << j);
1305  }
1306  }
1307
1308  // quadratic tetrahedral shape functions
1309  case TET10:
1310  {
1311  // The area coordinates are the same as used for the
1312  // shape() and shape_deriv() functions.
1313  // const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
1314  // const Real zeta1 = p(0);
1315  // const Real zeta2 = p(1);
1316  // const Real zeta3 = p(2);
1317  static const Real dzetadxi[4][3] =
1318  {
1319  {-1., -1., -1.},
1320  {1., 0., 0.},
1321  {0., 1., 0.},
1322  {0., 0., 1.}
1323  };
1324
1325  // Convert from j -> (j,k) indices for independent variable
1326  // (0=xi, 1=eta, 2=zeta)
1327  static const unsigned short int independent_var_indices[6][2] =
1328  {
1329  {0, 0}, // d^2 phi / dxi^2
1330  {0, 1}, // d^2 phi / dxi deta
1331  {1, 1}, // d^2 phi / deta^2
1332  {0, 2}, // d^2 phi / dxi dzeta
1333  {1, 2}, // d^2 phi / deta dzeta
1334  {2, 2} // d^2 phi / dzeta^2
1335  };
1336
1337  // Convert from i -> zeta indices. Each quadratic shape
1338  // function for the Tet10 depends on up to two of the zeta
1339  // area coordinate functions (see the shape() function above).
1340  // This table just tells which two area coords it uses.
1341  static const unsigned short int zeta_indices[10][2] =
1342  {
1343  {0, 0},
1344  {1, 1},
1345  {2, 2},
1346  {3, 3},
1347  {0, 1},
1348  {1, 2},
1349  {2, 0},
1350  {0, 3},
1351  {1, 3},
1352  {2, 3},
1353  };
1354
1355  // Look up the independent variable indices for this value of j.
1356  const unsigned int my_j = independent_var_indices[j][0];
1357  const unsigned int my_k = independent_var_indices[j][1];
1358
1359  if (i<4)
1360  {
1362  }
1363
1364  else if (i<10)
1365  {
1366  const unsigned short int my_m = zeta_indices[i][0];
1367  const unsigned short int my_n = zeta_indices[i][1];
1368
1371  }
1372  else
1373  libmesh_error_msg("Invalid shape function index " << i);
1374  }
1375
1376
1377  // quadratic prism shape functions
1378  case PRISM18:
1379  {
1380  libmesh_assert_less (i, 18);
1381
1382  // Compute prism shape functions as a tensor-product
1383  // of a triangle and an edge
1384
1385  Point p2d(p(0),p(1));
1386  Point p1d(p(2));
1387
1388  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
1389  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
1390  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};
1391
1392  switch (j)
1393  {
1394  // d^2()/dxi^2
1395  case 0:
1396  return (FE<2,L2_LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 0, p2d)*
1397  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
1398
1399  // d^2()/dxideta
1400  case 1:
1401  return (FE<2,L2_LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 1, p2d)*
1402  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
1403
1404  // d^2()/deta^2
1405  case 2:
1406  return (FE<2,L2_LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 2, p2d)*
1407  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
1408
1409  // d^2()/dxidzeta
1410  case 3:
1411  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)*
1412  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d));
1413
1415  case 4:
1416  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)*
1417  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d));
1418
1419  // d^2()/dzeta^2
1420  case 5:
1421  return (FE<2,L2_LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
1423
1424  default:
1425  libmesh_error_msg("Invalid shape function derivative j = " << j);
1426  }
1427  }
1428
1429
1430
1431  default:
1432  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
1433  }
1434  }
1435
1436
1437  // unsupported order
1438  default:
1439  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order);
1440  }
1441
1442 #endif
1443 }
1444
1445
1446
1447 template <>
1449  const Order order,
1450  const unsigned int i,
1451  const unsigned int j,
1452  const Point & p)
1453 {
1454  libmesh_assert(elem);
1455
1456  // call the orientation-independent shape function derivatives
1457  return FE<3,L2_LAGRANGE>::shape_second_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p);
1458 }
1459
1460 } // namespace libMesh
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
The base class for all geometric element types.
Definition: elem.h:100
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
unsigned int p_level() const
Definition: elem.h:2555
Template class which generates the different FE families and orders.
Definition: fe.h:89
OStreamProxy err(std::cerr)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual ElemType type() const =0
A geometric point in (x,y,z) 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)