15 #include "../coefficient.hpp"
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.;
146 const double x = ip.
x, y = ip.
y;
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.;
170 const double GaussBiLinear2DFiniteElement::p[] =
171 { 0.2113248654051871177454256, 0.7886751345948128822545744 };
189 const double x = ip.
x, y = ip.
y;
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]);
200 const double x = ip.
x, y = ip.
y;
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 double 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;
299 double x = ip.
x, y = ip.
y;
300 double 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;
313 double x = ip.
x, y = ip.
y;
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;
379 const double GaussQuad2DFiniteElement::p[] =
380 { 0.0915762135097707434595714634022015, 0.445948490915964886318329253883051 };
398 for (
int i = 0; i < 6; i++)
415 const double x = ip.
x, y = ip.
y;
429 const double x = ip.
x, y = ip.
y;
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;
467 double x = ip.
x, y = ip.
y;
468 double 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;
491 double x = ip.
x, y = ip.
y;
492 double l1x, l2x, l3x, l1y, l2y, l3y;
493 double 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 double p1 = 0.5*(1.-sqrt(3./5.));
584 const double a = sqrt(5./3.);
585 const double p1 = 0.5*(1.-sqrt(3./5.));
587 double x = a*(ip.
x-p1), y = a*(ip.
y-p1);
588 double 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;
611 const double a = sqrt(5./3.);
612 const double p1 = 0.5*(1.-sqrt(3./5.));
614 double x = a*(ip.
x-p1), y = a*(ip.
y-p1);
615 double l1x, l2x, l3x, l1y, l2y, l3y;
616 double 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;
700 double x = ip.
x, y = ip.
y;
702 double w1x, w2x, w3x, w1y, w2y, w3y;
703 double 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;
739 double x = ip.
x, y = ip.
y;
741 double w1x, w2x, w3x, w1y, w2y, w3y;
742 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
743 double 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;
789 double x = ip.
x, y = ip.
y;
791 double w1x, w2x, w3x, w1y, w2y, w3y;
792 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
793 double d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
794 double 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);
912 double x = ip.
x, y = ip.
y;
913 double 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;
932 double x = ip.
x, y = ip.
y;
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);
960 double x = ip.
x, y = ip.
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);
1072 double x = ip.
x, y = ip.
y, z = ip.
z;
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.;
1100 double x = ip.
x, y = ip.
y, z = ip.
z;
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 double *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];
1366 double x = ip.
x, y = ip.
y, z = ip.
z;
1367 double ox = 1.-x-z, oy = 1.-y-z, oz = 1.-z;
1384 double ozi = 1. / oz;
1386 shape(0) = ox * oy * ozi;
1387 shape(1) = x * oy * ozi;
1388 shape(2) = x * y * ozi;
1389 shape(3) = ox * y * ozi;
1396 double x = ip.
x, y = ip.
y, z = ip.
z;
1397 double 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;
1431 double ozi = 1. / oz;
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];
1503 double L0, L1, L2, L3;
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;
1583 double x = ip.
x, y = ip.
y, z = ip.
z;
1584 double 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;
1599 double x = ip.
x, y = ip.
y, z = ip.
z;
1600 double 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 double l1 = ip.
x+ip.
y-0.5, l2 = 1.-l1, l3 = ip.
x-ip.
y+0.5, l4 = 1.-l3;
1710 const double 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;
1733 double x = ip.
x, y = ip.
y;
1736 shape(0,1) = y - 1.;
1739 shape(2,0) = x - 1.;
1751 const double 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] ));
1843 double x = ip.
x, y = ip.
y;
1846 shape(0,1) = y - 1.;
1851 shape(3,0) = x - 1.;
1864 const double 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] ));
1963 double x = ip.
x, y = ip.
y;
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);
1986 double x = ip.
x, y = ip.
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);
1998 const double 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] ));
2117 double x = ip.
x, y = ip.
y;
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);
2154 double x = ip.
x, y = ip.
y;
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);
2170 const double 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] ));
2257 const double 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 double p = 0.11270166537925831148;
2356 double x = ip.
x, y = ip.
y;
2358 double Bx[15] = {1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y, 0., x*x*x,
2361 double 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 double 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];
2381 double x = ip.
x, y = ip.
y;
2383 double DivB[15] = {0., 0., 1., 0., 0., 1., 2.*x, 0., y, x, 0., 2.*y,
2384 4.*x*x, 4.*x*y, 4.*y*y
2387 for (
int i = 0; i < 15; i++)
2390 for (
int j = 0; j < 15; j++)
2392 div += M[i][j] * DivB[j];
2398 const double RT2QuadFiniteElement::pt[4] = {0.,1./3.,2./3.,1.};
2400 const double RT2QuadFiniteElement::dpt[3] = {0.25,0.5,0.75};
2443 double x = ip.
x, y = ip.
y;
2445 double ax0 = pt[0] - x;
2446 double ax1 = pt[1] - x;
2447 double ax2 = pt[2] - x;
2448 double ax3 = pt[3] - x;
2450 double by0 = dpt[0] - y;
2451 double by1 = dpt[1] - y;
2452 double by2 = dpt[2] - y;
2454 double ay0 = pt[0] - y;
2455 double ay1 = pt[1] - y;
2456 double ay2 = pt[2] - y;
2457 double ay3 = pt[3] - y;
2459 double bx0 = dpt[0] - x;
2460 double bx1 = dpt[1] - x;
2461 double bx2 = dpt[2] - x;
2463 double A01 = pt[0] - pt[1];
2464 double A02 = pt[0] - pt[2];
2465 double A12 = pt[1] - pt[2];
2466 double A03 = pt[0] - pt[3];
2467 double A13 = pt[1] - pt[3];
2468 double A23 = pt[2] - pt[3];
2470 double B01 = dpt[0] - dpt[1];
2471 double B02 = dpt[0] - dpt[2];
2472 double B12 = dpt[1] - dpt[2];
2474 double tx0 = (bx1*bx2)/(B01*B02);
2475 double tx1 = -(bx0*bx2)/(B01*B12);
2476 double tx2 = (bx0*bx1)/(B02*B12);
2478 double ty0 = (by1*by2)/(B01*B02);
2479 double ty1 = -(by0*by2)/(B01*B12);
2480 double ty2 = (by0*by1)/(B02*B12);
2484 shape(0, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx0;
2486 shape(1, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx1;
2488 shape(2, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx2;
2490 shape(3, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty0;
2492 shape(4, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty1;
2494 shape(5, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty2;
2498 shape(6, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx2;
2500 shape(7, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx1;
2502 shape(8, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx0;
2504 shape(9, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty2;
2506 shape(10, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty1;
2508 shape(11, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty0;
2511 shape(12, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty0;
2513 shape(13, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty1;
2515 shape(14, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty2;
2518 shape(15, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty0;
2520 shape(16, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty1;
2522 shape(17, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty2;
2526 shape(18, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx0;
2528 shape(19, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx1;
2530 shape(20, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx2;
2533 shape(21, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx0;
2535 shape(22, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx1;
2537 shape(23, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx2;
2543 double x = ip.
x, y = ip.
y;
2545 double a01 = pt[0]*pt[1];
2546 double a02 = pt[0]*pt[2];
2547 double a12 = pt[1]*pt[2];
2548 double a03 = pt[0]*pt[3];
2549 double a13 = pt[1]*pt[3];
2550 double a23 = pt[2]*pt[3];
2552 double bx0 = dpt[0] - x;
2553 double bx1 = dpt[1] - x;
2554 double bx2 = dpt[2] - x;
2556 double by0 = dpt[0] - y;
2557 double by1 = dpt[1] - y;
2558 double by2 = dpt[2] - y;
2560 double A01 = pt[0] - pt[1];
2561 double A02 = pt[0] - pt[2];
2562 double A12 = pt[1] - pt[2];
2563 double A03 = pt[0] - pt[3];
2564 double A13 = pt[1] - pt[3];
2565 double A23 = pt[2] - pt[3];
2567 double A012 = pt[0] + pt[1] + pt[2];
2568 double A013 = pt[0] + pt[1] + pt[3];
2569 double A023 = pt[0] + pt[2] + pt[3];
2570 double A123 = pt[1] + pt[2] + pt[3];
2572 double B01 = dpt[0] - dpt[1];
2573 double B02 = dpt[0] - dpt[2];
2574 double B12 = dpt[1] - dpt[2];
2576 double tx0 = (bx1*bx2)/(B01*B02);
2577 double tx1 = -(bx0*bx2)/(B01*B12);
2578 double tx2 = (bx0*bx1)/(B02*B12);
2580 double ty0 = (by1*by2)/(B01*B02);
2581 double ty1 = -(by0*by2)/(B01*B12);
2582 double ty2 = (by0*by1)/(B02*B12);
2585 divshape(0) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx0;
2586 divshape(1) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx1;
2587 divshape(2) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx2;
2589 divshape(3) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty0;
2590 divshape(4) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty1;
2591 divshape(5) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty2;
2593 divshape(6) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx2;
2594 divshape(7) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx1;
2595 divshape(8) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx0;
2597 divshape(9) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty2;
2598 divshape(10) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty1;
2599 divshape(11) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty0;
2601 divshape(12) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty0;
2602 divshape(13) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty1;
2603 divshape(14) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty2;
2605 divshape(15) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty0;
2606 divshape(16) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty1;
2607 divshape(17) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty2;
2609 divshape(18) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx0;
2610 divshape(19) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx1;
2611 divshape(20) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx2;
2613 divshape(21) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx0;
2614 divshape(22) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx1;
2615 divshape(23) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx2;
2618 const double RT2QuadFiniteElement::nk[24][2] =
2621 {0,-1}, {0,-1}, {0,-1},
2623 {1, 0}, {1, 0}, {1, 0},
2625 {0, 1}, {0, 1}, {0, 1},
2627 {-1,0}, {-1,0}, {-1,0},
2629 {1, 0}, {1, 0}, {1, 0},
2631 {1, 0}, {1, 0}, {1, 0},
2633 {0, 1}, {0, 1}, {0, 1},
2635 {0, 1}, {0, 1}, {0, 1}
2642 #ifdef MFEM_THREAD_SAFE
2647 for (k = 0; k < 24; k++)
2650 for (j = 0; j < 24; j++)
2653 if (j == k) { d -= 1.0; }
2654 if (fabs(d) > 1.0e-12)
2656 mfem::err <<
"RT2QuadFiniteElement::GetLocalInterpolation (...)\n"
2657 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2674 for (k = 0; k < 24; k++)
2677 ip.
x = vk[0]; ip.
y = vk[1];
2680 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
2681 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
2682 for (j = 0; j < 24; j++)
2683 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2696 for (
int k = 0; k < 24; k++)
2704 dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
2705 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
2721 shape(0) = 2. - 3. * x;
2722 shape(1) = 3. * x - 1.;
2736 const double p = 0.11270166537925831148;
2746 const double p = 0.11270166537925831148;
2747 const double w = 1./((1-2*
p)*(1-2*p));
2750 shape(0) = (2*x-1)*(x-1+p)*w;
2751 shape(1) = 4*(x-1+
p)*(p-x)*w;
2752 shape(2) = (2*x-1)*(x-p)*w;
2758 const double p = 0.11270166537925831148;
2759 const double w = 1./((1-2*
p)*(1-2*p));
2762 dshape(0,0) = (-3+4*x+2*
p)*w;
2763 dshape(1,0) = (4-8*x)*w;
2764 dshape(2,0) = (-1+4*x-2*
p)*w;
2775 for (i = 1; i < m; i++)
2781 #ifndef MFEM_THREAD_SAFE
2786 for (i = 1; i <= m; i++)
2788 rwk(i) = rwk(i-1) * ( (double)(m) / (double)(i) );
2790 for (i = 0; i < m/2+1; i++)
2792 rwk(m-i) = ( rwk(i) *= rwk(m-i) );
2794 for (i = m-1; i >= 0; i -= 2)
2803 double w, wk, x = ip.
x;
2806 #ifdef MFEM_THREAD_SAFE
2810 k = (int) floor ( m * x + 0.5 );
2811 k = k > m ? m : k < 0 ? 0 : k;
2814 for (i = 0; i <= m; i++)
2817 wk *= ( rxxk(i) = x - (double)(i) / m );
2819 w = wk * ( rxxk(k) = x - (double)(k) / m );
2823 shape(0) = w * rwk(0) / rxxk(0);
2827 shape(0) = wk * rwk(0);
2831 shape(1) = w * rwk(m) / rxxk(m);
2835 shape(1) = wk * rwk(k);
2837 for (i = 1; i < m; i++)
2840 shape(i+1) = w * rwk(i) / rxxk(i);
2844 shape(k+1) = wk * rwk(k);
2851 double s, srx, w, wk, x = ip.
x;
2854 #ifdef MFEM_THREAD_SAFE
2858 k = (int) floor ( m * x + 0.5 );
2859 k = k > m ? m : k < 0 ? 0 : k;
2862 for (i = 0; i <= m; i++)
2865 wk *= ( rxxk(i) = x - (double)(i) / m );
2867 w = wk * ( rxxk(k) = x - (double)(k) / m );
2869 for (i = 0; i <= m; i++)
2871 rxxk(i) = 1.0 / rxxk(i);
2874 for (i = 0; i <= m; i++)
2883 dshape(0,0) = (s - w * rxxk(0)) * rwk(0) * rxxk(0);
2887 dshape(0,0) = wk * srx * rwk(0);
2891 dshape(1,0) = (s - w * rxxk(m)) * rwk(m) * rxxk(m);
2895 dshape(1,0) = wk * srx * rwk(k);
2897 for (i = 1; i < m; i++)
2900 dshape(i+1,0) = (s - w * rxxk(i)) * rwk(i) * rxxk(i);
2904 dshape(k+1,0) = wk * srx * rwk(k);
2933 double L0, L1, L2, L3;
2935 L1 = ip.
x; L2 = ip.
y; L3 = ip.
z; L0 = 1.0 - L1 - L2 - L3;
2936 shape(0) = 1.0 - 3.0 * L0;
2937 shape(1) = 1.0 - 3.0 * L1;
2938 shape(2) = 1.0 - 3.0 * L2;
2939 shape(3) = 1.0 - 3.0 * L3;
2945 dshape(0,0) = 3.0; dshape(0,1) = 3.0; dshape(0,2) = 3.0;
2946 dshape(1,0) = -3.0; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
2947 dshape(2,0) = 0.0; dshape(2,1) = -3.0; dshape(2,2) = 0.0;
2948 dshape(3,0) = 0.0; dshape(3,1) = 0.0; dshape(3,2) = -3.0;
2969 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
2990 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3011 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3032 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3046 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3047 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3048 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3049 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3050 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3051 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3052 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3053 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3055 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3056 I[ 9] = 1; J[ 9] = 2; K[ 9] = 0;
3057 I[10] = 2; J[10] = 1; K[10] = 0;
3058 I[11] = 0; J[11] = 2; K[11] = 0;
3059 I[12] = 2; J[12] = 0; K[12] = 1;
3060 I[13] = 1; J[13] = 2; K[13] = 1;
3061 I[14] = 2; J[14] = 1; K[14] = 1;
3062 I[15] = 0; J[15] = 2; K[15] = 1;
3063 I[16] = 0; J[16] = 0; K[16] = 2;
3064 I[17] = 1; J[17] = 0; K[17] = 2;
3065 I[18] = 1; J[18] = 1; K[18] = 2;
3066 I[19] = 0; J[19] = 1; K[19] = 2;
3068 I[20] = 2; J[20] = 2; K[20] = 0;
3069 I[21] = 2; J[21] = 0; K[21] = 2;
3070 I[22] = 1; J[22] = 2; K[22] = 2;
3071 I[23] = 2; J[23] = 1; K[23] = 2;
3072 I[24] = 0; J[24] = 2; K[24] = 2;
3073 I[25] = 2; J[25] = 2; K[25] = 1;
3075 I[26] = 2; J[26] = 2; K[26] = 2;
3077 else if (degree == 3)
3083 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3084 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3085 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3086 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3087 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3088 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3089 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3090 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3092 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3093 I[ 9] = 3; J[ 9] = 0; K[ 9] = 0;
3094 I[10] = 1; J[10] = 2; K[10] = 0;
3095 I[11] = 1; J[11] = 3; K[11] = 0;
3096 I[12] = 2; J[12] = 1; K[12] = 0;
3097 I[13] = 3; J[13] = 1; K[13] = 0;
3098 I[14] = 0; J[14] = 2; K[14] = 0;
3099 I[15] = 0; J[15] = 3; K[15] = 0;
3100 I[16] = 2; J[16] = 0; K[16] = 1;
3101 I[17] = 3; J[17] = 0; K[17] = 1;
3102 I[18] = 1; J[18] = 2; K[18] = 1;
3103 I[19] = 1; J[19] = 3; K[19] = 1;
3104 I[20] = 2; J[20] = 1; K[20] = 1;
3105 I[21] = 3; J[21] = 1; K[21] = 1;
3106 I[22] = 0; J[22] = 2; K[22] = 1;
3107 I[23] = 0; J[23] = 3; K[23] = 1;
3108 I[24] = 0; J[24] = 0; K[24] = 2;
3109 I[25] = 0; J[25] = 0; K[25] = 3;
3110 I[26] = 1; J[26] = 0; K[26] = 2;
3111 I[27] = 1; J[27] = 0; K[27] = 3;
3112 I[28] = 1; J[28] = 1; K[28] = 2;
3113 I[29] = 1; J[29] = 1; K[29] = 3;
3114 I[30] = 0; J[30] = 1; K[30] = 2;
3115 I[31] = 0; J[31] = 1; K[31] = 3;
3117 I[32] = 2; J[32] = 3; K[32] = 0;
3118 I[33] = 3; J[33] = 3; K[33] = 0;
3119 I[34] = 2; J[34] = 2; K[34] = 0;
3120 I[35] = 3; J[35] = 2; K[35] = 0;
3121 I[36] = 2; J[36] = 0; K[36] = 2;
3122 I[37] = 3; J[37] = 0; K[37] = 2;
3123 I[38] = 2; J[38] = 0; K[38] = 3;
3124 I[39] = 3; J[39] = 0; K[39] = 3;
3125 I[40] = 1; J[40] = 2; K[40] = 2;
3126 I[41] = 1; J[41] = 3; K[41] = 2;
3127 I[42] = 1; J[42] = 2; K[42] = 3;
3128 I[43] = 1; J[43] = 3; K[43] = 3;
3129 I[44] = 3; J[44] = 1; K[44] = 2;
3130 I[45] = 2; J[45] = 1; K[45] = 2;
3131 I[46] = 3; J[46] = 1; K[46] = 3;
3132 I[47] = 2; J[47] = 1; K[47] = 3;
3133 I[48] = 0; J[48] = 3; K[48] = 2;
3134 I[49] = 0; J[49] = 2; K[49] = 2;
3135 I[50] = 0; J[50] = 3; K[50] = 3;
3136 I[51] = 0; J[51] = 2; K[51] = 3;
3137 I[52] = 2; J[52] = 2; K[52] = 1;
3138 I[53] = 3; J[53] = 2; K[53] = 1;
3139 I[54] = 2; J[54] = 3; K[54] = 1;
3140 I[55] = 3; J[55] = 3; K[55] = 1;
3142 I[56] = 2; J[56] = 2; K[56] = 2;
3143 I[57] = 3; J[57] = 2; K[57] = 2;
3144 I[58] = 3; J[58] = 3; K[58] = 2;
3145 I[59] = 2; J[59] = 3; K[59] = 2;
3146 I[60] = 2; J[60] = 2; K[60] = 3;
3147 I[61] = 3; J[61] = 2; K[61] = 3;
3148 I[62] = 3; J[62] = 3; K[62] = 3;
3149 I[63] = 2; J[63] = 3; K[63] = 3;
3153 mfem_error (
"LagrangeHexFiniteElement::LagrangeHexFiniteElement");
3157 dof1d = fe1d ->
GetDof();
3159 #ifndef MFEM_THREAD_SAFE
3169 for (
int n = 0; n <
dof; n++)
3184 #ifdef MFEM_THREAD_SAFE
3185 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
3192 for (
int n = 0; n <
dof; n++)
3194 shape(n) = shape1dx(I[n]) * shape1dy(J[n]) * shape1dz(K[n]);
3205 #ifdef MFEM_THREAD_SAFE
3206 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
3207 DenseMatrix dshape1dx(dof1d,1), dshape1dy(dof1d,1), dshape1dz(dof1d,1);
3218 for (
int n = 0; n <
dof; n++)
3220 dshape(n,0) = dshape1dx(I[n],0) * shape1dy(J[n]) * shape1dz(K[n]);
3221 dshape(n,1) = shape1dx(I[n]) * dshape1dy(J[n],0) * shape1dz(K[n]);
3222 dshape(n,2) = shape1dx(I[n]) * shape1dy(J[n]) * dshape1dz(K[n],0);
3251 shape(0) = 1.0 - 2.0 * x;
3258 shape(1) = 2.0 * x - 1.0;
3259 shape(2) = 2.0 - 2.0 * x;
3270 dshape(0,0) = - 2.0;
3278 dshape(2,0) = - 2.0;
3305 L0 = 2.0 * ( 1. - ip.
x - ip.
y );
3306 L1 = 2.0 * ( ip.
x );
3307 L2 = 2.0 * ( ip.
y );
3316 for (i = 0; i < 6; i++)
3323 shape(0) = L0 - 1.0;
3330 shape(1) = L1 - 1.0;
3337 shape(2) = L2 - 1.0;
3341 shape(3) = 1.0 - L2;
3342 shape(4) = 1.0 - L0;
3343 shape(5) = 1.0 - L1;
3353 L0 = 2.0 * ( 1. - ip.
x - ip.
y );
3354 L1 = 2.0 * ( ip.
x );
3355 L2 = 2.0 * ( ip.
y );
3357 double DL0[2], DL1[2], DL2[2];
3358 DL0[0] = -2.0; DL0[1] = -2.0;
3359 DL1[0] = 2.0; DL1[1] = 0.0;
3360 DL2[0] = 0.0; DL2[1] = 2.0;
3362 for (i = 0; i < 6; i++)
3363 for (j = 0; j < 2; j++)
3370 for (j = 0; j < 2; j++)
3372 dshape(0,j) = DL0[j];
3373 dshape(3,j) = DL1[j];
3374 dshape(5,j) = DL2[j];
3379 for (j = 0; j < 2; j++)
3381 dshape(3,j) = DL0[j];
3382 dshape(1,j) = DL1[j];
3383 dshape(4,j) = DL2[j];
3388 for (j = 0; j < 2; j++)
3390 dshape(5,j) = DL0[j];
3391 dshape(4,j) = DL1[j];
3392 dshape(2,j) = DL2[j];
3397 for (j = 0; j < 2; j++)
3399 dshape(3,j) = - DL2[j];
3400 dshape(4,j) = - DL0[j];
3401 dshape(5,j) = - DL1[j];
3446 double L0, L1, L2, L3, L4, L5;
3447 L0 = 2.0 * ( 1. - ip.
x - ip.
y - ip.
z );
3448 L1 = 2.0 * ( ip.
x );
3449 L2 = 2.0 * ( ip.
y );
3450 L3 = 2.0 * ( ip.
z );
3451 L4 = 2.0 * ( ip.
x + ip.
y );
3452 L5 = 2.0 * ( ip.
y + ip.
z );
3465 for (i = 0; i < 10; i++)
3472 shape(0) = L0 - 1.0;
3480 shape(1) = L1 - 1.0;
3488 shape(2) = L2 - 1.0;
3496 shape(3) = L3 - 1.0;
3498 else if ((L4 <= 1.0) && (L5 <= 1.0))
3500 shape(4) = 1.0 - L5;
3502 shape(6) = 1.0 - L4;
3503 shape(8) = 1.0 - L0;
3505 else if ((L4 >= 1.0) && (L5 <= 1.0))
3507 shape(4) = 1.0 - L5;
3508 shape(5) = 1.0 - L1;
3509 shape(7) = L4 - 1.0;
3512 else if ((L4 <= 1.0) && (L5 >= 1.0))
3514 shape(5) = 1.0 - L3;
3515 shape(6) = 1.0 - L4;
3517 shape(9) = L5 - 1.0;
3519 else if ((L4 >= 1.0) && (L5 >= 1.0))
3522 shape(7) = L4 - 1.0;
3523 shape(8) = 1.0 - L2;
3524 shape(9) = L5 - 1.0;
3533 double L0, L1, L2, L3, L4, L5;
3534 L0 = 2.0 * ( 1. - ip.
x - ip.
y - ip.
z );
3535 L1 = 2.0 * ( ip.
x );
3536 L2 = 2.0 * ( ip.
y );
3537 L3 = 2.0 * ( ip.
z );
3538 L4 = 2.0 * ( ip.
x + ip.
y );
3539 L5 = 2.0 * ( ip.
y + ip.
z );
3541 double DL0[3], DL1[3], DL2[3], DL3[3], DL4[3], DL5[3];
3542 DL0[0] = -2.0; DL0[1] = -2.0; DL0[2] = -2.0;
3543 DL1[0] = 2.0; DL1[1] = 0.0; DL1[2] = 0.0;
3544 DL2[0] = 0.0; DL2[1] = 2.0; DL2[2] = 0.0;
3545 DL3[0] = 0.0; DL3[1] = 0.0; DL3[2] = 2.0;
3546 DL4[0] = 2.0; DL4[1] = 2.0; DL4[2] = 0.0;
3547 DL5[0] = 0.0; DL5[1] = 2.0; DL5[2] = 2.0;
3549 for (i = 0; i < 10; i++)
3550 for (j = 0; j < 3; j++)
3557 for (j = 0; j < 3; j++)
3559 dshape(0,j) = DL0[j];
3560 dshape(4,j) = DL1[j];
3561 dshape(5,j) = DL2[j];
3562 dshape(6,j) = DL3[j];
3567 for (j = 0; j < 3; j++)
3569 dshape(4,j) = DL0[j];
3570 dshape(1,j) = DL1[j];
3571 dshape(7,j) = DL2[j];
3572 dshape(8,j) = DL3[j];
3577 for (j = 0; j < 3; j++)
3579 dshape(5,j) = DL0[j];
3580 dshape(7,j) = DL1[j];
3581 dshape(2,j) = DL2[j];
3582 dshape(9,j) = DL3[j];
3587 for (j = 0; j < 3; j++)
3589 dshape(6,j) = DL0[j];
3590 dshape(8,j) = DL1[j];
3591 dshape(9,j) = DL2[j];
3592 dshape(3,j) = DL3[j];
3595 else if ((L4 <= 1.0) && (L5 <= 1.0))
3597 for (j = 0; j < 3; j++)
3599 dshape(4,j) = - DL5[j];
3600 dshape(5,j) = DL2[j];
3601 dshape(6,j) = - DL4[j];
3602 dshape(8,j) = - DL0[j];
3605 else if ((L4 >= 1.0) && (L5 <= 1.0))
3607 for (j = 0; j < 3; j++)
3609 dshape(4,j) = - DL5[j];
3610 dshape(5,j) = - DL1[j];
3611 dshape(7,j) = DL4[j];
3612 dshape(8,j) = DL3[j];
3615 else if ((L4 <= 1.0) && (L5 >= 1.0))
3617 for (j = 0; j < 3; j++)
3619 dshape(5,j) = - DL3[j];
3620 dshape(6,j) = - DL4[j];
3621 dshape(8,j) = DL1[j];
3622 dshape(9,j) = DL5[j];
3625 else if ((L4 >= 1.0) && (L5 >= 1.0))
3627 for (j = 0; j < 3; j++)
3629 dshape(5,j) = DL0[j];
3630 dshape(7,j) = DL4[j];
3631 dshape(8,j) = - DL2[j];
3632 dshape(9,j) = DL5[j];
3665 double x = ip.
x, y = ip.
y;
3667 Lx = 2.0 * ( 1. - x );
3668 Ly = 2.0 * ( 1. - y );
3677 for (i = 0; i < 9; i++)
3682 if ((x <= 0.5) && (y <= 0.5))
3684 shape(0) = (Lx - 1.0) * (Ly - 1.0);
3685 shape(4) = (2.0 - Lx) * (Ly - 1.0);
3686 shape(8) = (2.0 - Lx) * (2.0 - Ly);
3687 shape(7) = (Lx - 1.0) * (2.0 - Ly);
3689 else if ((x >= 0.5) && (y <= 0.5))
3691 shape(4) = Lx * (Ly - 1.0);
3692 shape(1) = (1.0 - Lx) * (Ly - 1.0);
3693 shape(5) = (1.0 - Lx) * (2.0 - Ly);
3694 shape(8) = Lx * (2.0 - Ly);
3696 else if ((x >= 0.5) && (y >= 0.5))
3698 shape(8) = Lx * Ly ;
3699 shape(5) = (1.0 - Lx) * Ly ;
3700 shape(2) = (1.0 - Lx) * (1.0 - Ly);
3701 shape(6) = Lx * (1.0 - Ly);
3703 else if ((x <= 0.5) && (y >= 0.5))
3705 shape(7) = (Lx - 1.0) * Ly ;
3706 shape(8) = (2.0 - Lx) * Ly ;
3707 shape(6) = (2.0 - Lx) * (1.0 - Ly);
3708 shape(3) = (Lx - 1.0) * (1.0 - Ly);
3716 double x = ip.
x, y = ip.
y;
3718 Lx = 2.0 * ( 1. - x );
3719 Ly = 2.0 * ( 1. - y );
3721 for (i = 0; i < 9; i++)
3722 for (j = 0; j < 2; j++)
3727 if ((x <= 0.5) && (y <= 0.5))
3729 dshape(0,0) = 2.0 * (1.0 - Ly);
3730 dshape(0,1) = 2.0 * (1.0 - Lx);
3732 dshape(4,0) = 2.0 * (Ly - 1.0);
3733 dshape(4,1) = -2.0 * (2.0 - Lx);
3735 dshape(8,0) = 2.0 * (2.0 - Ly);
3736 dshape(8,1) = 2.0 * (2.0 - Lx);
3738 dshape(7,0) = -2.0 * (2.0 - Ly);
3739 dshape(7,0) = 2.0 * (Lx - 1.0);
3741 else if ((x >= 0.5) && (y <= 0.5))
3743 dshape(4,0) = -2.0 * (Ly - 1.0);
3744 dshape(4,1) = -2.0 * Lx;
3746 dshape(1,0) = 2.0 * (Ly - 1.0);
3747 dshape(1,1) = -2.0 * (1.0 - Lx);
3749 dshape(5,0) = 2.0 * (2.0 - Ly);
3750 dshape(5,1) = 2.0 * (1.0 - Lx);
3752 dshape(8,0) = -2.0 * (2.0 - Ly);
3753 dshape(8,1) = 2.0 * Lx;
3755 else if ((x >= 0.5) && (y >= 0.5))
3757 dshape(8,0) = -2.0 * Ly;
3758 dshape(8,1) = -2.0 * Lx;
3760 dshape(5,0) = 2.0 * Ly;
3761 dshape(5,1) = -2.0 * (1.0 - Lx);
3763 dshape(2,0) = 2.0 * (1.0 - Ly);
3764 dshape(2,1) = 2.0 * (1.0 - Lx);
3766 dshape(6,0) = -2.0 * (1.0 - Ly);
3767 dshape(6,1) = 2.0 * Lx;
3769 else if ((x <= 0.5) && (y >= 0.5))
3771 dshape(7,0) = -2.0 * Ly;
3772 dshape(7,1) = -2.0 * (Lx - 1.0);
3774 dshape(8,0) = 2.0 * Ly ;
3775 dshape(8,1) = -2.0 * (2.0 - Lx);
3777 dshape(6,0) = 2.0 * (1.0 - Ly);
3778 dshape(6,1) = 2.0 * (2.0 - Lx);
3780 dshape(3,0) = -2.0 * (1.0 - Ly);
3781 dshape(3,1) = 2.0 * (Lx - 1.0);
3792 I[ 0] = 0.0; J[ 0] = 0.0; K[ 0] = 0.0;
3793 I[ 1] = 1.0; J[ 1] = 0.0; K[ 1] = 0.0;
3794 I[ 2] = 1.0; J[ 2] = 1.0; K[ 2] = 0.0;
3795 I[ 3] = 0.0; J[ 3] = 1.0; K[ 3] = 0.0;
3796 I[ 4] = 0.0; J[ 4] = 0.0; K[ 4] = 1.0;
3797 I[ 5] = 1.0; J[ 5] = 0.0; K[ 5] = 1.0;
3798 I[ 6] = 1.0; J[ 6] = 1.0; K[ 6] = 1.0;
3799 I[ 7] = 0.0; J[ 7] = 1.0; K[ 7] = 1.0;
3801 I[ 8] = 0.5; J[ 8] = 0.0; K[ 8] = 0.0;
3802 I[ 9] = 1.0; J[ 9] = 0.5; K[ 9] = 0.0;
3803 I[10] = 0.5; J[10] = 1.0; K[10] = 0.0;
3804 I[11] = 0.0; J[11] = 0.5; K[11] = 0.0;
3805 I[12] = 0.5; J[12] = 0.0; K[12] = 1.0;
3806 I[13] = 1.0; J[13] = 0.5; K[13] = 1.0;
3807 I[14] = 0.5; J[14] = 1.0; K[14] = 1.0;
3808 I[15] = 0.0; J[15] = 0.5; K[15] = 1.0;
3809 I[16] = 0.0; J[16] = 0.0; K[16] = 0.5;
3810 I[17] = 1.0; J[17] = 0.0; K[17] = 0.5;
3811 I[18] = 1.0; J[18] = 1.0; K[18] = 0.5;
3812 I[19] = 0.0; J[19] = 1.0; K[19] = 0.5;
3814 I[20] = 0.5; J[20] = 0.5; K[20] = 0.0;
3815 I[21] = 0.5; J[21] = 0.0; K[21] = 0.5;
3816 I[22] = 1.0; J[22] = 0.5; K[22] = 0.5;
3817 I[23] = 0.5; J[23] = 1.0; K[23] = 0.5;
3818 I[24] = 0.0; J[24] = 0.5; K[24] = 0.5;
3819 I[25] = 0.5; J[25] = 0.5; K[25] = 1.0;
3821 I[26] = 0.5; J[26] = 0.5; K[26] = 0.5;
3823 for (
int n = 0; n < 27; n++)
3836 double x = ip.
x, y = ip.
y, z = ip.
z;
3838 for (i = 0; i < 27; i++)
3843 if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5))
3858 else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5))
3873 else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5))
3888 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5))
3903 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5))
3918 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5))
3933 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5))
3964 shape(N[0]) = Lx * Ly * Lz;
3965 shape(N[1]) = (1 - Lx) * Ly * Lz;
3966 shape(N[2]) = (1 - Lx) * (1 - Ly) * Lz;
3967 shape(N[3]) = Lx * (1 - Ly) * Lz;
3968 shape(N[4]) = Lx * Ly * (1 - Lz);
3969 shape(N[5]) = (1 - Lx) * Ly * (1 - Lz);
3970 shape(N[6]) = (1 - Lx) * (1 - Ly) * (1 - Lz);
3971 shape(N[7]) = Lx * (1 - Ly) * (1 - Lz);
3979 double x = ip.
x, y = ip.
y, z = ip.
z;
3981 for (i = 0; i < 27; i++)
3982 for (j = 0; j < 3; j++)
3987 if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5))
4002 else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5))
4017 else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5))
4032 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5))
4047 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5))
4062 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5))
4077 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5))
4108 dshape(N[0],0) = -2.0 * Ly * Lz ;
4109 dshape(N[0],1) = -2.0 * Lx * Lz ;
4110 dshape(N[0],2) = -2.0 * Lx * Ly ;
4112 dshape(N[1],0) = 2.0 * Ly * Lz ;
4113 dshape(N[1],1) = -2.0 * (1 - Lx) * Lz ;
4114 dshape(N[1],2) = -2.0 * (1 - Lx) * Ly ;
4116 dshape(N[2],0) = 2.0 * (1 - Ly) * Lz ;
4117 dshape(N[2],1) = 2.0 * (1 - Lx) * Lz ;
4118 dshape(N[2],2) = -2.0 * (1 - Lx) * (1 - Ly);
4120 dshape(N[3],0) = -2.0 * (1 - Ly) * Lz ;
4121 dshape(N[3],1) = 2.0 * Lx * Lz ;
4122 dshape(N[3],2) = -2.0 * Lx * (1 - Ly);
4124 dshape(N[4],0) = -2.0 * Ly * (1 - Lz);
4125 dshape(N[4],1) = -2.0 * Lx * (1 - Lz);
4126 dshape(N[4],2) = 2.0 * Lx * Ly ;
4128 dshape(N[5],0) = 2.0 * Ly * (1 - Lz);
4129 dshape(N[5],1) = -2.0 * (1 - Lx) * (1 - Lz);
4130 dshape(N[5],2) = 2.0 * (1 - Lx) * Ly ;
4132 dshape(N[6],0) = 2.0 * (1 - Ly) * (1 - Lz);
4133 dshape(N[6],1) = 2.0 * (1 - Lx) * (1 - Lz);
4134 dshape(N[6],2) = 2.0 * (1 - Lx) * (1 - Ly);
4136 dshape(N[7],0) = -2.0 * (1 - Ly) * (1 - Lz);
4137 dshape(N[7],1) = 2.0 * Lx * (1 - Lz);
4138 dshape(N[7],2) = 2.0 * Lx * (1 - Ly);
4199 double x = ip.
x, y = ip.
y, z = ip.
z;
4201 shape(0,0) = (1. - y) * (1. - z);
4205 shape(2,0) = y * (1. - z);
4209 shape(4,0) = z * (1. - y);
4218 shape(1,1) = x * (1. - z);
4222 shape(3,1) = (1. - x) * (1. - z);
4230 shape(7,1) = (1. - x) * z;
4235 shape(8,2) = (1. - x) * (1. - y);
4239 shape(9,2) = x * (1. - y);
4243 shape(10,2) = x * y;
4247 shape(11,2) = y * (1. - x);
4255 double x = ip.
x, y = ip.
y, z = ip.
z;
4257 curl_shape(0,0) = 0.;
4258 curl_shape(0,1) = y - 1.;
4259 curl_shape(0,2) = 1. - z;
4261 curl_shape(2,0) = 0.;
4262 curl_shape(2,1) = -y;
4263 curl_shape(2,2) = z - 1.;
4265 curl_shape(4,0) = 0;
4266 curl_shape(4,1) = 1. - y;
4267 curl_shape(4,2) = z;
4269 curl_shape(6,0) = 0.;
4270 curl_shape(6,1) = y;
4271 curl_shape(6,2) = -z;
4273 curl_shape(1,0) = x;
4274 curl_shape(1,1) = 0.;
4275 curl_shape(1,2) = 1. - z;
4277 curl_shape(3,0) = 1. - x;
4278 curl_shape(3,1) = 0.;
4279 curl_shape(3,2) = z - 1.;
4281 curl_shape(5,0) = -x;
4282 curl_shape(5,1) = 0.;
4283 curl_shape(5,2) = z;
4285 curl_shape(7,0) = x - 1.;
4286 curl_shape(7,1) = 0.;
4287 curl_shape(7,2) = -z;
4289 curl_shape(8,0) = x - 1.;
4290 curl_shape(8,1) = 1. - y;
4291 curl_shape(8,2) = 0.;
4293 curl_shape(9,0) = -x;
4294 curl_shape(9,1) = y - 1.;
4295 curl_shape(9,2) = 0;
4297 curl_shape(10,0) = x;
4298 curl_shape(10,1) = -y;
4299 curl_shape(10,2) = 0.;
4301 curl_shape(11,0) = 1. - x;
4302 curl_shape(11,1) = y;
4303 curl_shape(11,2) = 0.;
4306 const double Nedelec1HexFiniteElement::tk[12][3] =
4308 {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
4309 {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
4310 {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}
4317 #ifdef MFEM_THREAD_SAFE
4322 for (k = 0; k <
dof; k++)
4325 for (j = 0; j <
dof; j++)
4327 double d = (
vshape(j,0)*tk[k][0] +
vshape(j,1)*tk[k][1] +
4329 if (j == k) { d -= 1.0; }
4330 if (fabs(d) > 1.0e-12)
4332 mfem::err <<
"Nedelec1HexFiniteElement::GetLocalInterpolation (...)\n"
4333 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
4341 ip.
x = ip.
y = ip.
z = 0.0;
4348 for (k = 0; k <
dof; k++)
4351 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
4354 vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
4355 vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
4356 vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
4357 for (j = 0; j <
dof; j++)
4358 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
4359 vshape(j,2)*vk[2])) < 1.0e-12)
4373 for (
int k = 0; k <
dof; k++)
4381 vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
4382 vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
4383 vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
4395 for (
int k = 0; k <
dof; k++)
4398 dshape.Mult(tk[k], grad_k);
4399 for (
int j = 0; j < grad_k.Size(); j++)
4401 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
4439 double x = ip.
x, y = ip.
y, z = ip.
z;
4441 shape(0,0) = 1. - y - z;
4446 shape(1,1) = 1. - x - z;
4451 shape(2,2) = 1. - x - y;
4470 curl_shape(0,0) = 0.;
4471 curl_shape(0,1) = -2.;
4472 curl_shape(0,2) = 2.;
4474 curl_shape(1,0) = 2.;
4475 curl_shape(1,1) = 0.;
4476 curl_shape(1,2) = -2.;
4478 curl_shape(2,0) = -2.;
4479 curl_shape(2,1) = 2.;
4480 curl_shape(2,2) = 0.;
4482 curl_shape(3,0) = 0.;
4483 curl_shape(3,1) = 0.;
4484 curl_shape(3,2) = 2.;
4486 curl_shape(4,0) = 0.;
4487 curl_shape(4,1) = -2.;
4488 curl_shape(4,2) = 0.;
4490 curl_shape(5,0) = 2.;
4491 curl_shape(5,1) = 0.;
4492 curl_shape(5,2) = 0.;
4495 const double Nedelec1TetFiniteElement::tk[6][3] =
4496 {{1,0,0}, {0,1,0}, {0,0,1}, {-1,1,0}, {-1,0,1}, {0,-1,1}};
4502 #ifdef MFEM_THREAD_SAFE
4507 for (k = 0; k <
dof; k++)
4510 for (j = 0; j <
dof; j++)
4512 double d = (
vshape(j,0)*tk[k][0] +
vshape(j,1)*tk[k][1] +
4514 if (j == k) { d -= 1.0; }
4515 if (fabs(d) > 1.0e-12)
4517 mfem::err <<
"Nedelec1TetFiniteElement::GetLocalInterpolation (...)\n"
4518 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
4526 ip.
x = ip.
y = ip.
z = 0.0;
4533 for (k = 0; k <
dof; k++)
4536 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
4539 vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
4540 vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
4541 vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
4542 for (j = 0; j <
dof; j++)
4543 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
4544 vshape(j,2)*vk[2])) < 1.0e-12)
4558 for (
int k = 0; k <
dof; k++)
4566 vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
4567 vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
4568 vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
4580 for (
int k = 0; k <
dof; k++)
4583 dshape.Mult(tk[k], grad_k);
4584 for (
int j = 0; j < grad_k.Size(); j++)
4586 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
4636 double x = ip.
x, y = ip.
y, z = ip.
z;
4638 shape(0,0) = (1. - y) * (1. - z);
4639 shape(0,1) = x * (1. - z);
4642 shape(1,0) = - y * (1. - z);
4643 shape(1,1) = x * (1. - z);
4646 shape(2,0) = - y * (1. - z);
4647 shape(2,1) = - (1. - x) * (1. - z);
4650 shape(3,0) = (1. - y) * z;
4654 shape(4,0) = - y * z;
4658 shape(5,0) = - y * z;
4659 shape(5,1) = - (1. - x) * z;
4664 shape(6,2) = 1. - x - y;
4679 double x = ip.
x, y = ip.
y, z2 = 2. * ip.
z;
4681 curl_shape(0,0) = x;
4682 curl_shape(0,1) = - 1. + y;
4683 curl_shape(0,2) = 2. - z2;
4685 curl_shape(1,0) = x;
4686 curl_shape(1,1) = y;
4687 curl_shape(1,2) = 2. - z2;
4689 curl_shape(2,0) = - 1. + x;
4690 curl_shape(2,1) = y;
4691 curl_shape(2,2) = 2. - z2;
4693 curl_shape(3,0) = - x;
4694 curl_shape(3,1) = 1. - y;
4695 curl_shape(3,2) = z2;
4697 curl_shape(4,0) = - x;
4698 curl_shape(4,1) = - y;
4699 curl_shape(4,2) = z2;
4701 curl_shape(5,0) = 1. - x;
4702 curl_shape(5,1) = - y;
4703 curl_shape(5,2) = z2;
4705 curl_shape(6,0) = - 1.;
4706 curl_shape(6,1) = 1.;
4707 curl_shape(6,2) = 0.;
4709 curl_shape(7,0) = 0.;
4710 curl_shape(7,1) = - 1.;
4711 curl_shape(7,2) = 0.;
4713 curl_shape(8,0) = 1.;
4714 curl_shape(8,1) = 0.;
4715 curl_shape(8,2) = 0.;
4718 const double Nedelec1WdgFiniteElement::tk[9][3] =
4720 {1,0,0}, {-1,1,0}, {0,-1,0}, {1,0,0}, {-1,1,0}, {0,-1,0},
4721 {0,0,1}, {0,0,1}, {0,0,1}
4728 #ifdef MFEM_THREAD_SAFE
4733 for (k = 0; k <
dof; k++)
4736 for (j = 0; j <
dof; j++)
4738 double d = (
vshape(j,0)*tk[k][0] +
vshape(j,1)*tk[k][1] +
4740 if (j == k) { d -= 1.0; }
4741 if (fabs(d) > 1.0e-12)
4743 mfem::err <<
"Nedelec1WdgFiniteElement::GetLocalInterpolation (...)\n"
4744 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
4752 ip.
x = ip.
y = ip.
z = 0.0;
4759 for (k = 0; k <
dof; k++)
4762 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
4765 vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
4766 vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
4767 vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
4768 for (j = 0; j <
dof; j++)
4769 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
4770 vshape(j,2)*vk[2])) < 1.0e-12)
4784 for (
int k = 0; k <
dof; k++)
4792 vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
4793 vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
4794 vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
4806 for (
int k = 0; k <
dof; k++)
4809 dshape.Mult(tk[k], grad_k);
4810 for (
int j = 0; j < grad_k.Size(); j++)
4812 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
4858 double x = ip.
x, y = ip.
y, z = ip.
z, z2 = 2. * ip.
z;
4859 double ox = 1. - x - z, oy = 1. - y - z, oz = 1. - z;
4903 double ozi = 1.0 / oz;
4907 shape(0,2) = x * oy * ozi;
4911 shape(1,2) = x * y * ozi;
4915 shape(2,2) = x * y * ozi;
4919 shape(3,2) = ox * y * ozi;
4921 shape(4,0) = oy * z * ozi;
4922 shape(4,1) = ox * z * ozi;
4923 shape(4,2) = 1. - x - y + x * y * (1. - z2) * ozi * ozi;
4925 shape(5,0) = - oy * z * ozi;
4926 shape(5,1) = x * z * ozi;
4927 shape(5,2) = x * (1. - y * (1. - z2) * ozi * ozi);
4929 shape(6,0) = - y * z * ozi;
4930 shape(6,1) = - x * z * ozi;
4931 shape(6,2) = x * y * (1. - z2) * ozi * ozi;
4933 shape(7,0) = y * z * ozi;
4934 shape(7,1) = - ox * z * ozi;
4935 shape(7,2) = y * (1. - x * (1. - z2) * ozi * ozi);
4942 double x = ip.
x, y = ip.
y, z = ip.
z, z2 = 2. * z;
4943 double ox = 1. - x - z, oy = 1. - y - z, oz = 1. - z;
4953 curl_shape(0,0) = 0.;
4954 curl_shape(0,1) = - 2.;
4955 curl_shape(0,2) = 1.;
4957 curl_shape(1,0) = 0.;
4958 curl_shape(1,1) = 0.;
4959 curl_shape(1,2) = 1.;
4961 curl_shape(2,0) = 0.;
4962 curl_shape(2,1) = 0.;
4963 curl_shape(2,2) = - 1.;
4965 curl_shape(3,0) = 2.;
4966 curl_shape(3,1) = 0.;
4967 curl_shape(3,2) = - 1.;
4969 curl_shape(4,0) = - 2.;
4970 curl_shape(4,1) = 2.;
4971 curl_shape(4,2) = 0.;
4973 curl_shape(5,0) = 0.;
4974 curl_shape(5,1) = - 2.;
4975 curl_shape(5,2) = 0.;
4977 curl_shape(6,0) = 0.;
4978 curl_shape(6,1) = 0.;
4979 curl_shape(6,2) = 0.;
4981 curl_shape(7,0) = 2.;
4982 curl_shape(7,1) = 0.;
4983 curl_shape(7,2) = 0.;
4988 double ozi = 1. / oz;
4990 curl_shape(0,0) = - x * ozi;
4991 curl_shape(0,1) = - 2. + y * ozi;
4992 curl_shape(0,2) = 1.;
4994 curl_shape(1,0) = x * ozi;
4995 curl_shape(1,1) = - y * ozi;
4996 curl_shape(1,2) = 1.;
4998 curl_shape(2,0) = x * ozi;
4999 curl_shape(2,1) = - y * ozi;
5000 curl_shape(2,2) = - 1.;
5002 curl_shape(3,0) = (2. - x - z2) * ozi;
5003 curl_shape(3,1) = y * ozi;
5004 curl_shape(3,2) = - 1.;
5006 curl_shape(4,0) = - 2. * ox * ozi;
5007 curl_shape(4,1) = 2. * oy * ozi;
5008 curl_shape(4,2) = 0.;
5010 curl_shape(5,0) = - 2. * x * ozi;
5011 curl_shape(5,1) = - 2. * oy * ozi;
5012 curl_shape(5,2) = 0.;
5014 curl_shape(6,0) = 2. * x * ozi;
5015 curl_shape(6,1) = - 2. * y * ozi;
5016 curl_shape(6,2) = 0.;
5018 curl_shape(7,0) = 2. * ox * ozi;
5019 curl_shape(7,1) = 2. * y * ozi;
5020 curl_shape(7,2) = 0.;
5023 const double Nedelec1PyrFiniteElement::tk[8][3] =
5024 {{1,0,0}, {0,1,0}, {1,0,0}, {0,1,0}, {0,0,1}, {-1,0,1}, {-1,-1,1}, {0,-1,1}};
5030 #ifdef MFEM_THREAD_SAFE
5035 for (k = 0; k <
dof; k++)
5038 for (j = 0; j <
dof; j++)
5040 double d = (
vshape(j,0)*tk[k][0] +
vshape(j,1)*tk[k][1] +
5042 if (j == k) { d -= 1.0; }
5043 if (fabs(d) > 1.0e-12)
5045 mfem::err <<
"Nedelec1PyrFiniteElement::GetLocalInterpolation (...)\n"
5046 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5054 ip.
x = ip.
y = ip.
z = 0.0;
5061 for (k = 0; k <
dof; k++)
5064 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5067 vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
5068 vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
5069 vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
5070 for (j = 0; j <
dof; j++)
5071 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5072 vshape(j,2)*vk[2])) < 1.0e-12)
5086 for (
int k = 0; k <
dof; k++)
5094 vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
5095 vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
5096 vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
5108 for (
int k = 0; k <
dof; k++)
5111 dshape.Mult(tk[k], grad_k);
5112 for (
int j = 0; j < grad_k.Size(); j++)
5114 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
5153 double x = ip.
x, y = ip.
y, z = ip.
z;
5157 shape(0,2) = z - 1.;
5160 shape(1,1) = y - 1.;
5171 shape(4,0) = x - 1.;
5191 const double RT0HexFiniteElement::nk[6][3] =
5192 {{0,0,-1}, {0,-1,0}, {1,0,0}, {0,1,0}, {-1,0,0}, {0,0,1}};
5198 #ifdef MFEM_THREAD_SAFE
5203 for (k = 0; k < 6; k++)
5206 for (j = 0; j < 6; j++)
5208 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5210 if (j == k) { d -= 1.0; }
5211 if (fabs(d) > 1.0e-12)
5213 mfem::err <<
"RT0HexFiniteElement::GetLocalInterpolation (...)\n"
5214 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5222 ip.
x = ip.
y = ip.
z = 0.0;
5231 for (k = 0; k < 6; k++)
5234 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5237 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1]+Jinv(0,2)*nk[k][2];
5238 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1]+Jinv(1,2)*nk[k][2];
5239 vk[2] = Jinv(2,0)*nk[k][0]+Jinv(2,1)*nk[k][1]+Jinv(2,2)*nk[k][2];
5240 for (j = 0; j < 6; j++)
5241 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5242 vshape(j,2)*vk[2])) < 1.0e-12)
5256 for (
int k = 0; k < 6; k++)
5265 vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1]+Jinv(0,2)*nk[k][2] ) +
5266 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1]+Jinv(1,2)*nk[k][2] ) +
5267 vk[2] * ( Jinv(2,0)*nk[k][0]+Jinv(2,1)*nk[k][1]+Jinv(2,2)*nk[k][2] );
5397 double x = ip.
x, y = ip.
y, z = ip.
z;
5401 shape(2,2) = -(1. - 3.*z + 2.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5404 shape(3,2) = -(1. - 3.*z + 2.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5407 shape(0,2) = -(1. - 3.*z + 2.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5410 shape(1,2) = -(1. - 3.*z + 2.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5413 shape(4,1) = -(1. - 3.*y + 2.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5416 shape(5,1) = -(1. - 3.*y + 2.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5419 shape(6,1) = -(1. - 3.*y + 2.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5422 shape(7,1) = -(1. - 3.*y + 2.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5425 shape(8,0) = (-x + 2.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5428 shape(9,0) = (-x + 2.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5431 shape(10,0) = (-x + 2.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5434 shape(11,0) = (-x + 2.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5439 shape(13,1) = (-y + 2.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5442 shape(12,1) = (-y + 2.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5445 shape(15,1) = (-y + 2.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5448 shape(14,1) = (-y + 2.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5451 shape(17,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5454 shape(16,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5457 shape(19,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5460 shape(18,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5466 shape(20,2) = (-z + 2.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5469 shape(21,2) = (-z + 2.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5472 shape(22,2) = (-z + 2.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5475 shape(23,2) = (-z + 2.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5477 shape(24,0) = (4.*x - 4.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5480 shape(25,0) = (4.*x - 4.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5483 shape(26,0) = (4.*x - 4.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5486 shape(27,0) = (4.*x - 4.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5491 shape(28,1) = (4.*y - 4.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5494 shape(29,1) = (4.*y - 4.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5497 shape(30,1) = (4.*y - 4.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5500 shape(31,1) = (4.*y - 4.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5505 shape(32,2) = (4.*z - 4.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5508 shape(33,2) = (4.*z - 4.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5511 shape(34,2) = (4.*z - 4.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5514 shape(35,2) = (4.*z - 4.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5520 double x = ip.
x, y = ip.
y, z = ip.
z;
5522 divshape(2) = -(-3. + 4.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5523 divshape(3) = -(-3. + 4.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5524 divshape(0) = -(-3. + 4.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5525 divshape(1) = -(-3. + 4.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5527 divshape(4) = -(-3. + 4.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5528 divshape(5) = -(-3. + 4.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5529 divshape(6) = -(-3. + 4.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5530 divshape(7) = -(-3. + 4.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5532 divshape(8) = (-1. + 4.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5533 divshape(9) = (-1. + 4.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5534 divshape(10) = (-1. + 4.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5535 divshape(11) = (-1. + 4.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5537 divshape(13) = (-1. + 4.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5538 divshape(12) = (-1. + 4.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5539 divshape(15) = (-1. + 4.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5540 divshape(14) = (-1. + 4.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5542 divshape(17) = -(-3. + 4.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5543 divshape(16) = -(-3. + 4.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5544 divshape(19) = -(-3. + 4.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5545 divshape(18) = -(-3. + 4.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5547 divshape(20) = (-1. + 4.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5548 divshape(21) = (-1. + 4.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5549 divshape(22) = (-1. + 4.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5550 divshape(23) = (-1. + 4.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5552 divshape(24) = ( 4. - 8.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5553 divshape(25) = ( 4. - 8.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5554 divshape(26) = ( 4. - 8.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5555 divshape(27) = ( 4. - 8.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5557 divshape(28) = ( 4. - 8.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5558 divshape(29) = ( 4. - 8.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5559 divshape(30) = ( 4. - 8.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5560 divshape(31) = ( 4. - 8.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5562 divshape(32) = ( 4. - 8.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5563 divshape(33) = ( 4. - 8.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5564 divshape(34) = ( 4. - 8.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5565 divshape(35) = ( 4. - 8.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5568 const double RT1HexFiniteElement::nk[36][3] =
5570 {0, 0,-1}, {0, 0,-1}, {0, 0,-1}, {0, 0,-1},
5571 {0,-1, 0}, {0,-1, 0}, {0,-1, 0}, {0,-1, 0},
5572 {1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0},
5573 {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0},
5574 {-1,0, 0}, {-1,0, 0}, {-1,0, 0}, {-1,0, 0},
5575 {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1},
5576 {1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0},
5577 {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0},
5578 {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1}
5585 #ifdef MFEM_THREAD_SAFE
5590 for (k = 0; k < 36; k++)
5593 for (j = 0; j < 36; j++)
5595 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5597 if (j == k) { d -= 1.0; }
5598 if (fabs(d) > 1.0e-12)
5600 mfem::err <<
"RT0HexFiniteElement::GetLocalInterpolation (...)\n"
5601 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5609 ip.
x = ip.
y = ip.
z = 0.0;
5618 for (k = 0; k < 36; k++)
5621 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5624 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1]+Jinv(0,2)*nk[k][2];
5625 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1]+Jinv(1,2)*nk[k][2];
5626 vk[2] = Jinv(2,0)*nk[k][0]+Jinv(2,1)*nk[k][1]+Jinv(2,2)*nk[k][2];
5627 for (j = 0; j < 36; j++)
5628 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5629 vshape(j,2)*vk[2])) < 1.0e-12)
5643 for (
int k = 0; k < 36; k++)
5652 vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1]+Jinv(0,2)*nk[k][2] ) +
5653 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1]+Jinv(1,2)*nk[k][2] ) +
5654 vk[2] * ( Jinv(2,0)*nk[k][0]+Jinv(2,1)*nk[k][1]+Jinv(2,2)*nk[k][2] );
5682 double x2 = 2.0*ip.
x, y2 = 2.0*ip.
y, z2 = 2.0*ip.
z;
5688 shape(1,0) = x2 - 2.0;
5693 shape(2,1) = y2 - 2.0;
5698 shape(3,2) = z2 - 2.0;
5710 const double RT0TetFiniteElement::nk[4][3] =
5711 {{.5,.5,.5}, {-.5,0,0}, {0,-.5,0}, {0,0,-.5}};
5717 #ifdef MFEM_THREAD_SAFE
5722 for (k = 0; k < 4; k++)
5725 for (j = 0; j < 4; j++)
5727 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5729 if (j == k) { d -= 1.0; }
5730 if (fabs(d) > 1.0e-12)
5732 mfem::err <<
"RT0TetFiniteElement::GetLocalInterpolation (...)\n"
5733 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5741 ip.
x = ip.
y = ip.
z = 0.0;
5750 for (k = 0; k < 4; k++)
5753 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5756 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1]+Jinv(0,2)*nk[k][2];
5757 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1]+Jinv(1,2)*nk[k][2];
5758 vk[2] = Jinv(2,0)*nk[k][0]+Jinv(2,1)*nk[k][1]+Jinv(2,2)*nk[k][2];
5759 for (j = 0; j < 4; j++)
5760 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5761 vshape(j,2)*vk[2])) < 1.0e-12)
5775 for (
int k = 0; k < 4; k++)
5784 vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1]+Jinv(0,2)*nk[k][2] ) +
5785 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1]+Jinv(1,2)*nk[k][2] ) +
5786 vk[2] * ( Jinv(2,0)*nk[k][0]+Jinv(2,1)*nk[k][1]+Jinv(2,2)*nk[k][2] );
5818 double x = ip.
x, y = ip.
y, z2 = 2.0*ip.
z;
5822 shape(0,2) = z2 - 2.0;
5829 shape(2,1) = y - 1.0;
5836 shape(4,0) = x - 1.0;
5851 const double RT0WdgFiniteElement::nk[5][3] =
5852 {{0.,0.,-.5}, {0.,0.,.5}, {0,-1.,0}, {1.,1.,0}, {-1.,0,0}};
5858 #ifdef MFEM_THREAD_SAFE
5863 for (k = 0; k <
dof; k++)
5866 for (j = 0; j <
dof; j++)
5868 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5870 if (j == k) { d -= 1.0; }
5871 if (fabs(d) > 1.0e-12)
5873 mfem::err <<
"RT0WdgFiniteElement::GetLocalInterpolation (...)\n"
5874 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5882 ip.
x = ip.
y = ip.
z = 0.0;
5891 for (k = 0; k <
dof; k++)
5894 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5897 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1]+Jinv(0,2)*nk[k][2];
5898 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1]+Jinv(1,2)*nk[k][2];
5899 vk[2] = Jinv(2,0)*nk[k][0]+Jinv(2,1)*nk[k][1]+Jinv(2,2)*nk[k][2];
5900 for (j = 0; j <
dof; j++)
5901 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5902 vshape(j,2)*vk[2])) < 1.0e-12)
5916 for (
int k = 0; k < 5; k++)
5925 vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1]+Jinv(0,2)*nk[k][2] ) +
5926 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1]+Jinv(1,2)*nk[k][2] ) +
5927 vk[2] * ( Jinv(2,0)*nk[k][0]+Jinv(2,1)*nk[k][1]+Jinv(2,2)*nk[k][2] );
5939 for (
int k = 0; k <
dof; k++)
5942 curl_shape.Mult(nk[k], curl_k);
5943 for (
int j = 0; j < curl_k.Size(); j++)
5945 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
5978 double x = ip.
x, y = ip.
y, z = ip.
z, oz = 1.0 - z;
5979 double x2 = 2.0*ip.
x, y2 = 2.0*ip.
y, z2 = 2.0*ip.
z;
6013 for (
int i=1; i<5; i++)
6014 for (
int j=0; j<3; j++)
6023 double ozi = 1.0 / oz;
6027 shape(0,2) = z - 1.;
6029 shape(1,0) = - x * z * ozi;
6030 shape(1,1) = (y2 + z2 - y * z - 2.0) * ozi;
6033 shape(2,0) = x * (2.0 - z) * ozi;
6034 shape(2,1) = - y * z * ozi;;
6037 shape(3,0) = - x * z * ozi;
6038 shape(3,1) = y * (2.0 - z) * ozi;
6041 shape(4,0) = (x2 + z2 - x * z - 2.0) * ozi;
6042 shape(4,1) = - y * z * ozi;
6047 for (
int i=1; i<5; i++)
6048 for (
int j=0; j<3; j++)
6066 for (
int i=1; i<5; i++)
6073 const double RT0PyrFiniteElement::nk[5][3] =
6074 {{0.,0.,-1.}, {0,-.5,0}, {.5,0,.5}, {0,.5,.5}, {-.5,0,0}};
6080 #ifdef MFEM_THREAD_SAFE
6085 for (k = 0; k <
dof; k++)
6088 for (j = 0; j <
dof; j++)
6090 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
6092 if (j == k) { d -= 1.0; }
6093 if (fabs(d) > 1.0e-12)
6095 mfem::err <<
"RT0PyrFiniteElement::GetLocalInterpolation (...)\n"
6096 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
6104 ip.
x = ip.
y = ip.
z = 0.0;
6113 for (k = 0; k <
dof; k++)
6116 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
6119 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1]+Jinv(0,2)*nk[k][2];
6120 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1]+Jinv(1,2)*nk[k][2];
6121 vk[2] = Jinv(2,0)*nk[k][0]+Jinv(2,1)*nk[k][1]+Jinv(2,2)*nk[k][2];
6122 for (j = 0; j <
dof; j++)
6123 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
6124 vshape(j,2)*vk[2])) < 1.0e-12)
6138 for (
int k = 0; k <
dof; k++)
6147 vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1]+Jinv(0,2)*nk[k][2] ) +
6148 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1]+Jinv(1,2)*nk[k][2] ) +
6149 vk[2] * ( Jinv(2,0)*nk[k][0]+Jinv(2,1)*nk[k][1]+Jinv(2,2)*nk[k][2] );
6150 if (!rt0 && k > 0) { dofs[k] *= 2.0; }
6162 for (
int k = 0; k <
dof; k++)
6165 curl_shape.Mult(nk[k], curl_k);
6166 if (!rt0 && k > 0) { curl_k *= 2.0; }
6167 for (
int j = 0; j < curl_k.Size(); j++)
6169 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
6205 double x = 2. * ip.
x - 1.;
6206 double y = 2. * ip.
y - 1.;
6207 double z = 2. * ip.
z - 1.;
6208 double f5 = x * x - y * y;
6209 double f6 = y * y - z * z;
6211 shape(0) = (1./6.) * (1. - 3. * z - f5 - 2. * f6);
6212 shape(1) = (1./6.) * (1. - 3. * y - f5 + f6);
6213 shape(2) = (1./6.) * (1. + 3. * x + 2. * f5 + f6);
6214 shape(3) = (1./6.) * (1. + 3. * y - f5 + f6);
6215 shape(4) = (1./6.) * (1. - 3. * x + 2. * f5 + f6);
6216 shape(5) = (1./6.) * (1. + 3. * z - f5 - 2. * f6);
6222 const double a = 2./3.;
6224 double xt = a * (1. - 2. * ip.
x);
6225 double yt = a * (1. - 2. * ip.
y);
6226 double zt = a * (1. - 2. * ip.
z);
6230 dshape(0,2) = -1. - 2. * zt;
6233 dshape(1,1) = -1. - 2. * yt;
6236 dshape(2,0) = 1. - 2. * xt;
6241 dshape(3,1) = 1. - 2. * yt;
6244 dshape(4,0) = -1. - 2. * xt;
6250 dshape(5,2) = 1. - 2. * zt;
Abstract class for all finite elements.
RefinedLinear3DFiniteElement()
Construct the RefinedLinear3DFiniteElement.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
RT0TriangleFiniteElement()
Construct the RT0TriangleFiniteElement.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual function which evaluates the values of all shape functions at a given point ip and stores the...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
RT2TriangleFiniteElement()
Construct the RT2TriangleFiniteElement.
GaussQuad2DFiniteElement()
Construct the GaussQuad2DFiniteElement.
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Array< int > lex_ordering
CrouzeixRaviartQuadFiniteElement()
Construct the CrouzeixRaviartQuadFiniteElement.
P2SegmentFiniteElement()
Construct the P2SegmentFiniteElement.
Linear3DFiniteElement()
Construct the Linear3DFiniteElement.
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
RT1TriangleFiniteElement()
Construct the RT1TriangleFiniteElement.
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
RefinedBiLinear2DFiniteElement()
Construct the RefinedBiLinear2DFiniteElement.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Base class for vector Coefficients that optionally depend on time and space.
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Compute the discrete gradient matrix from the given FiniteElement onto 'this' FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual function which evaluates the values of all shape functions at a given point ip and stores the...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void ProjectCurl(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Compute the discrete curl matrix from the given FiniteElement onto 'this' FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
virtual void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Compute the discrete gradient matrix from the given FiniteElement onto 'this' FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
PointFiniteElement()
Construct the PointFiniteElement.
P0WdgFiniteElement()
Construct the P0WdgFiniteElement.
void SetSize(int s)
Resize the vector to size s.
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &h) const
Compute the Hessian of second order partial derivatives at ip.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
RT0WdgFiniteElement()
Construct the RT0WdgFiniteElement.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
evaluate derivatives of shape function - constant 0
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
int dim
Dimension of reference space.
LagrangeHexFiniteElement(int degree)
Construct the LagrangeHexFiniteElement with the provided degree.
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Data type dense matrix using column-major storage.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
GaussBiQuad2DFiniteElement()
Construct the GaussBiQuad2DFiniteElement.
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual function which evaluates the values of all partial derivatives of all shape functions at a gi...
virtual void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Compute the discrete gradient matrix from the given FiniteElement onto 'this' FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
RefinedTriLinear3DFiniteElement()
Construct the RefinedTriLinear3DFiniteElement.
P0QuadFiniteElement()
Construct the P0QuadFiniteElement.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Quadratic3DFiniteElement()
Construct the Quadratic3DFiniteElement.
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Intermediate class for finite elements whose basis functions return vector values.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
RT1HexFiniteElement()
Construct the RT1HexFiniteElement.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
RT0QuadFiniteElement()
Construct the RT0QuadFiniteElement.
RT1QuadFiniteElement()
Construct the RT1QuadFiniteElement.
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &h) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Linear2DFiniteElement()
Construct the Linear2DFiniteElement.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Cubic2DFiniteElement()
Construct the Cubic2DFiniteElement.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
P0PyrFiniteElement()
Construct the P0PyrFiniteElement.
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
evaluate shape function - constant 1
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &h) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Cubic3DFiniteElement()
Construct the Cubic3DFiniteElement.
Linear1DFiniteElement()
Construct the Linear1DFiniteElement.
A 1D element with uniform nodes.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
RT0HexFiniteElement()
Construct the RT0HexFiniteElement.
Class for standard nodal finite elements.
CrouzeixRaviartFiniteElement()
Construct the CrouzeixRaviartFiniteElement.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Nedelec1PyrFiniteElement()
Construct the Nedelec1PyrFiniteElement.
virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const
Get the dofs associated with the given face. *dofs is set to an internal array of the local dofc on t...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
RefinedLinear1DFiniteElement()
Construct the RefinedLinear1DFiniteElement.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
Cubic1DFiniteElement()
Construct the Cubic1DFiniteElement.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Nedelec1TetFiniteElement()
Construct the Nedelec1TetFiniteElement.
RT2QuadFiniteElement()
Construct the RT2QuadFiniteElement.
virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const
Get the dofs associated with the given face. *dofs is set to an internal array of the local dofc on t...
GaussLinear2DFiniteElement()
Construct the GaussLinear2DFiniteElement.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
void Invert()
Replaces the current matrix with its inverse.
virtual void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Compute the discrete gradient matrix from the given FiniteElement onto 'this' FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void ProjectCurl(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Compute the discrete curl matrix from the given FiniteElement onto 'this' FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
P1TetNonConfFiniteElement()
Construct the P1TetNonConfFiniteElement.
P0SegmentFiniteElement(int Ord=0)
Construct the P0SegmentFiniteElement with dummy order Ord.
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
Nedelec1HexFiniteElement()
Construct the Nedelec1HexFiniteElement.
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual function which evaluates the values of all partial derivatives of all shape functions at a gi...
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &h) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
RT0TetFiniteElement()
Construct the RT0TetFiniteElement.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
LinearPyramidFiniteElement()
Construct the LinearPyramidFiniteElement.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
P0HexFiniteElement()
Construct the P0HexFiniteElement.
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
P1SegmentFiniteElement()
Construct the P1SegmentFiniteElement.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
double p(const Vector &x, double t)
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
BiCubic2DFiniteElement()
Construct the BiCubic2DFiniteElement.
~LagrangeHexFiniteElement()
int GetDof() const
Returns the number of degrees of freedom in the finite element.
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
P0TetFiniteElement()
Construct the P0TetFiniteElement.
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Quad1DFiniteElement()
Construct the Quad1DFiniteElement.
Nedelec1WdgFiniteElement()
Construct the Nedelec1WdgFiniteElement.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual function which evaluates the values of all partial derivatives of all shape functions at a gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Class for integration point with weight.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
int dof
Number of degrees of freedom.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
RT0PyrFiniteElement(bool rt0tets=true)
Construct the RT0PyrFiniteElement.
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
RefinedLinear2DFiniteElement()
Construct the RefinedLinear2DFiniteElement.
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
BiQuad2DFiniteElement()
Construct the BiQuad2DFiniteElement.
virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const
Get the dofs associated with the given face. *dofs is set to an internal array of the local dofc on t...
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
P0TriangleFiniteElement()
Construct the P0TriangleFiniteElement.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
void Mult(const double *x, double *y) const
Matrix vector multiplication.
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Quad2DFiniteElement()
Construct the Quad2DFiniteElement.
Describes the function space on each element.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
GaussBiLinear2DFiniteElement()
Construct the FiniteElement.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual function which evaluates the values of all shape functions at a given point ip and stores the...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
P1OnQuadFiniteElement()
Construct the P1OnQuadFiniteElement.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
RotTriLinearHexFiniteElement()
Construct the RotTriLinearHexFiniteElement.
LinearWedgeFiniteElement()
Construct the LinearWedgeFiniteElement.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
TriLinear3DFiniteElement()
Construct the TriLinear3DFiniteElement.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
BiLinear2DFiniteElement()
Construct the BiLinear2DFiniteElement.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Lagrange1DFiniteElement(int degree)
Construct the Lagrange1DFiniteElement with the provided degree.
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...