23const real_t ND_HexahedronElement::tk[18] =
24{ 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,0.,0., 0.,-1.,0., 0.,0.,-1. };
27 const int cb_type,
const int ob_type)
30 dof2tk(dof), cp(
poly1d.ClosedPoints(
p, cb_type))
37 const int dof3 =
dof/3;
39#ifndef MFEM_THREAD_SAFE
53 for (
int i = 0; i <
p; i++)
55 dof_map[0*dof3 + i + (0 + 0*(
p + 1))*
p] = o++;
57 for (
int i = 0; i <
p; i++)
61 for (
int i = 0; i <
p; i++)
65 for (
int i = 0; i <
p; i++)
67 dof_map[1*dof3 + 0 + (i + 0*
p)*(
p + 1)] = o++;
69 for (
int i = 0; i <
p; i++)
73 for (
int i = 0; i <
p; i++)
77 for (
int i = 0; i <
p; i++)
81 for (
int i = 0; i <
p; i++)
85 for (
int i = 0; i <
p; i++)
87 dof_map[2*dof3 + 0 + (0 + i*(
p + 1))*(
p + 1)] = o++;
89 for (
int i = 0; i <
p; i++)
91 dof_map[2*dof3 +
p + (0 + i*(
p + 1))*(
p + 1)] = o++;
93 for (
int i = 0; i <
p; i++)
95 dof_map[2*dof3 +
p + (
p + i*(
p + 1))*(
p + 1)] = o++;
97 for (
int i = 0; i <
p; i++)
99 dof_map[2*dof3 + 0 + (
p + i*(
p + 1))*(
p + 1)] = o++;
104 for (
int j = 1; j <
p; j++)
105 for (
int i = 0; i <
p; i++)
107 dof_map[0*dof3 + i + ((
p - j) + 0*(
p + 1))*
p] = o++;
109 for (
int j = 0; j <
p; j++)
110 for (
int i = 1; i <
p; i++)
112 dof_map[1*dof3 + i + ((
p - 1 - j) + 0*
p)*(
p + 1)] = -1 - (o++);
115 for (
int k = 1; k <
p; k++)
116 for (
int i = 0; i <
p; i++)
118 dof_map[0*dof3 + i + (0 + k*(
p + 1))*
p] = o++;
120 for (
int k = 0; k <
p; k++)
121 for (
int i = 1; i <
p; i++ )
123 dof_map[2*dof3 + i + (0 + k*(
p + 1))*(
p + 1)] = o++;
126 for (
int k = 1; k <
p; k++)
127 for (
int j = 0; j <
p; j++)
129 dof_map[1*dof3 +
p + (j + k*
p)*(
p + 1)] = o++;
131 for (
int k = 0; k <
p; k++)
132 for (
int j = 1; j <
p; j++)
134 dof_map[2*dof3 +
p + (j + k*(
p + 1))*(
p + 1)] = o++;
137 for (
int k = 1; k <
p; k++)
138 for (
int i = 0; i <
p; i++)
140 dof_map[0*dof3 + (
p - 1 - i) + (
p + k*(
p + 1))*
p] = -1 - (o++);
142 for (
int k = 0; k <
p; k++)
143 for (
int i = 1; i <
p; i++)
145 dof_map[2*dof3 + (
p - i) + (
p + k*(
p + 1))*(
p + 1)] = o++;
148 for (
int k = 1; k <
p; k++)
149 for (
int j = 0; j <
p; j++)
151 dof_map[1*dof3 + 0 + ((
p - 1 - j) + k*
p)*(
p + 1)] = -1 - (o++);
153 for (
int k = 0; k <
p; k++)
154 for (
int j = 1; j <
p; j++)
156 dof_map[2*dof3 + 0 + ((
p - j) + k*(
p + 1))*(
p + 1)] = o++;
159 for (
int j = 1; j <
p; j++)
160 for (
int i = 0; i <
p; i++)
162 dof_map[0*dof3 + i + (j +
p*(
p + 1))*
p] = o++;
164 for (
int j = 0; j <
p; j++)
165 for (
int i = 1; i <
p; i++)
167 dof_map[1*dof3 + i + (j +
p*
p)*(
p + 1)] = o++;
172 for (
int k = 1; k <
p; k++)
173 for (
int j = 1; j <
p; j++)
174 for (
int i = 0; i <
p; i++)
176 dof_map[0*dof3 + i + (j + k*(
p + 1))*
p] = o++;
179 for (
int k = 1; k <
p; k++)
180 for (
int j = 0; j <
p; j++)
181 for (
int i = 1; i <
p; i++)
183 dof_map[1*dof3 + i + (j + k*
p)*(
p + 1)] = o++;
186 for (
int k = 0; k <
p; k++)
187 for (
int j = 1; j <
p; j++)
188 for (
int i = 1; i <
p; i++)
190 dof_map[2*dof3 + i + (j + k*(
p + 1))*(
p + 1)] = o++;
196 for (
int k = 0; k <=
p; k++)
197 for (
int j = 0; j <=
p; j++)
198 for (
int i = 0; i <
p; i++)
203 dof2tk[idx = -1 - idx] = 3;
212 for (
int k = 0; k <=
p; k++)
213 for (
int j = 0; j <
p; j++)
214 for (
int i = 0; i <=
p; i++)
219 dof2tk[idx = -1 - idx] = 4;
228 for (
int k = 0; k <
p; k++)
229 for (
int j = 0; j <=
p; j++)
230 for (
int i = 0; i <=
p; i++)
235 dof2tk[idx = -1 - idx] = 5;
259 for (
int c = 0; c < 3; ++c)
265 for (
int k = 0; k <= km; k++)
266 for (
int j = 0; j <= jm; j++)
267 for (
int i = 0; i <= im; i++)
275 const int id1 = c == 0 ? i : (c == 1 ? j : k);
276 const real_t h = cp[id1+1] - cp[id1];
280 for (
int q = 0; q < nqpt; q++)
286 ip3d.
Set3(cp[i] + (h*ip1d.
x), cp[j], cp[k]);
290 ip3d.
Set3(cp[i], cp[j] + (h*ip1d.
x), cp[k]);
294 ip3d.
Set3(cp[i], cp[j], cp[k] + (h*ip1d.
x));
298 vc.
Eval(xk, Trans, ip3d);
302 val += ip1d.
weight * ipval;
315#ifdef MFEM_THREAD_SAFE
316 Vector shape_cx(
p + 1), shape_ox(
p), shape_cy(
p + 1), shape_oy(
p);
317 Vector shape_cz(
p + 1), shape_oz(
p);
322#ifdef MFEM_THREAD_SAFE
323 Vector dshape_cx(
p + 1), dshape_cy(
p + 1), dshape_cz(
p + 1);
345 for (
int k = 0; k <=
p; k++)
346 for (
int j = 0; j <=
p; j++)
347 for (
int i = 0; i <
p; i++)
352 idx = -1 - idx, s = -1;
358 shape(idx,0) = s*shape_ox(i)*shape_cy(j)*shape_cz(k);
363 for (
int k = 0; k <=
p; k++)
364 for (
int j = 0; j <
p; j++)
365 for (
int i = 0; i <=
p; i++)
370 idx = -1 - idx, s = -1;
377 shape(idx,1) = s*shape_cx(i)*shape_oy(j)*shape_cz(k);
381 for (
int k = 0; k <
p; k++)
382 for (
int j = 0; j <=
p; j++)
383 for (
int i = 0; i <=
p; i++)
388 idx = -1 - idx, s = -1;
396 shape(idx,2) = s*shape_cx(i)*shape_cy(j)*shape_oz(k);
405#ifdef MFEM_THREAD_SAFE
406 Vector shape_cx(
p + 1), shape_ox(
p), shape_cy(
p + 1), shape_oy(
p);
407 Vector shape_cz(
p + 1), shape_oz(
p);
408 Vector dshape_cx(
p + 1), dshape_cy(
p + 1), dshape_cz(
p + 1);
430 for (
int k = 0; k <=
p; k++)
431 for (
int j = 0; j <=
p; j++)
432 for (
int i = 0; i <
p; i++)
437 idx = -1 - idx, s = -1;
443 curl_shape(idx,0) = 0.;
444 curl_shape(idx,1) = s*shape_ox(i)* shape_cy(j)*dshape_cz(k);
445 curl_shape(idx,2) = -s*shape_ox(i)*dshape_cy(j)* shape_cz(k);
448 for (
int k = 0; k <=
p; k++)
449 for (
int j = 0; j <
p; j++)
450 for (
int i = 0; i <=
p; i++)
455 idx = -1 - idx, s = -1;
461 curl_shape(idx,0) = -s* shape_cx(i)*shape_oy(j)*dshape_cz(k);
462 curl_shape(idx,1) = 0.;
463 curl_shape(idx,2) = s*dshape_cx(i)*shape_oy(j)* shape_cz(k);
466 for (
int k = 0; k <
p; k++)
467 for (
int j = 0; j <=
p; j++)
468 for (
int i = 0; i <=
p; i++)
473 idx = -1 - idx, s = -1;
479 curl_shape(idx,0) = s* shape_cx(i)*dshape_cy(j)*shape_oz(k);
480 curl_shape(idx,1) = -s*dshape_cx(i)* shape_cy(j)*shape_oz(k);
481 curl_shape(idx,2) = 0.;
489 const int pp1 =
p + 1;
490 const int n_face_dofs_per_component =
p*pp1;
491 const int n_dof_per_dim =
p*pp1*pp1;
493 std::vector<int> n_dofs = {
p, pp1, pp1,
p};
494 std::vector<int> offsets, strides;
496 const auto f = internal::GetFaceNormal3D(face_id);
497 const int face_normal =
f.first, level =
f.second;
498 if (face_normal == 0)
502 n_dof_per_dim + (level ? pp1 - 1 : 0),
503 2*n_dof_per_dim + (level ? pp1 - 1 : 0)
505 strides = {pp1,
p*pp1, pp1, pp1*pp1};
507 else if (face_normal == 1)
511 level ?
p*(pp1 - 1) : 0,
512 2*n_dof_per_dim + (level ? pp1*(pp1 - 1) : 0)
514 strides = {1,
p*pp1, 1, pp1*pp1};
516 else if (face_normal == 2)
520 level ?
p*pp1*(pp1 - 1) : 0,
521 n_dof_per_dim + (level ?
p*pp1*(pp1 - 1) : 0)
523 strides = {1,
p, 1, pp1};
526 internal::FillFaceMap(n_face_dofs_per_component, offsets, strides, n_dofs,
530const real_t ND_QuadrilateralElement::tk[8] =
531{ 1.,0., 0.,1., -1.,0., 0.,-1. };
539 cp(
poly1d.ClosedPoints(
p, cb_type))
546 const int dof2 =
dof/2;
548#ifndef MFEM_THREAD_SAFE
559 for (
int i = 0; i <
p; i++)
563 for (
int j = 0; j <
p; j++)
567 for (
int i = 0; i <
p; i++)
569 dof_map[0*dof2 + (
p - 1 - i) +
p*
p] = -1 - (o++);
571 for (
int j = 0; j <
p; j++)
573 dof_map[1*dof2 + 0 + (
p - 1 - j)*(
p + 1)] = -1 - (o++);
578 for (
int j = 1; j <
p; j++)
579 for (
int i = 0; i <
p; i++)
584 for (
int j = 0; j <
p; j++)
585 for (
int i = 1; i <
p; i++)
587 dof_map[1*dof2 + i + j*(
p + 1)] = o++;
593 for (
int j = 0; j <=
p; j++)
594 for (
int i = 0; i <
p; i++)
599 dof2tk[idx = -1 - idx] = 2;
608 for (
int j = 0; j <
p; j++)
609 for (
int i = 0; i <=
p; i++)
614 dof2tk[idx = -1 - idx] = 3;
639 for (
int j = 0; j <=
order; j++)
640 for (
int i = 0; i <
order; i++)
648 const real_t h = cp[i+1] - cp[i];
652 for (
int k = 0; k < nqpt; k++)
656 ip2d.
Set2(cp[i] + (h*ip1d.
x), cp[j]);
659 vc.
Eval(xk, Trans, ip2d);
663 val += ip1d.
weight * ipval;
669 for (
int j = 0; j <
order; j++)
670 for (
int i = 0; i <=
order; i++)
678 const real_t h = cp[j+1] - cp[j];
682 for (
int k = 0; k < nqpt; k++)
686 ip2d.
Set2(cp[i], cp[j] + (h*ip1d.
x));
689 vc.
Eval(xk, Trans, ip2d);
693 val += ip1d.
weight * ipval;
705#ifdef MFEM_THREAD_SAFE
706 Vector shape_cx(
p + 1), shape_ox(
p), shape_cy(
p + 1), shape_oy(
p);
711#ifdef MFEM_THREAD_SAFE
712 Vector dshape_cx(
p + 1), dshape_cy(
p + 1);
730 for (
int j = 0; j <=
p; j++)
731 for (
int i = 0; i <
p; i++)
736 idx = -1 - idx, s = -1;
742 shape(idx,0) = s*shape_ox(i)*shape_cy(j);
746 for (
int j = 0; j <
p; j++)
747 for (
int i = 0; i <=
p; i++)
752 idx = -1 - idx, s = -1;
759 shape(idx,1) = s*shape_cx(i)*shape_oy(j);
768#ifdef MFEM_THREAD_SAFE
769 Vector shape_cx(
p + 1), shape_ox(
p), shape_cy(
p + 1), shape_oy(
p);
770 Vector dshape_cx(
p + 1), dshape_cy(
p + 1);
789 for (
int j = 0; j <=
p; j++)
790 for (
int i = 0; i <
p; i++)
795 idx = -1 - idx, s = -1;
801 curl_shape(idx,0) = -s*shape_ox(i)*dshape_cy(j);
804 for (
int j = 0; j <
p; j++)
805 for (
int i = 0; i <=
p; i++)
810 idx = -1 - idx, s = -1;
816 curl_shape(idx,0) = s*dshape_cx(i)*shape_oy(j);
824 const int pp1 =
order + 1;
825 const int n_face_dofs_per_component =
p;
826 std::vector<int> strides = {(face_id == 0 || face_id == 2) ? 1 : pp1};
827 std::vector<int> n_dofs = {
p};
828 std::vector<int> offsets;
831 case 0: offsets = {0};
break;
832 case 1: offsets = {
p*pp1 + pp1 - 1};
break;
833 case 2: offsets = {
p*(pp1 - 1)};
break;
834 case 3: offsets = {
p*pp1};
break;
836 internal::FillFaceMap(n_face_dofs_per_component, offsets, strides, n_dofs,
841const real_t ND_TetrahedronElement::tk[18] =
842{ 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,1.,0., -1.,0.,1., 0.,-1.,1. };
844const real_t ND_TetrahedronElement::c = 1./4.;
854 const int pm1 =
p - 1, pm2 =
p - 2, pm3 =
p - 3;
856#ifndef MFEM_THREAD_SAFE
867 Vector shape_x(
p), shape_y(
p), shape_z(
p), shape_l(
p);
872 for (
int i = 0; i <
p; i++)
877 for (
int i = 0; i <
p; i++)
882 for (
int i = 0; i <
p; i++)
887 for (
int i = 0; i <
p; i++)
892 for (
int i = 0; i <
p; i++)
897 for (
int i = 0; i <
p; i++)
904 for (
int j = 0; j <= pm2; j++)
905 for (
int i = 0; i + j <= pm2; i++)
907 real_t w = fop[i] + fop[j] + fop[pm2-i-j];
913 for (
int j = 0; j <= pm2; j++)
914 for (
int i = 0; i + j <= pm2; i++)
916 real_t w = fop[i] + fop[j] + fop[pm2-i-j];
922 for (
int j = 0; j <= pm2; j++)
923 for (
int i = 0; i + j <= pm2; i++)
925 real_t w = fop[i] + fop[j] + fop[pm2-i-j];
931 for (
int j = 0; j <= pm2; j++)
932 for (
int i = 0; i + j <= pm2; i++)
934 real_t w = fop[i] + fop[j] + fop[pm2-i-j];
942 for (
int k = 0; k <= pm3; k++)
943 for (
int j = 0; j + k <= pm3; j++)
944 for (
int i = 0; i + j + k <= pm3; i++)
946 real_t w = iop[i] + iop[j] + iop[k] + iop[pm3-i-j-k];
956 for (
int m = 0; m <
dof; m++)
959 const real_t *tm = tk + 3*dof2tk[m];
967 for (
int k = 0; k <= pm1; k++)
968 for (
int j = 0; j + k <= pm1; j++)
969 for (
int i = 0; i + j + k <= pm1; i++)
971 real_t s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
972 T(o++, m) = s * tm[0];
973 T(o++, m) = s * tm[1];
974 T(o++, m) = s * tm[2];
976 for (
int k = 0; k <= pm1; k++)
977 for (
int j = 0; j + k <= pm1; j++)
979 real_t s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
980 T(o++, m) = s*((ip.
y - c)*tm[0] - (ip.
x - c)*tm[1]);
981 T(o++, m) = s*((ip.
z - c)*tm[0] - (ip.
x - c)*tm[2]);
983 for (
int k = 0; k <= pm1; k++)
986 shape_y(pm1-k)*shape_z(k)*((ip.
z - c)*tm[1] - (ip.
y - c)*tm[2]);
997 const int pm1 =
order - 1;
999#ifdef MFEM_THREAD_SAFE
1001 Vector shape_x(
p), shape_y(
p), shape_z(
p), shape_l(
p);
1011 for (
int k = 0; k <= pm1; k++)
1012 for (
int j = 0; j + k <= pm1; j++)
1013 for (
int i = 0; i + j + k <= pm1; i++)
1015 real_t s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
1016 u(n,0) = s;
u(n,1) = 0.;
u(n,2) = 0.; n++;
1017 u(n,0) = 0.;
u(n,1) = s;
u(n,2) = 0.; n++;
1018 u(n,0) = 0.;
u(n,1) = 0.;
u(n,2) = s; n++;
1020 for (
int k = 0; k <= pm1; k++)
1021 for (
int j = 0; j + k <= pm1; j++)
1023 real_t s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
1024 u(n,0) = s*(ip.
y - c);
u(n,1) = -s*(ip.
x - c);
u(n,2) = 0.; n++;
1025 u(n,0) = s*(ip.
z - c);
u(n,1) = 0.;
u(n,2) = -s*(ip.
x - c); n++;
1027 for (
int k = 0; k <= pm1; k++)
1029 real_t s = shape_y(pm1-k)*shape_z(k);
1030 u(n,0) = 0.;
u(n,1) = s*(ip.
z - c);
u(n,2) = -s*(ip.
y - c); n++;
1039 const int pm1 =
order - 1;
1041#ifdef MFEM_THREAD_SAFE
1043 Vector shape_x(
p), shape_y(
p), shape_z(
p), shape_l(
p);
1044 Vector dshape_x(
p), dshape_y(
p), dshape_z(
p), dshape_l(
p);
1054 for (
int k = 0; k <= pm1; k++)
1055 for (
int j = 0; j + k <= pm1; j++)
1056 for (
int i = 0; i + j + k <= pm1; i++)
1059 const real_t dx = (dshape_x(i)*shape_l(l) -
1060 shape_x(i)*dshape_l(l))*shape_y(j)*shape_z(k);
1061 const real_t dy = (dshape_y(j)*shape_l(l) -
1062 shape_y(j)*dshape_l(l))*shape_x(i)*shape_z(k);
1063 const real_t dz = (dshape_z(k)*shape_l(l) -
1064 shape_z(k)*dshape_l(l))*shape_x(i)*shape_y(j);
1066 u(n,0) = 0.;
u(n,1) = dz;
u(n,2) = -dy; n++;
1067 u(n,0) = -dz;
u(n,1) = 0.;
u(n,2) = dx; n++;
1068 u(n,0) = dy;
u(n,1) = -dx;
u(n,2) = 0.; n++;
1070 for (
int k = 0; k <= pm1; k++)
1071 for (
int j = 0; j + k <= pm1; j++)
1073 int i = pm1 - j - k;
1076 u(n,0) = shape_x(i)*(ip.
x - c)*shape_y(j)*dshape_z(k);
1077 u(n,1) = shape_x(i)*shape_y(j)*(ip.
y - c)*dshape_z(k);
1079 -((dshape_x(i)*(ip.
x - c) + shape_x(i))*shape_y(j)*shape_z(k) +
1080 (dshape_y(j)*(ip.
y - c) + shape_y(j))*shape_x(i)*shape_z(k));
1083 u(n,0) = -shape_x(i)*(ip.
x - c)*dshape_y(j)*shape_z(k);
1084 u(n,1) = (shape_x(i)*shape_y(j)*(dshape_z(k)*(ip.
z - c) + shape_z(k)) +
1085 (dshape_x(i)*(ip.
x - c) + shape_x(i))*shape_y(j)*shape_z(k));
1086 u(n,2) = -shape_x(i)*dshape_y(j)*shape_z(k)*(ip.
z - c);
1089 for (
int k = 0; k <= pm1; k++)
1093 u(n,0) = -((dshape_y(j)*(ip.
y - c) + shape_y(j))*shape_z(k) +
1094 shape_y(j)*(dshape_z(k)*(ip.
z - c) + shape_z(k)));
1099 Ti.
Mult(
u, curl_shape);
1103const real_t ND_TriangleElement::tk[8] =
1104{ 1.,0., -1.,1., 0.,-1., 0.,1. };
1106const real_t ND_TriangleElement::c = 1./3.;
1111 dof2tk(dof), doftrans(
p)
1116 const int pm1 =
p - 1, pm2 =
p - 2;
1118#ifndef MFEM_THREAD_SAFE
1128 Vector shape_x(
p), shape_y(
p), shape_l(
p);
1133 for (
int i = 0; i <
p; i++)
1138 for (
int i = 0; i <
p; i++)
1143 for (
int i = 0; i <
p; i++)
1150 for (
int j = 0; j <= pm2; j++)
1151 for (
int i = 0; i + j <= pm2; i++)
1153 real_t w = iop[i] + iop[j] + iop[pm2-i-j];
1161 for (
int m = 0; m <
dof; m++)
1164 const real_t *tm = tk + 2*dof2tk[m];
1171 for (
int j = 0; j <= pm1; j++)
1172 for (
int i = 0; i + j <= pm1; i++)
1174 real_t s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
1175 T(n++, m) = s * tm[0];
1176 T(n++, m) = s * tm[1];
1178 for (
int j = 0; j <= pm1; j++)
1181 shape_x(pm1-j)*shape_y(j)*((ip.
y - c)*tm[0] - (ip.
x - c)*tm[1]);
1192 const int pm1 =
order - 1;
1194#ifdef MFEM_THREAD_SAFE
1196 Vector shape_x(
p), shape_y(
p), shape_l(
p);
1205 for (
int j = 0; j <= pm1; j++)
1206 for (
int i = 0; i + j <= pm1; i++)
1208 real_t s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
1209 u(n,0) = s;
u(n,1) = 0; n++;
1210 u(n,0) = 0;
u(n,1) = s; n++;
1212 for (
int j = 0; j <= pm1; j++)
1214 real_t s = shape_x(pm1-j)*shape_y(j);
1215 u(n,0) = s*(ip.
y - c);
1216 u(n,1) = -s*(ip.
x - c);
1226 const int pm1 =
order - 1;
1228#ifdef MFEM_THREAD_SAFE
1230 Vector shape_x(
p), shape_y(
p), shape_l(
p);
1231 Vector dshape_x(
p), dshape_y(
p), dshape_l(
p);
1240 for (
int j = 0; j <= pm1; j++)
1241 for (
int i = 0; i + j <= pm1; i++)
1244 const real_t dx = (dshape_x(i)*shape_l(l) -
1245 shape_x(i)*dshape_l(l)) * shape_y(j);
1246 const real_t dy = (dshape_y(j)*shape_l(l) -
1247 shape_y(j)*dshape_l(l)) * shape_x(i);
1253 for (
int j = 0; j <= pm1; j++)
1257 curlu(n++) = -((dshape_x(i)*(ip.
x - c) + shape_x(i)) * shape_y(j) +
1258 (dshape_y(j)*(ip.
y - c) + shape_y(j)) * shape_x(i));
1262 Ti.
Mult(curlu, curl2d);
1266const real_t ND_SegmentElement::tk[1] = { 1. };
1278 for (
int i = 0; i <
p; i++)
1293const real_t ND_WedgeElement::tk[15] =
1294{ 1.,0.,0., -1.,1.,0., 0.,-1.,0., 0.,0.,1., 0.,1.,0. };
1300 3 *
p * ((
p + 1) * (
p + 2))/2,
p,
1306 H1TriangleFE(
p, cb_type),
1308 H1SegmentFE(
p, cb_type),
1309 NDSegmentFE(
p, ob_type)
1311 MFEM_ASSERT(H1TriangleFE.
GetDof() * NDSegmentFE.
GetDof() +
1313 "Mismatch in number of degrees of freedom "
1314 "when building ND_WedgeElement!");
1316#ifndef MFEM_THREAD_SAFE
1326 const int pm1 =
p - 1, pm2 =
p - 2;
1335 for (
int i = 0; i <
p; i++)
1337 t_dof[o] = i; s_dof[o] = 0; dof2tk[o] = 0;
1342 for (
int i = 0; i <
p; i++)
1344 t_dof[o] =
p + i; s_dof[o] = 0; dof2tk[o] = 1;
1349 for (
int i = 0; i <
p; i++)
1351 t_dof[o] = 2 *
p + i; s_dof[o] = 0; dof2tk[o] = 2;
1356 for (
int i = 0; i <
p; i++)
1358 t_dof[o] = i; s_dof[o] = 1; dof2tk[o] = 0;
1363 for (
int i = 0; i <
p; i++)
1365 t_dof[o] =
p + i; s_dof[o] = 1; dof2tk[o] = 1;
1370 for (
int i = 0; i <
p; i++)
1372 t_dof[o] = 2 *
p + i; s_dof[o] = 1; dof2tk[o] = 2;
1377 for (
int i = 0; i <
p; i++)
1379 t_dof[o] = 0; s_dof[o] = i; dof2tk[o] = 3;
1384 for (
int i = 0; i <
p; i++)
1386 t_dof[o] = 1; s_dof[o] = i; dof2tk[o] = 3;
1391 for (
int i = 0; i <
p; i++)
1393 t_dof[o] = 2; s_dof[o] = i; dof2tk[o] = 3;
1402 for (
int j = 0; j <= pm2; j++)
1403 for (
int i = 0; i + j <= pm2; i++)
1405 l = j + ( 2 *
p - 1 - i) * i / 2;
1406 t_dof[o] = 3 *
p + 2*l+1; s_dof[o] = 0; dof2tk[o] = 4;
1410 t_dof[o] = 3 *
p + 2*l; s_dof[o] = 0; dof2tk[o] = 0;
1417 for (
int j = 0; j <= pm2; j++)
1418 for (
int i = 0; i + j <= pm2; i++)
1420 t_dof[o] = 3 *
p + m; s_dof[o] = 1; dof2tk[o] = 0; m++;
1424 t_dof[o] = 3 *
p + m; s_dof[o] = 1; dof2tk[o] = 4; m++;
1430 for (
int j = 2; j <=
p; j++)
1431 for (
int i = 0; i <
p; i++)
1433 t_dof[o] = i; s_dof[o] = j; dof2tk[o] = 0;
1438 for (
int j = 0; j <
p; j++)
1439 for (
int i = 0; i < pm1; i++)
1441 t_dof[o] = 3 + i; s_dof[o] = j; dof2tk[o] = 3;
1447 for (
int j = 2; j <=
p; j++)
1448 for (
int i = 0; i <
p; i++)
1450 t_dof[o] =
p + i; s_dof[o] = j; dof2tk[o] = 1;
1455 for (
int j = 0; j <
p; j++)
1456 for (
int i = 0; i < pm1; i++)
1458 t_dof[o] =
p + 2 + i; s_dof[o] = j; dof2tk[o] = 3;
1464 for (
int j = 2; j <=
p; j++)
1465 for (
int i = 0; i <
p; i++)
1467 t_dof[o] = 2 *
p + i; s_dof[o] = j; dof2tk[o] = 2;
1472 for (
int j = 0; j <
p; j++)
1473 for (
int i = 0; i < pm1; i++)
1475 t_dof[o] = 2 *
p + 1 + i; s_dof[o] = j; dof2tk[o] = 3;
1482 for (
int k = 2; k <=
p; k++)
1485 for (
int j = 0; j <= pm2; j++)
1486 for (
int i = 0; i + j <= pm2; i++)
1488 t_dof[o] = 3 *
p + l; s_dof[o] = k; dof2tk[o] = 0; l++;
1492 t_dof[o] = 3 *
p + l; s_dof[o] = k; dof2tk[o] = 4; l++;
1498 for (
int k = 0; k <
p; k++)
1501 for (
int j = 0; j < pm2; j++)
1502 for (
int i = 0; i + j < pm2; i++)
1504 t_dof[o] = 3 *
p + l; s_dof[o] = k; dof2tk[o] = 3; l++;
1515#ifdef MFEM_THREAD_SAFE
1529 for (
int i=0; i<
dof; i++)
1531 if ( dof2tk[i] != 3 )
1533 shape(i, 0) = tn_shape(t_dof[i], 0) * s1_shape[s_dof[i]];
1534 shape(i, 1) = tn_shape(t_dof[i], 1) * s1_shape[s_dof[i]];
1541 shape(i, 2) = t1_shape[t_dof[i]] * sn_shape(s_dof[i], 0);
1549#ifdef MFEM_THREAD_SAFE
1567 for (
int i=0; i<
dof; i++)
1569 if ( dof2tk[i] != 3 )
1571 curl_shape(i, 0) = -tn_shape(t_dof[i], 1) * s1_dshape(s_dof[i], 0);
1572 curl_shape(i, 1) = tn_shape(t_dof[i], 0) * s1_dshape(s_dof[i], 0);
1573 curl_shape(i, 2) = tn_dshape(t_dof[i], 0) * s1_shape[s_dof[i]];
1577 curl_shape(i, 0) = t1_dshape(t_dof[i], 1) * sn_shape(s_dof[i], 0);
1578 curl_shape(i, 1) = -t1_dshape(t_dof[i], 0) * sn_shape(s_dof[i], 0);
1579 curl_shape(i, 2) = 0.0;
1584const real_t ND_FuentesPyramidElement::tk[27] =
1586 1., 0., 0., 0., 1., 0., 0., 0., 1.,
1587 -1., 0., 1., -1.,-1., 1., 0.,-1., 1.,
1588 -1., 0., 0., 0.,-1., 0., -M_SQRT1_2,-M_SQRT1_2,M_SQRT2
1596 dof2tk(dof), doftrans(
p)
1605 const int pm2 =
p - 2;
1607#ifndef MFEM_THREAD_SAFE
1643 for (
int i = 0; i <
p; i++)
1648 for (
int i = 0; i <
p; i++)
1653 for (
int i = 0; i <
p; i++)
1658 for (
int i = 0; i <
p; i++)
1663 for (
int i = 0; i <
p; i++)
1668 for (
int i = 0; i <
p; i++)
1673 for (
int i = 0; i <
p; i++)
1678 for (
int i = 0; i <
p; i++)
1686 for (
int j = 1; j <
p; j++)
1687 for (
int i = 0; i <
p; i++)
1694 for (
int j = 0; j <
p; j++)
1695 for (
int i = 1; i <
p; i++)
1702 for (
int j = 0; j <= pm2; j++)
1703 for (
int i = 0; i + j <= pm2; i++)
1705 real_t w = top[i] + top[j] + top[pm2-i-j];
1711 for (
int j = 0; j <= pm2; j++)
1712 for (
int i = 0; i + j <= pm2; i++)
1714 real_t w = top[i] + top[j] + top[pm2-i-j];
1720 for (
int j = 0; j <= pm2; j++)
1721 for (
int i = 0; i + j <= pm2; i++)
1723 real_t w = top[i] + top[j] + top[pm2-i-j];
1731 for (
int j = 0; j <= pm2; j++)
1732 for (
int i = 0; i + j <= pm2; i++)
1734 real_t w = top[i] + top[j] + top[pm2-i-j];
1743 for (
int k = 1; k <
p; k++)
1744 for (
int j = 1; j <
p; j++)
1745 for (
int i = 0; i <
p; i++)
1752 for (
int k = 1; k <
p; k++)
1753 for (
int j = 0; j <
p; j++)
1754 for (
int i = 1; i <
p; i++)
1761 for (
int k = 0; k <
p; k++)
1762 for (
int j = 1; j <
p; j++)
1763 for (
int i = 1; i <
p; i++)
1772 for (
int m = 0; m <
dof; m++)
1775 calcBasis(
p, ip, tmp_E_E_ij, tmp_E_Q1_ijk, tmp_E_Q2_ijk, tmp_E_T_ijk,
1776 tmp_phi_Q1_ij, tmp_dphi_Q1_ij, tmp_phi_Q2_ij,
1777 tmp_phi_E_i, tmp_dphi_E_i,
u);
1779 const Vector tm({tk[3*dof2tk[m]], tk[3*dof2tk[m]+1], tk[3*dof2tk[m]+2]});
1791#ifdef MFEM_THREAD_SAFE
1804 calcBasis(
p, ip, tmp_E_E_ij, tmp_E_Q1_ijk, tmp_E_Q2_ijk, tmp_E_T_ijk,
1805 tmp_phi_Q1_ij, tmp_dphi_Q1_ij, tmp_phi_Q2_ij,
1806 tmp_phi_E_i, tmp_dphi_E_i,
u);
1816#ifdef MFEM_THREAD_SAFE
1832 calcCurlBasis(
p, ip, tmp_E_E_ij, tmp_dE_E_ij, tmp_E_Q1_ijk, tmp_dE_Q1_ijk,
1833 tmp_E_Q2_ijk, tmp_dE_Q2_ijk, tmp_E_T_ijk, tmp_dE_T_ijk,
1834 tmp_phi_Q2_ij, tmp_dphi_Q2_ij, tmp_phi_E_i, tmp_dphi_E_i,
1837 Ti.
Mult(curlu, curl_shape);
1845#ifdef MFEM_THREAD_SAFE
1857 calcBasis(
p, ip, tmp_E_E_ij, tmp_E_Q1_ijk, tmp_E_Q2_ijk, tmp_E_T_ijk,
1858 tmp_phi_Q1_ij, tmp_dphi_Q1_ij, tmp_phi_Q2_ij,
1859 tmp_phi_E_i, tmp_dphi_E_i, shape);
1867#ifdef MFEM_THREAD_SAFE
1882 calcCurlBasis(
p, ip, tmp_E_E_ij, tmp_dE_E_ij, tmp_E_Q1_ijk, tmp_dE_Q1_ijk,
1883 tmp_E_Q2_ijk, tmp_dE_Q2_ijk, tmp_E_T_ijk, tmp_dE_T_ijk,
1884 tmp_phi_Q2_ij, tmp_dphi_Q2_ij, tmp_phi_E_i, tmp_dphi_E_i,
1888void ND_FuentesPyramidElement::calcBasis(
const int p,
1904 Vector xy({x,y}), dmu(3);
1910 y = 0.5 * (1.0 - z);
1911 x = 0.5 * (1.0 - z);
1912 xy(0) = x; xy(1) = y;
1914 zmax = std::max(z, zmax);
1926 for (
int i=0; i<
p; i++, o++)
1927 for (
int k=0; k<3; k++)
1929 W(o, k) =
mu * E_E_ik(i, k);
1934 for (
int i=0; i<
p; i++, o++)
1935 for (
int k=0; k<3; k++)
1937 W(o, k) =
mu * E_E_ik(i, k);
1943 for (
int i=0; i<
p; i++, o++)
1944 for (
int k=0; k<3; k++)
1946 W(o, k) =
mu * E_E_ik(i, k);
1951 for (
int i=0; i<
p; i++, o++)
1952 for (
int k=0; k<3; k++)
1954 W(o, k) =
mu * E_E_ik(i, k);
1962 for (
int i=0; i<
p; i++, o++)
1963 for (
int k=0; k<3; k++)
1965 W(o, k) = E_E_ik(i, k);
1969 for (
int i=0; i<
p; i++, o++)
1970 for (
int k=0; k<3; k++)
1972 W(o, k) = E_E_ik(i, k);
1976 for (
int i=0; i<
p; i++, o++)
1977 for (
int k=0; k<3; k++)
1979 W(o, k) = E_E_ik(i, k);
1983 for (
int i=0; i<
p; i++, o++)
1984 for (
int k=0; k<3; k++)
1986 W(o, k) = E_E_ik(i, k);
1991 if (z < 1.0 && p >= 2)
1997 E_Q(
p,
mu01(z, xy, 1),
mu01_grad_mu01(z, xy, 1),
mu01(z, xy, 2),
1999 for (
int j=2; j<=
p; j++)
2000 for (
int i=0; i<
p; i++, o++)
2001 for (
int k=0; k<3; k++)
2003 W(o, k) = mu2 * E_Q1_ijk(i, j, k);
2007 E_Q(
p,
mu01(z, xy, 2),
mu01_grad_mu01(z, xy, 2),
mu01(z, xy, 1),
2009 for (
int j=2; j<=
p; j++)
2010 for (
int i=0; i<
p; i++, o++)
2011 for (
int k=0; k<3; k++)
2013 W(o, k) = mu2 * E_Q2_ijk(i, j, k);
2018 if (z < 1.0 && p >= 2)
2024 for (
int j=1; j<
p; j++)
2025 for (
int i=0; i+j<
p; i++, o++)
2026 for (
int k=0; k<3; k++)
2028 W(o, k) =
mu * E_T_ijk(i, j, k);
2033 for (
int j=1; j<
p; j++)
2034 for (
int i=0; i+j<
p; i++, o++)
2035 for (
int k=0; k<3; k++)
2037 W(o, k) =
mu * E_T_ijk(i, j, k);
2043 for (
int j=1; j<
p; j++)
2044 for (
int i=0; i+j<
p; i++, o++)
2045 for (
int k=0; k<3; k++)
2047 W(o, k) =
mu * E_T_ijk(i, j, k);
2052 for (
int j=1; j<
p; j++)
2053 for (
int i=0; i+j<
p; i++, o++)
2054 for (
int k=0; k<3; k++)
2056 W(o, k) =
mu * E_T_ijk(i, j, k);
2063 for (
int j=1; j<
p; j++)
2064 for (
int i=0; i+j<
p; i++, o++)
2065 for (
int k=0; k<3; k++)
2067 W(o, k) =
mu * E_T_ijk(i, j, k);
2072 for (
int j=1; j<
p; j++)
2073 for (
int i=0; i+j<
p; i++, o++)
2074 for (
int k=0; k<3; k++)
2076 W(o, k) =
mu * E_T_ijk(i, j, k);
2082 for (
int j=1; j<
p; j++)
2083 for (
int i=0; i+j<
p; i++, o++)
2084 for (
int k=0; k<3; k++)
2086 W(o, k) =
mu * E_T_ijk(i, j, k);
2091 for (
int j=1; j<
p; j++)
2092 for (
int i=0; i+j<
p; i++, o++)
2093 for (
int k=0; k<3; k++)
2095 W(o, k) =
mu * E_T_ijk(i, j, k);
2100 if (z < 1.0 && p >= 2)
2103 phi_Q(
p,
mu01(z, xy, 1),
grad_mu01(z, xy, 1),
mu01(z, xy, 2),
2104 grad_mu01(z, xy, 2), phi_Q1_ij, dphi_Q1_ij);
2106 for (
int k=2; k<=
p; k++)
2107 for (
int j=2; j<=
p; j++)
2108 for (
int i=2; i<=
p; i++, o++)
2109 for (
int l=0; l<3; l++)
2110 W(o, l) = dphi_Q1_ij(i, j, l) * phi_E_k(k) +
2111 phi_Q1_ij(i, j) * dphi_E_k(k, l);
2115 for (
int k=2; k<=
p; k++)
2116 for (
int j=2; j<=
p; j++)
2117 for (
int i=0; i<
p; i++, o++)
2118 for (
int l=0; l<3; l++)
2120 W(o, l) =
mu * E_Q1_ijk(i, j, l) * phi_E_k(k);
2124 for (
int k=2; k<=
p; k++)
2125 for (
int j=2; j<=
p; j++)
2126 for (
int i=0; i<
p; i++, o++)
2127 for (
int l=0; l<3; l++)
2129 W(o, l) =
mu * E_Q2_ijk(i, j, l) * phi_E_k(k);
2136 for (
int j=2; j<=
p; j++)
2137 for (
int i=2; i<=
p; i++, o++)
2139 const int n = std::max(i,j);
2140 const real_t nmu = n * pow(
mu, n-1);
2141 for (
int l=0; l<3; l++)
2143 W(o, l) = nmu * phi_Q2_ij(i, j) * dmu(l);
2149void ND_FuentesPyramidElement::calcCurlBasis(
const int p,
2150 const IntegrationPoint &ip,
2151 DenseMatrix & E_E_ik,
2152 DenseMatrix & dE_E_ik,
2153 DenseTensor & E_Q1_ijk,
2154 DenseTensor & dE_Q1_ijk,
2155 DenseTensor & E_Q2_ijk,
2156 DenseTensor & dE_Q2_ijk,
2157 DenseTensor & E_T_ijk,
2158 DenseTensor & dE_T_ijk,
2159 DenseMatrix & phi_Q2_ij,
2160 DenseTensor & dphi_Q2_ij,
2162 DenseMatrix & dphi_E_k,
2163 DenseMatrix & dW)
const
2168 Vector xy({x,y}), dmu(3);
2169 Vector dmuxE(3), E(3), dphi(3), muphi(3);
2176 y = 0.5 * (1.0 - z);
2177 x = 0.5 * (1.0 - z);
2178 xy(0) = x; xy(1) = y;
2180 zmax = std::max(z, zmax);
2193 for (
int i=0; i<
p; i++, o++)
2195 E(0) = E_E_ik(i, 0); E(1) = E_E_ik(i, 1); E(2) = E_E_ik(i, 2);
2196 dmu.cross3D(E, dmuxE);
2197 for (
int k=0; k<3; k++)
2199 dW(o, k) =
mu * dE_E_ik(i, k) + dmuxE(k);
2206 for (
int i=0; i<
p; i++, o++)
2208 E(0) = E_E_ik(i, 0); E(1) = E_E_ik(i, 1); E(2) = E_E_ik(i, 2);
2209 dmu.cross3D(E, dmuxE);
2210 for (
int k=0; k<3; k++)
2212 dW(o, k) =
mu * dE_E_ik(i, k) + dmuxE(k);
2220 for (
int i=0; i<
p; i++, o++)
2222 E(0) = E_E_ik(i, 0); E(1) = E_E_ik(i, 1); E(2) = E_E_ik(i, 2);
2223 dmu.cross3D(E, dmuxE);
2224 for (
int k=0; k<3; k++)
2226 dW(o, k) =
mu * dE_E_ik(i, k) + dmuxE(k);
2233 for (
int i=0; i<
p; i++, o++)
2235 E(0) = E_E_ik(i, 0); E(1) = E_E_ik(i, 1); E(2) = E_E_ik(i, 2);
2236 dmu.cross3D(E, dmuxE);
2237 for (
int k=0; k<3; k++)
2239 dW(o, k) =
mu * dE_E_ik(i, k) + dmuxE(k);
2248 for (
int i=0; i<
p; i++, o++)
2249 for (
int k=0; k<3; k++)
2251 dW(o, k) = dE_E_ik(i, k);
2255 for (
int i=0; i<
p; i++, o++)
2256 for (
int k=0; k<3; k++)
2258 dW(o, k) = dE_E_ik(i, k);
2262 for (
int i=0; i<
p; i++, o++)
2263 for (
int k=0; k<3; k++)
2265 dW(o, k) = dE_E_ik(i, k);
2269 for (
int i=0; i<
p; i++, o++)
2270 for (
int k=0; k<3; k++)
2272 dW(o, k) = dE_E_ik(i, k);
2277 if (z < 1.0 && p >= 2)
2286 for (
int j=2; j<=
p; j++)
2287 for (
int i=0; i<
p; i++, o++)
2289 E(0) = E_Q1_ijk(i, j, 0);
2290 E(1) = E_Q1_ijk(i, j, 1);
2291 E(2) = E_Q1_ijk(i, j, 2);
2292 dmu.cross3D(E, dmuxE);
2293 for (
int k=0; k<3; k++)
2295 dW(o, k) = mu2 * dE_Q1_ijk(i, j, k) + 2.0 *
mu * dmuxE(k);
2302 for (
int j=2; j<=
p; j++)
2303 for (
int i=0; i<
p; i++, o++)
2305 E(0) = E_Q2_ijk(i, j, 0);
2306 E(1) = E_Q2_ijk(i, j, 1);
2307 E(2) = E_Q2_ijk(i, j, 2);
2308 dmu.cross3D(E, dmuxE);
2309 for (
int k=0; k<3; k++)
2311 dW(o, k) = mu2 * dE_Q2_ijk(i, j, k) + 2.0 *
mu * dmuxE(k);
2317 if (z < 1.0 && p >= 2)
2324 for (
int j=1; j<
p; j++)
2325 for (
int i=0; i+j<
p; i++, o++)
2327 E(0) = E_T_ijk(i, j, 0);
2328 E(1) = E_T_ijk(i, j, 1);
2329 E(2) = E_T_ijk(i, j, 2);
2330 dmu.cross3D(E, dmuxE);
2331 for (
int k=0; k<3; k++)
2333 dW(o, k) =
mu * dE_T_ijk(i, j, k) + dmuxE(k);
2340 for (
int j=1; j<
p; j++)
2341 for (
int i=0; i+j<
p; i++, o++)
2343 E(0) = E_T_ijk(i, j, 0);
2344 E(1) = E_T_ijk(i, j, 1);
2345 E(2) = E_T_ijk(i, j, 2);
2346 dmu.cross3D(E, dmuxE);
2347 for (
int k=0; k<3; k++)
2349 dW(o, k) =
mu * dE_T_ijk(i, j, k) + dmuxE(k);
2357 for (
int j=1; j<
p; j++)
2358 for (
int i=0; i+j<
p; i++, o++)
2360 E(0) = E_T_ijk(i, j, 0);
2361 E(1) = E_T_ijk(i, j, 1);
2362 E(2) = E_T_ijk(i, j, 2);
2363 dmu.cross3D(E, dmuxE);
2364 for (
int k=0; k<3; k++)
2366 dW(o, k) =
mu * dE_T_ijk(i, j, k) + dmuxE(k);
2373 for (
int j=1; j<
p; j++)
2374 for (
int i=0; i+j<
p; i++, o++)
2376 E(0) = E_T_ijk(i, j, 0);
2377 E(1) = E_T_ijk(i, j, 1);
2378 E(2) = E_T_ijk(i, j, 2);
2379 dmu.cross3D(E, dmuxE);
2380 for (
int k=0; k<3; k++)
2382 dW(o, k) =
mu * dE_T_ijk(i, j, k) + dmuxE(k);
2391 for (
int j=1; j<
p; j++)
2392 for (
int i=0; i+j<
p; i++, o++)
2394 E(0) = E_T_ijk(i, j, 0);
2395 E(1) = E_T_ijk(i, j, 1);
2396 E(2) = E_T_ijk(i, j, 2);
2397 dmu.cross3D(E, dmuxE);
2398 for (
int k=0; k<3; k++)
2400 dW(o, k) =
mu * dE_T_ijk(i, j, k) + dmuxE(k);
2407 for (
int j=1; j<
p; j++)
2408 for (
int i=0; i+j<
p; i++, o++)
2410 E(0) = E_T_ijk(i, j, 0);
2411 E(1) = E_T_ijk(i, j, 1);
2412 E(2) = E_T_ijk(i, j, 2);
2413 dmu.cross3D(E, dmuxE);
2414 for (
int k=0; k<3; k++)
2416 dW(o, k) =
mu * dE_T_ijk(i, j, k) + dmuxE(k);
2424 for (
int j=1; j<
p; j++)
2425 for (
int i=0; i+j<
p; i++, o++)
2427 E(0) = E_T_ijk(i, j, 0);
2428 E(1) = E_T_ijk(i, j, 1);
2429 E(2) = E_T_ijk(i, j, 2);
2430 dmu.cross3D(E, dmuxE);
2431 for (
int k=0; k<3; k++)
2433 dW(o, k) =
mu * dE_T_ijk(i, j, k) + dmuxE(k);
2440 for (
int j=1; j<
p; j++)
2441 for (
int i=0; i+j<
p; i++, o++)
2443 E(0) = E_T_ijk(i, j, 0);
2444 E(1) = E_T_ijk(i, j, 1);
2445 E(2) = E_T_ijk(i, j, 2);
2446 dmu.cross3D(E, dmuxE);
2447 for (
int k=0; k<3; k++)
2449 dW(o, k) =
mu * dE_T_ijk(i, j, k) + dmuxE(k);
2455 if (z < 1.0 && p >= 2)
2459 o += (
p - 1) * (
p - 1) * (
p - 1);
2465 for (
int k=2; k<=
p; k++)
2467 dphi(0) = dphi_E_k(k, 0);
2468 dphi(1) = dphi_E_k(k, 1);
2469 dphi(2) = dphi_E_k(k, 2);
2470 add(
mu, dphi, phi_E_k(k), dmu, muphi);
2472 for (
int j=2; j<=
p; j++)
2473 for (
int i=0; i<
p; i++, o++)
2475 E(0) = E_Q1_ijk(i, j, 0);
2476 E(1) = E_Q1_ijk(i, j, 1);
2477 E(2) = E_Q1_ijk(i, j, 2);
2478 muphi.cross3D(E, dmuxE);
2479 for (
int l=0; l<3; l++)
2481 dW(o, l) =
mu * dE_Q1_ijk(i, j, l) * phi_E_k(k) + dmuxE(l);
2487 for (
int k=2; k<=
p; k++)
2489 dphi(0) = dphi_E_k(k, 0);
2490 dphi(1) = dphi_E_k(k, 1);
2491 dphi(2) = dphi_E_k(k, 2);
2492 add(
mu, dphi, phi_E_k(k), dmu, muphi);
2494 for (
int j=2; j<=
p; j++)
2495 for (
int i=0; i<
p; i++, o++)
2497 E(0) = E_Q2_ijk(i, j, 0);
2498 E(1) = E_Q2_ijk(i, j, 1);
2499 E(2) = E_Q2_ijk(i, j, 2);
2500 muphi.cross3D(E, dmuxE);
2501 for (
int l=0; l<3; l++)
2503 dW(o, l) =
mu * dE_Q2_ijk(i, j, l) * phi_E_k(k) + dmuxE(l);
2511 phi_Q(
p,
mu01(z, xy, 2),
grad_mu01(z, xy, 2),
mu01(z, xy, 1),
2512 grad_mu01(z, xy, 1), phi_Q2_ij, dphi_Q2_ij);
2513 for (
int j=2; j<=
p; j++)
2514 for (
int i=2; i<=
p; i++, o++)
2516 const int n = std::max(i,j);
2517 const real_t nmu = n * pow(
mu, n-1);
2519 dphi(0) = dphi_Q2_ij(i, j, 0);
2520 dphi(1) = dphi_Q2_ij(i, j, 1);
2521 dphi(2) = dphi_Q2_ij(i, j, 2);
2522 dphi.cross3D(dmu, muphi);
2524 for (
int l=0; l<3; l++)
2526 dW(o, l) = nmu * muphi(l);
2559const real_t ND_R1D_SegmentElement::tk[9] = { 1.,0.,0., 0.,1.,0., 0.,0.,1. };
2567 cbasis1d(
poly1d.GetBasis(
p, VerifyClosed(cb_type))),
2568 obasis1d(
poly1d.GetBasis(
p - 1, VerifyOpen(ob_type)))
2582#ifndef MFEM_THREAD_SAFE
2594 dof_map[
p] = o; dof2tk[o++] = 1;
2596 dof_map[2*
p+1] = o; dof2tk[o++] = 2;
2600 dof_map[2*
p] = o; dof2tk[o++] = 1;
2602 dof_map[3*
p+1] = o; dof2tk[o++] = 2;
2606 for (
int i = 0; i <
p; i++)
2609 dof_map[i] = o; dof2tk[o++] = 0;
2612 for (
int i = 1; i <
p; i++)
2615 dof_map[
p+i] = o; dof2tk[o++] = 1;
2618 for (
int i = 1; i <
p; i++)
2621 dof_map[2*
p+1+i] = o; dof2tk[o++] = 2;
2630#ifdef MFEM_THREAD_SAFE
2631 Vector shape_cx(
p + 1), shape_ox(
p);
2634 cbasis1d.
Eval(ip.
x, shape_cx);
2635 obasis1d.
Eval(ip.
x, shape_ox);
2639 for (
int i = 0; i <
p; i++)
2641 int idx = dof_map[o++];
2642 shape(idx,0) = shape_ox(i);
2647 for (
int i = 0; i <=
p; i++)
2649 int idx = dof_map[o++];
2651 shape(idx,1) = shape_cx(i);
2655 for (
int i = 0; i <=
p; i++)
2657 int idx = dof_map[o++];
2660 shape(idx,2) = shape_cx(i);
2670 "ND_R1D_SegmentElement cannot be embedded in "
2671 "2 or 3 dimensional spaces");
2672 for (
int i=0; i<
dof; i++)
2674 shape(i, 0) *= JI(0,0);
2683#ifdef MFEM_THREAD_SAFE
2684 Vector shape_cx(
p + 1), shape_ox(
p);
2688 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
2689 obasis1d.
Eval(ip.
x, shape_ox);
2693 for (
int i = 0; i <
p; i++)
2695 int idx = dof_map[o++];
2696 curl_shape(idx,0) = 0.;
2697 curl_shape(idx,1) = 0.;
2698 curl_shape(idx,2) = 0.;
2701 for (
int i = 0; i <=
p; i++)
2703 int idx = dof_map[o++];
2704 curl_shape(idx,0) = 0.;
2705 curl_shape(idx,1) = 0.;
2706 curl_shape(idx,2) = dshape_cx(i);
2709 for (
int i = 0; i <=
p; i++)
2711 int idx = dof_map[o++];
2712 curl_shape(idx,0) = 0.;
2713 curl_shape(idx,1) = -dshape_cx(i);
2714 curl_shape(idx,2) = 0.;
2724 "ND_R1D_SegmentElement cannot be embedded in "
2725 "2 or 3 dimensional spaces");
2726 curl_shape *= (1.0 / J.
Weight());
2736 for (
int k = 0; k <
dof; k++)
2743 dofs(k) = Trans.
Jacobian()(0,0) * t(0) * vk(0) +
2744 t(1) * vk(1) + t(2) * vk(2);
2761 for (
int k = 0; k <
dof; k++)
2765 Vector t1(&tk_ptr[dof2tk[k] * 3], 1);
2766 Vector t3(&tk_ptr[dof2tk[k] * 3], 3);
2778 for (
int d = 0; d <
vdim; d++)
2784 for (
int j = 0; j < shape.
Size(); j++)
2787 if (fabs(s) < 1e-12)
2793 for (
int d = 0; d <
vdim; d++)
2795 I(k, j + d*shape.
Size()) = s*vk[d];
2808 for (
int k = 0; k <
dof; k++)
2812 Vector t1(&tk_ptr[dof2tk[k] * 3], 1);
2813 Vector t3(&tk_ptr[dof2tk[k] * 3], 3);
2825 I(k, j) +=
vshape(j, 0) * vk[0];
2828 I(k, j) +=
vshape(j, 1) * t3(1);
2829 I(k, j) +=
vshape(j, 2) * t3(2);
2836const real_t ND_R2D_SegmentElement::tk[4] = { 1.,0., 0.,1. };
2844 cbasis1d(
poly1d.GetBasis(
p, VerifyClosed(cb_type))),
2845 obasis1d(
poly1d.GetBasis(
p - 1, VerifyOpen(ob_type)))
2854#ifndef MFEM_THREAD_SAFE
2866 dof_map[
p] = o; dof2tk[o++] = 1;
2870 dof_map[2*
p] = o; dof2tk[o++] = 1;
2874 for (
int i = 0; i <
p; i++)
2877 dof_map[i] = o; dof2tk[o++] = 0;
2880 for (
int i = 1; i <
p; i++)
2883 dof_map[
p+i] = o; dof2tk[o++] = 1;
2892#ifdef MFEM_THREAD_SAFE
2893 Vector shape_cx(
p + 1), shape_ox(
p);
2896 cbasis1d.
Eval(ip.
x, shape_cx);
2897 obasis1d.
Eval(ip.
x, shape_ox);
2901 for (
int i = 0; i <
p; i++)
2903 int idx = dof_map[o++];
2904 shape(idx,0) = shape_ox(i);
2908 for (
int i = 0; i <=
p; i++)
2910 int idx = dof_map[o++];
2912 shape(idx,1) = shape_cx(i);
2922 "ND_R2D_SegmentElement cannot be embedded in "
2923 "2 or 3 dimensional spaces");
2924 for (
int i=0; i<
dof; i++)
2926 shape(i, 0) *= JI(0,0);
2935#ifdef MFEM_THREAD_SAFE
2936 Vector shape_cx(
p + 1), shape_ox(
p);
2940 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
2941 obasis1d.
Eval(ip.
x, shape_ox);
2945 for (
int i = 0; i <
p; i++)
2947 int idx = dof_map[o++];
2948 curl_shape(idx,0) = 0.;
2951 for (
int i = 0; i <=
p; i++)
2953 int idx = dof_map[o++];
2954 curl_shape(idx,0) = -dshape_cx(i);
2974 for (
int k = 0; k <
dof; k++)
2976 Vector t1(&tk_ptr[dof2tk[k] * 2], 1);
2977 Vector t2(&tk_ptr[dof2tk[k] * 2], 2);
2988 for (
int i = 0; i <
dim; i++)
2990 Ikj +=
vshape(j, i) * vk[i];
2992 Ikj +=
vshape(j, 1) * t2(1);
2993 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
3009 for (
int k = 0; k <
dof; k++)
3015 Vector t1(&tk_ptr[dof2tk[k] * 2], 1);
3016 Vector t2(&tk_ptr[dof2tk[k] * 2], 2);
3047 "ND_R2D_FiniteElement cannot be embedded in "
3048 "3 dimensional spaces");
3049 for (
int i=0; i<
dof; i++)
3053 shape(i, 0) = sx * JI(0, 0) + sy * JI(1, 0);
3054 shape(i, 1) = sx * JI(0, 1) + sy * JI(1, 1);
3064 "ND_R2D_FiniteElement cannot be embedded in "
3065 "3 dimensional spaces");
3066 for (
int i=0; i<
dof; i++)
3068 real_t sx = curl_shape(i, 0);
3069 real_t sy = curl_shape(i, 1);
3070 curl_shape(i, 0) = sx * J(0, 0) + sy * J(0, 1);
3071 curl_shape(i, 1) = sx * J(1, 0) + sy * J(1, 1);
3073 curl_shape *= (1.0 / Trans.
Weight());
3076void ND_R2D_FiniteElement::LocalInterpolation(
3084#ifdef MFEM_THREAD_SAFE
3097 for (
int k = 0; k <
dof; k++)
3111 for (
int i = 0; i <
dim; i++)
3113 Ikj +=
vshape(j, i) * vk[i];
3115 Ikj +=
vshape(j, 2) * t3(2);
3116 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
3128#ifdef MFEM_THREAD_SAFE
3136 for (
int j = 0; j <
dof; j++)
3146 Jinv.
Mult(t2, pt_data);
3147 for (
int k = 0; k <
dof; k++)
3150 for (
int d = 0; d <
dim; d++)
3152 R_jk +=
vshape(k,d)*pt_data[d];
3154 R_jk +=
vshape(k, 2) * t3(2);
3177 for (
int k = 0; k <
dof; k++)
3203 for (
int k = 0; k <
dof; k++)
3219 for (
int d = 0; d <
vdim; d++)
3225 for (
int j = 0; j < shape.
Size(); j++)
3228 if (fabs(s) < 1e-12)
3234 for (
int d = 0; d <
vdim; d++)
3236 I(k, j + d*shape.
Size()) = s*vk[d];
3249 for (
int k = 0; k <
dof; k++)
3266 for (
int i=0; i<2; i++)
3268 I(k, j) +=
vshape(j, i) * vk[i];
3272 I(k, j) +=
vshape(j, 2) * t3(2);
3289 for (
int k = 0; k <
dof; k++)
3293 for (
int j = 0; j < grad_k.
Size(); j++)
3295 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
3300const real_t ND_R2D_TriangleElement::tk_t[15] =
3301{ 1.,0.,0., -1.,1.,0., 0.,-1.,0., 0.,1.,0., 0.,0.,1. };
3309 int pm1 =
p - 1, pm2 =
p - 2;
3311#ifndef MFEM_THREAD_SAFE
3327 for (
int e=0; e<3; e++)
3330 for (
int i=0; i<
p; i++)
3335 for (
int i=0; i<pm1; i++)
3342 for (
int j = 0; j <= pm2; j++)
3343 for (
int i = 0; i + j <= pm2; i++)
3350 for (
int j = 0; j < pm1; j++)
3351 for (
int i = 0; i + j < pm2; i++)
3356 MFEM_VERIFY(n == ND_FE.
GetDof(),
3357 "ND_R2D_Triangle incorrect number of ND dofs.");
3358 MFEM_VERIFY(h == H1_FE.
GetDof(),
3359 "ND_R2D_Triangle incorrect number of H1 dofs.");
3360 MFEM_VERIFY(o ==
GetDof(),
3361 "ND_R2D_Triangle incorrect number of dofs.");
3366 for (
int i=0; i<
dof; i++)
3385#ifdef MFEM_THREAD_SAFE
3393 for (
int i=0; i<
dof; i++)
3398 shape(i, 0) = nd_shape(idx, 0);
3399 shape(i, 1) = nd_shape(idx, 1);
3406 shape(i, 2) = h1_shape(-idx-1);
3414#ifdef MFEM_THREAD_SAFE
3422 for (
int i=0; i<
dof; i++)
3427 curl_shape(i, 0) = 0.0;
3428 curl_shape(i, 1) = 0.0;
3429 curl_shape(i, 2) = nd_dshape(idx, 0);
3433 curl_shape(i, 0) = h1_dshape(-idx-1, 1);
3434 curl_shape(i, 1) = -h1_dshape(-idx-1, 0);
3435 curl_shape(i, 2) = 0.0;
3441const real_t ND_R2D_QuadrilateralElement::tk_q[15] =
3442{ 1.,0.,0., 0.,1.,0., -1.,0.,0., 0.,-1.,0., 0.,0.,1. };
3448 cbasis1d(
poly1d.GetBasis(
p, VerifyClosed(cb_type))),
3449 obasis1d(
poly1d.GetBasis(
p - 1, VerifyOpen(ob_type)))
3453 const int dofx =
p*(
p+1);
3454 const int dofy =
p*(
p+1);
3455 const int dofxy = dofx+dofy;
3457#ifndef MFEM_THREAD_SAFE
3476 for (
int i = 0; i <
p; i++)
3480 for (
int i = 1; i <
p; i++)
3482 dof_map[dofxy + i + 0*(
p+1)] = o++;
3484 for (
int j = 0; j <
p; j++)
3488 for (
int j = 1; j <
p; j++)
3492 for (
int i = 0; i <
p; i++)
3496 for (
int i = 1; i <
p; i++)
3500 for (
int j = 0; j <
p; j++)
3502 dof_map[dofx + 0 + (
p - 1 - j)*(
p + 1)] = -1 - (o++);
3504 for (
int j = 1; j <
p; j++)
3511 for (
int j = 1; j <
p; j++)
3512 for (
int i = 0; i <
p; i++)
3517 for (
int j = 0; j <
p; j++)
3518 for (
int i = 1; i <
p; i++)
3520 dof_map[dofx + i + j*(
p + 1)] = o++;
3523 for (
int j = 1; j <
p; j++)
3524 for (
int i = 1; i <
p; i++)
3526 dof_map[dofxy + i + j*(
p + 1)] = o++;
3532 for (
int j = 0; j <=
p; j++)
3533 for (
int i = 0; i <
p; i++)
3538 dof2tk[idx = -1 - idx] = 2;
3547 for (
int j = 0; j <
p; j++)
3548 for (
int i = 0; i <=
p; i++)
3553 dof2tk[idx = -1 - idx] = 3;
3562 for (
int j = 0; j <=
p; j++)
3563 for (
int i = 0; i <=
p; i++)
3576#ifdef MFEM_THREAD_SAFE
3577 Vector shape_cx(
p + 1), shape_ox(
p), shape_cy(
p + 1), shape_oy(
p);
3580 cbasis1d.
Eval(ip.
x, shape_cx);
3581 obasis1d.
Eval(ip.
x, shape_ox);
3582 cbasis1d.
Eval(ip.
y, shape_cy);
3583 obasis1d.
Eval(ip.
y, shape_oy);
3587 for (
int j = 0; j <=
p; j++)
3588 for (
int i = 0; i <
p; i++)
3593 idx = -1 - idx, s = -1;
3599 shape(idx,0) = s*shape_ox(i)*shape_cy(j);
3604 for (
int j = 0; j <
p; j++)
3605 for (
int i = 0; i <=
p; i++)
3610 idx = -1 - idx, s = -1;
3617 shape(idx,1) = s*shape_cx(i)*shape_oy(j);
3621 for (
int j = 0; j <=
p; j++)
3622 for (
int i = 0; i <=
p; i++)
3627 shape(idx,2) = shape_cx(i)*shape_cy(j);
3636#ifdef MFEM_THREAD_SAFE
3637 Vector shape_cx(
p + 1), shape_ox(
p), shape_cy(
p + 1), shape_oy(
p);
3638 Vector dshape_cx(
p + 1), dshape_cy(
p + 1);
3641 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
3642 obasis1d.
Eval(ip.
x, shape_ox);
3643 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
3644 obasis1d.
Eval(ip.
y, shape_oy);
3648 for (
int j = 0; j <=
p; j++)
3649 for (
int i = 0; i <
p; i++)
3654 idx = -1 - idx, s = -1;
3660 curl_shape(idx,0) = 0.;
3661 curl_shape(idx,1) = 0.;
3662 curl_shape(idx,2) = -s*shape_ox(i)*dshape_cy(j);
3665 for (
int j = 0; j <
p; j++)
3666 for (
int i = 0; i <=
p; i++)
3671 idx = -1 - idx, s = -1;
3677 curl_shape(idx,0) = 0.;
3678 curl_shape(idx,1) = 0.;
3679 curl_shape(idx,2) = s*dshape_cx(i)*shape_oy(j);
3682 for (
int j = 0; j <=
p; j++)
3683 for (
int i = 0; i <=
p; i++)
3686 curl_shape(idx,0) = shape_cx(i)*dshape_cy(j);
3687 curl_shape(idx,1) = -dshape_cx(i)*shape_cy(j);
3688 curl_shape(idx,2) = 0.;
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
void Factor()
Factor the current DenseMatrix, *a.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication with the inverse of dense matrix.
Data type dense matrix using column-major storage.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
void Threshold(real_t eps)
Replace small entries, abs(a_ij) <= eps, with zero.
void SetRow(int r, const real_t *row)
real_t InnerProduct(const real_t *x, const real_t *y) const
Compute y^t A x.
real_t * Data() const
Returns the matrix data array.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
void GetColumn(int c, Vector &col) const
Rank 3 tensor (array of matrices)
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
Abstract class for all finite elements.
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...
int dof
Number of degrees of freedom.
int GetRangeDim() const
Returns the vector dimension for vector-valued finite elements, which is also the dimension of the in...
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...
int GetDim() const
Returns the reference space dimension for the finite element.
int vdim
Vector dimension of vector-valued basis functions.
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
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...
int cdim
Dimension of curl for vector-valued basis functions.
@ CURL
Implements CalcCurlShape methods.
Geometry::Type geom_type
Geometry::Type of the reference element.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
int order
Order/degree of the shape functions.
int dim
Dimension of reference space.
static DenseMatrix grad_mu01(real_t z)
static Vector lam45_grad_lam45(real_t x, real_t y, real_t z)
void phi_Q(int p, Vector s, Vector t, DenseMatrix &u) const
static constexpr real_t apex_tol
static Vector lam45(real_t x, real_t y, real_t z)
static void phi_E(int p, real_t s0, real_t s1, real_t *u)
static Vector nu012(real_t z, Vector xy, unsigned int ab)
static DenseMatrix grad_nu120(real_t z, Vector xy, unsigned int ab)
static DenseMatrix grad_lam25(real_t x, real_t y, real_t z)
void E_E(int p, Vector s, Vector sds, DenseMatrix &u) const
static Vector nu01_grad_nu01(real_t z, Vector xy, unsigned int ab)
static Vector lam35(real_t x, real_t y, real_t z)
static Vector lam25(real_t x, real_t y, real_t z)
static Vector nu01(real_t z, Vector xy, unsigned int ab)
static Vector lam25_grad_lam25(real_t x, real_t y, real_t z)
static Vector mu01(real_t z)
static DenseMatrix grad_nu012(real_t z, Vector xy, unsigned int ab)
static Vector grad_mu1(real_t z)
void E_T(int p, Vector s, Vector sds, DenseTensor &u) const
static real_t mu0(real_t z)
static Vector nu120(real_t z, Vector xy, unsigned int ab)
static Vector lam15_grad_lam15(real_t x, real_t y, real_t z)
Computes .
static Vector lam35_grad_lam35(real_t x, real_t y, real_t z)
static Vector grad_mu0(real_t z)
static DenseMatrix grad_lam35(real_t x, real_t y, real_t z)
static DenseMatrix grad_lam15(real_t x, real_t y, real_t z)
Gradients of the above two component vectors.
void E_Q(int p, Vector s, Vector ds, Vector t, DenseTensor &u) const
static Vector lam15(real_t x, real_t y, real_t z)
Two component vectors associated with edges touching the apex.
static Vector mu01_grad_mu01(real_t z, Vector xy, unsigned int ab)
static real_t mu1(real_t z)
static DenseMatrix grad_lam45(real_t x, real_t y, real_t z)
static DenseMatrix grad_nu01(real_t z, Vector xy, unsigned int ab)
static Vector nu12_grad_nu12(real_t z, Vector xy, unsigned int ab)
Describes the function space on each element.
const IntegrationPoint & GetCenter(int GeomType) const
Return the center of the given Geometry::Type, GeomType.
static bool CheckPoint(int GeomType, const IntegrationPoint &ip)
Check if the given point is inside the given reference element.
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Class for integration point with weight.
void Set(const real_t *p, const int dim)
void Set2(const real_t x1, const real_t x2)
void Set3(const real_t x1, const real_t x2, const real_t x3)
Class for an integration rule - an Array of IntegrationPoint.
int GetNPoints() const
Returns the number of the points in the integration rule.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
void CalcRawVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
void CalcRawCurlShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
ND_FuentesPyramidElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
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 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...
ND_HexahedronElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Construct the ND_HexahedronElement of order p and closed and open BasisType cb_type and ob_type.
void GetFaceMap(const int face_id, Array< int > &face_map) const override
Return the mapping from lexicographic face DOFs to lexicographic element DOFs for the given local fac...
void ProjectIntegrated(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
ND_QuadrilateralElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Construct the ND_QuadrilateralElement of order p and closed and open BasisType cb_type and ob_type.
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
void GetFaceMap(const int face_id, Array< int > &face_map) const override
Return the mapping from lexicographic face DOFs to lexicographic element DOFs for the given local fac...
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
void ProjectIntegrated(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
ND_R1D_PointElement(int p)
Construct the ND_R1D_PointElement.
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
ND_R1D_SegmentElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Construct the ND_R1D_SegmentElement of order p and closed and open BasisType cb_type and ob_type.
void CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in physical space at the point de...
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const override
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
ND_R2D_FiniteElement(int p, Geometry::Type G, int Do, const real_t *tk_fe)
void CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in physical space at the point de...
void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const override
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
void CalcVShape(ElementTransformation &Trans, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in physical space at the point ...
void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const override
Compute the discrete gradient matrix from the given FiniteElement onto 'this' FiniteElement....
void GetLocalRestriction(ElementTransformation &Trans, DenseMatrix &R) const override
Return a local restriction matrix R (Dof x Dof) mapping fine dofs to coarse dofs.
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
ND_R2D_QuadrilateralElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Construct the ND_R2D_QuadrilateralElement of order p and closed and open BasisType cb_type and ob_typ...
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const override
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
ND_R2D_SegmentElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Construct the ND_R2D_SegmentElement of order p and closed and open BasisType cb_type and ob_type.
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
ND_R2D_TriangleElement(const int p, const int cb_type=BasisType::GaussLobatto)
Construct the ND_R2D_TriangleElement of order p.
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
ND_SegmentElement(const int p, const int ob_type=BasisType::GaussLegendre)
Construct the ND_SegmentElement of order p and open BasisType ob_type.
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
ND_TetrahedronElement(const int p)
Construct the ND_TetrahedronElement of order p.
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
ND_TriangleElement(const int p)
Construct the ND_TriangleElement of order p.
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
ND_WedgeElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
void Eval(const real_t x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
bool IsIntegratedType() const
Returns true if the basis is "integrated", false otherwise.
void ScaleIntegrated(bool scale_integrated_)
Set whether the "integrated" basis should be scaled by the subcell sizes. Has no effect for non-integ...
void EvalIntegrated(const Vector &d, Vector &i) const
Evaluate the "integrated" basis type using pre-computed closed basis derivatives.
const real_t * ClosedPoints(const int p, const int btype=BasisType::GaussLobatto, bool on_device=false)
Get coordinates of a closed (GaussLegendre) set of points if degree p.
const real_t * OpenPoints(const int p, const int btype=BasisType::GaussLegendre, bool on_device=false)
Get coordinates of an open (GaussLegendre) set of points if degree p.
static void CalcBasis(const int p, const real_t x, real_t *u)
Evaluate the values of a hierarchical 1D basis at point x hierarchical = k-th basis function is degre...
Base class for vector Coefficients that optionally depend on time and space.
int GetVDim()
Returns dimension of the vector.
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 ...
Intermediate class for finite elements whose basis functions return vector values.
Poly_1D::Basis & obasis1d
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
real_t u(const Vector &xvec)
void add(const Vector &v1, const Vector &v2, Vector &v)
std::function< real_t(const Vector &)> f(real_t mass_coeff)
void InvertLinearTrans(ElementTransformation &trans, const IntegrationPoint &pt, Vector &x)
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
real_t p(const Vector &x, real_t t)