libMesh
fe_szabab_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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 // Local includes
19 #include "libmesh/libmesh_config.h"
20 
21 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
22 
23 // libmesh includes
24 #include "libmesh/fe.h"
25 #include "libmesh/elem.h"
26 #include "libmesh/utility.h"
27 #include "libmesh/enum_to_string.h"
28 
29 // C++ includes
30 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
31 #include <cmath> // for std::sqrt
32 
33 // Anonymous namespace to hold static std::sqrt values
34 namespace
35 {
36 using namespace libMesh;
37 
38 static const Real sqrt2 = std::sqrt(2.);
39 static const Real sqrt6 = std::sqrt(6.);
40 static const Real sqrt10 = std::sqrt(10.);
41 static const Real sqrt14 = std::sqrt(14.);
42 static const Real sqrt22 = std::sqrt(22.);
43 static const Real sqrt26 = std::sqrt(26.);
44 
45 
46 // Take care of edge orientation, this is needed at edge shapes with
47 // (y=0)-asymmetric 1D shapes
48 Real quad_flip(const Elem * elem,
49  const unsigned int totalorder,
50  const unsigned int i)
51 {
52  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
53 
54  // vertex and interior shape functions don't flip
55  if (i < 4 || i >= 4*totalorder)
56  return 1;
57 
58  const int edge = (i-4) / (totalorder-1);
59  libmesh_assert_less(edge, 4);
60  libmesh_assert_greater_equal(edge, 0);
61 
62  const int edge_end = (edge+1)%4;
63 
64  const int edge_i = i - 4 - edge*(totalorder-1);
65 
66  // The "natural" orientation is p1>p0 on edge 0,
67  // p2>p1 on e1, p2>p3 on e2, p3>p0 on e3
68  if (edge_i%2 &&
69  ((elem->point(edge) > elem->point(edge_end)) ==
70  (edge < 2)))
71  return -1;
72 
73  return 1;
74 }
75 
76 } // anonymous namespace
77 
78 
79 namespace libMesh
80 {
81 
82 
84 
85 
86 template <>
87 Real FE<2,SZABAB>::shape(const Elem * elem,
88  const Order order,
89  const unsigned int i,
90  const Point & p,
91  const bool add_p_level)
92 {
93  libmesh_assert(elem);
94 
95  const ElemType type = elem->type();
96 
97  const Order totalorder = static_cast<Order>(order + add_p_level * elem->p_level());
98 
99  // Declare that we are using our own special power function
100  // from the Utility namespace. This saves typing later.
101  using Utility::pow;
102 
103  switch (totalorder)
104  {
105  // 1st & 2nd-order Szabo-Babuska.
106  case FIRST:
107  case SECOND:
108  {
109  switch (type)
110  {
111 
112  // Szabo-Babuska shape functions on the triangle.
113  case TRI3:
114  case TRI6:
115  case TRI7:
116  {
117  const Real l1 = 1-p(0)-p(1);
118  const Real l2 = p(0);
119  const Real l3 = p(1);
120 
121  libmesh_assert_less (i, 6);
122 
123  switch (i)
124  {
125  case 0: return l1;
126  case 1: return l2;
127  case 2: return l3;
128 
129  case 3: return l1*l2*(-4.*sqrt6);
130  case 4: return l2*l3*(-4.*sqrt6);
131  case 5: return l3*l1*(-4.*sqrt6);
132 
133  default:
134  libmesh_error_msg("Invalid i = " << i);
135  }
136  }
137 
138 
139  // Szabo-Babuska shape functions on the quadrilateral.
140  case QUAD4:
141  case QUAD8:
142  case QUAD9:
143  {
144  // Compute quad shape functions as a tensor-product
145  const Real xi = p(0);
146  const Real eta = p(1);
147 
148  libmesh_assert_less (i, 9);
149 
150  // 0 1 2 3 4 5 6 7 8
151  static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
152  static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
153 
154  return (FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
155  FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
156 
157  }
158 
159  default:
160  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
161  }
162  }
163 
164 
165  // 3rd-order Szabo-Babuska.
166  case THIRD:
167  {
168  switch (type)
169  {
170 
171  // Szabo-Babuska shape functions on the triangle.
172  case TRI6:
173  case TRI7:
174  {
175  Real l1 = 1-p(0)-p(1);
176  Real l2 = p(0);
177  Real l3 = p(1);
178 
179  Real f=1;
180 
181  libmesh_assert_less (i, 10);
182 
183 
184  if (i==4 && (elem->point(0) > elem->point(1)))f=-1;
185  if (i==6 && (elem->point(1) > elem->point(2)))f=-1;
186  if (i==8 && (elem->point(2) > elem->point(0)))f=-1;
187 
188 
189  switch (i)
190  {
191  //nodal modes
192  case 0: return l1;
193  case 1: return l2;
194  case 2: return l3;
195 
196  //side modes
197  case 3: return l1*l2*(-4.*sqrt6);
198  case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
199 
200  case 5: return l2*l3*(-4.*sqrt6);
201  case 6: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
202 
203  case 7: return l3*l1*(-4.*sqrt6);
204  case 8: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
205 
206  //internal modes
207  case 9: return l1*l2*l3;
208 
209  default:
210  libmesh_error_msg("Invalid i = " << i);
211  }
212  }
213 
214 
215  // Szabo-Babuska shape functions on the quadrilateral.
216  case QUAD8:
217  case QUAD9:
218  {
219  // Compute quad shape functions as a tensor-product
220  const Real xi = p(0);
221  const Real eta = p(1);
222 
223  libmesh_assert_less (i, 16);
224 
225  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
226  static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3};
227  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3};
228 
229  const Real f = quad_flip(elem, totalorder, i);
230 
231  return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
232  FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
233  }
234 
235  default:
236  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
237  }
238  }
239 
240 
241 
242 
243  // 4th-order Szabo-Babuska.
244  case FOURTH:
245  {
246  switch (type)
247  {
248  // Szabo-Babuska shape functions on the triangle.
249  case TRI6:
250  case TRI7:
251  {
252  Real l1 = 1-p(0)-p(1);
253  Real l2 = p(0);
254  Real l3 = p(1);
255 
256  Real f=1;
257 
258  libmesh_assert_less (i, 15);
259 
260 
261  if (i== 4 && (elem->point(0) > elem->point(1)))f=-1;
262  if (i== 7 && (elem->point(1) > elem->point(2)))f=-1;
263  if (i==10 && (elem->point(2) > elem->point(0)))f=-1;
264 
265 
266  switch (i)
267  {
268  //nodal modes
269  case 0: return l1;
270  case 1: return l2;
271  case 2: return l3;
272 
273  //side modes
274  case 3: return l1*l2*(-4.*sqrt6);
275  case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
276  case 5: return l1*l2*(-sqrt14)*(5.*pow<2>(l2-l1)-1);
277 
278  case 6: return l2*l3*(-4.*sqrt6);
279  case 7: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
280  case 8: return l2*l3*(-sqrt14)*(5.*pow<2>(l3-l2)-1);
281 
282  case 9: return l3*l1*(-4.*sqrt6);
283  case 10: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
284  case 11: return l3*l1*(-sqrt14)*(5.*pow<2>(l1-l3)-1);
285 
286  //internal modes
287  case 12: return l1*l2*l3;
288 
289  case 13: return l1*l2*l3*(l2-l1);
290  case 14: return l1*l2*l3*(2*l3-1);
291 
292  default:
293  libmesh_error_msg("Invalid i = " << i);
294  }
295  }
296 
297 
298  // Szabo-Babuska shape functions on the quadrilateral.
299  case QUAD8:
300  case QUAD9:
301  {
302  // Compute quad shape functions as a tensor-product
303  const Real xi = p(0);
304  const Real eta = p(1);
305 
306  libmesh_assert_less (i, 25);
307 
308  // 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
309  static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4};
310  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4};
311 
312  const Real f = quad_flip(elem, totalorder, i);
313 
314  return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
315  FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
316  }
317 
318  default:
319  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
320  }
321  }
322 
323 
324 
325 
326  // 5th-order Szabo-Babuska.
327  case FIFTH:
328  {
329  switch (type)
330  {
331  // Szabo-Babuska shape functions on the triangle.
332  case TRI6:
333  case TRI7:
334  {
335  Real l1 = 1-p(0)-p(1);
336  Real l2 = p(0);
337  Real l3 = p(1);
338 
339  const Real x=l2-l1;
340  const Real y=2.*l3-1;
341 
342  Real f=1;
343 
344  libmesh_assert_less (i, 21);
345 
346 
347  if ((i== 4||i== 6) && (elem->point(0) > elem->point(1)))f=-1;
348  if ((i== 8||i==10) && (elem->point(1) > elem->point(2)))f=-1;
349  if ((i==12||i==14) && (elem->point(2) > elem->point(0)))f=-1;
350 
351 
352  switch (i)
353  {
354  //nodal modes
355  case 0: return l1;
356  case 1: return l2;
357  case 2: return l3;
358 
359  //side modes
360  case 3: return l1*l2*(-4.*sqrt6);
361  case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
362  case 5: return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
363  case 6: return f*(-sqrt2)*l1*l2*((9.-21.*l1*l1)*l1+(-9.+63.*l1*l1+(-63.*l1+21.*l2)*l2)*l2);
364 
365  case 7: return l2*l3*(-4.*sqrt6);
366  case 8: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
367  case 9: return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
368  case 10: return -f*sqrt2*l2*l3*((-9.0+21.0*l3*l3)*l3+(-63.0*l3*l3+9.0+(63.0*l3-21.0*l2)*l2)*l2);
369 
370  case 11: return l3*l1*(-4.*sqrt6);
371  case 12: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
372  case 13: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
373  case 14: return f*(-sqrt2)*l3*l1*((9.0-21.0*l3*l3)*l3+(-9.0+63.0*l3*l3+(-63.0*l3+21.0*l1)*l1)*l1);
374 
375  //internal modes
376  case 15: return l1*l2*l3;
377 
378  case 16: return l1*l2*l3*x;
379  case 17: return l1*l2*l3*y;
380 
381  case 18: return l1*l2*l3*(1.5*l1*l1-.5+(-3.0*l1+1.5*l2)*l2);
382  case 19: return l1*l2*l3*x*y;
383  case 20: return l1*l2*l3*(1.0+(-6.0+6.0*l3)*l3);
384 
385  default:
386  libmesh_error_msg("Invalid i = " << i);
387  }
388  } // case TRI6/TRI7
389 
390  // Szabo-Babuska shape functions on the quadrilateral.
391  case QUAD8:
392  case QUAD9:
393  {
394  // Compute quad shape functions as a tensor-product
395  const Real xi = p(0);
396  const Real eta = p(1);
397 
398  libmesh_assert_less (i, 36);
399 
400  // 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 28 29 30 31 32 33 34 35
401  static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5};
402  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5};
403 
404  const Real f = quad_flip(elem, totalorder, i);
405 
406  return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
407  FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
408 
409  } // case QUAD8/QUAD9
410 
411  default:
412  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
413 
414  } // switch type
415 
416  } // case FIFTH
417 
418  // 6th-order Szabo-Babuska.
419  case SIXTH:
420  {
421  switch (type)
422  {
423  // Szabo-Babuska shape functions on the triangle.
424  case TRI6:
425  case TRI7:
426  {
427  Real l1 = 1-p(0)-p(1);
428  Real l2 = p(0);
429  Real l3 = p(1);
430 
431  const Real x=l2-l1;
432  const Real y=2.*l3-1;
433 
434  Real f=1;
435 
436  libmesh_assert_less (i, 28);
437 
438 
439  if ((i== 4||i== 6) && (elem->point(0) > elem->point(1)))f=-1;
440  if ((i== 9||i==11) && (elem->point(1) > elem->point(2)))f=-1;
441  if ((i==14||i==16) && (elem->point(2) > elem->point(0)))f=-1;
442 
443 
444  switch (i)
445  {
446  //nodal modes
447  case 0: return l1;
448  case 1: return l2;
449  case 2: return l3;
450 
451  //side modes
452  case 3: return l1*l2*(-4.*sqrt6);
453  case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
454  case 5: return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
455  case 6: return f*(-sqrt2)*l1*l2*((9.0-21.0*l1*l1)*l1+(-9.0+63.0*l1*l1+(-63.0*l1+21.0*l2)*l2)*l2);
456  case 7: return -sqrt22*l1*l2*(0.5+(-7.0+0.105E2*l1*l1)*l1*l1+((14.0-0.42E2*l1*l1)*l1+(-7.0+0.63E2*l1*l1+(-0.42E2*l1+0.105E2*l2)*l2)*l2)*l2);
457 
458  case 8: return l2*l3*(-4.*sqrt6);
459  case 9: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
460  case 10: return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
461  case 11: return f*(-sqrt2)*l2*l3*((-9.0+21.0*l3*l3)*l3+(-63.0*l3*l3+9.0+(63.0*l3-21.0*l2)*l2)*l2);
462  case 12: return -sqrt22*l2*l3*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l2)*l2)*l2)*l2);
463 
464  case 13: return l3*l1*(-4.*sqrt6);
465  case 14: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
466  case 15: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
467  case 16: return f*(-sqrt2)*l3*l1*((9.0-21.0*l3*l3)*l3+(-9.0+63.0*l3*l3+(-63.0*l3+21.0*l1)*l1)*l1);
468  case 17: return -sqrt22*l3*l1*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l1)*l1)*l1)*l1);
469 
470 
471 
472  //internal modes
473  case 18: return l1*l2*l3;
474 
475  case 19: return l1*l2*l3*x;
476  case 20: return l1*l2*l3*y;
477 
478  case 21: return 0.5*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2);
479  case 22: return l1*l2*l3*(l2-l1)*(2.0*l3-1.0);
480  case 23: return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3);
481  case 24: return 0.5*l1*l2*l3*((3.0-5.0*l1*l1)*l1+(-3.0+15.0*l1*l1+(-15.0*l1+5.0*l2)*l2)*l2);
482  case 25: return 0.5*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2)*(2.0*l3-1.0);
483  case 26: return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3)*(l2-l1);
484  case 27: return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3);
485 
486 
487  default:
488  libmesh_error_msg("Invalid i = " << i);
489  }
490  } // case TRI6/TRI7
491 
492  // Szabo-Babuska shape functions on the quadrilateral.
493  case QUAD8:
494  case QUAD9:
495  {
496  // Compute quad shape functions as a tensor-product
497  const Real xi = p(0);
498  const Real eta = p(1);
499 
500  libmesh_assert_less (i, 49);
501 
502  // 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 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
503  static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6};
504  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6};
505 
506  const Real f = quad_flip(elem, totalorder, i);
507 
508  return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
509  FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
510 
511  } // case QUAD8/QUAD9
512 
513  default:
514  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
515 
516  } // switch type
517 
518  } // case SIXTH
519 
520 
521  // 7th-order Szabo-Babuska.
522  case SEVENTH:
523  {
524  switch (type)
525  {
526  // Szabo-Babuska shape functions on the triangle.
527  case TRI6:
528  case TRI7:
529  {
530 
531  Real l1 = 1-p(0)-p(1);
532  Real l2 = p(0);
533  Real l3 = p(1);
534 
535  const Real x=l2-l1;
536  const Real y=2.*l3-1.;
537 
538  Real f=1;
539 
540  libmesh_assert_less (i, 36);
541 
542 
543  if ((i>= 4&&i<= 8) && (elem->point(0) > elem->point(1)))f=-1;
544  if ((i>=10&&i<=14) && (elem->point(1) > elem->point(2)))f=-1;
545  if ((i>=16&&i<=20) && (elem->point(2) > elem->point(0)))f=-1;
546 
547 
548  switch (i)
549  {
550  //nodal modes
551  case 0: return l1;
552  case 1: return l2;
553  case 2: return l3;
554 
555  //side modes
556  case 3: return l1*l2*(-4.*sqrt6);
557  case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
558 
559  case 5: return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
560  case 6: return f*-sqrt2*l1*l2*((9.0-21.0*l1*l1)*l1+(-9.0+63.0*l1*l1+(-63.0*l1+21.0*l2)*l2)*l2);
561  case 7: return -sqrt22*l1*l2*(0.5+(-7.0+0.105E2*l1*l1)*l1*l1+((14.0-0.42E2*l1*l1)*l1+(-7.0+0.63E2*l1*l1+(-0.42E2*l1+0.105E2*l2)*l2)*l2)*l2);
562  case 8: return f*-sqrt26*l1*l2*((-0.25E1+(15.0-0.165E2*l1*l1)*l1*l1)*l1+(0.25E1+(-45.0+0.825E2*l1*l1)*l1*l1+((45.0-0.165E3*l1*l1)*l1+(-15.0+0.165E3*l1*l1+(-0.825E2*l1+0.165E2*l2)*l2)*l2)*l2)*l2);
563 
564  case 9: return l2*l3*(-4.*sqrt6);
565  case 10: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
566 
567  case 11: return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
568  case 12: return f*-sqrt2*l2*l3*((-9.0+21.0*l3*l3)*l3+(-63.0*l3*l3+9.0+(63.0*l3-21.0*l2)*l2)*l2);
569  case 13: return -sqrt22*l2*l3*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l2)*l2)*l2)*l2);
570  case 14: return f*-sqrt26*l2*l3*((0.25E1+(-15.0+0.165E2*l3*l3)*l3*l3)*l3+(-0.25E1+(45.0-0.825E2*l3*l3)*l3*l3+((-45.0+0.165E3*l3*l3)*l3+(15.0-0.165E3*l3*l3+(0.825E2*l3-0.165E2*l2)*l2)*l2)*l2)*l2);
571 
572  case 15: return l3*l1*(-4.*sqrt6);
573  case 16: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
574 
575  case 17: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
576  case 18: return -f*sqrt2*l3*l1*((9.-21.*l3*l3)*l3+(-9.+63.*l3*l3+(-63.*l3+21.*l1)*l1)*l1);
577  case 19: return -sqrt22*l3*l1*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l1)*l1)*l1)*l1);
578  case 20: return f*-sqrt26*l3*l1*((-0.25E1+(15.0-0.165E2*l3*l3)*l3*l3)*l3+(0.25E1+(-45.0+0.825E2*l3*l3)*l3*l3+((45.0-0.165E3*l3*l3)*l3+(-15.0+0.165E3*l3*l3+(-0.825E2*l3+0.165E2*l1)*l1)*l1)*l1)*l1);
579 
580 
581  //internal modes
582  case 21: return l1*l2*l3;
583 
584  case 22: return l1*l2*l3*x;
585  case 23: return l1*l2*l3*y;
586 
587  case 24: return l1*l2*l3*0.5*(3.*pow<2>(x)-1.);
588  case 25: return l1*l2*l3*x*y;
589  case 26: return l1*l2*l3*0.5*(3.*pow<2>(y)-1.);
590 
591  case 27: return 0.5*l1*l2*l3*((3.0-5.0*l1*l1)*l1+(-3.0+15.0*l1*l1+(-15.0*l1+5.0*l2)*l2)*l2);
592  case 28: return 0.5*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2)*(2.0*l3-1.0);
593  case 29: return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3)*(l2-l1);
594  case 30: return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3);
595  case 31: return 0.125*l1*l2*l3*((-15.0+(70.0-63.0*l1*l1)*l1*l1)*l1+(15.0+(-210.0+315.0*l1*l1)*l1*l1+((210.0-630.0*l1*l1)*l1+(-70.0+630.0*l1*l1+(-315.0*l1+63.0*l2)*l2)*l2)*l2)*l2);
596  case 32: return 0.5*l1*l2*l3*((3.0-5.0*l1*l1)*l1+(-3.0+15.0*l1*l1+(-15.0*l1+5.0*l2)*l2)*l2)*(2.0*l3-1.0);
597  case 33: return 0.25*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2)*(2.0+(-12.0+12.0*l3)*l3);
598  case 34: return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3)*(l2-l1);
599  case 35: return 0.125*l1*l2*l3*(-8.0+(240.0+(-1680.0+(4480.0+(-5040.0+2016.0*l3)*l3)*l3)*l3)*l3);
600 
601  default:
602  libmesh_error_msg("Invalid i = " << i);
603  }
604  } // case TRI6/TRI7
605 
606  // Szabo-Babuska shape functions on the quadrilateral.
607  case QUAD8:
608  case QUAD9:
609  {
610  // Compute quad shape functions as a tensor-product
611  const Real xi = p(0);
612  const Real eta = p(1);
613 
614  libmesh_assert_less (i, 64);
615 
616  // 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 28 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 61 62 63
617  static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7};
618  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7};
619 
620  const Real f = quad_flip(elem, totalorder, i);
621 
622  return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
623  FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
624 
625  } // case QUAD8/QUAD9
626 
627  default:
628  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
629 
630  } // switch type
631 
632  } // case SEVENTH
633 
634 
635  // by default throw an error
636  default:
637  libmesh_error_msg("ERROR: Unsupported polynomial order!");
638  } // switch order
639 }
640 
641 
642 
643 
644 
645 template <>
647  const Order,
648  const unsigned int,
649  const Point &)
650 {
651  libmesh_error_msg("Szabo-Babuska polynomials require the element type \nbecause edge orientation is needed.");
652  return 0.;
653 }
654 
655 
656 
657 template <>
659  const Elem * elem,
660  const unsigned int i,
661  const Point & p,
662  const bool add_p_level)
663 {
664  return FE<2,SZABAB>::shape(elem, fet.order, i, p, add_p_level);
665 }
666 
667 
668 
669 
670 
671 template <>
673  const Order order,
674  const unsigned int i,
675  const unsigned int j,
676  const Point & p,
677  const bool add_p_level)
678 {
679  libmesh_assert(elem);
680 
681  const ElemType type = elem->type();
682 
683  const Order totalorder = static_cast<Order>(order + add_p_level * elem->p_level());
684 
685  switch (totalorder)
686  {
687 
688  // 1st & 2nd-order Szabo-Babuska.
689  case FIRST:
690  case SECOND:
691  {
692  switch (type)
693  {
694 
695  // Szabo-Babuska shape functions on the triangle.
696  case TRI3:
697  case TRI6:
698  case TRI7:
699  {
700  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
701  }
702 
703 
704  // Szabo-Babuska shape functions on the quadrilateral.
705  case QUAD4:
706  case QUAD8:
707  case QUAD9:
708  {
709  // Compute quad shape functions as a tensor-product
710  const Real xi = p(0);
711  const Real eta = p(1);
712 
713  libmesh_assert_less (i, 9);
714 
715  // 0 1 2 3 4 5 6 7 8
716  static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
717  static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
718 
719  switch (j)
720  {
721  // d()/dxi
722  case 0:
723  return (FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
724  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
725 
726  // d()/deta
727  case 1:
728  return (FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
729  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
730 
731  default:
732  libmesh_error_msg("Invalid j = " << j);
733  }
734  }
735 
736  default:
737  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
738  }
739  }
740 
741 
742 
743  // 3rd-order Szabo-Babuska.
744  case THIRD:
745  {
746  switch (type)
747  {
748  // Szabo-Babuska shape functions on the triangle.
749  case TRI6:
750  case TRI7:
751  {
752  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
753  }
754 
755 
756  // Szabo-Babuska shape functions on the quadrilateral.
757  case QUAD8:
758  case QUAD9:
759  {
760  // Compute quad shape functions as a tensor-product
761  const Real xi = p(0);
762  const Real eta = p(1);
763 
764  libmesh_assert_less (i, 16);
765 
766  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
767  static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3};
768  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3};
769 
770  const Real f = quad_flip(elem, totalorder, i);
771 
772  switch (j)
773  {
774  // d()/dxi
775  case 0:
776  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
777  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
778 
779  // d()/deta
780  case 1:
781  return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
782  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
783 
784  default:
785  libmesh_error_msg("Invalid j = " << j);
786  }
787  }
788 
789  default:
790  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
791  }
792  }
793 
794 
795 
796 
797  // 4th-order Szabo-Babuska.
798  case FOURTH:
799  {
800  switch (type)
801  {
802 
803  // Szabo-Babuska shape functions on the triangle.
804  case TRI6:
805  case TRI7:
806  {
807  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
808  }
809 
810 
811  // Szabo-Babuska shape functions on the quadrilateral.
812  case QUAD8:
813  case QUAD9:
814  {
815  // Compute quad shape functions as a tensor-product
816  const Real xi = p(0);
817  const Real eta = p(1);
818 
819  libmesh_assert_less (i, 25);
820 
821  // 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
822  static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4};
823  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4};
824 
825  const Real f = quad_flip(elem, totalorder, i);
826 
827  switch (j)
828  {
829  // d()/dxi
830  case 0:
831  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
832  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
833 
834  // d()/deta
835  case 1:
836  return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
837  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
838 
839  default:
840  libmesh_error_msg("Invalid j = " << j);
841  }
842  }
843 
844  default:
845  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
846  }
847  }
848 
849 
850 
851 
852  // 5th-order Szabo-Babuska.
853  case FIFTH:
854  {
855  // Szabo-Babuska shape functions on the quadrilateral.
856  switch (type)
857  {
858 
859  // Szabo-Babuska shape functions on the triangle.
860  case TRI6:
861  case TRI7:
862  {
863  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
864  }
865 
866 
867  case QUAD8:
868  case QUAD9:
869  {
870  // Compute quad shape functions as a tensor-product
871  const Real xi = p(0);
872  const Real eta = p(1);
873 
874  libmesh_assert_less (i, 36);
875 
876  // 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 28 29 30 31 32 33 34 35
877  static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5};
878  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5};
879 
880  const Real f = quad_flip(elem, totalorder, i);
881 
882  switch (j)
883  {
884  // d()/dxi
885  case 0:
886  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
887  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
888 
889  // d()/deta
890  case 1:
891  return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
892  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
893 
894  default:
895  libmesh_error_msg("Invalid j = " << j);
896  }
897  }
898 
899  default:
900  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
901  }
902  }
903 
904 
905  // 6th-order Szabo-Babuska.
906  case SIXTH:
907  {
908  // Szabo-Babuska shape functions on the quadrilateral.
909  switch (type)
910  {
911 
912  // Szabo-Babuska shape functions on the triangle.
913  case TRI6:
914  case TRI7:
915  {
916  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
917  }
918 
919 
920  case QUAD8:
921  case QUAD9:
922  {
923  // Compute quad shape functions as a tensor-product
924  const Real xi = p(0);
925  const Real eta = p(1);
926 
927  libmesh_assert_less (i, 49);
928 
929  // 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 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
930  static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6};
931  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6};
932 
933  const Real f = quad_flip(elem, totalorder, i);
934 
935  switch (j)
936  {
937  // d()/dxi
938  case 0:
939  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
940  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
941 
942  // d()/deta
943  case 1:
944  return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
945  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
946 
947  default:
948  libmesh_error_msg("Invalid j = " << j);
949  }
950  }
951 
952  default:
953  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
954  }
955  }
956 
957 
958  // 7th-order Szabo-Babuska.
959  case SEVENTH:
960  {
961  // Szabo-Babuska shape functions on the quadrilateral.
962  switch (type)
963  {
964 
965  // Szabo-Babuska shape functions on the triangle.
966  case TRI6:
967  case TRI7:
968  {
969  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
970  }
971 
972 
973  case QUAD8:
974  case QUAD9:
975  {
976  // Compute quad shape functions as a tensor-product
977  const Real xi = p(0);
978  const Real eta = p(1);
979 
980  libmesh_assert_less (i, 64);
981 
982  // 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 28 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 61 62 63
983  static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7};
984  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7};
985 
986  const Real f = quad_flip(elem, totalorder, i);
987 
988  switch (j)
989  {
990  // d()/dxi
991  case 0:
992  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
993  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
994 
995  // d()/deta
996  case 1:
997  return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
998  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
999 
1000  default:
1001  libmesh_error_msg("Invalid j = " << j);
1002  }
1003  }
1004 
1005  default:
1006  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1007  }
1008  }
1009 
1010 
1011 
1012  // by default throw an error;call the orientation-independent shape functions
1013  default:
1014  libmesh_error_msg("ERROR: Unsupported polynomial order!");
1015  }
1016 }
1017 
1018 
1019 
1020 
1021 
1022 template <>
1024  const Order,
1025  const unsigned int,
1026  const unsigned int,
1027  const Point &)
1028 {
1029  libmesh_error_msg("Szabo-Babuska polynomials require the element type \nbecause edge orientation is needed.");
1030  return 0.;
1031 }
1032 
1033 
1034 template <>
1036  const Elem * elem,
1037  const unsigned int i,
1038  const unsigned int j,
1039  const Point & p,
1040  const bool add_p_level)
1041 {
1042  return FE<2,SZABAB>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
1043 }
1044 
1045 
1046 
1047 
1048 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1049 
1050 template <>
1052  const Order,
1053  const unsigned int,
1054  const unsigned int,
1055  const Point &)
1056 {
1057  static bool warning_given = false;
1058 
1059  if (!warning_given)
1060  libMesh::err << "Second derivatives for Szabab elements "
1061  << " are not yet implemented!"
1062  << std::endl;
1063 
1064  warning_given = true;
1065  return 0.;
1066 }
1067 
1068 
1069 
1070 template <>
1072  const Order,
1073  const unsigned int,
1074  const unsigned int,
1075  const Point &,
1076  const bool)
1077 {
1078  static bool warning_given = false;
1079 
1080  if (!warning_given)
1081  libMesh::err << "Second derivatives for Szabab elements "
1082  << " are not yet implemented!"
1083  << std::endl;
1084 
1085  warning_given = true;
1086  return 0.;
1087 }
1088 
1089 
1090 template <>
1092  const Elem *,
1093  const unsigned int,
1094  const unsigned int,
1095  const Point &,
1096  const bool)
1097 {
1098  static bool warning_given = false;
1099 
1100  if (!warning_given)
1101  libMesh::err << "Second derivatives for Szabab elements "
1102  << " are not yet implemented!"
1103  << std::endl;
1104 
1105  warning_given = true;
1106  return 0.;
1107 }
1108 
1109 #endif
1110 
1111 } // namespace libMesh
1112 
1113 #endif// LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
OStreamProxy err
ElemType
Defines an enum for geometric element types.
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
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:2945
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:201
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:53
The libMesh namespace provides an interface to certain functionality in the library.
OutputShape fe_fdm_deriv(const Elem *elem, const Order order, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level, OutputShape(*shape_func)(const Elem *, const Order, const unsigned int, const Point &, const bool))
Helper functions for finite differenced derivatives in cases where analytical calculations haven&#39;t be...
Definition: fe.C:744
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
A specific instantiation of the FEBase class.
Definition: fe.h:127
T pow(const T &x)
Definition: utility.h:328
libmesh_assert(ctx)
std::string enum_to_string(const T e)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
const Point & point(const unsigned int i) const
Definition: elem.h:2277
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)