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