77 shape(0) = 1. - ip.
x - ip.
y;
85 dshape(0,0) = -1.; dshape(0,1) = -1.;
86 dshape(1,0) = 1.; dshape(1,1) = 0.;
87 dshape(2,0) = 0.; dshape(2,1) = 1.;
107 shape(0) = (1. - ip.
x) * (1. - ip.
y) ;
108 shape(1) = ip.
x * (1. - ip.
y) ;
109 shape(2) = ip.
x * ip.
y ;
110 shape(3) = (1. - ip.
x) * ip.
y ;
116 dshape(0,0) = -1. + ip.
y; dshape(0,1) = -1. + ip.
x ;
117 dshape(1,0) = 1. - ip.
y; dshape(1,1) = -ip.
x ;
118 dshape(2,0) = ip.
y ; dshape(2,1) = ip.
x ;
119 dshape(3,0) = -ip.
y ; dshape(3,1) = 1. - ip.
x ;
125 h(0,0) = 0.; h(0,1) = 1.; h(0,2) = 0.;
126 h(1,0) = 0.; h(1,1) = -1.; h(1,2) = 0.;
127 h(2,0) = 0.; h(2,1) = 1.; h(2,2) = 0.;
128 h(3,0) = 0.; h(3,1) = -1.; h(3,2) = 0.;
148 shape(0) = 5./3. - 2. * (x + y);
149 shape(1) = 2. * (x - 1./6.);
150 shape(2) = 2. * (y - 1./6.);
156 dshape(0,0) = -2.; dshape(0,1) = -2.;
157 dshape(1,0) = 2.; dshape(1,1) = 0.;
158 dshape(2,0) = 0.; dshape(2,1) = 2.;
163 dofs(vertex) = 2./3.;
164 dofs((vertex+1)%3) = 1./6.;
165 dofs((vertex+2)%3) = 1./6.;
170const real_t GaussBiLinear2DFiniteElement::p[] =
171{ 0.2113248654051871177454256, 0.7886751345948128822545744 };
191 shape(0) = 3. * (p[1] - x) * (p[1] - y);
192 shape(1) = 3. * (x - p[0]) * (p[1] - y);
193 shape(2) = 3. * (x - p[0]) * (y - p[0]);
194 shape(3) = 3. * (p[1] - x) * (y - p[0]);
202 dshape(0,0) = 3. * (y - p[1]); dshape(0,1) = 3. * (x - p[1]);
203 dshape(1,0) = 3. * (p[1] - y); dshape(1,1) = 3. * (p[0] - x);
204 dshape(2,0) = 3. * (y - p[0]); dshape(2,1) = 3. * (x - p[0]);
205 dshape(3,0) = 3. * (p[0] - y); dshape(3,1) = 3. * (p[1] - x);
211 dofs(vertex) = p[1]*p[1];
212 dofs((vertex+1)%4) =
p[0]*
p[1];
213 dofs((vertex+2)%4) =
p[0]*
p[0];
214 dofs((vertex+3)%4) =
p[0]*
p[1];
235 shape(0) = 1. - ip.
x - ip.
y;
243 dshape(0,0) = -1.; dshape(0,1) = -1.;
244 dshape(1,0) = 1.; dshape(1,1) = 0.;
245 dshape(2,0) = 0.; dshape(2,1) = 1.;
261 real_t l1 = 1.0 - x, l2 = x, l3 = 2. * x - 1.;
263 shape(0) = l1 * (-l3);
265 shape(2) = 4. * l1 * l2;
273 dshape(0,0) = 4. * x - 3.;
274 dshape(1,0) = 4. * x - 1.;
275 dshape(2,0) = 4. - 8. * x;
300 real_t l1 = 1.-x-y, l2 = x, l3 = y;
302 shape(0) = l1 * (2. * l1 - 1.);
303 shape(1) = l2 * (2. * l2 - 1.);
304 shape(2) = l3 * (2. * l3 - 1.);
305 shape(3) = 4. * l1 * l2;
306 shape(4) = 4. * l2 * l3;
307 shape(5) = 4. * l3 * l1;
316 dshape(0,1) = 4. * (x + y) - 3.;
318 dshape(1,0) = 4. * x - 1.;
322 dshape(2,1) = 4. * y - 1.;
324 dshape(3,0) = -4. * (2. * x + y - 1.);
325 dshape(3,1) = -4. * x;
327 dshape(4,0) = 4. * y;
328 dshape(4,1) = 4. * x;
330 dshape(5,0) = -4. * y;
331 dshape(5,1) = -4. * (x + 2. * y - 1.);
371 case 0: dofs(3) = 0.25; dofs(5) = 0.25;
break;
372 case 1: dofs(3) = 0.25; dofs(4) = 0.25;
break;
373 case 2: dofs(4) = 0.25; dofs(5) = 0.25;
break;
379const real_t GaussQuad2DFiniteElement::p[] =
380{ 0.0915762135097707434595714634022015, 0.445948490915964886318329253883051 };
398 for (
int i = 0; i < 6; i++)
430 D(0,0) = 0.; D(0,1) = 0.;
431 D(1,0) = 1.; D(1,1) = 0.;
432 D(2,0) = 0.; D(2,1) = 1.;
433 D(3,0) = 2. * x; D(3,1) = 0.;
434 D(4,0) = y; D(4,1) = x;
435 D(5,0) = 0.; D(5,1) = 2. * y;
468 real_t l1x, l2x, l3x, l1y, l2y, l3y;
470 l1x = (x - 1.) * (2. * x - 1);
471 l2x = 4. * x * (1. - x);
472 l3x = x * (2. * x - 1.);
473 l1y = (y - 1.) * (2. * y - 1);
474 l2y = 4. * y * (1. - y);
475 l3y = y * (2. * y - 1.);
477 shape(0) = l1x * l1y;
478 shape(4) = l2x * l1y;
479 shape(1) = l3x * l1y;
480 shape(7) = l1x * l2y;
481 shape(8) = l2x * l2y;
482 shape(5) = l3x * l2y;
483 shape(3) = l1x * l3y;
484 shape(6) = l2x * l3y;
485 shape(2) = l3x * l3y;
492 real_t l1x, l2x, l3x, l1y, l2y, l3y;
493 real_t d1x, d2x, d3x, d1y, d2y, d3y;
495 l1x = (x - 1.) * (2. * x - 1);
496 l2x = 4. * x * (1. - x);
497 l3x = x * (2. * x - 1.);
498 l1y = (y - 1.) * (2. * y - 1);
499 l2y = 4. * y * (1. - y);
500 l3y = y * (2. * y - 1.);
509 dshape(0,0) = d1x * l1y;
510 dshape(0,1) = l1x * d1y;
512 dshape(4,0) = d2x * l1y;
513 dshape(4,1) = l2x * d1y;
515 dshape(1,0) = d3x * l1y;
516 dshape(1,1) = l3x * d1y;
518 dshape(7,0) = d1x * l2y;
519 dshape(7,1) = l1x * d2y;
521 dshape(8,0) = d2x * l2y;
522 dshape(8,1) = l2x * d2y;
524 dshape(5,0) = d3x * l2y;
525 dshape(5,1) = l3x * d2y;
527 dshape(3,0) = d1x * l3y;
528 dshape(3,1) = l1x * d3y;
530 dshape(6,0) = d2x * l3y;
531 dshape(6,1) = l2x * d3y;
533 dshape(2,0) = d3x * l3y;
534 dshape(2,1) = l3x * d3y;
546 case 0: dofs(4) = 0.25; dofs(7) = 0.25;
break;
547 case 1: dofs(4) = 0.25; dofs(5) = 0.25;
break;
548 case 2: dofs(5) = 0.25; dofs(6) = 0.25;
break;
549 case 3: dofs(6) = 0.25; dofs(7) = 0.25;
break;
559 const real_t p1 = 0.5*(1.-sqrt(3./5.));
585 const real_t p1 = 0.5*(1.-sqrt(3./5.));
588 real_t l1x, l2x, l3x, l1y, l2y, l3y;
590 l1x = (x - 1.) * (2. * x - 1);
591 l2x = 4. * x * (1. - x);
592 l3x = x * (2. * x - 1.);
593 l1y = (y - 1.) * (2. * y - 1);
594 l2y = 4. * y * (1. - y);
595 l3y = y * (2. * y - 1.);
597 shape(0) = l1x * l1y;
598 shape(4) = l2x * l1y;
599 shape(1) = l3x * l1y;
600 shape(7) = l1x * l2y;
601 shape(8) = l2x * l2y;
602 shape(5) = l3x * l2y;
603 shape(3) = l1x * l3y;
604 shape(6) = l2x * l3y;
605 shape(2) = l3x * l3y;
612 const real_t p1 = 0.5*(1.-sqrt(3./5.));
615 real_t l1x, l2x, l3x, l1y, l2y, l3y;
616 real_t d1x, d2x, d3x, d1y, d2y, d3y;
618 l1x = (x - 1.) * (2. * x - 1);
619 l2x = 4. * x * (1. - x);
620 l3x = x * (2. * x - 1.);
621 l1y = (y - 1.) * (2. * y - 1);
622 l2y = 4. * y * (1. - y);
623 l3y = y * (2. * y - 1.);
625 d1x =
a * (4. * x - 3.);
626 d2x =
a * (4. - 8. * x);
627 d3x =
a * (4. * x - 1.);
628 d1y =
a * (4. * y - 3.);
629 d2y =
a * (4. - 8. * y);
630 d3y =
a * (4. * y - 1.);
632 dshape(0,0) = d1x * l1y;
633 dshape(0,1) = l1x * d1y;
635 dshape(4,0) = d2x * l1y;
636 dshape(4,1) = l2x * d1y;
638 dshape(1,0) = d3x * l1y;
639 dshape(1,1) = l3x * d1y;
641 dshape(7,0) = d1x * l2y;
642 dshape(7,1) = l1x * d2y;
644 dshape(8,0) = d2x * l2y;
645 dshape(8,1) = l2x * d2y;
647 dshape(5,0) = d3x * l2y;
648 dshape(5,1) = l3x * d2y;
650 dshape(3,0) = d1x * l3y;
651 dshape(3,1) = l1x * d3y;
653 dshape(6,0) = d2x * l3y;
654 dshape(6,1) = l2x * d3y;
656 dshape(2,0) = d3x * l3y;
657 dshape(2,1) = l3x * d3y;
702 real_t w1x, w2x, w3x, w1y, w2y, w3y;
703 real_t l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
705 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
706 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
708 l0x = (- 4.5) * w1x * w2x * w3x;
709 l1x = ( 13.5) * x * w2x * w3x;
710 l2x = (-13.5) * x * w1x * w3x;
711 l3x = ( 4.5) * x * w1x * w2x;
713 l0y = (- 4.5) * w1y * w2y * w3y;
714 l1y = ( 13.5) * y * w2y * w3y;
715 l2y = (-13.5) * y * w1y * w3y;
716 l3y = ( 4.5) * y * w1y * w2y;
718 shape(0) = l0x * l0y;
719 shape(1) = l3x * l0y;
720 shape(2) = l3x * l3y;
721 shape(3) = l0x * l3y;
722 shape(4) = l1x * l0y;
723 shape(5) = l2x * l0y;
724 shape(6) = l3x * l1y;
725 shape(7) = l3x * l2y;
726 shape(8) = l2x * l3y;
727 shape(9) = l1x * l3y;
728 shape(10) = l0x * l2y;
729 shape(11) = l0x * l1y;
730 shape(12) = l1x * l1y;
731 shape(13) = l2x * l1y;
732 shape(14) = l1x * l2y;
733 shape(15) = l2x * l2y;
741 real_t w1x, w2x, w3x, w1y, w2y, w3y;
742 real_t l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
743 real_t d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
745 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
746 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
748 l0x = (- 4.5) * w1x * w2x * w3x;
749 l1x = ( 13.5) * x * w2x * w3x;
750 l2x = (-13.5) * x * w1x * w3x;
751 l3x = ( 4.5) * x * w1x * w2x;
753 l0y = (- 4.5) * w1y * w2y * w3y;
754 l1y = ( 13.5) * y * w2y * w3y;
755 l2y = (-13.5) * y * w1y * w3y;
756 l3y = ( 4.5) * y * w1y * w2y;
758 d0x = -5.5 + ( 18. - 13.5 * x) * x;
759 d1x = 9. + (-45. + 40.5 * x) * x;
760 d2x = -4.5 + ( 36. - 40.5 * x) * x;
761 d3x = 1. + (- 9. + 13.5 * x) * x;
763 d0y = -5.5 + ( 18. - 13.5 * y) * y;
764 d1y = 9. + (-45. + 40.5 * y) * y;
765 d2y = -4.5 + ( 36. - 40.5 * y) * y;
766 d3y = 1. + (- 9. + 13.5 * y) * y;
768 dshape( 0,0) = d0x * l0y; dshape( 0,1) = l0x * d0y;
769 dshape( 1,0) = d3x * l0y; dshape( 1,1) = l3x * d0y;
770 dshape( 2,0) = d3x * l3y; dshape( 2,1) = l3x * d3y;
771 dshape( 3,0) = d0x * l3y; dshape( 3,1) = l0x * d3y;
772 dshape( 4,0) = d1x * l0y; dshape( 4,1) = l1x * d0y;
773 dshape( 5,0) = d2x * l0y; dshape( 5,1) = l2x * d0y;
774 dshape( 6,0) = d3x * l1y; dshape( 6,1) = l3x * d1y;
775 dshape( 7,0) = d3x * l2y; dshape( 7,1) = l3x * d2y;
776 dshape( 8,0) = d2x * l3y; dshape( 8,1) = l2x * d3y;
777 dshape( 9,0) = d1x * l3y; dshape( 9,1) = l1x * d3y;
778 dshape(10,0) = d0x * l2y; dshape(10,1) = l0x * d2y;
779 dshape(11,0) = d0x * l1y; dshape(11,1) = l0x * d1y;
780 dshape(12,0) = d1x * l1y; dshape(12,1) = l1x * d1y;
781 dshape(13,0) = d2x * l1y; dshape(13,1) = l2x * d1y;
782 dshape(14,0) = d1x * l2y; dshape(14,1) = l1x * d2y;
783 dshape(15,0) = d2x * l2y; dshape(15,1) = l2x * d2y;
791 real_t w1x, w2x, w3x, w1y, w2y, w3y;
792 real_t l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
793 real_t d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
794 real_t h0x, h1x, h2x, h3x, h0y, h1y, h2y, h3y;
796 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
797 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
799 l0x = (- 4.5) * w1x * w2x * w3x;
800 l1x = ( 13.5) * x * w2x * w3x;
801 l2x = (-13.5) * x * w1x * w3x;
802 l3x = ( 4.5) * x * w1x * w2x;
804 l0y = (- 4.5) * w1y * w2y * w3y;
805 l1y = ( 13.5) * y * w2y * w3y;
806 l2y = (-13.5) * y * w1y * w3y;
807 l3y = ( 4.5) * y * w1y * w2y;
809 d0x = -5.5 + ( 18. - 13.5 * x) * x;
810 d1x = 9. + (-45. + 40.5 * x) * x;
811 d2x = -4.5 + ( 36. - 40.5 * x) * x;
812 d3x = 1. + (- 9. + 13.5 * x) * x;
814 d0y = -5.5 + ( 18. - 13.5 * y) * y;
815 d1y = 9. + (-45. + 40.5 * y) * y;
816 d2y = -4.5 + ( 36. - 40.5 * y) * y;
817 d3y = 1. + (- 9. + 13.5 * y) * y;
819 h0x = -27. * x + 18.;
821 h2x = -81. * x + 36.;
824 h0y = -27. * y + 18.;
826 h2y = -81. * y + 36.;
829 h( 0,0) = h0x * l0y; h( 0,1) = d0x * d0y; h( 0,2) = l0x * h0y;
830 h( 1,0) = h3x * l0y; h( 1,1) = d3x * d0y; h( 1,2) = l3x * h0y;
831 h( 2,0) = h3x * l3y; h( 2,1) = d3x * d3y; h( 2,2) = l3x * h3y;
832 h( 3,0) = h0x * l3y; h( 3,1) = d0x * d3y; h( 3,2) = l0x * h3y;
833 h( 4,0) = h1x * l0y; h( 4,1) = d1x * d0y; h( 4,2) = l1x * h0y;
834 h( 5,0) = h2x * l0y; h( 5,1) = d2x * d0y; h( 5,2) = l2x * h0y;
835 h( 6,0) = h3x * l1y; h( 6,1) = d3x * d1y; h( 6,2) = l3x * h1y;
836 h( 7,0) = h3x * l2y; h( 7,1) = d3x * d2y; h( 7,2) = l3x * h2y;
837 h( 8,0) = h2x * l3y; h( 8,1) = d2x * d3y; h( 8,2) = l2x * h3y;
838 h( 9,0) = h1x * l3y; h( 9,1) = d1x * d3y; h( 9,2) = l1x * h3y;
839 h(10,0) = h0x * l2y; h(10,1) = d0x * d2y; h(10,2) = l0x * h2y;
840 h(11,0) = h0x * l1y; h(11,1) = d0x * d1y; h(11,2) = l0x * h1y;
841 h(12,0) = h1x * l1y; h(12,1) = d1x * d1y; h(12,2) = l1x * h1y;
842 h(13,0) = h2x * l1y; h(13,1) = d2x * d1y; h(13,2) = l2x * h1y;
843 h(14,0) = h1x * l2y; h(14,1) = d1x * d2y; h(14,2) = l1x * h2y;
844 h(15,0) = h2x * l2y; h(15,1) = d2x * d2y; h(15,2) = l2x * h2y;
863 l3 = (0.33333333333333333333-x),
864 l4 = (0.66666666666666666667-x);
866 shape(0) = 4.5 * l2 * l3 * l4;
867 shape(1) = 4.5 * l1 * l3 * l4;
868 shape(2) = 13.5 * l1 * l2 * l4;
869 shape(3) = -13.5 * l1 * l2 * l3;
877 dshape(0,0) = -5.5 + x * (18. - 13.5 * x);
878 dshape(1,0) = 1. - x * (9. - 13.5 * x);
879 dshape(2,0) = 9. - x * (45. - 40.5 * x);
880 dshape(3,0) = -4.5 + x * (36. - 40.5 * x);
913 real_t l1 = (-1. + x + y),
917 shape(0) = -0.5*l1*(3.*l1 + 1.)*(3.*l1 + 2.);
918 shape(1) = 0.5*x*(lx - 1.)*lx;
919 shape(2) = 0.5*y*(-1. + ly)*ly;
920 shape(3) = 4.5*x*l1*(3.*l1 + 1.);
921 shape(4) = -4.5*x*lx*l1;
922 shape(5) = 4.5*x*lx*y;
923 shape(6) = 4.5*x*y*ly;
924 shape(7) = -4.5*y*l1*ly;
925 shape(8) = 4.5*y*l1*(1. + 3.*l1);
926 shape(9) = -27.*x*y*l1;
934 dshape(0,0) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
935 dshape(1,0) = 1. + 4.5*x*(-2. + 3.*x);
937 dshape(3,0) = 4.5*(2. + 9.*x*x - 5.*y + 3.*y*y + 2.*x*(-5. + 6.*y));
938 dshape(4,0) = -4.5*(1. - 1.*y + x*(-8. + 9.*x + 6.*y));
939 dshape(5,0) = 4.5*(-1. + 6.*x)*y;
940 dshape(6,0) = 4.5*y*(-1. + 3.*y);
941 dshape(7,0) = 4.5*(1. - 3.*y)*y;
942 dshape(8,0) = 4.5*y*(-5. + 6.*x + 6.*y);
943 dshape(9,0) = -27.*y*(-1. + 2.*x + y);
945 dshape(0,1) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
947 dshape(2,1) = 1. + 4.5*y*(-2. + 3.*y);
948 dshape(3,1) = 4.5*x*(-5. + 6.*x + 6.*y);
949 dshape(4,1) = 4.5*(1. - 3.*x)*x;
950 dshape(5,1) = 4.5*x*(-1. + 3.*x);
951 dshape(6,1) = 4.5*x*(-1. + 6.*y);
952 dshape(7,1) = -4.5*(1. + x*(-1. + 6.*y) + y*(-8. + 9.*y));
953 dshape(8,1) = 4.5*(2. + 3.*x*x + y*(-10. + 9.*y) + x*(-5. + 12.*y));
954 dshape(9,1) = -27.*x*(-1. + x + 2.*y);
962 h(0,0) = 18.-27.*(x+y);
963 h(0,1) = 18.-27.*(x+y);
964 h(0,2) = 18.-27.*(x+y);
974 h(3,0) = -45.+81.*x+54.*y;
975 h(3,1) = -22.5+54.*x+27.*y;
978 h(4,0) = 36.-81.*x-27.*y;
992 h(7,2) = 36.-27.*x-81.*y;
995 h(8,1) = -22.5+27.*x+54.*y;
996 h(8,2) = -45.+54.*x+81.*y;
999 h(9,1) = 27.-54.*(x+y);
1074 shape(0) = -((-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z)*
1075 (-1 + 3*x + 3*y + 3*z))/2.;
1076 shape(4) = (9*x*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
1077 shape(5) = (-9*x*(-1 + 3*x)*(-1 + x + y + z))/2.;
1078 shape(1) = (x*(2 + 9*(-1 + x)*x))/2.;
1079 shape(6) = (9*y*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
1080 shape(19) = -27*x*y*(-1 + x + y + z);
1081 shape(10) = (9*x*(-1 + 3*x)*y)/2.;
1082 shape(7) = (-9*y*(-1 + 3*y)*(-1 + x + y + z))/2.;
1083 shape(11) = (9*x*y*(-1 + 3*y))/2.;
1084 shape(2) = (y*(2 + 9*(-1 + y)*y))/2.;
1085 shape(8) = (9*z*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
1086 shape(18) = -27*x*z*(-1 + x + y + z);
1087 shape(12) = (9*x*(-1 + 3*x)*z)/2.;
1088 shape(17) = -27*y*z*(-1 + x + y + z);
1089 shape(16) = 27*x*y*z;
1090 shape(14) = (9*y*(-1 + 3*y)*z)/2.;
1091 shape(9) = (-9*z*(-1 + x + y + z)*(-1 + 3*z))/2.;
1092 shape(13) = (9*x*z*(-1 + 3*z))/2.;
1093 shape(15) = (9*y*z*(-1 + 3*z))/2.;
1094 shape(3) = (z*(2 + 9*(-1 + z)*z))/2.;
1102 dshape(0,0) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
1103 x*(-4 + 6*y + 6*z)))/2.;
1104 dshape(0,1) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
1105 x*(-4 + 6*y + 6*z)))/2.;
1106 dshape(0,2) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
1107 x*(-4 + 6*y + 6*z)))/2.;
1108 dshape(4,0) = (9*(9*pow(x,2) + (-1 + y + z)*(-2 + 3*y + 3*z) +
1109 2*x*(-5 + 6*y + 6*z)))/2.;
1110 dshape(4,1) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
1111 dshape(4,2) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
1112 dshape(5,0) = (-9*(1 - y - z + x*(-8 + 9*x + 6*y + 6*z)))/2.;
1113 dshape(5,1) = (9*(1 - 3*x)*x)/2.;
1114 dshape(5,2) = (9*(1 - 3*x)*x)/2.;
1115 dshape(1,0) = 1 + (9*x*(-2 + 3*x))/2.;
1118 dshape(6,0) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
1119 dshape(6,1) = (9*(2 + 3*pow(x,2) - 10*y - 5*z + 3*(y + z)*(3*y + z) +
1120 x*(-5 + 12*y + 6*z)))/2.;
1121 dshape(6,2) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
1122 dshape(19,0) = -27*y*(-1 + 2*x + y + z);
1123 dshape(19,1) = -27*x*(-1 + x + 2*y + z);
1124 dshape(19,2) = -27*x*y;
1125 dshape(10,0) = (9*(-1 + 6*x)*y)/2.;
1126 dshape(10,1) = (9*x*(-1 + 3*x))/2.;
1128 dshape(7,0) = (9*(1 - 3*y)*y)/2.;
1129 dshape(7,1) = (-9*(1 + x*(-1 + 6*y) - z + y*(-8 + 9*y + 6*z)))/2.;
1130 dshape(7,2) = (9*(1 - 3*y)*y)/2.;
1131 dshape(11,0) = (9*y*(-1 + 3*y))/2.;
1132 dshape(11,1) = (9*x*(-1 + 6*y))/2.;
1135 dshape(2,1) = 1 + (9*y*(-2 + 3*y))/2.;
1137 dshape(8,0) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
1138 dshape(8,1) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
1139 dshape(8,2) = (9*(2 + 3*pow(x,2) - 5*y - 10*z + 3*(y + z)*(y + 3*z) +
1140 x*(-5 + 6*y + 12*z)))/2.;
1141 dshape(18,0) = -27*z*(-1 + 2*x + y + z);
1142 dshape(18,1) = -27*x*z;
1143 dshape(18,2) = -27*x*(-1 + x + y + 2*z);
1144 dshape(12,0) = (9*(-1 + 6*x)*z)/2.;
1146 dshape(12,2) = (9*x*(-1 + 3*x))/2.;
1147 dshape(17,0) = -27*y*z;
1148 dshape(17,1) = -27*z*(-1 + x + 2*y + z);
1149 dshape(17,2) = -27*y*(-1 + x + y + 2*z);
1150 dshape(16,0) = 27*y*z;
1151 dshape(16,1) = 27*x*z;
1152 dshape(16,2) = 27*x*y;
1154 dshape(14,1) = (9*(-1 + 6*y)*z)/2.;
1155 dshape(14,2) = (9*y*(-1 + 3*y))/2.;
1156 dshape(9,0) = (9*(1 - 3*z)*z)/2.;
1157 dshape(9,1) = (9*(1 - 3*z)*z)/2.;
1158 dshape(9,2) = (9*(-1 + x + y + 8*z - 6*(x + y)*z - 9*pow(z,2)))/2.;
1159 dshape(13,0) = (9*z*(-1 + 3*z))/2.;
1161 dshape(13,2) = (9*x*(-1 + 6*z))/2.;
1163 dshape(15,1) = (9*z*(-1 + 3*z))/2.;
1164 dshape(15,2) = (9*y*(-1 + 6*z))/2.;
1167 dshape(3,2) = 1 + (9*z*(-2 + 3*z))/2.;
1233 shape(0) = 1. - ip.
x - ip.
y - ip.
z;
1242 if (dshape.
Height() == 4)
1244 real_t *A = &dshape(0,0);
1245 A[0] = -1.; A[4] = -1.; A[8] = -1.;
1246 A[1] = 1.; A[5] = 0.; A[9] = 0.;
1247 A[2] = 0.; A[6] = 1.; A[10] = 0.;
1248 A[3] = 0.; A[7] = 0.; A[11] = 1.;
1252 dshape(0,0) = -1.; dshape(0,1) = -1.; dshape(0,2) = -1.;
1253 dshape(1,0) = 1.; dshape(1,1) = 0.; dshape(1,2) = 0.;
1254 dshape(2,0) = 0.; dshape(2,1) = 1.; dshape(2,2) = 0.;
1255 dshape(3,0) = 0.; dshape(3,1) = 0.; dshape(3,2) = 1.;
1262 static int face_dofs[4][3] = {{1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {0, 1, 2}};
1265 *dofs = face_dofs[face];
1296 shape(0) = (1. - ip.
x - ip.
y) * (1. - ip.
z);
1297 shape(1) = ip.
x * (1. - ip.
z);
1298 shape(2) = ip.
y * (1. - ip.
z);
1299 shape(3) = (1. - ip.
x - ip.
y) * ip.
z;
1300 shape(4) = ip.
x * ip.
z;
1301 shape(5) = ip.
y * ip.
z;
1307 dshape(0,0) = -1. + ip.
z;
1308 dshape(0,1) = -1. + ip.
z;
1309 dshape(0,2) = -1. + ip.
x + ip.
y;
1311 dshape(1,0) = 1. - ip.
z;
1313 dshape(1,2) = -ip.
x;
1316 dshape(2,1) = 1. - ip.
z;
1317 dshape(2,2) = -ip.
y;
1319 dshape(3,0) = -ip.
z;
1320 dshape(3,1) = -ip.
z;
1321 dshape(3,2) = 1. - ip.
x - ip.
y;
1335 static int face_dofs[5][4] =
1336 {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {1, 2, 5, 4}, {2, 0, 3, 5}};
1338 *ndofs = (face < 2) ? 3 : 4;
1339 *dofs = face_dofs[face];
1367 real_t ox = 1.-x-z, oy = 1.-y-z, oz = 1.-z;
1386 shape(0) = ox * oy * ozi;
1387 shape(1) = x * oy * ozi;
1388 shape(2) = x * y * ozi;
1389 shape(3) = ox * y * ozi;
1397 real_t ox = 1.-x-z, oy = 1.-y-z, oz = 1.-z;
1408 dshape(0,0) = - 0.5;
1409 dshape(0,1) = - 0.5;
1410 dshape(0,2) = - 0.75;
1413 dshape(1,1) = - 0.5;
1414 dshape(1,2) = - 0.25;
1420 dshape(3,0) = - 0.5;
1422 dshape(3,2) = - 0.25;
1433 dshape(0,0) = - oy * ozi;
1434 dshape(0,1) = - ox * ozi;
1435 dshape(0,2) = x * y * ozi * ozi - 1.;
1437 dshape(1,0) = oy * ozi;
1438 dshape(1,1) = - x * ozi;
1439 dshape(1,2) = - x * y * ozi * ozi;
1441 dshape(2,0) = y * ozi;
1442 dshape(2,1) = x * ozi;
1443 dshape(2,2) = x * y * ozi * ozi;
1445 dshape(3,0) = - y * ozi;
1446 dshape(3,1) = ox * ozi;
1447 dshape(3,2) = - x * y * ozi * ozi;
1457 static int face_dofs[5][4] =
1458 {{3, 2, 1, 0}, {0, 1, 4, -1}, {1, 2, 4, -1}, {2, 3, 4, -1}, {3, 0, 4, -1}};
1460 *ndofs = (face < 1) ? 4 : 3;
1461 *dofs = face_dofs[face];
1505 L0 = 1. - ip.
x - ip.
y - ip.
z;
1510 shape(0) = L0 * ( 2.0 * L0 - 1.0 );
1511 shape(1) = L1 * ( 2.0 * L1 - 1.0 );
1512 shape(2) = L2 * ( 2.0 * L2 - 1.0 );
1513 shape(3) = L3 * ( 2.0 * L3 - 1.0 );
1514 shape(4) = 4.0 * L0 * L1;
1515 shape(5) = 4.0 * L0 * L2;
1516 shape(6) = 4.0 * L0 * L3;
1517 shape(7) = 4.0 * L1 * L2;
1518 shape(8) = 4.0 * L1 * L3;
1519 shape(9) = 4.0 * L2 * L3;
1530 L0 = 1.0 - x - y - z;
1532 dshape(0,0) = dshape(0,1) = dshape(0,2) = 1.0 - 4.0 * L0;
1533 dshape(1,0) = -1.0 + 4.0 * x; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
1534 dshape(2,0) = 0.0; dshape(2,1) = -1.0 + 4.0 * y; dshape(2,2) = 0.0;
1535 dshape(3,0) = dshape(3,1) = 0.0; dshape(3,2) = -1.0 + 4.0 * z;
1536 dshape(4,0) = 4.0 * (L0 - x); dshape(4,1) = dshape(4,2) = -4.0 * x;
1537 dshape(5,0) = dshape(5,2) = -4.0 * y; dshape(5,1) = 4.0 * (L0 - y);
1538 dshape(6,0) = dshape(6,1) = -4.0 * z; dshape(6,2) = 4.0 * (L0 - z);
1539 dshape(7,0) = 4.0 * y; dshape(7,1) = 4.0 * x; dshape(7,2) = 0.0;
1540 dshape(8,0) = 4.0 * z; dshape(8,1) = 0.0; dshape(8,2) = 4.0 * x;
1541 dshape(9,0) = 0.0; dshape(9,1) = 4.0 * z; dshape(9,2) = 4.0 * y;
1584 real_t ox = 1.-x, oy = 1.-y, oz = 1.-z;
1586 shape(0) = ox * oy * oz;
1587 shape(1) = x * oy * oz;
1588 shape(2) = x * y * oz;
1589 shape(3) = ox * y * oz;
1590 shape(4) = ox * oy * z;
1591 shape(5) = x * oy * z;
1592 shape(6) = x * y * z;
1593 shape(7) = ox * y * z;
1600 real_t ox = 1.-x, oy = 1.-y, oz = 1.-z;
1602 dshape(0,0) = - oy * oz;
1603 dshape(0,1) = - ox * oz;
1604 dshape(0,2) = - ox * oy;
1606 dshape(1,0) = oy * oz;
1607 dshape(1,1) = - x * oz;
1608 dshape(1,2) = - x * oy;
1610 dshape(2,0) = y * oz;
1611 dshape(2,1) = x * oz;
1612 dshape(2,2) = - x * y;
1614 dshape(3,0) = - y * oz;
1615 dshape(3,1) = ox * oz;
1616 dshape(3,2) = - ox * y;
1618 dshape(4,0) = - oy * z;
1619 dshape(4,1) = - ox * z;
1620 dshape(4,2) = ox * oy;
1622 dshape(5,0) = oy * z;
1623 dshape(5,1) = - x * z;
1624 dshape(5,2) = x * oy;
1626 dshape(6,0) = y * z;
1627 dshape(6,1) = x * z;
1628 dshape(6,2) = x * y;
1630 dshape(7,0) = - y * z;
1631 dshape(7,1) = ox * z;
1632 dshape(7,2) = ox * y;
1668 shape(0) = 1.0 - 2.0 * ip.
y;
1669 shape(1) = -1.0 + 2.0 * ( ip.
x + ip.
y );
1670 shape(2) = 1.0 - 2.0 * ip.
x;
1676 dshape(0,0) = 0.0; dshape(0,1) = -2.0;
1677 dshape(1,0) = 2.0; dshape(1,1) = 2.0;
1678 dshape(2,0) = -2.0; dshape(2,1) = 0.0;
1699 const real_t l1 = ip.
x+ip.
y-0.5, l2 = 1.-l1, l3 = ip.
x-ip.
y+0.5, l4 = 1.-l3;
1710 const real_t x2 = 2.*ip.
x, y2 = 2.*ip.
y;
1712 dshape(0,0) = 1. - x2; dshape(0,1) = -2. + y2;
1713 dshape(1,0) = x2; dshape(1,1) = 1. - y2;
1714 dshape(2,0) = 1. - x2; dshape(2,1) = y2;
1715 dshape(3,0) = -2. + x2; dshape(3,1) = 1. - y2;
1736 shape(0,1) = y - 1.;
1739 shape(2,0) = x - 1.;
1751const real_t RT0TriangleFiniteElement::nk[3][2] =
1752{ {0, -1}, {1, 1}, {-1, 0} };
1758#ifdef MFEM_THREAD_SAFE
1763 for (k = 0; k < 3; k++)
1766 for (j = 0; j < 3; j++)
1769 if (j == k) { d -= 1.0; }
1770 if (fabs(d) > 1.0e-12)
1772 mfem::err <<
"RT0TriangleFiniteElement::GetLocalInterpolation (...)\n"
1773 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
1790 for (k = 0; k < 3; k++)
1793 ip.
x = vk[0]; ip.
y = vk[1];
1796 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
1797 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
1798 for (j = 0; j < 3; j++)
1799 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
1813 for (
int k = 0; k < 3; k++)
1821 dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
1822 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
1846 shape(0,1) = y - 1.;
1851 shape(3,0) = x - 1.;
1864const real_t RT0QuadFiniteElement::nk[4][2] =
1865{ {0, -1}, {1, 0}, {0, 1}, {-1, 0} };
1871#ifdef MFEM_THREAD_SAFE
1876 for (k = 0; k < 4; k++)
1879 for (j = 0; j < 4; j++)
1882 if (j == k) { d -= 1.0; }
1883 if (fabs(d) > 1.0e-12)
1885 mfem::err <<
"RT0QuadFiniteElement::GetLocalInterpolation (...)\n"
1886 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
1903 for (k = 0; k < 4; k++)
1906 ip.
x = vk[0]; ip.
y = vk[1];
1909 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
1910 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
1911 for (j = 0; j < 4; j++)
1912 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
1926 for (
int k = 0; k < 4; k++)
1934 dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
1935 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
1965 shape(0,0) = -2 * x * (-1 + x + 2 * y);
1966 shape(0,1) = -2 * (-1 + y) * (-1 + x + 2 * y);
1967 shape(1,0) = 2 * x * (x - y);
1968 shape(1,1) = 2 * (x - y) * (-1 + y);
1969 shape(2,0) = 2 * x * (-1 + 2 * x + y);
1970 shape(2,1) = 2 * y * (-1 + 2 * x + y);
1971 shape(3,0) = 2 * x * (-1 + x + 2 * y);
1972 shape(3,1) = 2 * y * (-1 + x + 2 * y);
1973 shape(4,0) = -2 * (-1 + x) * (x - y);
1974 shape(4,1) = 2 * y * (-x + y);
1975 shape(5,0) = -2 * (-1 + x) * (-1 + 2 * x + y);
1976 shape(5,1) = -2 * y * (-1 + 2 * x + y);
1977 shape(6,0) = -3 * x * (-2 + 2 * x + y);
1978 shape(6,1) = -3 * y * (-1 + 2 * x + y);
1979 shape(7,0) = -3 * x * (-1 + x + 2 * y);
1980 shape(7,1) = -3 * y * (-2 + x + 2 * y);
1988 divshape(0) = -2 * (-4 + 3 * x + 6 * y);
1989 divshape(1) = 2 + 6 * x - 6 * y;
1990 divshape(2) = -4 + 12 * x + 6 * y;
1991 divshape(3) = -4 + 6 * x + 12 * y;
1992 divshape(4) = 2 - 6 * x + 6 * y;
1993 divshape(5) = -2 * (-4 + 6 * x + 3 * y);
1994 divshape(6) = -9 * (-1 + 2 * x + y);
1995 divshape(7) = -9 * (-1 + x + 2 * y);
1998const real_t RT1TriangleFiniteElement::nk[8][2] =
2010#ifdef MFEM_THREAD_SAFE
2015 for (k = 0; k < 8; k++)
2018 for (j = 0; j < 8; j++)
2021 if (j == k) { d -= 1.0; }
2022 if (fabs(d) > 1.0e-12)
2024 mfem::err <<
"RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
2025 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2042 for (k = 0; k < 8; k++)
2045 ip.
x = vk[0]; ip.
y = vk[1];
2048 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
2049 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
2050 for (j = 0; j < 8; j++)
2051 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2064 for (
int k = 0; k < 8; k++)
2072 dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
2073 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
2121 shape(0,1) = -( 1. - 3.*y + 2.*y*y)*( 2. - 3.*x);
2123 shape(1,1) = -( 1. - 3.*y + 2.*y*y)*(-1. + 3.*x);
2125 shape(2,0) = (-x + 2.*x*x)*( 2. - 3.*y);
2127 shape(3,0) = (-x + 2.*x*x)*(-1. + 3.*y);
2131 shape(4,1) = (-y + 2.*y*y)*(-1. + 3.*x);
2133 shape(5,1) = (-y + 2.*y*y)*( 2. - 3.*x);
2135 shape(6,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y);
2137 shape(7,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y);
2140 shape(8,0) = (4.*x - 4.*x*x)*( 2. - 3.*y);
2142 shape(9,0) = (4.*x - 4.*x*x)*(-1. + 3.*y);
2146 shape(10,1) = (4.*y - 4.*y*y)*( 2. - 3.*x);
2148 shape(11,1) = (4.*y - 4.*y*y)*(-1. + 3.*x);
2156 divshape(0) = -(-3. + 4.*y)*( 2. - 3.*x);
2157 divshape(1) = -(-3. + 4.*y)*(-1. + 3.*x);
2158 divshape(2) = (-1. + 4.*x)*( 2. - 3.*y);
2159 divshape(3) = (-1. + 4.*x)*(-1. + 3.*y);
2160 divshape(4) = (-1. + 4.*y)*(-1. + 3.*x);
2161 divshape(5) = (-1. + 4.*y)*( 2. - 3.*x);
2162 divshape(6) = -(-3. + 4.*x)*(-1. + 3.*y);
2163 divshape(7) = -(-3. + 4.*x)*( 2. - 3.*y);
2164 divshape(8) = ( 4. - 8.*x)*( 2. - 3.*y);
2165 divshape(9) = ( 4. - 8.*x)*(-1. + 3.*y);
2166 divshape(10) = ( 4. - 8.*y)*( 2. - 3.*x);
2167 divshape(11) = ( 4. - 8.*y)*(-1. + 3.*x);
2170const real_t RT1QuadFiniteElement::nk[12][2] =
2190#ifdef MFEM_THREAD_SAFE
2195 for (k = 0; k < 12; k++)
2198 for (j = 0; j < 12; j++)
2201 if (j == k) { d -= 1.0; }
2202 if (fabs(d) > 1.0e-12)
2204 mfem::err <<
"RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
2205 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2222 for (k = 0; k < 12; k++)
2225 ip.
x = vk[0]; ip.
y = vk[1];
2228 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
2229 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
2230 for (j = 0; j < 12; j++)
2231 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2244 for (
int k = 0; k < 12; k++)
2252 dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
2253 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
2257const real_t RT2TriangleFiniteElement::M[15][15] =
2261 0, -5.3237900077244501311, 5.3237900077244501311, 16.647580015448900262,
2262 0, 24.442740046346700787, -16.647580015448900262, -12.,
2263 -19.118950038622250656, -47.237900077244501311, 0, -34.414110069520051180,
2264 12., 30.590320061795601049, 15.295160030897800524
2267 0, 1.5, -1.5, -15., 0, 2.625, 15., 15., -4.125, 30., 0, -14.625, -15.,
2271 0, -0.67620999227554986889, 0.67620999227554986889, 7.3524199845510997378,
2272 0, -3.4427400463467007866, -7.3524199845510997378, -12.,
2273 4.1189500386222506555, -0.76209992275549868892, 0, 7.4141100695200511800,
2274 12., -6.5903200617956010489, -3.2951600308978005244
2277 0, 0, 1.5, 0, 0, 1.5, -11.471370023173350393, 0, 2.4713700231733503933,
2278 -11.471370023173350393, 0, 2.4713700231733503933, 15.295160030897800524,
2279 0, -3.2951600308978005244
2282 0, 0, 4.875, 0, 0, 4.875, -16.875, 0, -16.875, -16.875, 0, -16.875, 10.5,
2286 0, 0, 1.5, 0, 0, 1.5, 2.4713700231733503933, 0, -11.471370023173350393,
2287 2.4713700231733503933, 0, -11.471370023173350393, -3.2951600308978005244,
2288 0, 15.295160030897800524
2291 -0.67620999227554986889, 0, -3.4427400463467007866, 0,
2292 7.3524199845510997378, 0.67620999227554986889, 7.4141100695200511800, 0,
2293 -0.76209992275549868892, 4.1189500386222506555, -12.,
2294 -7.3524199845510997378, -3.2951600308978005244, -6.5903200617956010489,
2298 1.5, 0, 2.625, 0, -15., -1.5, -14.625, 0, 30., -4.125, 15., 15., 10.5,
2302 -5.3237900077244501311, 0, 24.442740046346700787, 0, 16.647580015448900262,
2303 5.3237900077244501311, -34.414110069520051180, 0, -47.237900077244501311,
2304 -19.118950038622250656, -12., -16.647580015448900262, 15.295160030897800524,
2305 30.590320061795601049, 12.
2307 { 0, 0, 18., 0, 0, 6., -42., 0, -30., -26., 0, -14., 24., 32., 8.},
2308 { 0, 0, 6., 0, 0, 18., -14., 0, -26., -30., 0, -42., 8., 32., 24.},
2309 { 0, 0, -6., 0, 0, -4., 30., 0, 4., 22., 0, 4., -24., -16., 0},
2310 { 0, 0, -4., 0, 0, -8., 20., 0, 8., 36., 0, 8., -16., -32., 0},
2311 { 0, 0, -8., 0, 0, -4., 8., 0, 36., 8., 0, 20., 0, -32., -16.},
2312 { 0, 0, -4., 0, 0, -6., 4., 0, 22., 4., 0, 30., 0, -16., -24.}
2319 const real_t p = 0.11270166537925831148;
2358 real_t Bx[15] = {1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y, 0., x*x*x,
2361 real_t By[15] = {0., 1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y,
2365 for (
int i = 0; i < 15; i++)
2367 real_t cx = 0.0, cy = 0.0;
2368 for (
int j = 0; j < 15; j++)
2370 cx += M[i][j] * Bx[j];
2371 cy += M[i][j] * By[j];
2382 constexpr real_t f2 = 2.0;
2383 constexpr real_t f4 = 4.0;
2385 real_t DivB[15] = {0., 0., 1., 0., 0., 1., f2*x, 0., y, x, 0., f2*y,
2386 f4*x*x, f4*x*y, f4*y*y
2389 for (
int i = 0; i < 15; i++)
2392 for (
int j = 0; j < 15; j++)
2394 div += M[i][j] * DivB[j];
2400const real_t RT2QuadFiniteElement::pt[4] = {0.,1./3.,2./3.,1.};
2402const real_t RT2QuadFiniteElement::dpt[3] = {0.25,0.5,0.75};
2465 real_t A01 = pt[0] - pt[1];
2466 real_t A02 = pt[0] - pt[2];
2467 real_t A12 = pt[1] - pt[2];
2468 real_t A03 = pt[0] - pt[3];
2469 real_t A13 = pt[1] - pt[3];
2470 real_t A23 = pt[2] - pt[3];
2472 real_t B01 = dpt[0] - dpt[1];
2473 real_t B02 = dpt[0] - dpt[2];
2474 real_t B12 = dpt[1] - dpt[2];
2476 real_t tx0 = (bx1*bx2)/(B01*B02);
2477 real_t tx1 = -(bx0*bx2)/(B01*B12);
2478 real_t tx2 = (bx0*bx1)/(B02*B12);
2480 real_t ty0 = (by1*by2)/(B01*B02);
2481 real_t ty1 = -(by0*by2)/(B01*B12);
2482 real_t ty2 = (by0*by1)/(B02*B12);
2486 shape(0, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx0;
2488 shape(1, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx1;
2490 shape(2, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx2;
2492 shape(3, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty0;
2494 shape(4, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty1;
2496 shape(5, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty2;
2500 shape(6, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx2;
2502 shape(7, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx1;
2504 shape(8, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx0;
2506 shape(9, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty2;
2508 shape(10, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty1;
2510 shape(11, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty0;
2513 shape(12, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty0;
2515 shape(13, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty1;
2517 shape(14, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty2;
2520 shape(15, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty0;
2522 shape(16, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty1;
2524 shape(17, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty2;
2528 shape(18, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx0;
2530 shape(19, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx1;
2532 shape(20, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx2;
2535 shape(21, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx0;
2537 shape(22, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx1;
2539 shape(23, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx2;
2547 real_t a01 = pt[0]*pt[1];
2548 real_t a02 = pt[0]*pt[2];
2549 real_t a12 = pt[1]*pt[2];
2550 real_t a03 = pt[0]*pt[3];
2551 real_t a13 = pt[1]*pt[3];
2552 real_t a23 = pt[2]*pt[3];
2562 real_t A01 = pt[0] - pt[1];
2563 real_t A02 = pt[0] - pt[2];
2564 real_t A12 = pt[1] - pt[2];
2565 real_t A03 = pt[0] - pt[3];
2566 real_t A13 = pt[1] - pt[3];
2567 real_t A23 = pt[2] - pt[3];
2569 real_t A012 = pt[0] + pt[1] + pt[2];
2570 real_t A013 = pt[0] + pt[1] + pt[3];
2571 real_t A023 = pt[0] + pt[2] + pt[3];
2572 real_t A123 = pt[1] + pt[2] + pt[3];
2574 real_t B01 = dpt[0] - dpt[1];
2575 real_t B02 = dpt[0] - dpt[2];
2576 real_t B12 = dpt[1] - dpt[2];
2578 real_t tx0 = (bx1*bx2)/(B01*B02);
2579 real_t tx1 = -(bx0*bx2)/(B01*B12);
2580 real_t tx2 = (bx0*bx1)/(B02*B12);
2582 real_t ty0 = (by1*by2)/(B01*B02);
2583 real_t ty1 = -(by0*by2)/(B01*B12);
2584 real_t ty2 = (by0*by1)/(B02*B12);
2587 divshape(0) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx0;
2588 divshape(1) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx1;
2589 divshape(2) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx2;
2591 divshape(3) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty0;
2592 divshape(4) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty1;
2593 divshape(5) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty2;
2595 divshape(6) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx2;
2596 divshape(7) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx1;
2597 divshape(8) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx0;
2599 divshape(9) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty2;
2600 divshape(10) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty1;
2601 divshape(11) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty0;
2603 divshape(12) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty0;
2604 divshape(13) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty1;
2605 divshape(14) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty2;
2607 divshape(15) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty0;
2608 divshape(16) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty1;
2609 divshape(17) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty2;
2611 divshape(18) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx0;
2612 divshape(19) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx1;
2613 divshape(20) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx2;
2615 divshape(21) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx0;
2616 divshape(22) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx1;
2617 divshape(23) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx2;
2620const real_t RT2QuadFiniteElement::nk[24][2] =
2623 {0,-1}, {0,-1}, {0,-1},
2625 {1, 0}, {1, 0}, {1, 0},
2627 {0, 1}, {0, 1}, {0, 1},
2629 {-1,0}, {-1,0}, {-1,0},
2631 {1, 0}, {1, 0}, {1, 0},
2633 {1, 0}, {1, 0}, {1, 0},
2635 {0, 1}, {0, 1}, {0, 1},
2637 {0, 1}, {0, 1}, {0, 1}
2644#ifdef MFEM_THREAD_SAFE
2649 for (k = 0; k < 24; k++)
2652 for (j = 0; j < 24; j++)
2655 if (j == k) { d -= 1.0; }
2656 if (fabs(d) > 1.0e-12)
2658 mfem::err <<
"RT2QuadFiniteElement::GetLocalInterpolation (...)\n"
2659 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2676 for (k = 0; k < 24; k++)
2679 ip.
x = vk[0]; ip.
y = vk[1];
2682 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
2683 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
2684 for (j = 0; j < 24; j++)
2685 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2698 for (
int k = 0; k < 24; k++)
2706 dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
2707 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
2723 shape(0) = 2. - 3. * x;
2724 shape(1) = 3. * x - 1.;
2738 const real_t p = 0.11270166537925831148;
2748 const real_t p = 0.11270166537925831148;
2749 const real_t w = 1./((1-2*
p)*(1-2*
p));
2752 shape(0) = (2*x-1)*(x-1+
p)*w;
2753 shape(1) = 4*(x-1+
p)*(
p-x)*w;
2754 shape(2) = (2*x-1)*(x-
p)*w;
2760 const real_t p = 0.11270166537925831148;
2761 const real_t w = 1./((1-2*
p)*(1-2*
p));
2764 dshape(0,0) = (-3+4*x+2*
p)*w;
2765 dshape(1,0) = (4-8*x)*w;
2766 dshape(2,0) = (-1+4*x-2*
p)*w;
2777 for (i = 1; i < m; i++)
2783#ifndef MFEM_THREAD_SAFE
2788 for (i = 1; i <= m; i++)
2792 for (i = 0; i < m/2+1; i++)
2794 rwk(m-i) = ( rwk(i) *= rwk(m-i) );
2796 for (i = m-1; i >= 0; i -= 2)
2808#ifdef MFEM_THREAD_SAFE
2812 k = (int) floor ( m * x + 0.5 );
2813 k = k > m ? m : k < 0 ? 0 : k;
2816 for (i = 0; i <= m; i++)
2819 wk *= ( rxxk(i) = x - (
real_t)(i) / m );
2821 w = wk * ( rxxk(k) = x - (
real_t)(k) / m );
2825 shape(0) = w * rwk(0) / rxxk(0);
2829 shape(0) = wk * rwk(0);
2833 shape(1) = w * rwk(m) / rxxk(m);
2837 shape(1) = wk * rwk(k);
2839 for (i = 1; i < m; i++)
2842 shape(i+1) = w * rwk(i) / rxxk(i);
2846 shape(k+1) = wk * rwk(k);
2856#ifdef MFEM_THREAD_SAFE
2860 k = (int) floor ( m * x + 0.5 );
2861 k = k > m ? m : k < 0 ? 0 : k;
2864 for (i = 0; i <= m; i++)
2867 wk *= ( rxxk(i) = x - (
real_t)(i) / m );
2869 w = wk * ( rxxk(k) = x - (
real_t)(k) / m );
2871 for (i = 0; i <= m; i++)
2873 rxxk(i) = 1.0 / rxxk(i);
2876 for (i = 0; i <= m; i++)
2885 dshape(0,0) = (
s - w * rxxk(0)) * rwk(0) * rxxk(0);
2889 dshape(0,0) = wk * srx * rwk(0);
2893 dshape(1,0) = (
s - w * rxxk(m)) * rwk(m) * rxxk(m);
2897 dshape(1,0) = wk * srx * rwk(k);
2899 for (i = 1; i < m; i++)
2902 dshape(i+1,0) = (
s - w * rxxk(i)) * rwk(i) * rxxk(i);
2906 dshape(k+1,0) = wk * srx * rwk(k);
2937 L1 = ip.
x; L2 = ip.
y; L3 = ip.
z; L0 = 1.0 - L1 - L2 - L3;
2938 shape(0) = 1.0 - 3.0 * L0;
2939 shape(1) = 1.0 - 3.0 * L1;
2940 shape(2) = 1.0 - 3.0 * L2;
2941 shape(3) = 1.0 - 3.0 * L3;
2947 dshape(0,0) = 3.0; dshape(0,1) = 3.0; dshape(0,2) = 3.0;
2948 dshape(1,0) = -3.0; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
2949 dshape(2,0) = 0.0; dshape(2,1) = -3.0; dshape(2,2) = 0.0;
2950 dshape(3,0) = 0.0; dshape(3,1) = 0.0; dshape(3,2) = -3.0;
2971 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
2992 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3013 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3034 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3048 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3049 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3050 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3051 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3052 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3053 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3054 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3055 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3057 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3058 I[ 9] = 1; J[ 9] = 2; K[ 9] = 0;
3059 I[10] = 2; J[10] = 1; K[10] = 0;
3060 I[11] = 0; J[11] = 2; K[11] = 0;
3061 I[12] = 2; J[12] = 0; K[12] = 1;
3062 I[13] = 1; J[13] = 2; K[13] = 1;
3063 I[14] = 2; J[14] = 1; K[14] = 1;
3064 I[15] = 0; J[15] = 2; K[15] = 1;
3065 I[16] = 0; J[16] = 0; K[16] = 2;
3066 I[17] = 1; J[17] = 0; K[17] = 2;
3067 I[18] = 1; J[18] = 1; K[18] = 2;
3068 I[19] = 0; J[19] = 1; K[19] = 2;
3070 I[20] = 2; J[20] = 2; K[20] = 0;
3071 I[21] = 2; J[21] = 0; K[21] = 2;
3072 I[22] = 1; J[22] = 2; K[22] = 2;
3073 I[23] = 2; J[23] = 1; K[23] = 2;
3074 I[24] = 0; J[24] = 2; K[24] = 2;
3075 I[25] = 2; J[25] = 2; K[25] = 1;
3077 I[26] = 2; J[26] = 2; K[26] = 2;
3079 else if (degree == 3)
3085 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3086 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3087 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3088 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3089 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3090 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3091 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3092 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3094 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3095 I[ 9] = 3; J[ 9] = 0; K[ 9] = 0;
3096 I[10] = 1; J[10] = 2; K[10] = 0;
3097 I[11] = 1; J[11] = 3; K[11] = 0;
3098 I[12] = 2; J[12] = 1; K[12] = 0;
3099 I[13] = 3; J[13] = 1; K[13] = 0;
3100 I[14] = 0; J[14] = 2; K[14] = 0;
3101 I[15] = 0; J[15] = 3; K[15] = 0;
3102 I[16] = 2; J[16] = 0; K[16] = 1;
3103 I[17] = 3; J[17] = 0; K[17] = 1;
3104 I[18] = 1; J[18] = 2; K[18] = 1;
3105 I[19] = 1; J[19] = 3; K[19] = 1;
3106 I[20] = 2; J[20] = 1; K[20] = 1;
3107 I[21] = 3; J[21] = 1; K[21] = 1;
3108 I[22] = 0; J[22] = 2; K[22] = 1;
3109 I[23] = 0; J[23] = 3; K[23] = 1;
3110 I[24] = 0; J[24] = 0; K[24] = 2;
3111 I[25] = 0; J[25] = 0; K[25] = 3;
3112 I[26] = 1; J[26] = 0; K[26] = 2;
3113 I[27] = 1; J[27] = 0; K[27] = 3;
3114 I[28] = 1; J[28] = 1; K[28] = 2;
3115 I[29] = 1; J[29] = 1; K[29] = 3;
3116 I[30] = 0; J[30] = 1; K[30] = 2;
3117 I[31] = 0; J[31] = 1; K[31] = 3;
3119 I[32] = 2; J[32] = 3; K[32] = 0;
3120 I[33] = 3; J[33] = 3; K[33] = 0;
3121 I[34] = 2; J[34] = 2; K[34] = 0;
3122 I[35] = 3; J[35] = 2; K[35] = 0;
3123 I[36] = 2; J[36] = 0; K[36] = 2;
3124 I[37] = 3; J[37] = 0; K[37] = 2;
3125 I[38] = 2; J[38] = 0; K[38] = 3;
3126 I[39] = 3; J[39] = 0; K[39] = 3;
3127 I[40] = 1; J[40] = 2; K[40] = 2;
3128 I[41] = 1; J[41] = 3; K[41] = 2;
3129 I[42] = 1; J[42] = 2; K[42] = 3;
3130 I[43] = 1; J[43] = 3; K[43] = 3;
3131 I[44] = 3; J[44] = 1; K[44] = 2;
3132 I[45] = 2; J[45] = 1; K[45] = 2;
3133 I[46] = 3; J[46] = 1; K[46] = 3;
3134 I[47] = 2; J[47] = 1; K[47] = 3;
3135 I[48] = 0; J[48] = 3; K[48] = 2;
3136 I[49] = 0; J[49] = 2; K[49] = 2;
3137 I[50] = 0; J[50] = 3; K[50] = 3;
3138 I[51] = 0; J[51] = 2; K[51] = 3;
3139 I[52] = 2; J[52] = 2; K[52] = 1;
3140 I[53] = 3; J[53] = 2; K[53] = 1;
3141 I[54] = 2; J[54] = 3; K[54] = 1;
3142 I[55] = 3; J[55] = 3; K[55] = 1;
3144 I[56] = 2; J[56] = 2; K[56] = 2;
3145 I[57] = 3; J[57] = 2; K[57] = 2;
3146 I[58] = 3; J[58] = 3; K[58] = 2;
3147 I[59] = 2; J[59] = 3; K[59] = 2;
3148 I[60] = 2; J[60] = 2; K[60] = 3;
3149 I[61] = 3; J[61] = 2; K[61] = 3;
3150 I[62] = 3; J[62] = 3; K[62] = 3;
3151 I[63] = 2; J[63] = 3; K[63] = 3;
3155 mfem_error (
"LagrangeHexFiniteElement::LagrangeHexFiniteElement");
3159 dof1d = fe1d ->
GetDof();
3161#ifndef MFEM_THREAD_SAFE
3171 for (
int n = 0; n <
dof; n++)
3186#ifdef MFEM_THREAD_SAFE
3187 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
3194 for (
int n = 0; n <
dof; n++)
3196 shape(n) = shape1dx(I[n]) * shape1dy(J[n]) * shape1dz(K[n]);
3207#ifdef MFEM_THREAD_SAFE
3208 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
3209 DenseMatrix dshape1dx(dof1d,1), dshape1dy(dof1d,1), dshape1dz(dof1d,1);
3220 for (
int n = 0; n <
dof; n++)
3222 dshape(n,0) = dshape1dx(I[n],0) * shape1dy(J[n]) * shape1dz(K[n]);
3223 dshape(n,1) = shape1dx(I[n]) * dshape1dy(J[n],0) * shape1dz(K[n]);
3224 dshape(n,2) = shape1dx(I[n]) * shape1dy(J[n]) * dshape1dz(K[n],0);
3253 shape(0) = 1.0 - 2.0 * x;
3260 shape(1) = 2.0 * x - 1.0;
3261 shape(2) = 2.0 - 2.0 * x;
3272 dshape(0,0) = - 2.0;
3280 dshape(2,0) = - 2.0;
3307 L0 = 2.0 * ( 1. - ip.
x - ip.
y );
3308 L1 = 2.0 * ( ip.
x );
3309 L2 = 2.0 * ( ip.
y );
3318 for (i = 0; i < 6; i++)
3325 shape(0) = L0 - 1.0;
3332 shape(1) = L1 - 1.0;
3339 shape(2) = L2 - 1.0;
3343 shape(3) = 1.0 - L2;
3344 shape(4) = 1.0 - L0;
3345 shape(5) = 1.0 - L1;
3355 L0 = 2.0 * ( 1. - ip.
x - ip.
y );
3356 L1 = 2.0 * ( ip.
x );
3357 L2 = 2.0 * ( ip.
y );
3359 real_t DL0[2], DL1[2], DL2[2];
3360 DL0[0] = -2.0; DL0[1] = -2.0;
3361 DL1[0] = 2.0; DL1[1] = 0.0;
3362 DL2[0] = 0.0; DL2[1] = 2.0;
3364 for (i = 0; i < 6; i++)
3365 for (j = 0; j < 2; j++)
3372 for (j = 0; j < 2; j++)
3374 dshape(0,j) = DL0[j];
3375 dshape(3,j) = DL1[j];
3376 dshape(5,j) = DL2[j];
3381 for (j = 0; j < 2; j++)
3383 dshape(3,j) = DL0[j];
3384 dshape(1,j) = DL1[j];
3385 dshape(4,j) = DL2[j];
3390 for (j = 0; j < 2; j++)
3392 dshape(5,j) = DL0[j];
3393 dshape(4,j) = DL1[j];
3394 dshape(2,j) = DL2[j];
3399 for (j = 0; j < 2; j++)
3401 dshape(3,j) = - DL2[j];
3402 dshape(4,j) = - DL0[j];
3403 dshape(5,j) = - DL1[j];
3448 real_t L0, L1, L2, L3, L4, L5;
3449 L0 = 2.0 * ( 1. - ip.
x - ip.
y - ip.
z );
3450 L1 = 2.0 * ( ip.
x );
3451 L2 = 2.0 * ( ip.
y );
3452 L3 = 2.0 * ( ip.
z );
3453 L4 = 2.0 * ( ip.
x + ip.
y );
3454 L5 = 2.0 * ( ip.
y + ip.
z );
3467 for (i = 0; i < 10; i++)
3474 shape(0) = L0 - 1.0;
3482 shape(1) = L1 - 1.0;
3490 shape(2) = L2 - 1.0;
3498 shape(3) = L3 - 1.0;
3500 else if ((L4 <= 1.0) && (L5 <= 1.0))
3502 shape(4) = 1.0 - L5;
3504 shape(6) = 1.0 - L4;
3505 shape(8) = 1.0 - L0;
3507 else if ((L4 >= 1.0) && (L5 <= 1.0))
3509 shape(4) = 1.0 - L5;
3510 shape(5) = 1.0 - L1;
3511 shape(7) = L4 - 1.0;
3514 else if ((L4 <= 1.0) && (L5 >= 1.0))
3516 shape(5) = 1.0 - L3;
3517 shape(6) = 1.0 - L4;
3519 shape(9) = L5 - 1.0;
3521 else if ((L4 >= 1.0) && (L5 >= 1.0))
3524 shape(7) = L4 - 1.0;
3525 shape(8) = 1.0 - L2;
3526 shape(9) = L5 - 1.0;
3535 real_t L0, L1, L2, L3, L4, L5;
3536 L0 = 2.0 * ( 1. - ip.
x - ip.
y - ip.
z );
3537 L1 = 2.0 * ( ip.
x );
3538 L2 = 2.0 * ( ip.
y );
3539 L3 = 2.0 * ( ip.
z );
3540 L4 = 2.0 * ( ip.
x + ip.
y );
3541 L5 = 2.0 * ( ip.
y + ip.
z );
3543 real_t DL0[3], DL1[3], DL2[3], DL3[3], DL4[3], DL5[3];
3544 DL0[0] = -2.0; DL0[1] = -2.0; DL0[2] = -2.0;
3545 DL1[0] = 2.0; DL1[1] = 0.0; DL1[2] = 0.0;
3546 DL2[0] = 0.0; DL2[1] = 2.0; DL2[2] = 0.0;
3547 DL3[0] = 0.0; DL3[1] = 0.0; DL3[2] = 2.0;
3548 DL4[0] = 2.0; DL4[1] = 2.0; DL4[2] = 0.0;
3549 DL5[0] = 0.0; DL5[1] = 2.0; DL5[2] = 2.0;
3551 for (i = 0; i < 10; i++)
3552 for (j = 0; j < 3; j++)
3559 for (j = 0; j < 3; j++)
3561 dshape(0,j) = DL0[j];
3562 dshape(4,j) = DL1[j];
3563 dshape(5,j) = DL2[j];
3564 dshape(6,j) = DL3[j];
3569 for (j = 0; j < 3; j++)
3571 dshape(4,j) = DL0[j];
3572 dshape(1,j) = DL1[j];
3573 dshape(7,j) = DL2[j];
3574 dshape(8,j) = DL3[j];
3579 for (j = 0; j < 3; j++)
3581 dshape(5,j) = DL0[j];
3582 dshape(7,j) = DL1[j];
3583 dshape(2,j) = DL2[j];
3584 dshape(9,j) = DL3[j];
3589 for (j = 0; j < 3; j++)
3591 dshape(6,j) = DL0[j];
3592 dshape(8,j) = DL1[j];
3593 dshape(9,j) = DL2[j];
3594 dshape(3,j) = DL3[j];
3597 else if ((L4 <= 1.0) && (L5 <= 1.0))
3599 for (j = 0; j < 3; j++)
3601 dshape(4,j) = - DL5[j];
3602 dshape(5,j) = DL2[j];
3603 dshape(6,j) = - DL4[j];
3604 dshape(8,j) = - DL0[j];
3607 else if ((L4 >= 1.0) && (L5 <= 1.0))
3609 for (j = 0; j < 3; j++)
3611 dshape(4,j) = - DL5[j];
3612 dshape(5,j) = - DL1[j];
3613 dshape(7,j) = DL4[j];
3614 dshape(8,j) = DL3[j];
3617 else if ((L4 <= 1.0) && (L5 >= 1.0))
3619 for (j = 0; j < 3; j++)
3621 dshape(5,j) = - DL3[j];
3622 dshape(6,j) = - DL4[j];
3623 dshape(8,j) = DL1[j];
3624 dshape(9,j) = DL5[j];
3627 else if ((L4 >= 1.0) && (L5 >= 1.0))
3629 for (j = 0; j < 3; j++)
3631 dshape(5,j) = DL0[j];
3632 dshape(7,j) = DL4[j];
3633 dshape(8,j) = - DL2[j];
3634 dshape(9,j) = DL5[j];
3669 Lx = 2.0 * ( 1. - x );
3670 Ly = 2.0 * ( 1. - y );
3679 for (i = 0; i < 9; i++)
3684 if ((x <= 0.5) && (y <= 0.5))
3686 shape(0) = (Lx - 1.0) * (Ly - 1.0);
3687 shape(4) = (2.0 - Lx) * (Ly - 1.0);
3688 shape(8) = (2.0 - Lx) * (2.0 - Ly);
3689 shape(7) = (Lx - 1.0) * (2.0 - Ly);
3691 else if ((x >= 0.5) && (y <= 0.5))
3693 shape(4) = Lx * (Ly - 1.0);
3694 shape(1) = (1.0 - Lx) * (Ly - 1.0);
3695 shape(5) = (1.0 - Lx) * (2.0 - Ly);
3696 shape(8) = Lx * (2.0 - Ly);
3698 else if ((x >= 0.5) && (y >= 0.5))
3700 shape(8) = Lx * Ly ;
3701 shape(5) = (1.0 - Lx) * Ly ;
3702 shape(2) = (1.0 - Lx) * (1.0 - Ly);
3703 shape(6) = Lx * (1.0 - Ly);
3705 else if ((x <= 0.5) && (y >= 0.5))
3707 shape(7) = (Lx - 1.0) * Ly ;
3708 shape(8) = (2.0 - Lx) * Ly ;
3709 shape(6) = (2.0 - Lx) * (1.0 - Ly);
3710 shape(3) = (Lx - 1.0) * (1.0 - Ly);
3720 Lx = 2.0 * ( 1. - x );
3721 Ly = 2.0 * ( 1. - y );
3723 for (i = 0; i < 9; i++)
3724 for (j = 0; j < 2; j++)
3729 if ((x <= 0.5) && (y <= 0.5))
3731 dshape(0,0) = 2.0 * (1.0 - Ly);
3732 dshape(0,1) = 2.0 * (1.0 - Lx);
3734 dshape(4,0) = 2.0 * (Ly - 1.0);
3735 dshape(4,1) = -2.0 * (2.0 - Lx);
3737 dshape(8,0) = 2.0 * (2.0 - Ly);
3738 dshape(8,1) = 2.0 * (2.0 - Lx);
3740 dshape(7,0) = -2.0 * (2.0 - Ly);
3741 dshape(7,0) = 2.0 * (Lx - 1.0);
3743 else if ((x >= 0.5) && (y <= 0.5))
3745 dshape(4,0) = -2.0 * (Ly - 1.0);
3746 dshape(4,1) = -2.0 * Lx;
3748 dshape(1,0) = 2.0 * (Ly - 1.0);
3749 dshape(1,1) = -2.0 * (1.0 - Lx);
3751 dshape(5,0) = 2.0 * (2.0 - Ly);
3752 dshape(5,1) = 2.0 * (1.0 - Lx);
3754 dshape(8,0) = -2.0 * (2.0 - Ly);
3755 dshape(8,1) = 2.0 * Lx;
3757 else if ((x >= 0.5) && (y >= 0.5))
3759 dshape(8,0) = -2.0 * Ly;
3760 dshape(8,1) = -2.0 * Lx;
3762 dshape(5,0) = 2.0 * Ly;
3763 dshape(5,1) = -2.0 * (1.0 - Lx);
3765 dshape(2,0) = 2.0 * (1.0 - Ly);
3766 dshape(2,1) = 2.0 * (1.0 - Lx);
3768 dshape(6,0) = -2.0 * (1.0 - Ly);
3769 dshape(6,1) = 2.0 * Lx;
3771 else if ((x <= 0.5) && (y >= 0.5))
3773 dshape(7,0) = -2.0 * Ly;
3774 dshape(7,1) = -2.0 * (Lx - 1.0);
3776 dshape(8,0) = 2.0 * Ly ;
3777 dshape(8,1) = -2.0 * (2.0 - Lx);
3779 dshape(6,0) = 2.0 * (1.0 - Ly);
3780 dshape(6,1) = 2.0 * (2.0 - Lx);
3782 dshape(3,0) = -2.0 * (1.0 - Ly);
3783 dshape(3,1) = 2.0 * (Lx - 1.0);
3794 I[ 0] = 0.0; J[ 0] = 0.0; K[ 0] = 0.0;
3795 I[ 1] = 1.0; J[ 1] = 0.0; K[ 1] = 0.0;
3796 I[ 2] = 1.0; J[ 2] = 1.0; K[ 2] = 0.0;
3797 I[ 3] = 0.0; J[ 3] = 1.0; K[ 3] = 0.0;
3798 I[ 4] = 0.0; J[ 4] = 0.0; K[ 4] = 1.0;
3799 I[ 5] = 1.0; J[ 5] = 0.0; K[ 5] = 1.0;
3800 I[ 6] = 1.0; J[ 6] = 1.0; K[ 6] = 1.0;
3801 I[ 7] = 0.0; J[ 7] = 1.0; K[ 7] = 1.0;
3803 I[ 8] = 0.5; J[ 8] = 0.0; K[ 8] = 0.0;
3804 I[ 9] = 1.0; J[ 9] = 0.5; K[ 9] = 0.0;
3805 I[10] = 0.5; J[10] = 1.0; K[10] = 0.0;
3806 I[11] = 0.0; J[11] = 0.5; K[11] = 0.0;
3807 I[12] = 0.5; J[12] = 0.0; K[12] = 1.0;
3808 I[13] = 1.0; J[13] = 0.5; K[13] = 1.0;
3809 I[14] = 0.5; J[14] = 1.0; K[14] = 1.0;
3810 I[15] = 0.0; J[15] = 0.5; K[15] = 1.0;
3811 I[16] = 0.0; J[16] = 0.0; K[16] = 0.5;
3812 I[17] = 1.0; J[17] = 0.0; K[17] = 0.5;
3813 I[18] = 1.0; J[18] = 1.0; K[18] = 0.5;
3814 I[19] = 0.0; J[19] = 1.0; K[19] = 0.5;
3816 I[20] = 0.5; J[20] = 0.5; K[20] = 0.0;
3817 I[21] = 0.5; J[21] = 0.0; K[21] = 0.5;
3818 I[22] = 1.0; J[22] = 0.5; K[22] = 0.5;
3819 I[23] = 0.5; J[23] = 1.0; K[23] = 0.5;
3820 I[24] = 0.0; J[24] = 0.5; K[24] = 0.5;
3821 I[25] = 0.5; J[25] = 0.5; K[25] = 1.0;
3823 I[26] = 0.5; J[26] = 0.5; K[26] = 0.5;
3825 for (
int n = 0; n < 27; n++)
3840 for (i = 0; i < 27; i++)
3845 if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5))
3860 else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5))
3875 else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5))
3890 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5))
3905 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5))
3920 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5))
3935 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5))
3966 shape(N[0]) = Lx * Ly * Lz;
3967 shape(N[1]) = (1 - Lx) * Ly * Lz;
3968 shape(N[2]) = (1 - Lx) * (1 - Ly) * Lz;
3969 shape(N[3]) = Lx * (1 - Ly) * Lz;
3970 shape(N[4]) = Lx * Ly * (1 - Lz);
3971 shape(N[5]) = (1 - Lx) * Ly * (1 - Lz);
3972 shape(N[6]) = (1 - Lx) * (1 - Ly) * (1 - Lz);
3973 shape(N[7]) = Lx * (1 - Ly) * (1 - Lz);
3983 for (i = 0; i < 27; i++)
3984 for (j = 0; j < 3; j++)
3989 if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5))
4004 else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5))
4019 else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5))
4034 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5))
4049 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5))
4064 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5))
4079 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5))
4110 dshape(N[0],0) = -2.0 * Ly * Lz ;
4111 dshape(N[0],1) = -2.0 * Lx * Lz ;
4112 dshape(N[0],2) = -2.0 * Lx * Ly ;
4114 dshape(N[1],0) = 2.0 * Ly * Lz ;
4115 dshape(N[1],1) = -2.0 * (1 - Lx) * Lz ;
4116 dshape(N[1],2) = -2.0 * (1 - Lx) * Ly ;
4118 dshape(N[2],0) = 2.0 * (1 - Ly) * Lz ;
4119 dshape(N[2],1) = 2.0 * (1 - Lx) * Lz ;
4120 dshape(N[2],2) = -2.0 * (1 - Lx) * (1 - Ly);
4122 dshape(N[3],0) = -2.0 * (1 - Ly) * Lz ;
4123 dshape(N[3],1) = 2.0 * Lx * Lz ;
4124 dshape(N[3],2) = -2.0 * Lx * (1 - Ly);
4126 dshape(N[4],0) = -2.0 * Ly * (1 - Lz);
4127 dshape(N[4],1) = -2.0 * Lx * (1 - Lz);
4128 dshape(N[4],2) = 2.0 * Lx * Ly ;
4130 dshape(N[5],0) = 2.0 * Ly * (1 - Lz);
4131 dshape(N[5],1) = -2.0 * (1 - Lx) * (1 - Lz);
4132 dshape(N[5],2) = 2.0 * (1 - Lx) * Ly ;
4134 dshape(N[6],0) = 2.0 * (1 - Ly) * (1 - Lz);
4135 dshape(N[6],1) = 2.0 * (1 - Lx) * (1 - Lz);
4136 dshape(N[6],2) = 2.0 * (1 - Lx) * (1 - Ly);
4138 dshape(N[7],0) = -2.0 * (1 - Ly) * (1 - Lz);
4139 dshape(N[7],1) = 2.0 * Lx * (1 - Lz);
4140 dshape(N[7],2) = 2.0 * Lx * (1 - Ly);
4203 shape(0,0) = (1. - y) * (1. - z);
4207 shape(2,0) = y * (1. - z);
4211 shape(4,0) = z * (1. - y);
4220 shape(1,1) = x * (1. - z);
4224 shape(3,1) = (1. - x) * (1. - z);
4232 shape(7,1) = (1. - x) * z;
4237 shape(8,2) = (1. - x) * (1. - y);
4241 shape(9,2) = x * (1. - y);
4245 shape(10,2) = x * y;
4249 shape(11,2) = y * (1. - x);
4259 curl_shape(0,0) = 0.;
4260 curl_shape(0,1) = y - 1.;
4261 curl_shape(0,2) = 1. - z;
4263 curl_shape(2,0) = 0.;
4264 curl_shape(2,1) = -y;
4265 curl_shape(2,2) = z - 1.;
4267 curl_shape(4,0) = 0;
4268 curl_shape(4,1) = 1. - y;
4269 curl_shape(4,2) = z;
4271 curl_shape(6,0) = 0.;
4272 curl_shape(6,1) = y;
4273 curl_shape(6,2) = -z;
4275 curl_shape(1,0) = x;
4276 curl_shape(1,1) = 0.;
4277 curl_shape(1,2) = 1. - z;
4279 curl_shape(3,0) = 1. - x;
4280 curl_shape(3,1) = 0.;
4281 curl_shape(3,2) = z - 1.;
4283 curl_shape(5,0) = -x;
4284 curl_shape(5,1) = 0.;
4285 curl_shape(5,2) = z;
4287 curl_shape(7,0) = x - 1.;
4288 curl_shape(7,1) = 0.;
4289 curl_shape(7,2) = -z;
4291 curl_shape(8,0) = x - 1.;
4292 curl_shape(8,1) = 1. - y;
4293 curl_shape(8,2) = 0.;
4295 curl_shape(9,0) = -x;
4296 curl_shape(9,1) = y - 1.;
4297 curl_shape(9,2) = 0;
4299 curl_shape(10,0) = x;
4300 curl_shape(10,1) = -y;
4301 curl_shape(10,2) = 0.;
4303 curl_shape(11,0) = 1. - x;
4304 curl_shape(11,1) = y;
4305 curl_shape(11,2) = 0.;
4308const real_t Nedelec1HexFiniteElement::tk[12][3] =
4310 {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
4311 {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
4312 {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}
4319#ifdef MFEM_THREAD_SAFE
4324 for (k = 0; k <
dof; k++)
4327 for (j = 0; j <
dof; j++)
4331 if (j == k) { d -= 1.0; }
4332 if (fabs(d) > 1.0e-12)
4334 mfem::err <<
"Nedelec1HexFiniteElement::GetLocalInterpolation (...)\n"
4335 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
4343 ip.
x = ip.
y = ip.
z = 0.0;
4350 for (k = 0; k <
dof; k++)
4353 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
4356 vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
4357 vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
4358 vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
4359 for (j = 0; j <
dof; j++)
4360 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
4361 vshape(j,2)*vk[2])) < 1.0e-12)
4375 for (
int k = 0; k <
dof; k++)
4383 vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
4384 vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
4385 vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
4397 for (
int k = 0; k <
dof; k++)
4400 dshape.
Mult(tk[k], grad_k);
4401 for (
int j = 0; j < grad_k.
Size(); j++)
4403 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
4443 shape(0,0) = 1. - y - z;
4448 shape(1,1) = 1. - x - z;
4453 shape(2,2) = 1. - x - y;
4472 curl_shape(0,0) = 0.;
4473 curl_shape(0,1) = -2.;
4474 curl_shape(0,2) = 2.;
4476 curl_shape(1,0) = 2.;
4477 curl_shape(1,1) = 0.;
4478 curl_shape(1,2) = -2.;
4480 curl_shape(2,0) = -2.;
4481 curl_shape(2,1) = 2.;
4482 curl_shape(2,2) = 0.;
4484 curl_shape(3,0) = 0.;
4485 curl_shape(3,1) = 0.;
4486 curl_shape(3,2) = 2.;
4488 curl_shape(4,0) = 0.;
4489 curl_shape(4,1) = -2.;
4490 curl_shape(4,2) = 0.;
4492 curl_shape(5,0) = 2.;
4493 curl_shape(5,1) = 0.;
4494 curl_shape(5,2) = 0.;
4497const real_t Nedelec1TetFiniteElement::tk[6][3] =
4498{{1,0,0}, {0,1,0}, {0,0,1}, {-1,1,0}, {-1,0,1}, {0,-1,1}};
4504#ifdef MFEM_THREAD_SAFE
4509 for (k = 0; k <
dof; k++)
4512 for (j = 0; j <
dof; j++)
4516 if (j == k) { d -= 1.0; }
4517 if (fabs(d) > 1.0e-12)
4519 mfem::err <<
"Nedelec1TetFiniteElement::GetLocalInterpolation (...)\n"
4520 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
4528 ip.
x = ip.
y = ip.
z = 0.0;
4535 for (k = 0; k <
dof; k++)
4538 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
4541 vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
4542 vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
4543 vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
4544 for (j = 0; j <
dof; j++)
4545 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
4546 vshape(j,2)*vk[2])) < 1.0e-12)
4560 for (
int k = 0; k <
dof; k++)
4568 vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
4569 vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
4570 vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
4582 for (
int k = 0; k <
dof; k++)
4585 dshape.
Mult(tk[k], grad_k);
4586 for (
int j = 0; j < grad_k.
Size(); j++)
4588 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
4640 shape(0,0) = (1. - y) * (1. - z);
4641 shape(0,1) = x * (1. - z);
4644 shape(1,0) = - y * (1. - z);
4645 shape(1,1) = x * (1. - z);
4648 shape(2,0) = - y * (1. - z);
4649 shape(2,1) = - (1. - x) * (1. - z);
4652 shape(3,0) = (1. - y) * z;
4656 shape(4,0) = - y * z;
4660 shape(5,0) = - y * z;
4661 shape(5,1) = - (1. - x) * z;
4666 shape(6,2) = 1. - x - y;
4681 real_t x = ip.
x, y = ip.
y, z2 = 2. * ip.
z;
4683 curl_shape(0,0) = x;
4684 curl_shape(0,1) = - 1. + y;
4685 curl_shape(0,2) = 2. - z2;
4687 curl_shape(1,0) = x;
4688 curl_shape(1,1) = y;
4689 curl_shape(1,2) = 2. - z2;
4691 curl_shape(2,0) = - 1. + x;
4692 curl_shape(2,1) = y;
4693 curl_shape(2,2) = 2. - z2;
4695 curl_shape(3,0) = - x;
4696 curl_shape(3,1) = 1. - y;
4697 curl_shape(3,2) = z2;
4699 curl_shape(4,0) = - x;
4700 curl_shape(4,1) = - y;
4701 curl_shape(4,2) = z2;
4703 curl_shape(5,0) = 1. - x;
4704 curl_shape(5,1) = - y;
4705 curl_shape(5,2) = z2;
4707 curl_shape(6,0) = - 1.;
4708 curl_shape(6,1) = 1.;
4709 curl_shape(6,2) = 0.;
4711 curl_shape(7,0) = 0.;
4712 curl_shape(7,1) = - 1.;
4713 curl_shape(7,2) = 0.;
4715 curl_shape(8,0) = 1.;
4716 curl_shape(8,1) = 0.;
4717 curl_shape(8,2) = 0.;
4720const real_t Nedelec1WdgFiniteElement::tk[9][3] =
4722 {1,0,0}, {-1,1,0}, {0,-1,0}, {1,0,0}, {-1,1,0}, {0,-1,0},
4723 {0,0,1}, {0,0,1}, {0,0,1}
4730#ifdef MFEM_THREAD_SAFE
4735 for (k = 0; k <
dof; k++)
4738 for (j = 0; j <
dof; j++)
4742 if (j == k) { d -= 1.0; }
4743 if (fabs(d) > 1.0e-12)
4745 mfem::err <<
"Nedelec1WdgFiniteElement::GetLocalInterpolation (...)\n"
4746 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
4754 ip.
x = ip.
y = ip.
z = 0.0;
4761 for (k = 0; k <
dof; k++)
4764 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
4767 vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
4768 vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
4769 vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
4770 for (j = 0; j <
dof; j++)
4771 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
4772 vshape(j,2)*vk[2])) < 1.0e-12)
4786 for (
int k = 0; k <
dof; k++)
4794 vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
4795 vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
4796 vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
4808 for (
int k = 0; k <
dof; k++)
4811 dshape.
Mult(tk[k], grad_k);
4812 for (
int j = 0; j < grad_k.
Size(); j++)
4814 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
4860 real_t x = ip.
x, y = ip.
y, z = ip.
z, z2 = 2. * ip.
z;
4861 real_t ox = 1. - x - z, oy = 1. - y - z, oz = 1. - z;
4909 shape(0,2) = x * oy * ozi;
4913 shape(1,2) = x * y * ozi;
4917 shape(2,2) = x * y * ozi;
4921 shape(3,2) = ox * y * ozi;
4923 shape(4,0) = oy * z * ozi;
4924 shape(4,1) = ox * z * ozi;
4925 shape(4,2) = 1. - x - y + x * y * (1. - z2) * ozi * ozi;
4927 shape(5,0) = - oy * z * ozi;
4928 shape(5,1) = x * z * ozi;
4929 shape(5,2) = x * (1. - y * (1. - z2) * ozi * ozi);
4931 shape(6,0) = - y * z * ozi;
4932 shape(6,1) = - x * z * ozi;
4933 shape(6,2) = x * y * (1. - z2) * ozi * ozi;
4935 shape(7,0) = y * z * ozi;
4936 shape(7,1) = - ox * z * ozi;
4937 shape(7,2) = y * (1. - x * (1. - z2) * ozi * ozi);
4944 real_t x = ip.
x, y = ip.
y, z = ip.
z, z2 = 2. * z;
4945 real_t ox = 1. - x - z, oy = 1. - y - z, oz = 1. - z;
4955 curl_shape(0,0) = 0.;
4956 curl_shape(0,1) = - 2.;
4957 curl_shape(0,2) = 1.;
4959 curl_shape(1,0) = 0.;
4960 curl_shape(1,1) = 0.;
4961 curl_shape(1,2) = 1.;
4963 curl_shape(2,0) = 0.;
4964 curl_shape(2,1) = 0.;
4965 curl_shape(2,2) = - 1.;
4967 curl_shape(3,0) = 2.;
4968 curl_shape(3,1) = 0.;
4969 curl_shape(3,2) = - 1.;
4971 curl_shape(4,0) = - 2.;
4972 curl_shape(4,1) = 2.;
4973 curl_shape(4,2) = 0.;
4975 curl_shape(5,0) = 0.;
4976 curl_shape(5,1) = - 2.;
4977 curl_shape(5,2) = 0.;
4979 curl_shape(6,0) = 0.;
4980 curl_shape(6,1) = 0.;
4981 curl_shape(6,2) = 0.;
4983 curl_shape(7,0) = 2.;
4984 curl_shape(7,1) = 0.;
4985 curl_shape(7,2) = 0.;
4992 curl_shape(0,0) = - x * ozi;
4993 curl_shape(0,1) = - 2. + y * ozi;
4994 curl_shape(0,2) = 1.;
4996 curl_shape(1,0) = x * ozi;
4997 curl_shape(1,1) = - y * ozi;
4998 curl_shape(1,2) = 1.;
5000 curl_shape(2,0) = x * ozi;
5001 curl_shape(2,1) = - y * ozi;
5002 curl_shape(2,2) = - 1.;
5004 curl_shape(3,0) = (2. - x - z2) * ozi;
5005 curl_shape(3,1) = y * ozi;
5006 curl_shape(3,2) = - 1.;
5008 curl_shape(4,0) = - 2. * ox * ozi;
5009 curl_shape(4,1) = 2. * oy * ozi;
5010 curl_shape(4,2) = 0.;
5012 curl_shape(5,0) = - 2. * x * ozi;
5013 curl_shape(5,1) = - 2. * oy * ozi;
5014 curl_shape(5,2) = 0.;
5016 curl_shape(6,0) = 2. * x * ozi;
5017 curl_shape(6,1) = - 2. * y * ozi;
5018 curl_shape(6,2) = 0.;
5020 curl_shape(7,0) = 2. * ox * ozi;
5021 curl_shape(7,1) = 2. * y * ozi;
5022 curl_shape(7,2) = 0.;
5025const real_t Nedelec1PyrFiniteElement::tk[8][3] =
5026{{1,0,0}, {0,1,0}, {1,0,0}, {0,1,0}, {0,0,1}, {-1,0,1}, {-1,-1,1}, {0,-1,1}};
5032#ifdef MFEM_THREAD_SAFE
5037 for (k = 0; k <
dof; k++)
5040 for (j = 0; j <
dof; j++)
5044 if (j == k) { d -= 1.0; }
5045 if (fabs(d) > 1.0e-12)
5047 mfem::err <<
"Nedelec1PyrFiniteElement::GetLocalInterpolation (...)\n"
5048 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5056 ip.
x = ip.
y = ip.
z = 0.0;
5063 for (k = 0; k <
dof; k++)
5066 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5069 vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
5070 vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
5071 vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
5072 for (j = 0; j <
dof; j++)
5073 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5074 vshape(j,2)*vk[2])) < 1.0e-12)
5088 for (
int k = 0; k <
dof; k++)
5096 vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
5097 vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
5098 vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
5110 for (
int k = 0; k <
dof; k++)
5113 dshape.
Mult(tk[k], grad_k);
5114 for (
int j = 0; j < grad_k.
Size(); j++)
5116 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);