33 mfem_error (
"FiniteElement::CalcVShape (...)\n"
34 " is not implemented for this class!");
40 mfem_error (
"FiniteElement::CalcVShape 2 (...)\n"
41 " is not implemented for this class!");
47 mfem_error (
"FiniteElement::CalcDivShape (...)\n"
48 " is not implemented for this class!");
54 mfem_error (
"FiniteElement::CalcCurlShape (...)\n"
55 " is not implemented for this class!");
60 mfem_error (
"FiniteElement::GetFaceDofs (...)");
66 mfem_error (
"FiniteElement::CalcHessian (...) is not overloaded !");
72 mfem_error (
"GetLocalInterpolation (...) is not overloaded !");
78 mfem_error (
"FiniteElement::Project (...) is not overloaded !");
84 mfem_error (
"FiniteElement::Project (...) (vector) is not overloaded !");
89 mfem_error(
"FiniteElement::ProjectDelta(...) is not implemented for "
96 mfem_error(
"FiniteElement::Project(...) (fe version) is not implemented "
104 mfem_error(
"FiniteElement::ProjectGrad(...) is not implemented for "
112 mfem_error(
"FiniteElement::ProjectCurl(...) is not implemented for "
120 mfem_error(
"FiniteElement::ProjectDiv(...) is not implemented for "
133 #ifdef MFEM_THREAD_SAFE
139 for (
int i = 0; i < fine_fe.
Dof; i++)
144 for (
int j = 0; j <
Dof; j++)
145 if (fabs (I (i,j) = c_shape (j)) < 1.0e-12)
159 for (
int i = 0; i <
Dof; i++)
165 dofs(i) = coeff.
Eval (Trans, ip);
167 dofs(i) *= Trans.
Weight();
174 MFEM_ASSERT(vc.
GetVDim() <= 3,
"");
179 for (
int i = 0; i <
Dof; i++)
183 vc.
Eval (x, Trans, ip);
186 for (
int j = 0; j < x.Size(); j++)
187 dofs(Dof*j+i) = v[j];
201 for (
int k = 0; k <
Dof; k++)
204 for (
int j = 0; j < shape.Size(); j++)
205 I(k,j) = (fabs(shape(j)) < 1e-12) ? 0.0 : shape(j);
213 for (
int k = 0; k <
Dof; k++)
219 for (
int j = 0; j < vshape.Height(); j++)
220 for (
int d = 0; d < vshape.Width(); d++)
221 I(k+d*Dof,j) = vshape(j,d);
235 for (
int k = 0; k <
Dof; k++)
241 Mult(dshape, Jinv, grad_k);
244 for (
int j = 0; j < grad_k.Height(); j++)
245 for (
int d = 0; d <
Dim; d++)
246 grad(k+d*Dof,j) = grad_k(j,d);
258 for (
int k = 0; k <
Dof; k++)
266 for (
int j = 0; j < div_shape.Size(); j++)
267 div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j)/detJ;
271 for (
int j = 0; j < div_shape.Size(); j++)
272 div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j);
281 for (
int i = 0; i <
Dof; i++)
285 dofs(i) = coeff.
Eval(Trans, ip);
311 Mult(pos_mass, mixed_mass, I);
316 void VectorFiniteElement::CalcShape (
319 mfem_error (
"Error: Cannot use scalar CalcShape(...) function with\n"
320 " VectorFiniteElements!");
323 void VectorFiniteElement::CalcDShape (
324 const IntegrationPoint &ip, DenseMatrix &dshape )
const
326 mfem_error (
"Error: Cannot use scalar CalcDShape(...) function with\n"
327 " VectorFiniteElements!");
333 #ifdef MFEM_THREAD_SAFE
340 shape *= (1.0 / Trans.
Weight());
348 #ifdef MFEM_THREAD_SAFE
368 #ifdef MFEM_THREAD_SAFE
372 for (
int k = 0; k <
Dof; k++)
392 #ifdef MFEM_THREAD_SAFE
397 for (
int k = 0; k <
Dof; k++)
404 Jinv.Mult(nk + d2n[k]*
Dim, vk);
407 double w = 1.0/Trans.
Weight();
408 for (
int d = 0; d <
Dim; d++)
412 for (
int j = 0; j < shape.Size(); j++)
417 for (
int d = 0; d <
Dim; d++)
418 I(k,j+d*shape.Size()) = s*vk[d];
424 mfem_error(
"VectorFiniteElement::Project_RT (fe version)");
433 mfem_error(
"VectorFiniteElement::ProjectGrad_RT works only in 2D!");
440 for (
int k = 0; k <
Dof; k++)
443 tk[0] = nk[d2n[k]*
Dim+1];
444 tk[1] = -nk[d2n[k]*
Dim];
445 dshape.Mult(tk, grad_k);
446 for (
int j = 0; j < grad_k.Size(); j++)
447 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
459 for (
int k = 0; k <
Dof; k++)
462 curl_shape.Mult(nk + d2n[k]*
Dim, curl_k);
463 for (
int j = 0; j < curl_k.Size(); j++)
464 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
475 for (
int k = 0; k <
Dof; k++)
495 for (
int k = 0; k <
Dof; k++)
504 double w = 1.0/Trans.
Weight();
505 for (
int d = 0; d <
Dim; d++)
509 for (
int j = 0; j < shape.Size(); j++)
514 for (
int d = 0; d <
Dim; d++)
515 I(k, j + d*shape.Size()) = s*vk[d];
521 mfem_error(
"VectorFiniteElement::Project_ND (fe version)");
535 for (
int k = 0; k <
Dof; k++)
538 dshape.Mult(tk + d2t[k]*
Dim, grad_k);
539 for (
int j = 0; j < grad_k.Size(); j++)
540 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
551 #ifdef MFEM_THREAD_SAFE
560 for (
int k = 0; k <
Dof; k++)
566 Jinv.
Mult(nk + d2n[k]*
Dim, vk);
568 for (
int j = 0; j <
Dof; j++)
571 for (
int i = 0; i <
Dim; i++)
572 Ikj +=
vshape(j, i) * vk[i];
573 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
585 #ifdef MFEM_THREAD_SAFE
592 for (
int k = 0; k <
Dof; k++)
600 for (
int j = 0; j <
Dof; j++)
603 for (
int i = 0; i <
Dim; i++)
604 Ikj +=
vshape(j, i) * vk[i];
605 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
639 shape(0) = 1. - ip.
x;
664 shape(0) = 1. - ip.
x - ip.
y;
672 dshape(0,0) = -1.; dshape(0,1) = -1.;
673 dshape(1,0) = 1.; dshape(1,1) = 0.;
674 dshape(2,0) = 0.; dshape(2,1) = 1.;
693 shape(0) = (1. - ip.
x) * (1. - ip.
y) ;
694 shape(1) = ip.
x * (1. - ip.
y) ;
695 shape(2) = ip.
x * ip.
y ;
696 shape(3) = (1. - ip.
x) * ip.
y ;
702 dshape(0,0) = -1. + ip.
y; dshape(0,1) = -1. + ip.
x ;
703 dshape(1,0) = 1. - ip.
y; dshape(1,1) = -ip.
x ;
704 dshape(2,0) = ip.
y ; dshape(2,1) = ip.
x ;
705 dshape(3,0) = -ip.
y ; dshape(3,1) = 1. - ip.
x ;
711 h( 0,0) = 0.; h( 0,1) = 1.; h( 0,2) = 0.;
712 h( 1,0) = 0.; h( 1,1) = -1.; h( 1,2) = 0.;
713 h( 2,0) = 0.; h( 2,1) = 1.; h( 2,2) = 0.;
714 h( 3,0) = 0.; h( 3,1) = -1.; h( 3,2) = 0.;
732 const double x = ip.
x, y = ip.
y;
734 shape(0) = 5./3. - 2. * (x + y);
735 shape(1) = 2. * (x - 1./6.);
736 shape(2) = 2. * (y - 1./6.);
742 dshape(0,0) = -2.; dshape(0,1) = -2.;
743 dshape(1,0) = 2.; dshape(1,1) = 0.;
744 dshape(2,0) = 0.; dshape(2,1) = 2.;
749 dofs(vertex) = 2./3.;
750 dofs((vertex+1)%3) = 1./6.;
751 dofs((vertex+2)%3) = 1./6.;
756 const double GaussBiLinear2DFiniteElement::p[] =
757 { 0.2113248654051871177454256, 0.7886751345948128822545744 };
775 const double x = ip.
x, y = ip.
y;
777 shape(0) = 3. * (p[1] - x) * (p[1] - y);
778 shape(1) = 3. * (x - p[0]) * (p[1] - y);
779 shape(2) = 3. * (x - p[0]) * (y - p[0]);
780 shape(3) = 3. * (p[1] - x) * (y - p[0]);
786 const double x = ip.
x, y = ip.
y;
788 dshape(0,0) = 3. * (y - p[1]); dshape(0,1) = 3. * (x - p[1]);
789 dshape(1,0) = 3. * (p[1] - y); dshape(1,1) = 3. * (p[0] - x);
790 dshape(2,0) = 3. * (y - p[0]); dshape(2,1) = 3. * (x - p[0]);
791 dshape(3,0) = 3. * (p[0] - y); dshape(3,1) = 3. * (p[1] - x);
797 dofs(vertex) = p[1]*p[1];
798 dofs((vertex+1)%4) = p[0]*p[1];
799 dofs((vertex+2)%4) = p[0]*p[0];
800 dofs((vertex+3)%4) = p[0]*p[1];
821 shape(0) = 1. - ip.
x - ip.
y;
829 dshape(0,0) = -1.; dshape(0,1) = -1.;
830 dshape(1,0) = 1.; dshape(1,1) = 0.;
831 dshape(2,0) = 0.; dshape(2,1) = 1.;
847 double l1 = 1.0 - x, l2 = x, l3 = 2. * x - 1.;
849 shape(0) = l1 * (-l3);
851 shape(2) = 4. * l1 * l2;
859 dshape(0,0) = 4. * x - 3.;
860 dshape(1,0) = 4. * x - 1.;
861 dshape(2,0) = 4. - 8. * x;
876 const double x = ip.
x, x1 = 1. - x;
880 shape(2) = 2. * x * x1;
886 const double x = ip.
x;
888 dshape(0,0) = 2. * x - 2.;
889 dshape(1,0) = 2. * x;
890 dshape(2,0) = 2. - 4. * x;
913 double x = ip.
x, y = ip.
y;
914 double l1 = 1.-x-y, l2 = x, l3 = y;
916 shape(0) = l1 * (2. * l1 - 1.);
917 shape(1) = l2 * (2. * l2 - 1.);
918 shape(2) = l3 * (2. * l3 - 1.);
919 shape(3) = 4. * l1 * l2;
920 shape(4) = 4. * l2 * l3;
921 shape(5) = 4. * l3 * l1;
927 double x = ip.
x, y = ip.
y;
930 dshape(0,1) = 4. * (x + y) - 3.;
932 dshape(1,0) = 4. * x - 1.;
936 dshape(2,1) = 4. * y - 1.;
938 dshape(3,0) = -4. * (2. * x + y - 1.);
939 dshape(3,1) = -4. * x;
941 dshape(4,0) = 4. * y;
942 dshape(4,1) = 4. * x;
944 dshape(5,0) = -4. * y;
945 dshape(5,1) = -4. * (x + 2. * y - 1.);
985 case 0: dofs(3) = 0.25; dofs(5) = 0.25;
break;
986 case 1: dofs(3) = 0.25; dofs(4) = 0.25;
break;
987 case 2: dofs(4) = 0.25; dofs(5) = 0.25;
break;
993 const double GaussQuad2DFiniteElement::p[] =
994 { 0.0915762135097707434595714634022015, 0.445948490915964886318329253883051 };
1012 for (
int i = 0; i < 6; i++)
1029 const double x = ip.
x, y = ip.
y;
1043 const double x = ip.
x, y = ip.
y;
1044 D(0,0) = 0.; D(0,1) = 0.;
1045 D(1,0) = 1.; D(1,1) = 0.;
1046 D(2,0) = 0.; D(2,1) = 1.;
1047 D(3,0) = 2. * x; D(3,1) = 0.;
1048 D(4,0) = y; D(4,1) = x;
1049 D(5,0) = 0.; D(5,1) = 2. * y;
1081 double x = ip.
x, y = ip.
y;
1082 double l1x, l2x, l3x, l1y, l2y, l3y;
1084 l1x = (x - 1.) * (2. * x - 1);
1085 l2x = 4. * x * (1. - x);
1086 l3x = x * (2. * x - 1.);
1087 l1y = (y - 1.) * (2. * y - 1);
1088 l2y = 4. * y * (1. - y);
1089 l3y = y * (2. * y - 1.);
1091 shape(0) = l1x * l1y;
1092 shape(4) = l2x * l1y;
1093 shape(1) = l3x * l1y;
1094 shape(7) = l1x * l2y;
1095 shape(8) = l2x * l2y;
1096 shape(5) = l3x * l2y;
1097 shape(3) = l1x * l3y;
1098 shape(6) = l2x * l3y;
1099 shape(2) = l3x * l3y;
1105 double x = ip.
x, y = ip.
y;
1106 double l1x, l2x, l3x, l1y, l2y, l3y;
1107 double d1x, d2x, d3x, d1y, d2y, d3y;
1109 l1x = (x - 1.) * (2. * x - 1);
1110 l2x = 4. * x * (1. - x);
1111 l3x = x * (2. * x - 1.);
1112 l1y = (y - 1.) * (2. * y - 1);
1113 l2y = 4. * y * (1. - y);
1114 l3y = y * (2. * y - 1.);
1123 dshape(0,0) = d1x * l1y;
1124 dshape(0,1) = l1x * d1y;
1126 dshape(4,0) = d2x * l1y;
1127 dshape(4,1) = l2x * d1y;
1129 dshape(1,0) = d3x * l1y;
1130 dshape(1,1) = l3x * d1y;
1132 dshape(7,0) = d1x * l2y;
1133 dshape(7,1) = l1x * d2y;
1135 dshape(8,0) = d2x * l2y;
1136 dshape(8,1) = l2x * d2y;
1138 dshape(5,0) = d3x * l2y;
1139 dshape(5,1) = l3x * d2y;
1141 dshape(3,0) = d1x * l3y;
1142 dshape(3,1) = l1x * d3y;
1144 dshape(6,0) = d2x * l3y;
1145 dshape(6,1) = l2x * d3y;
1147 dshape(2,0) = d3x * l3y;
1148 dshape(2,1) = l3x * d3y;
1160 case 0: dofs(4) = 0.25; dofs(7) = 0.25;
break;
1161 case 1: dofs(4) = 0.25; dofs(5) = 0.25;
break;
1162 case 2: dofs(5) = 0.25; dofs(6) = 0.25;
break;
1163 case 3: dofs(6) = 0.25; dofs(7) = 0.25;
break;
1195 double x = ip.
x, y = ip.
y;
1196 double l1x, l2x, l3x, l1y, l2y, l3y;
1198 l1x = (1. - x) * (1. - x);
1199 l2x = 2. * x * (1. - x);
1201 l1y = (1. - y) * (1. - y);
1202 l2y = 2. * y * (1. - y);
1205 shape(0) = l1x * l1y;
1206 shape(4) = l2x * l1y;
1207 shape(1) = l3x * l1y;
1208 shape(7) = l1x * l2y;
1209 shape(8) = l2x * l2y;
1210 shape(5) = l3x * l2y;
1211 shape(3) = l1x * l3y;
1212 shape(6) = l2x * l3y;
1213 shape(2) = l3x * l3y;
1219 double x = ip.
x, y = ip.
y;
1220 double l1x, l2x, l3x, l1y, l2y, l3y;
1221 double d1x, d2x, d3x, d1y, d2y, d3y;
1223 l1x = (1. - x) * (1. - x);
1224 l2x = 2. * x * (1. - x);
1226 l1y = (1. - y) * (1. - y);
1227 l2y = 2. * y * (1. - y);
1237 dshape(0,0) = d1x * l1y;
1238 dshape(0,1) = l1x * d1y;
1240 dshape(4,0) = d2x * l1y;
1241 dshape(4,1) = l2x * d1y;
1243 dshape(1,0) = d3x * l1y;
1244 dshape(1,1) = l3x * d1y;
1246 dshape(7,0) = d1x * l2y;
1247 dshape(7,1) = l1x * d2y;
1249 dshape(8,0) = d2x * l2y;
1250 dshape(8,1) = l2x * d2y;
1252 dshape(5,0) = d3x * l2y;
1253 dshape(5,1) = l3x * d2y;
1255 dshape(3,0) = d1x * l3y;
1256 dshape(3,1) = l1x * d3y;
1258 dshape(6,0) = d2x * l3y;
1259 dshape(6,1) = l2x * d3y;
1261 dshape(2,0) = d3x * l3y;
1262 dshape(2,1) = l3x * d3y;
1270 Vector xx(&tr_ip.
x, 2), shape(s, 9);
1272 for (
int i = 0; i < 9; i++)
1276 for (
int j = 0; j < 9; j++)
1277 if (fabs(I(i,j) = s[j]) < 1.0e-12)
1280 for (
int i = 0; i < 9; i++)
1282 double *d = &I(0,i);
1283 d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
1284 d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
1285 d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
1286 d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
1287 d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
1288 0.25 * (d[0] + d[1] + d[2] + d[3]);
1297 for (
int i = 0; i < 9; i++)
1301 d[i] = coeff.
Eval(Trans, ip);
1303 d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
1304 d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
1305 d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
1306 d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
1307 d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
1308 0.25 * (d[0] + d[1] + d[2] + d[3]);
1318 for (
int i = 0; i < 9; i++)
1322 vc.
Eval (x, Trans, ip);
1323 for (
int j = 0; j < x.Size(); j++)
1326 for (
int j = 0; j < x.Size(); j++)
1328 double *d = &dofs(9*j);
1330 d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
1331 d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
1332 d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
1333 d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
1334 d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
1335 0.25 * (d[0] + d[1] + d[2] + d[3]);
1343 const double p1 = 0.5*(1.-sqrt(3./5.));
1368 const double a = sqrt(5./3.);
1369 const double p1 = 0.5*(1.-sqrt(3./5.));
1371 double x = a*(ip.
x-p1), y = a*(ip.
y-p1);
1372 double l1x, l2x, l3x, l1y, l2y, l3y;
1374 l1x = (x - 1.) * (2. * x - 1);
1375 l2x = 4. * x * (1. - x);
1376 l3x = x * (2. * x - 1.);
1377 l1y = (y - 1.) * (2. * y - 1);
1378 l2y = 4. * y * (1. - y);
1379 l3y = y * (2. * y - 1.);
1381 shape(0) = l1x * l1y;
1382 shape(4) = l2x * l1y;
1383 shape(1) = l3x * l1y;
1384 shape(7) = l1x * l2y;
1385 shape(8) = l2x * l2y;
1386 shape(5) = l3x * l2y;
1387 shape(3) = l1x * l3y;
1388 shape(6) = l2x * l3y;
1389 shape(2) = l3x * l3y;
1395 const double a = sqrt(5./3.);
1396 const double p1 = 0.5*(1.-sqrt(3./5.));
1398 double x = a*(ip.
x-p1), y = a*(ip.
y-p1);
1399 double l1x, l2x, l3x, l1y, l2y, l3y;
1400 double d1x, d2x, d3x, d1y, d2y, d3y;
1402 l1x = (x - 1.) * (2. * x - 1);
1403 l2x = 4. * x * (1. - x);
1404 l3x = x * (2. * x - 1.);
1405 l1y = (y - 1.) * (2. * y - 1);
1406 l2y = 4. * y * (1. - y);
1407 l3y = y * (2. * y - 1.);
1409 d1x = a * (4. * x - 3.);
1410 d2x = a * (4. - 8. * x);
1411 d3x = a * (4. * x - 1.);
1412 d1y = a * (4. * y - 3.);
1413 d2y = a * (4. - 8. * y);
1414 d3y = a * (4. * y - 1.);
1416 dshape(0,0) = d1x * l1y;
1417 dshape(0,1) = l1x * d1y;
1419 dshape(4,0) = d2x * l1y;
1420 dshape(4,1) = l2x * d1y;
1422 dshape(1,0) = d3x * l1y;
1423 dshape(1,1) = l3x * d1y;
1425 dshape(7,0) = d1x * l2y;
1426 dshape(7,1) = l1x * d2y;
1428 dshape(8,0) = d2x * l2y;
1429 dshape(8,1) = l2x * d2y;
1431 dshape(5,0) = d3x * l2y;
1432 dshape(5,1) = l3x * d2y;
1434 dshape(3,0) = d1x * l3y;
1435 dshape(3,1) = l1x * d3y;
1437 dshape(6,0) = d2x * l3y;
1438 dshape(6,1) = l2x * d3y;
1440 dshape(2,0) = d3x * l3y;
1441 dshape(2,1) = l3x * d3y;
1484 double x = ip.
x, y = ip.
y;
1486 double w1x, w2x, w3x, w1y, w2y, w3y;
1487 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
1489 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
1490 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
1492 l0x = (- 4.5) * w1x * w2x * w3x;
1493 l1x = ( 13.5) * x * w2x * w3x;
1494 l2x = (-13.5) * x * w1x * w3x;
1495 l3x = ( 4.5) * x * w1x * w2x;
1497 l0y = (- 4.5) * w1y * w2y * w3y;
1498 l1y = ( 13.5) * y * w2y * w3y;
1499 l2y = (-13.5) * y * w1y * w3y;
1500 l3y = ( 4.5) * y * w1y * w2y;
1502 shape(0) = l0x * l0y;
1503 shape(1) = l3x * l0y;
1504 shape(2) = l3x * l3y;
1505 shape(3) = l0x * l3y;
1506 shape(4) = l1x * l0y;
1507 shape(5) = l2x * l0y;
1508 shape(6) = l3x * l1y;
1509 shape(7) = l3x * l2y;
1510 shape(8) = l2x * l3y;
1511 shape(9) = l1x * l3y;
1512 shape(10) = l0x * l2y;
1513 shape(11) = l0x * l1y;
1514 shape(12) = l1x * l1y;
1515 shape(13) = l2x * l1y;
1516 shape(14) = l1x * l2y;
1517 shape(15) = l2x * l2y;
1523 double x = ip.
x, y = ip.
y;
1525 double w1x, w2x, w3x, w1y, w2y, w3y;
1526 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
1527 double d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
1529 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
1530 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
1532 l0x = (- 4.5) * w1x * w2x * w3x;
1533 l1x = ( 13.5) * x * w2x * w3x;
1534 l2x = (-13.5) * x * w1x * w3x;
1535 l3x = ( 4.5) * x * w1x * w2x;
1537 l0y = (- 4.5) * w1y * w2y * w3y;
1538 l1y = ( 13.5) * y * w2y * w3y;
1539 l2y = (-13.5) * y * w1y * w3y;
1540 l3y = ( 4.5) * y * w1y * w2y;
1542 d0x = -5.5 + ( 18. - 13.5 * x) * x;
1543 d1x = 9. + (-45. + 40.5 * x) * x;
1544 d2x = -4.5 + ( 36. - 40.5 * x) * x;
1545 d3x = 1. + (- 9. + 13.5 * x) * x;
1547 d0y = -5.5 + ( 18. - 13.5 * y) * y;
1548 d1y = 9. + (-45. + 40.5 * y) * y;
1549 d2y = -4.5 + ( 36. - 40.5 * y) * y;
1550 d3y = 1. + (- 9. + 13.5 * y) * y;
1552 dshape( 0,0) = d0x * l0y; dshape( 0,1) = l0x * d0y;
1553 dshape( 1,0) = d3x * l0y; dshape( 1,1) = l3x * d0y;
1554 dshape( 2,0) = d3x * l3y; dshape( 2,1) = l3x * d3y;
1555 dshape( 3,0) = d0x * l3y; dshape( 3,1) = l0x * d3y;
1556 dshape( 4,0) = d1x * l0y; dshape( 4,1) = l1x * d0y;
1557 dshape( 5,0) = d2x * l0y; dshape( 5,1) = l2x * d0y;
1558 dshape( 6,0) = d3x * l1y; dshape( 6,1) = l3x * d1y;
1559 dshape( 7,0) = d3x * l2y; dshape( 7,1) = l3x * d2y;
1560 dshape( 8,0) = d2x * l3y; dshape( 8,1) = l2x * d3y;
1561 dshape( 9,0) = d1x * l3y; dshape( 9,1) = l1x * d3y;
1562 dshape(10,0) = d0x * l2y; dshape(10,1) = l0x * d2y;
1563 dshape(11,0) = d0x * l1y; dshape(11,1) = l0x * d1y;
1564 dshape(12,0) = d1x * l1y; dshape(12,1) = l1x * d1y;
1565 dshape(13,0) = d2x * l1y; dshape(13,1) = l2x * d1y;
1566 dshape(14,0) = d1x * l2y; dshape(14,1) = l1x * d2y;
1567 dshape(15,0) = d2x * l2y; dshape(15,1) = l2x * d2y;
1573 double x = ip.
x, y = ip.
y;
1575 double w1x, w2x, w3x, w1y, w2y, w3y;
1576 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
1577 double d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
1578 double h0x, h1x, h2x, h3x, h0y, h1y, h2y, h3y;
1580 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
1581 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
1583 l0x = (- 4.5) * w1x * w2x * w3x;
1584 l1x = ( 13.5) * x * w2x * w3x;
1585 l2x = (-13.5) * x * w1x * w3x;
1586 l3x = ( 4.5) * x * w1x * w2x;
1588 l0y = (- 4.5) * w1y * w2y * w3y;
1589 l1y = ( 13.5) * y * w2y * w3y;
1590 l2y = (-13.5) * y * w1y * w3y;
1591 l3y = ( 4.5) * y * w1y * w2y;
1593 d0x = -5.5 + ( 18. - 13.5 * x) * x;
1594 d1x = 9. + (-45. + 40.5 * x) * x;
1595 d2x = -4.5 + ( 36. - 40.5 * x) * x;
1596 d3x = 1. + (- 9. + 13.5 * x) * x;
1598 d0y = -5.5 + ( 18. - 13.5 * y) * y;
1599 d1y = 9. + (-45. + 40.5 * y) * y;
1600 d2y = -4.5 + ( 36. - 40.5 * y) * y;
1601 d3y = 1. + (- 9. + 13.5 * y) * y;
1603 h0x = -27. * x + 18.;
1604 h1x = 81. * x - 45.;
1605 h2x = -81. * x + 36.;
1608 h0y = -27. * y + 18.;
1609 h1y = 81. * y - 45.;
1610 h2y = -81. * y + 36.;
1613 h( 0,0) = h0x * l0y; h( 0,1) = d0x * d0y; h( 0,2) = l0x * h0y;
1614 h( 1,0) = h3x * l0y; h( 1,1) = d3x * d0y; h( 1,2) = l3x * h0y;
1615 h( 2,0) = h3x * l3y; h( 2,1) = d3x * d3y; h( 2,2) = l3x * h3y;
1616 h( 3,0) = h0x * l3y; h( 3,1) = d0x * d3y; h( 3,2) = l0x * h3y;
1617 h( 4,0) = h1x * l0y; h( 4,1) = d1x * d0y; h( 4,2) = l1x * h0y;
1618 h( 5,0) = h2x * l0y; h( 5,1) = d2x * d0y; h( 5,2) = l2x * h0y;
1619 h( 6,0) = h3x * l1y; h( 6,1) = d3x * d1y; h( 6,2) = l3x * h1y;
1620 h( 7,0) = h3x * l2y; h( 7,1) = d3x * d2y; h( 7,2) = l3x * h2y;
1621 h( 8,0) = h2x * l3y; h( 8,1) = d2x * d3y; h( 8,2) = l2x * h3y;
1622 h( 9,0) = h1x * l3y; h( 9,1) = d1x * d3y; h( 9,2) = l1x * h3y;
1623 h(10,0) = h0x * l2y; h(10,1) = d0x * d2y; h(10,2) = l0x * h2y;
1624 h(11,0) = h0x * l1y; h(11,1) = d0x * d1y; h(11,2) = l0x * h1y;
1625 h(12,0) = h1x * l1y; h(12,1) = d1x * d1y; h(12,2) = l1x * h1y;
1626 h(13,0) = h2x * l1y; h(13,1) = d2x * d1y; h(13,2) = l2x * h1y;
1627 h(14,0) = h1x * l2y; h(14,1) = d1x * d2y; h(14,2) = l1x * h2y;
1628 h(15,0) = h2x * l2y; h(15,1) = d2x * d2y; h(15,2) = l2x * h2y;
1647 l3 = (0.33333333333333333333-x),
1648 l4 = (0.66666666666666666667-x);
1650 shape(0) = 4.5 * l2 * l3 * l4;
1651 shape(1) = 4.5 * l1 * l3 * l4;
1652 shape(2) = 13.5 * l1 * l2 * l4;
1653 shape(3) = -13.5 * l1 * l2 * l3;
1661 dshape(0,0) = -5.5 + x * (18. - 13.5 * x);
1662 dshape(1,0) = 1. - x * (9. - 13.5 * x);
1663 dshape(2,0) = 9. - x * (45. - 40.5 * x);
1664 dshape(3,0) = -4.5 + x * (36. - 40.5 * x);
1696 double x = ip.
x, y = ip.
y;
1697 double l1 = (-1. + x + y),
1701 shape(0) = -0.5*l1*(3.*l1 + 1.)*(3.*l1 + 2.);
1702 shape(1) = 0.5*x*(lx - 1.)*lx;
1703 shape(2) = 0.5*y*(-1. + ly)*ly;
1704 shape(3) = 4.5*x*l1*(3.*l1 + 1.);
1705 shape(4) = -4.5*x*lx*l1;
1706 shape(5) = 4.5*x*lx*y;
1707 shape(6) = 4.5*x*y*ly;
1708 shape(7) = -4.5*y*l1*ly;
1709 shape(8) = 4.5*y*l1*(1. + 3.*l1);
1710 shape(9) = -27.*x*y*l1;
1716 double x = ip.
x, y = ip.
y;
1718 dshape(0,0) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
1719 dshape(1,0) = 1. + 4.5*x*(-2. + 3.*x);
1721 dshape(3,0) = 4.5*(2. + 9.*x*x - 5.*y + 3.*y*y + 2.*x*(-5. + 6.*y));
1722 dshape(4,0) = -4.5*(1. - 1.*y + x*(-8. + 9.*x + 6.*y));
1723 dshape(5,0) = 4.5*(-1. + 6.*x)*y;
1724 dshape(6,0) = 4.5*y*(-1. + 3.*y);
1725 dshape(7,0) = 4.5*(1. - 3.*y)*y;
1726 dshape(8,0) = 4.5*y*(-5. + 6.*x + 6.*y);
1727 dshape(9,0) = -27.*y*(-1. + 2.*x + y);
1729 dshape(0,1) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
1731 dshape(2,1) = 1. + 4.5*y*(-2. + 3.*y);
1732 dshape(3,1) = 4.5*x*(-5. + 6.*x + 6.*y);
1733 dshape(4,1) = 4.5*(1. - 3.*x)*x;
1734 dshape(5,1) = 4.5*x*(-1. + 3.*x);
1735 dshape(6,1) = 4.5*x*(-1. + 6.*y);
1736 dshape(7,1) = -4.5*(1. + x*(-1. + 6.*y) + y*(-8. + 9.*y));
1737 dshape(8,1) = 4.5*(2. + 3.*x*x + y*(-10. + 9.*y) + x*(-5. + 12.*y));
1738 dshape(9,1) = -27.*x*(-1. + x + 2.*y);
1744 double x = ip.
x, y = ip.
y;
1746 h(0,0) = 18.-27.*(x+y);
1747 h(0,1) = 18.-27.*(x+y);
1748 h(0,2) = 18.-27.*(x+y);
1758 h(3,0) = -45.+81.*x+54.*y;
1759 h(3,1) = -22.5+54.*x+27.*y;
1762 h(4,0) = 36.-81.*x-27.*y;
1767 h(5,1) = -4.5+27.*x;
1771 h(6,1) = -4.5+27.*y;
1776 h(7,2) = 36.-27.*x-81.*y;
1779 h(8,1) = -22.5+27.*x+54.*y;
1780 h(8,2) = -45.+54.*x+81.*y;
1783 h(9,1) = 27.-54.*(x+y);
1856 double x = ip.
x, y = ip.
y, z = ip.
z;
1858 shape(0) = -((-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z)*
1859 (-1 + 3*x + 3*y + 3*z))/2.;
1860 shape(4) = (9*x*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
1861 shape(5) = (-9*x*(-1 + 3*x)*(-1 + x + y + z))/2.;
1862 shape(1) = (x*(2 + 9*(-1 + x)*x))/2.;
1863 shape(6) = (9*y*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
1864 shape(19) = -27*x*y*(-1 + x + y + z);
1865 shape(10) = (9*x*(-1 + 3*x)*y)/2.;
1866 shape(7) = (-9*y*(-1 + 3*y)*(-1 + x + y + z))/2.;
1867 shape(11) = (9*x*y*(-1 + 3*y))/2.;
1868 shape(2) = (y*(2 + 9*(-1 + y)*y))/2.;
1869 shape(8) = (9*z*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
1870 shape(18) = -27*x*z*(-1 + x + y + z);
1871 shape(12) = (9*x*(-1 + 3*x)*z)/2.;
1872 shape(17) = -27*y*z*(-1 + x + y + z);
1873 shape(16) = 27*x*y*z;
1874 shape(14) = (9*y*(-1 + 3*y)*z)/2.;
1875 shape(9) = (-9*z*(-1 + x + y + z)*(-1 + 3*z))/2.;
1876 shape(13) = (9*x*z*(-1 + 3*z))/2.;
1877 shape(15) = (9*y*z*(-1 + 3*z))/2.;
1878 shape(3) = (z*(2 + 9*(-1 + z)*z))/2.;
1884 double x = ip.
x, y = ip.
y, z = ip.
z;
1886 dshape(0,0) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
1887 x*(-4 + 6*y + 6*z)))/2.;
1888 dshape(0,1) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
1889 x*(-4 + 6*y + 6*z)))/2.;
1890 dshape(0,2) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
1891 x*(-4 + 6*y + 6*z)))/2.;
1892 dshape(4,0) = (9*(9*pow(x,2) + (-1 + y + z)*(-2 + 3*y + 3*z) +
1893 2*x*(-5 + 6*y + 6*z)))/2.;
1894 dshape(4,1) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
1895 dshape(4,2) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
1896 dshape(5,0) = (-9*(1 - y - z + x*(-8 + 9*x + 6*y + 6*z)))/2.;
1897 dshape(5,1) = (9*(1 - 3*x)*x)/2.;
1898 dshape(5,2) = (9*(1 - 3*x)*x)/2.;
1899 dshape(1,0) = 1 + (9*x*(-2 + 3*x))/2.;
1902 dshape(6,0) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
1903 dshape(6,1) = (9*(2 + 3*pow(x,2) - 10*y - 5*z + 3*(y + z)*(3*y + z) +
1904 x*(-5 + 12*y + 6*z)))/2.;
1905 dshape(6,2) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
1906 dshape(19,0) = -27*y*(-1 + 2*x + y + z);
1907 dshape(19,1) = -27*x*(-1 + x + 2*y + z);
1908 dshape(19,2) = -27*x*y;
1909 dshape(10,0) = (9*(-1 + 6*x)*y)/2.;
1910 dshape(10,1) = (9*x*(-1 + 3*x))/2.;
1912 dshape(7,0) = (9*(1 - 3*y)*y)/2.;
1913 dshape(7,1) = (-9*(1 + x*(-1 + 6*y) - z + y*(-8 + 9*y + 6*z)))/2.;
1914 dshape(7,2) = (9*(1 - 3*y)*y)/2.;
1915 dshape(11,0) = (9*y*(-1 + 3*y))/2.;
1916 dshape(11,1) = (9*x*(-1 + 6*y))/2.;
1919 dshape(2,1) = 1 + (9*y*(-2 + 3*y))/2.;
1921 dshape(8,0) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
1922 dshape(8,1) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
1923 dshape(8,2) = (9*(2 + 3*pow(x,2) - 5*y - 10*z + 3*(y + z)*(y + 3*z) +
1924 x*(-5 + 6*y + 12*z)))/2.;
1925 dshape(18,0) = -27*z*(-1 + 2*x + y + z);
1926 dshape(18,1) = -27*x*z;
1927 dshape(18,2) = -27*x*(-1 + x + y + 2*z);
1928 dshape(12,0) = (9*(-1 + 6*x)*z)/2.;
1930 dshape(12,2) = (9*x*(-1 + 3*x))/2.;
1931 dshape(17,0) = -27*y*z;
1932 dshape(17,1) = -27*z*(-1 + x + 2*y + z);
1933 dshape(17,2) = -27*y*(-1 + x + y + 2*z);
1934 dshape(16,0) = 27*y*z;
1935 dshape(16,1) = 27*x*z;
1936 dshape(16,2) = 27*x*y;
1938 dshape(14,1) = (9*(-1 + 6*y)*z)/2.;
1939 dshape(14,2) = (9*y*(-1 + 3*y))/2.;
1940 dshape(9,0) = (9*(1 - 3*z)*z)/2.;
1941 dshape(9,1) = (9*(1 - 3*z)*z)/2.;
1942 dshape(9,2) = (9*(-1 + x + y + 8*z - 6*(x + y)*z - 9*pow(z,2)))/2.;
1943 dshape(13,0) = (9*z*(-1 + 3*z))/2.;
1945 dshape(13,2) = (9*x*(-1 + 6*z))/2.;
1947 dshape(15,1) = (9*z*(-1 + 3*z))/2.;
1948 dshape(15,2) = (9*y*(-1 + 6*z))/2.;
1951 dshape(3,2) = 1 + (9*z*(-2 + 3*z))/2.;
2017 shape(0) = 1. - ip.
x - ip.
y - ip.
z;
2026 if (dshape.
Height() == 4)
2028 double *A = &dshape(0,0);
2029 A[0] = -1.; A[4] = -1.; A[8] = -1.;
2030 A[1] = 1.; A[5] = 0.; A[9] = 0.;
2031 A[2] = 0.; A[6] = 1.; A[10] = 0.;
2032 A[3] = 0.; A[7] = 0.; A[11] = 1.;
2036 dshape(0,0) = -1.; dshape(0,1) = -1.; dshape(0,2) = -1.;
2037 dshape(1,0) = 1.; dshape(1,1) = 0.; dshape(1,2) = 0.;
2038 dshape(2,0) = 0.; dshape(2,1) = 1.; dshape(2,2) = 0.;
2039 dshape(3,0) = 0.; dshape(3,1) = 0.; dshape(3,2) = 1.;
2046 static int face_dofs[4][3] = {{1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {0, 1, 2}};
2049 *dofs = face_dofs[face];
2091 double L0, L1, L2, L3;
2093 L0 = 1. - ip.
x - ip.
y - ip.
z;
2098 shape(0) = L0 * ( 2.0 * L0 - 1.0 );
2099 shape(1) = L1 * ( 2.0 * L1 - 1.0 );
2100 shape(2) = L2 * ( 2.0 * L2 - 1.0 );
2101 shape(3) = L3 * ( 2.0 * L3 - 1.0 );
2102 shape(4) = 4.0 * L0 * L1;
2103 shape(5) = 4.0 * L0 * L2;
2104 shape(6) = 4.0 * L0 * L3;
2105 shape(7) = 4.0 * L1 * L2;
2106 shape(8) = 4.0 * L1 * L3;
2107 shape(9) = 4.0 * L2 * L3;
2118 L0 = 1.0 - x - y - z;
2120 dshape(0,0) = dshape(0,1) = dshape(0,2) = 1.0 - 4.0 * L0;
2121 dshape(1,0) = -1.0 + 4.0 * x; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
2122 dshape(2,0) = 0.0; dshape(2,1) = -1.0 + 4.0 * y; dshape(2,2) = 0.0;
2123 dshape(3,0) = dshape(3,1) = 0.0; dshape(3,2) = -1.0 + 4.0 * z;
2124 dshape(4,0) = 4.0 * (L0 - x); dshape(4,1) = dshape(4,2) = -4.0 * x;
2125 dshape(5,0) = dshape(5,2) = -4.0 * y; dshape(5,1) = 4.0 * (L0 - y);
2126 dshape(6,0) = dshape(6,1) = -4.0 * z; dshape(6,2) = 4.0 * (L0 - z);
2127 dshape(7,0) = 4.0 * y; dshape(7,1) = 4.0 * x; dshape(7,2) = 0.0;
2128 dshape(8,0) = 4.0 * z; dshape(8,1) = 0.0; dshape(8,2) = 4.0 * x;
2129 dshape(9,0) = 0.0; dshape(9,1) = 4.0 * z; dshape(9,2) = 4.0 * y;
2171 double x = ip.
x, y = ip.
y, z = ip.
z;
2172 double ox = 1.-x, oy = 1.-y, oz = 1.-z;
2174 shape(0) = ox * oy * oz;
2175 shape(1) = x * oy * oz;
2176 shape(2) = x * y * oz;
2177 shape(3) = ox * y * oz;
2178 shape(4) = ox * oy * z;
2179 shape(5) = x * oy * z;
2180 shape(6) = x * y * z;
2181 shape(7) = ox * y * z;
2187 double x = ip.
x, y = ip.
y, z = ip.
z;
2188 double ox = 1.-x, oy = 1.-y, oz = 1.-z;
2190 dshape(0,0) = - oy * oz;
2191 dshape(0,1) = - ox * oz;
2192 dshape(0,2) = - ox * oy;
2194 dshape(1,0) = oy * oz;
2195 dshape(1,1) = - x * oz;
2196 dshape(1,2) = - x * oy;
2198 dshape(2,0) = y * oz;
2199 dshape(2,1) = x * oz;
2200 dshape(2,2) = - x * y;
2202 dshape(3,0) = - y * oz;
2203 dshape(3,1) = ox * oz;
2204 dshape(3,2) = - ox * y;
2206 dshape(4,0) = - oy * z;
2207 dshape(4,1) = - ox * z;
2208 dshape(4,2) = ox * oy;
2210 dshape(5,0) = oy * z;
2211 dshape(5,1) = - x * z;
2212 dshape(5,2) = x * oy;
2214 dshape(6,0) = y * z;
2215 dshape(6,1) = x * z;
2216 dshape(6,2) = x * y;
2218 dshape(7,0) = - y * z;
2219 dshape(7,1) = ox * z;
2220 dshape(7,2) = ox * y;
2255 shape(0) = 1.0 - 2.0 * ip.
y;
2256 shape(1) = -1.0 + 2.0 * ( ip.
x + ip.
y );
2257 shape(2) = 1.0 - 2.0 * ip.
x;
2263 dshape(0,0) = 0.0; dshape(0,1) = -2.0;
2264 dshape(1,0) = 2.0; dshape(1,1) = 2.0;
2265 dshape(2,0) = -2.0; dshape(2,1) = 0.0;
2287 const double l1 = ip.
x+ip.
y-0.5, l2 = 1.-l1, l3 = ip.
x-ip.
y+0.5, l4 = 1.-l3;
2298 const double x2 = 2.*ip.
x, y2 = 2.*ip.
y;
2300 dshape(0,0) = 1. - x2; dshape(0,1) = -2. + y2;
2301 dshape(1,0) = x2; dshape(1,1) = 1. - y2;
2302 dshape(2,0) = 1. - x2; dshape(2,1) = y2;
2303 dshape(3,0) = -2. + x2; dshape(3,1) = 1. - y2;
2321 double x = ip.
x, y = ip.
y;
2324 shape(0,1) = y - 1.;
2327 shape(2,0) = x - 1.;
2339 const double RT0TriangleFiniteElement::nk[3][2] =
2340 { {0, -1}, {1, 1}, {-1, 0} };
2346 #ifdef MFEM_THREAD_SAFE
2352 for (k = 0; k < 3; k++)
2355 for (j = 0; j < 3; j++)
2358 if (j == k) d -= 1.0;
2359 if (fabs(d) > 1.0e-12)
2361 cerr <<
"RT0TriangleFiniteElement::GetLocalInterpolation (...)\n"
2362 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2378 for (k = 0; k < 3; k++)
2381 ip.
x = vk[0]; ip.
y = vk[1];
2384 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2385 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2386 for (j = 0; j < 3; j++)
2387 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2398 #ifdef MFEM_THREAD_SAFE
2402 for (
int k = 0; k < 3; k++)
2410 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2411 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
2431 double x = ip.
x, y = ip.
y;
2434 shape(0,1) = y - 1.;
2439 shape(3,0) = x - 1.;
2452 const double RT0QuadFiniteElement::nk[4][2] =
2453 { {0, -1}, {1, 0}, {0, 1}, {-1, 0} };
2459 #ifdef MFEM_THREAD_SAFE
2465 for (k = 0; k < 4; k++)
2468 for (j = 0; j < 4; j++)
2471 if (j == k) d -= 1.0;
2472 if (fabs(d) > 1.0e-12)
2474 cerr <<
"RT0QuadFiniteElement::GetLocalInterpolation (...)\n"
2475 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2491 for (k = 0; k < 4; k++)
2494 ip.
x = vk[0]; ip.
y = vk[1];
2497 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2498 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2499 for (j = 0; j < 4; j++)
2500 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2511 #ifdef MFEM_THREAD_SAFE
2515 for (
int k = 0; k < 4; k++)
2523 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2524 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
2552 double x = ip.
x, y = ip.
y;
2554 shape(0,0) = -2 * x * (-1 + x + 2 * y);
2555 shape(0,1) = -2 * (-1 + y) * (-1 + x + 2 * y);
2556 shape(1,0) = 2 * x * (x - y);
2557 shape(1,1) = 2 * (x - y) * (-1 + y);
2558 shape(2,0) = 2 * x * (-1 + 2 * x + y);
2559 shape(2,1) = 2 * y * (-1 + 2 * x + y);
2560 shape(3,0) = 2 * x * (-1 + x + 2 * y);
2561 shape(3,1) = 2 * y * (-1 + x + 2 * y);
2562 shape(4,0) = -2 * (-1 + x) * (x - y);
2563 shape(4,1) = 2 * y * (-x + y);
2564 shape(5,0) = -2 * (-1 + x) * (-1 + 2 * x + y);
2565 shape(5,1) = -2 * y * (-1 + 2 * x + y);
2566 shape(6,0) = -3 * x * (-2 + 2 * x + y);
2567 shape(6,1) = -3 * y * (-1 + 2 * x + y);
2568 shape(7,0) = -3 * x * (-1 + x + 2 * y);
2569 shape(7,1) = -3 * y * (-2 + x + 2 * y);
2575 double x = ip.
x, y = ip.
y;
2577 divshape(0) = -2 * (-4 + 3 * x + 6 * y);
2578 divshape(1) = 2 + 6 * x - 6 * y;
2579 divshape(2) = -4 + 12 * x + 6 * y;
2580 divshape(3) = -4 + 6 * x + 12 * y;
2581 divshape(4) = 2 - 6 * x + 6 * y;
2582 divshape(5) = -2 * (-4 + 6 * x + 3 * y);
2583 divshape(6) = -9 * (-1 + 2 * x + y);
2584 divshape(7) = -9 * (-1 + x + 2 * y);
2587 const double RT1TriangleFiniteElement::nk[8][2] =
2599 #ifdef MFEM_THREAD_SAFE
2605 for (k = 0; k < 8; k++)
2608 for (j = 0; j < 8; j++)
2611 if (j == k) d -= 1.0;
2612 if (fabs(d) > 1.0e-12)
2614 cerr <<
"RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
2615 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2631 for (k = 0; k < 8; k++)
2634 ip.
x = vk[0]; ip.
y = vk[1];
2637 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2638 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2639 for (j = 0; j < 8; j++)
2640 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2650 #ifdef MFEM_THREAD_SAFE
2654 for (
int k = 0; k < 8; k++)
2662 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2663 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
2706 double x = ip.
x, y = ip.
y;
2710 shape(0,1) = -( 1. - 3.*y + 2.*y*y)*( 2. - 3.*x);
2712 shape(1,1) = -( 1. - 3.*y + 2.*y*y)*(-1. + 3.*x);
2714 shape(2,0) = (-x + 2.*x*x)*( 2. - 3.*y);
2716 shape(3,0) = (-x + 2.*x*x)*(-1. + 3.*y);
2720 shape(4,1) = (-y + 2.*y*y)*(-1. + 3.*x);
2722 shape(5,1) = (-y + 2.*y*y)*( 2. - 3.*x);
2724 shape(6,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y);
2726 shape(7,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y);
2729 shape(8,0) = (4.*x - 4.*x*x)*( 2. - 3.*y);
2731 shape(9,0) = (4.*x - 4.*x*x)*(-1. + 3.*y);
2735 shape(10,1) = (4.*y - 4.*y*y)*( 2. - 3.*x);
2737 shape(11,1) = (4.*y - 4.*y*y)*(-1. + 3.*x);
2743 double x = ip.
x, y = ip.
y;
2745 divshape(0) = -(-3. + 4.*y)*( 2. - 3.*x);
2746 divshape(1) = -(-3. + 4.*y)*(-1. + 3.*x);
2747 divshape(2) = (-1. + 4.*x)*( 2. - 3.*y);
2748 divshape(3) = (-1. + 4.*x)*(-1. + 3.*y);
2749 divshape(4) = (-1. + 4.*y)*(-1. + 3.*x);
2750 divshape(5) = (-1. + 4.*y)*( 2. - 3.*x);
2751 divshape(6) = -(-3. + 4.*x)*(-1. + 3.*y);
2752 divshape(7) = -(-3. + 4.*x)*( 2. - 3.*y);
2753 divshape(8) = ( 4. - 8.*x)*( 2. - 3.*y);
2754 divshape(9) = ( 4. - 8.*x)*(-1. + 3.*y);
2755 divshape(10) = ( 4. - 8.*y)*( 2. - 3.*x);
2756 divshape(11) = ( 4. - 8.*y)*(-1. + 3.*x);
2759 const double RT1QuadFiniteElement::nk[12][2] =
2779 #ifdef MFEM_THREAD_SAFE
2785 for (k = 0; k < 12; k++)
2788 for (j = 0; j < 12; j++)
2791 if (j == k) d -= 1.0;
2792 if (fabs(d) > 1.0e-12)
2794 cerr <<
"RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
2795 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2811 for (k = 0; k < 12; k++)
2814 ip.
x = vk[0]; ip.
y = vk[1];
2817 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2818 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2819 for (j = 0; j < 12; j++)
2820 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2830 #ifdef MFEM_THREAD_SAFE
2834 for (
int k = 0; k < 12; k++)
2842 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2843 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
2847 const double RT2TriangleFiniteElement::M[15][15] =
2848 {{ 0, -5.3237900077244501311, 5.3237900077244501311, 16.647580015448900262,
2849 0, 24.442740046346700787, -16.647580015448900262, -12.,
2850 -19.118950038622250656, -47.237900077244501311, 0, -34.414110069520051180,
2851 12., 30.590320061795601049, 15.295160030897800524},
2852 { 0, 1.5, -1.5, -15., 0, 2.625, 15., 15., -4.125, 30., 0, -14.625, -15.,
2854 { 0, -0.67620999227554986889, 0.67620999227554986889, 7.3524199845510997378,
2855 0, -3.4427400463467007866, -7.3524199845510997378, -12.,
2856 4.1189500386222506555, -0.76209992275549868892, 0, 7.4141100695200511800,
2857 12., -6.5903200617956010489, -3.2951600308978005244},
2858 { 0, 0, 1.5, 0, 0, 1.5, -11.471370023173350393, 0, 2.4713700231733503933,
2859 -11.471370023173350393, 0, 2.4713700231733503933, 15.295160030897800524,
2860 0, -3.2951600308978005244},
2861 { 0, 0, 4.875, 0, 0, 4.875, -16.875, 0, -16.875, -16.875, 0, -16.875, 10.5,
2863 { 0, 0, 1.5, 0, 0, 1.5, 2.4713700231733503933, 0, -11.471370023173350393,
2864 2.4713700231733503933, 0, -11.471370023173350393, -3.2951600308978005244,
2865 0, 15.295160030897800524},
2866 { -0.67620999227554986889, 0, -3.4427400463467007866, 0,
2867 7.3524199845510997378, 0.67620999227554986889, 7.4141100695200511800, 0,
2868 -0.76209992275549868892, 4.1189500386222506555, -12.,
2869 -7.3524199845510997378, -3.2951600308978005244, -6.5903200617956010489,
2871 { 1.5, 0, 2.625, 0, -15., -1.5, -14.625, 0, 30., -4.125, 15., 15., 10.5,
2873 { -5.3237900077244501311, 0, 24.442740046346700787, 0, 16.647580015448900262,
2874 5.3237900077244501311, -34.414110069520051180, 0, -47.237900077244501311,
2875 -19.118950038622250656, -12., -16.647580015448900262, 15.295160030897800524,
2876 30.590320061795601049, 12.},
2877 { 0, 0, 18., 0, 0, 6., -42., 0, -30., -26., 0, -14., 24., 32., 8.},
2878 { 0, 0, 6., 0, 0, 18., -14., 0, -26., -30., 0, -42., 8., 32., 24.},
2879 { 0, 0, -6., 0, 0, -4., 30., 0, 4., 22., 0, 4., -24., -16., 0},
2880 { 0, 0, -4., 0, 0, -8., 20., 0, 8., 36., 0, 8., -16., -32., 0},
2881 { 0, 0, -8., 0, 0, -4., 8., 0, 36., 8., 0, 20., 0, -32., -16.},
2882 { 0, 0, -4., 0, 0, -6., 4., 0, 22., 4., 0, 30., 0, -16., -24.}
2888 const double p = 0.11270166537925831148;
2925 double x = ip.
x, y = ip.
y;
2927 double Bx[15] = {1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y, 0., x*x*x,
2929 double By[15] = {0., 1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y,
2930 x*x*y, x*y*y, y*y*y};
2932 for (
int i = 0; i < 15; i++)
2934 double cx = 0.0, cy = 0.0;
2935 for (
int j = 0; j < 15; j++)
2937 cx += M[i][j] * Bx[j];
2938 cy += M[i][j] * By[j];
2948 double x = ip.
x, y = ip.
y;
2950 double DivB[15] = {0., 0., 1., 0., 0., 1., 2.*x, 0., y, x, 0., 2.*y,
2951 4.*x*x, 4.*x*y, 4.*y*y};
2953 for (
int i = 0; i < 15; i++)
2956 for (
int j = 0; j < 15; j++)
2957 div += M[i][j] * DivB[j];
2962 const double RT2QuadFiniteElement::pt[4] = {0.,1./3.,2./3.,1.};
2964 const double RT2QuadFiniteElement::dpt[3] = {0.25,0.5,0.75};
3006 double x = ip.
x, y = ip.
y;
3008 double ax0 = pt[0] - x;
3009 double ax1 = pt[1] - x;
3010 double ax2 = pt[2] - x;
3011 double ax3 = pt[3] - x;
3013 double by0 = dpt[0] - y;
3014 double by1 = dpt[1] - y;
3015 double by2 = dpt[2] - y;
3017 double ay0 = pt[0] - y;
3018 double ay1 = pt[1] - y;
3019 double ay2 = pt[2] - y;
3020 double ay3 = pt[3] - y;
3022 double bx0 = dpt[0] - x;
3023 double bx1 = dpt[1] - x;
3024 double bx2 = dpt[2] - x;
3026 double A01 = pt[0] - pt[1];
3027 double A02 = pt[0] - pt[2];
3028 double A12 = pt[1] - pt[2];
3029 double A03 = pt[0] - pt[3];
3030 double A13 = pt[1] - pt[3];
3031 double A23 = pt[2] - pt[3];
3033 double B01 = dpt[0] - dpt[1];
3034 double B02 = dpt[0] - dpt[2];
3035 double B12 = dpt[1] - dpt[2];
3037 double tx0 = (bx1*bx2)/(B01*B02);
3038 double tx1 = -(bx0*bx2)/(B01*B12);
3039 double tx2 = (bx0*bx1)/(B02*B12);
3041 double ty0 = (by1*by2)/(B01*B02);
3042 double ty1 = -(by0*by2)/(B01*B12);
3043 double ty2 = (by0*by1)/(B02*B12);
3047 shape(0, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx0;
3049 shape(1, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx1;
3051 shape(2, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx2;
3053 shape(3, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty0;
3055 shape(4, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty1;
3057 shape(5, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty2;
3061 shape(6, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx2;
3063 shape(7, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx1;
3065 shape(8, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx0;
3067 shape(9, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty2;
3069 shape(10, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty1;
3071 shape(11, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty0;
3074 shape(12, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty0;
3076 shape(13, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty1;
3078 shape(14, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty2;
3081 shape(15, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty0;
3083 shape(16, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty1;
3085 shape(17, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty2;
3089 shape(18, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx0;
3091 shape(19, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx1;
3093 shape(20, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx2;
3096 shape(21, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx0;
3098 shape(22, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx1;
3100 shape(23, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx2;
3106 double x = ip.
x, y = ip.
y;
3108 double a01 = pt[0]*pt[1];
3109 double a02 = pt[0]*pt[2];
3110 double a12 = pt[1]*pt[2];
3111 double a03 = pt[0]*pt[3];
3112 double a13 = pt[1]*pt[3];
3113 double a23 = pt[2]*pt[3];
3115 double bx0 = dpt[0] - x;
3116 double bx1 = dpt[1] - x;
3117 double bx2 = dpt[2] - x;
3119 double by0 = dpt[0] - y;
3120 double by1 = dpt[1] - y;
3121 double by2 = dpt[2] - y;
3123 double A01 = pt[0] - pt[1];
3124 double A02 = pt[0] - pt[2];
3125 double A12 = pt[1] - pt[2];
3126 double A03 = pt[0] - pt[3];
3127 double A13 = pt[1] - pt[3];
3128 double A23 = pt[2] - pt[3];
3130 double A012 = pt[0] + pt[1] + pt[2];
3131 double A013 = pt[0] + pt[1] + pt[3];
3132 double A023 = pt[0] + pt[2] + pt[3];
3133 double A123 = pt[1] + pt[2] + pt[3];
3135 double B01 = dpt[0] - dpt[1];
3136 double B02 = dpt[0] - dpt[2];
3137 double B12 = dpt[1] - dpt[2];
3139 double tx0 = (bx1*bx2)/(B01*B02);
3140 double tx1 = -(bx0*bx2)/(B01*B12);
3141 double tx2 = (bx0*bx1)/(B02*B12);
3143 double ty0 = (by1*by2)/(B01*B02);
3144 double ty1 = -(by0*by2)/(B01*B12);
3145 double ty2 = (by0*by1)/(B02*B12);
3148 divshape(0) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx0;
3149 divshape(1) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx1;
3150 divshape(2) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx2;
3152 divshape(3) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty0;
3153 divshape(4) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty1;
3154 divshape(5) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty2;
3156 divshape(6) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx2;
3157 divshape(7) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx1;
3158 divshape(8) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx0;
3160 divshape(9) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty2;
3161 divshape(10) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty1;
3162 divshape(11) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty0;
3164 divshape(12) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty0;
3165 divshape(13) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty1;
3166 divshape(14) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty2;
3168 divshape(15) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty0;
3169 divshape(16) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty1;
3170 divshape(17) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty2;
3172 divshape(18) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx0;
3173 divshape(19) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx1;
3174 divshape(20) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx2;
3176 divshape(21) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx0;
3177 divshape(22) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx1;
3178 divshape(23) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx2;
3181 const double RT2QuadFiniteElement::nk[24][2] =
3184 {0,-1}, {0,-1}, {0,-1},
3186 {1, 0}, {1, 0}, {1, 0},
3188 {0, 1}, {0, 1}, {0, 1},
3190 {-1,0}, {-1,0}, {-1,0},
3192 {1, 0}, {1, 0}, {1, 0},
3194 {1, 0}, {1, 0}, {1, 0},
3196 {0, 1}, {0, 1}, {0, 1},
3198 {0, 1}, {0, 1}, {0, 1}
3205 #ifdef MFEM_THREAD_SAFE
3211 for (k = 0; k < 24; k++)
3214 for (j = 0; j < 24; j++)
3217 if (j == k) d -= 1.0;
3218 if (fabs(d) > 1.0e-12)
3220 cerr <<
"RT2QuadFiniteElement::GetLocalInterpolation (...)\n"
3221 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
3237 for (k = 0; k < 24; k++)
3240 ip.
x = vk[0]; ip.
y = vk[1];
3243 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
3244 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
3245 for (j = 0; j < 24; j++)
3246 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
3256 #ifdef MFEM_THREAD_SAFE
3260 for (
int k = 0; k < 24; k++)
3268 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
3269 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
3285 shape(0) = 2. - 3. * x;
3286 shape(1) = 3. * x - 1.;
3300 const double p = 0.11270166537925831148;
3310 const double p = 0.11270166537925831148;
3311 const double w = 1./((1-2*p)*(1-2*p));
3314 shape(0) = (2*x-1)*(x-1+p)*w;
3315 shape(1) = 4*(x-1+p)*(p-x)*w;
3316 shape(2) = (2*x-1)*(x-p)*w;
3322 const double p = 0.11270166537925831148;
3323 const double w = 1./((1-2*p)*(1-2*p));
3326 dshape(0,0) = (-3+4*x+2*p)*w;
3327 dshape(1,0) = (4-8*x)*w;
3328 dshape(2,0) = (-1+4*x-2*p)*w;
3339 for (i = 1; i < m; i++)
3343 #ifndef MFEM_THREAD_SAFE
3348 for (i = 1; i <= m; i++)
3349 rwk(i) = rwk(i-1) * ( (double)(m) / (double)(i) );
3350 for (i = 0; i < m/2+1; i++)
3351 rwk(m-i) = ( rwk(i) *= rwk(m-i) );
3352 for (i = m-1; i >= 0; i -= 2)
3359 double w, wk, x = ip.
x;
3362 #ifdef MFEM_THREAD_SAFE
3366 k = (int) floor ( m * x + 0.5 );
3368 for (i = 0; i <= m; i++)
3370 wk *= ( rxxk(i) = x - (double)(i) / m );
3371 w = wk * ( rxxk(k) = x - (double)(k) / m );
3374 shape(0) = w * rwk(0) / rxxk(0);
3376 shape(0) = wk * rwk(0);
3378 shape(1) = w * rwk(m) / rxxk(m);
3380 shape(1) = wk * rwk(k);
3381 for (i = 1; i < m; i++)
3383 shape(i+1) = w * rwk(i) / rxxk(i);
3385 shape(k+1) = wk * rwk(k);
3391 double s, srx, w, wk, x = ip.
x;
3394 #ifdef MFEM_THREAD_SAFE
3398 k = (int) floor ( m * x + 0.5 );
3400 for (i = 0; i <= m; i++)
3402 wk *= ( rxxk(i) = x - (double)(i) / m );
3403 w = wk * ( rxxk(k) = x - (double)(k) / m );
3405 for (i = 0; i <= m; i++)
3406 rxxk(i) = 1.0 / rxxk(i);
3408 for (i = 0; i <= m; i++)
3414 dshape(0,0) = (s - w * rxxk(0)) * rwk(0) * rxxk(0);
3416 dshape(0,0) = wk * srx * rwk(0);
3418 dshape(1,0) = (s - w * rxxk(m)) * rwk(m) * rxxk(m);
3420 dshape(1,0) = wk * srx * rwk(k);
3421 for (i = 1; i < m; i++)
3423 dshape(i+1,0) = (s - w * rxxk(i)) * rwk(i) * rxxk(i);
3425 dshape(k+1,0) = wk * srx * rwk(k);
3453 double L0, L1, L2, L3;
3455 L1 = ip.
x; L2 = ip.
y; L3 = ip.
z; L0 = 1.0 - L1 - L2 - L3;
3456 shape(0) = 1.0 - 3.0 * L0;
3457 shape(1) = 1.0 - 3.0 * L1;
3458 shape(2) = 1.0 - 3.0 * L2;
3459 shape(3) = 1.0 - 3.0 * L3;
3465 dshape(0,0) = 3.0; dshape(0,1) = 3.0; dshape(0,2) = 3.0;
3466 dshape(1,0) = -3.0; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
3467 dshape(2,0) = 0.0; dshape(2,1) = -3.0; dshape(2,2) = 0.0;
3468 dshape(3,0) = 0.0; dshape(3,1) = 0.0; dshape(3,2) = -3.0;
3489 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3510 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3524 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3525 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3526 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3527 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3528 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3529 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3530 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3531 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3533 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3534 I[ 9] = 1; J[ 9] = 2; K[ 9] = 0;
3535 I[10] = 2; J[10] = 1; K[10] = 0;
3536 I[11] = 0; J[11] = 2; K[11] = 0;
3537 I[12] = 2; J[12] = 0; K[12] = 1;
3538 I[13] = 1; J[13] = 2; K[13] = 1;
3539 I[14] = 2; J[14] = 1; K[14] = 1;
3540 I[15] = 0; J[15] = 2; K[15] = 1;
3541 I[16] = 0; J[16] = 0; K[16] = 2;
3542 I[17] = 1; J[17] = 0; K[17] = 2;
3543 I[18] = 1; J[18] = 1; K[18] = 2;
3544 I[19] = 0; J[19] = 1; K[19] = 2;
3546 I[20] = 2; J[20] = 2; K[20] = 0;
3547 I[21] = 2; J[21] = 0; K[21] = 2;
3548 I[22] = 1; J[22] = 2; K[22] = 2;
3549 I[23] = 2; J[23] = 1; K[23] = 2;
3550 I[24] = 0; J[24] = 2; K[24] = 2;
3551 I[25] = 2; J[25] = 2; K[25] = 1;
3553 I[26] = 2; J[26] = 2; K[26] = 2;
3555 else if (degree == 3)
3561 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3562 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3563 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3564 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3565 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3566 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3567 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3568 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3570 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3571 I[ 9] = 3; J[ 9] = 0; K[ 9] = 0;
3572 I[10] = 1; J[10] = 2; K[10] = 0;
3573 I[11] = 1; J[11] = 3; K[11] = 0;
3574 I[12] = 2; J[12] = 1; K[12] = 0;
3575 I[13] = 3; J[13] = 1; K[13] = 0;
3576 I[14] = 0; J[14] = 2; K[14] = 0;
3577 I[15] = 0; J[15] = 3; K[15] = 0;
3578 I[16] = 2; J[16] = 0; K[16] = 1;
3579 I[17] = 3; J[17] = 0; K[17] = 1;
3580 I[18] = 1; J[18] = 2; K[18] = 1;
3581 I[19] = 1; J[19] = 3; K[19] = 1;
3582 I[20] = 2; J[20] = 1; K[20] = 1;
3583 I[21] = 3; J[21] = 1; K[21] = 1;
3584 I[22] = 0; J[22] = 2; K[22] = 1;
3585 I[23] = 0; J[23] = 3; K[23] = 1;
3586 I[24] = 0; J[24] = 0; K[24] = 2;
3587 I[25] = 0; J[25] = 0; K[25] = 3;
3588 I[26] = 1; J[26] = 0; K[26] = 2;
3589 I[27] = 1; J[27] = 0; K[27] = 3;
3590 I[28] = 1; J[28] = 1; K[28] = 2;
3591 I[29] = 1; J[29] = 1; K[29] = 3;
3592 I[30] = 0; J[30] = 1; K[30] = 2;
3593 I[31] = 0; J[31] = 1; K[31] = 3;
3595 I[32] = 2; J[32] = 3; K[32] = 0;
3596 I[33] = 3; J[33] = 3; K[33] = 0;
3597 I[34] = 2; J[34] = 2; K[34] = 0;
3598 I[35] = 3; J[35] = 2; K[35] = 0;
3599 I[36] = 2; J[36] = 0; K[36] = 2;
3600 I[37] = 3; J[37] = 0; K[37] = 2;
3601 I[38] = 2; J[38] = 0; K[38] = 3;
3602 I[39] = 3; J[39] = 0; K[39] = 3;
3603 I[40] = 1; J[40] = 2; K[40] = 2;
3604 I[41] = 1; J[41] = 3; K[41] = 2;
3605 I[42] = 1; J[42] = 2; K[42] = 3;
3606 I[43] = 1; J[43] = 3; K[43] = 3;
3607 I[44] = 3; J[44] = 1; K[44] = 2;
3608 I[45] = 2; J[45] = 1; K[45] = 2;
3609 I[46] = 3; J[46] = 1; K[46] = 3;
3610 I[47] = 2; J[47] = 1; K[47] = 3;
3611 I[48] = 0; J[48] = 3; K[48] = 2;
3612 I[49] = 0; J[49] = 2; K[49] = 2;
3613 I[50] = 0; J[50] = 3; K[50] = 3;
3614 I[51] = 0; J[51] = 2; K[51] = 3;
3615 I[52] = 2; J[52] = 2; K[52] = 1;
3616 I[53] = 3; J[53] = 2; K[53] = 1;
3617 I[54] = 2; J[54] = 3; K[54] = 1;
3618 I[55] = 3; J[55] = 3; K[55] = 1;
3620 I[56] = 2; J[56] = 2; K[56] = 2;
3621 I[57] = 3; J[57] = 2; K[57] = 2;
3622 I[58] = 3; J[58] = 3; K[58] = 2;
3623 I[59] = 2; J[59] = 3; K[59] = 2;
3624 I[60] = 2; J[60] = 2; K[60] = 3;
3625 I[61] = 3; J[61] = 2; K[61] = 3;
3626 I[62] = 3; J[62] = 3; K[62] = 3;
3627 I[63] = 2; J[63] = 3; K[63] = 3;
3631 mfem_error (
"LagrangeHexFiniteElement::LagrangeHexFiniteElement");
3635 dof1d = fe1d ->
GetDof();
3637 #ifndef MFEM_THREAD_SAFE
3647 for (
int n = 0; n <
Dof; n++)
3662 #ifdef MFEM_THREAD_SAFE
3663 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
3670 for (
int n = 0; n <
Dof; n++)
3671 shape(n) = shape1dx(I[n]) * shape1dy(J[n]) * shape1dz(K[n]);
3681 #ifdef MFEM_THREAD_SAFE
3682 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
3683 DenseMatrix dshape1dx(dof1d,1), dshape1dy(dof1d,1), dshape1dz(dof1d,1);
3694 for (
int n = 0; n <
Dof; n++) {
3695 dshape(n,0) = dshape1dx(I[n],0) * shape1dy(J[n]) * shape1dz(K[n]);
3696 dshape(n,1) = shape1dx(I[n]) * dshape1dy(J[n],0) * shape1dz(K[n]);
3697 dshape(n,2) = shape1dx(I[n]) * shape1dy(J[n]) * dshape1dz(K[n],0);
3725 shape(0) = 1.0 - 2.0 * x;
3730 shape(1) = 2.0 * x - 1.0;
3731 shape(2) = 2.0 - 2.0 * x;
3741 dshape(0,0) = - 2.0;
3747 dshape(2,0) = - 2.0;
3774 L0 = 2.0 * ( 1. - ip.
x - ip.
y );
3775 L1 = 2.0 * ( ip.
x );
3776 L2 = 2.0 * ( ip.
y );
3785 for (i = 0; i < 6; i++)
3789 shape(0) = L0 - 1.0;
3793 else if (L1 >= 1.0) {
3795 shape(1) = L1 - 1.0;
3798 else if (L2 >= 1.0) {
3801 shape(2) = L2 - 1.0;
3804 shape(3) = 1.0 - L2;
3805 shape(4) = 1.0 - L0;
3806 shape(5) = 1.0 - L1;
3816 L0 = 2.0 * ( 1. - ip.
x - ip.
y );
3817 L1 = 2.0 * ( ip.
x );
3818 L2 = 2.0 * ( ip.
y );
3820 double DL0[2], DL1[2], DL2[2];
3821 DL0[0] = -2.0; DL0[1] = -2.0;
3822 DL1[0] = 2.0; DL1[1] = 0.0;
3823 DL2[0] = 0.0; DL2[1] = 2.0;
3825 for (i = 0; i < 6; i++)
3826 for (j = 0; j < 2; j++)
3830 for (j = 0; j < 2; j++) {
3831 dshape(0,j) = DL0[j];
3832 dshape(3,j) = DL1[j];
3833 dshape(5,j) = DL2[j];
3836 else if (L1 >= 1.0) {
3837 for (j = 0; j < 2; j++) {
3838 dshape(3,j) = DL0[j];
3839 dshape(1,j) = DL1[j];
3840 dshape(4,j) = DL2[j];
3843 else if (L2 >= 1.0) {
3844 for (j = 0; j < 2; j++) {
3845 dshape(5,j) = DL0[j];
3846 dshape(4,j) = DL1[j];
3847 dshape(2,j) = DL2[j];
3851 for (j = 0; j < 2; j++) {
3852 dshape(3,j) = - DL2[j];
3853 dshape(4,j) = - DL0[j];
3854 dshape(5,j) = - DL1[j];
3899 double L0, L1, L2, L3, L4, L5;
3900 L0 = 2.0 * ( 1. - ip.
x - ip.
y - ip.
z );
3901 L1 = 2.0 * ( ip.
x );
3902 L2 = 2.0 * ( ip.
y );
3903 L3 = 2.0 * ( ip.
z );
3904 L4 = 2.0 * ( ip.
x + ip.
y );
3905 L5 = 2.0 * ( ip.
y + ip.
z );
3918 for (i = 0; i < 10; i++)
3922 shape(0) = L0 - 1.0;
3927 else if (L1 >= 1.0) {
3929 shape(1) = L1 - 1.0;
3933 else if (L2 >= 1.0) {
3936 shape(2) = L2 - 1.0;
3939 else if (L3 >= 1.0) {
3943 shape(3) = L3 - 1.0;
3945 else if ((L4 <= 1.0) && (L5 <= 1.0)) {
3946 shape(4) = 1.0 - L5;
3948 shape(6) = 1.0 - L4;
3949 shape(8) = 1.0 - L0;
3951 else if ((L4 >= 1.0) && (L5 <= 1.0)) {
3952 shape(4) = 1.0 - L5;
3953 shape(5) = 1.0 - L1;
3954 shape(7) = L4 - 1.0;
3957 else if ((L4 <= 1.0) && (L5 >= 1.0)) {
3958 shape(5) = 1.0 - L3;
3959 shape(6) = 1.0 - L4;
3961 shape(9) = L5 - 1.0;
3963 else if ((L4 >= 1.0) && (L5 >= 1.0)) {
3965 shape(7) = L4 - 1.0;
3966 shape(8) = 1.0 - L2;
3967 shape(9) = L5 - 1.0;
3976 double L0, L1, L2, L3, L4, L5;
3977 L0 = 2.0 * ( 1. - ip.
x - ip.
y - ip.
z );
3978 L1 = 2.0 * ( ip.
x );
3979 L2 = 2.0 * ( ip.
y );
3980 L3 = 2.0 * ( ip.
z );
3981 L4 = 2.0 * ( ip.
x + ip.
y );
3982 L5 = 2.0 * ( ip.
y + ip.
z );
3984 double DL0[3], DL1[3], DL2[3], DL3[3], DL4[3], DL5[3];
3985 DL0[0] = -2.0; DL0[1] = -2.0; DL0[2] = -2.0;
3986 DL1[0] = 2.0; DL1[1] = 0.0; DL1[2] = 0.0;
3987 DL2[0] = 0.0; DL2[1] = 2.0; DL2[2] = 0.0;
3988 DL3[0] = 0.0; DL3[1] = 0.0; DL3[2] = 2.0;
3989 DL4[0] = 2.0; DL4[1] = 2.0; DL4[2] = 0.0;
3990 DL5[0] = 0.0; DL5[1] = 2.0; DL5[2] = 2.0;
3992 for (i = 0; i < 10; i++)
3993 for (j = 0; j < 3; j++)
3997 for (j = 0; j < 3; j++) {
3998 dshape(0,j) = DL0[j];
3999 dshape(4,j) = DL1[j];
4000 dshape(5,j) = DL2[j];
4001 dshape(6,j) = DL3[j];
4004 else if (L1 >= 1.0) {
4005 for (j = 0; j < 3; j++) {
4006 dshape(4,j) = DL0[j];
4007 dshape(1,j) = DL1[j];
4008 dshape(7,j) = DL2[j];
4009 dshape(8,j) = DL3[j];
4012 else if (L2 >= 1.0) {
4013 for (j = 0; j < 3; j++) {
4014 dshape(5,j) = DL0[j];
4015 dshape(7,j) = DL1[j];
4016 dshape(2,j) = DL2[j];
4017 dshape(9,j) = DL3[j];
4020 else if (L3 >= 1.0) {
4021 for (j = 0; j < 3; j++) {
4022 dshape(6,j) = DL0[j];
4023 dshape(8,j) = DL1[j];
4024 dshape(9,j) = DL2[j];
4025 dshape(3,j) = DL3[j];
4028 else if ((L4 <= 1.0) && (L5 <= 1.0)) {
4029 for (j = 0; j < 3; j++) {
4030 dshape(4,j) = - DL5[j];
4031 dshape(5,j) = DL2[j];
4032 dshape(6,j) = - DL4[j];
4033 dshape(8,j) = - DL0[j];
4036 else if ((L4 >= 1.0) && (L5 <= 1.0)) {
4037 for (j = 0; j < 3; j++) {
4038 dshape(4,j) = - DL5[j];
4039 dshape(5,j) = - DL1[j];
4040 dshape(7,j) = DL4[j];
4041 dshape(8,j) = DL3[j];
4044 else if ((L4 <= 1.0) && (L5 >= 1.0)) {
4045 for (j = 0; j < 3; j++) {
4046 dshape(5,j) = - DL3[j];
4047 dshape(6,j) = - DL4[j];
4048 dshape(8,j) = DL1[j];
4049 dshape(9,j) = DL5[j];
4052 else if ((L4 >= 1.0) && (L5 >= 1.0)) {
4053 for (j = 0; j < 3; j++) {
4054 dshape(5,j) = DL0[j];
4055 dshape(7,j) = DL4[j];
4056 dshape(8,j) = - DL2[j];
4057 dshape(9,j) = DL5[j];
4090 double x = ip.
x, y = ip.
y;
4092 Lx = 2.0 * ( 1. - x );
4093 Ly = 2.0 * ( 1. - y );
4102 for (i = 0; i < 9; i++)
4105 if ((x <= 0.5) && (y <= 0.5)) {
4106 shape(0) = (Lx - 1.0) * (Ly - 1.0);
4107 shape(4) = (2.0 - Lx) * (Ly - 1.0);
4108 shape(8) = (2.0 - Lx) * (2.0 - Ly);
4109 shape(7) = (Lx - 1.0) * (2.0 - Ly);
4111 else if ((x >= 0.5) && (y <= 0.5)) {
4112 shape(4) = Lx * (Ly - 1.0);
4113 shape(1) = (1.0 - Lx) * (Ly - 1.0);
4114 shape(5) = (1.0 - Lx) * (2.0 - Ly);
4115 shape(8) = Lx * (2.0 - Ly);
4117 else if ((x >= 0.5) && (y >= 0.5)) {
4118 shape(8) = Lx * Ly ;
4119 shape(5) = (1.0 - Lx) * Ly ;
4120 shape(2) = (1.0 - Lx) * (1.0 - Ly);
4121 shape(6) = Lx * (1.0 - Ly);
4123 else if ((x <= 0.5) && (y >= 0.5)) {
4124 shape(7) = (Lx - 1.0) * Ly ;
4125 shape(8) = (2.0 - Lx) * Ly ;
4126 shape(6) = (2.0 - Lx) * (1.0 - Ly);
4127 shape(3) = (Lx - 1.0) * (1.0 - Ly);
4135 double x = ip.
x, y = ip.
y;
4137 Lx = 2.0 * ( 1. - x );
4138 Ly = 2.0 * ( 1. - y );
4140 for (i = 0; i < 9; i++)
4141 for (j = 0; j < 2; j++)
4144 if ((x <= 0.5) && (y <= 0.5)) {
4145 dshape(0,0) = 2.0 * (1.0 - Ly);
4146 dshape(0,1) = 2.0 * (1.0 - Lx);
4148 dshape(4,0) = 2.0 * (Ly - 1.0);
4149 dshape(4,1) = -2.0 * (2.0 - Lx);
4151 dshape(8,0) = 2.0 * (2.0 - Ly);
4152 dshape(8,1) = 2.0 * (2.0 - Lx);
4154 dshape(7,0) = -2.0 * (2.0 - Ly);
4155 dshape(7,0) = 2.0 * (Lx - 1.0);
4157 else if ((x >= 0.5) && (y <= 0.5)) {
4158 dshape(4,0) = -2.0 * (Ly - 1.0);
4159 dshape(4,1) = -2.0 * Lx;
4161 dshape(1,0) = 2.0 * (Ly - 1.0);
4162 dshape(1,1) = -2.0 * (1.0 - Lx);
4164 dshape(5,0) = 2.0 * (2.0 - Ly);
4165 dshape(5,1) = 2.0 * (1.0 - Lx);
4167 dshape(8,0) = -2.0 * (2.0 - Ly);
4168 dshape(8,1) = 2.0 * Lx;
4170 else if ((x >= 0.5) && (y >= 0.5)) {
4171 dshape(8,0) = -2.0 * Ly;
4172 dshape(8,1) = -2.0 * Lx;
4174 dshape(5,0) = 2.0 * Ly;
4175 dshape(5,1) = -2.0 * (1.0 - Lx);
4177 dshape(2,0) = 2.0 * (1.0 - Ly);
4178 dshape(2,1) = 2.0 * (1.0 - Lx);
4180 dshape(6,0) = -2.0 * (1.0 - Ly);
4181 dshape(6,1) = 2.0 * Lx;
4183 else if ((x <= 0.5) && (y >= 0.5)) {
4184 dshape(7,0) = -2.0 * Ly;
4185 dshape(7,1) = -2.0 * (Lx - 1.0);
4187 dshape(8,0) = 2.0 * Ly ;
4188 dshape(8,1) = -2.0 * (2.0 - Lx);
4190 dshape(6,0) = 2.0 * (1.0 - Ly);
4191 dshape(6,1) = 2.0 * (2.0 - Lx);
4193 dshape(3,0) = -2.0 * (1.0 - Ly);
4194 dshape(3,1) = 2.0 * (Lx - 1.0);
4205 I[ 0] = 0.0; J[ 0] = 0.0; K[ 0] = 0.0;
4206 I[ 1] = 1.0; J[ 1] = 0.0; K[ 1] = 0.0;
4207 I[ 2] = 1.0; J[ 2] = 1.0; K[ 2] = 0.0;
4208 I[ 3] = 0.0; J[ 3] = 1.0; K[ 3] = 0.0;
4209 I[ 4] = 0.0; J[ 4] = 0.0; K[ 4] = 1.0;
4210 I[ 5] = 1.0; J[ 5] = 0.0; K[ 5] = 1.0;
4211 I[ 6] = 1.0; J[ 6] = 1.0; K[ 6] = 1.0;
4212 I[ 7] = 0.0; J[ 7] = 1.0; K[ 7] = 1.0;
4214 I[ 8] = 0.5; J[ 8] = 0.0; K[ 8] = 0.0;
4215 I[ 9] = 1.0; J[ 9] = 0.5; K[ 9] = 0.0;
4216 I[10] = 0.5; J[10] = 1.0; K[10] = 0.0;
4217 I[11] = 0.0; J[11] = 0.5; K[11] = 0.0;
4218 I[12] = 0.5; J[12] = 0.0; K[12] = 1.0;
4219 I[13] = 1.0; J[13] = 0.5; K[13] = 1.0;
4220 I[14] = 0.5; J[14] = 1.0; K[14] = 1.0;
4221 I[15] = 0.0; J[15] = 0.5; K[15] = 1.0;
4222 I[16] = 0.0; J[16] = 0.0; K[16] = 0.5;
4223 I[17] = 1.0; J[17] = 0.0; K[17] = 0.5;
4224 I[18] = 1.0; J[18] = 1.0; K[18] = 0.5;
4225 I[19] = 0.0; J[19] = 1.0; K[19] = 0.5;
4227 I[20] = 0.5; J[20] = 0.5; K[20] = 0.0;
4228 I[21] = 0.5; J[21] = 0.0; K[21] = 0.5;
4229 I[22] = 1.0; J[22] = 0.5; K[22] = 0.5;
4230 I[23] = 0.5; J[23] = 1.0; K[23] = 0.5;
4231 I[24] = 0.0; J[24] = 0.5; K[24] = 0.5;
4232 I[25] = 0.5; J[25] = 0.5; K[25] = 1.0;
4234 I[26] = 0.5; J[26] = 0.5; K[26] = 0.5;
4236 for (
int n = 0; n < 27; n++) {
4248 double x = ip.
x, y = ip.
y, z = ip.
z;
4250 for (i = 0; i < 27; i++)
4253 if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5)) {
4267 else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5)) {
4281 else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5)) {
4295 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5)) {
4309 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5)) {
4323 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5)) {
4337 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5)) {
4366 shape(N[0]) = Lx * Ly * Lz;
4367 shape(N[1]) = (1 - Lx) * Ly * Lz;
4368 shape(N[2]) = (1 - Lx) * (1 - Ly) * Lz;
4369 shape(N[3]) = Lx * (1 - Ly) * Lz;
4370 shape(N[4]) = Lx * Ly * (1 - Lz);
4371 shape(N[5]) = (1 - Lx) * Ly * (1 - Lz);
4372 shape(N[6]) = (1 - Lx) * (1 - Ly) * (1 - Lz);
4373 shape(N[7]) = Lx * (1 - Ly) * (1 - Lz);
4381 double x = ip.
x, y = ip.
y, z = ip.
z;
4383 for (i = 0; i < 27; i++)
4384 for (j = 0; j < 3; j++)
4387 if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5)) {
4401 else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5)) {
4415 else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5)) {
4429 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5)) {
4443 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5)) {
4457 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5)) {
4471 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5)) {
4500 dshape(N[0],0) = -2.0 * Ly * Lz ;
4501 dshape(N[0],1) = -2.0 * Lx * Lz ;
4502 dshape(N[0],2) = -2.0 * Lx * Ly ;
4504 dshape(N[1],0) = 2.0 * Ly * Lz ;
4505 dshape(N[1],1) = -2.0 * (1 - Lx) * Lz ;
4506 dshape(N[1],2) = -2.0 * (1 - Lx) * Ly ;
4508 dshape(N[2],0) = 2.0 * (1 - Ly) * Lz ;
4509 dshape(N[2],1) = 2.0 * (1 - Lx) * Lz ;
4510 dshape(N[2],2) = -2.0 * (1 - Lx) * (1 - Ly);
4512 dshape(N[3],0) = -2.0 * (1 - Ly) * Lz ;
4513 dshape(N[3],1) = 2.0 * Lx * Lz ;
4514 dshape(N[3],2) = -2.0 * Lx * (1 - Ly);
4516 dshape(N[4],0) = -2.0 * Ly * (1 - Lz);
4517 dshape(N[4],1) = -2.0 * Lx * (1 - Lz);
4518 dshape(N[4],2) = 2.0 * Lx * Ly ;
4520 dshape(N[5],0) = 2.0 * Ly * (1 - Lz);
4521 dshape(N[5],1) = -2.0 * (1 - Lx) * (1 - Lz);
4522 dshape(N[5],2) = 2.0 * (1 - Lx) * Ly ;
4524 dshape(N[6],0) = 2.0 * (1 - Ly) * (1 - Lz);
4525 dshape(N[6],1) = 2.0 * (1 - Lx) * (1 - Lz);
4526 dshape(N[6],2) = 2.0 * (1 - Lx) * (1 - Ly);
4528 dshape(N[7],0) = -2.0 * (1 - Ly) * (1 - Lz);
4529 dshape(N[7],1) = 2.0 * Lx * (1 - Lz);
4530 dshape(N[7],2) = 2.0 * Lx * (1 - Ly);
4590 double x = ip.
x, y = ip.
y, z = ip.
z;
4592 shape(0,0) = (1. - y) * (1. - z);
4596 shape(2,0) = y * (1. - z);
4600 shape(4,0) = z * (1. - y);
4609 shape(1,1) = x * (1. - z);
4613 shape(3,1) = (1. - x) * (1. - z);
4621 shape(7,1) = (1. - x) * z;
4626 shape(8,2) = (1. - x) * (1. - y);
4630 shape(9,2) = x * (1. - y);
4634 shape(10,2) = x * y;
4638 shape(11,2) = y * (1. - x);
4646 double x = ip.
x, y = ip.
y, z = ip.
z;
4648 curl_shape(0,0) = 0.;
4649 curl_shape(0,1) = y - 1.;
4650 curl_shape(0,2) = 1. - z;
4652 curl_shape(2,0) = 0.;
4653 curl_shape(2,1) = -y;
4654 curl_shape(2,2) = z - 1.;
4656 curl_shape(4,0) = 0;
4657 curl_shape(4,1) = 1. - y;
4658 curl_shape(4,2) = z;
4660 curl_shape(6,0) = 0.;
4661 curl_shape(6,1) = y;
4662 curl_shape(6,2) = -z;
4664 curl_shape(1,0) = x;
4665 curl_shape(1,1) = 0.;
4666 curl_shape(1,2) = 1. - z;
4668 curl_shape(3,0) = 1. - x;
4669 curl_shape(3,1) = 0.;
4670 curl_shape(3,2) = z - 1.;
4672 curl_shape(5,0) = -x;
4673 curl_shape(5,1) = 0.;
4674 curl_shape(5,2) = z;
4676 curl_shape(7,0) = x - 1.;
4677 curl_shape(7,1) = 0.;
4678 curl_shape(7,2) = -z;
4680 curl_shape(8,0) = x - 1.;
4681 curl_shape(8,1) = 1. - y;
4682 curl_shape(8,2) = 0.;
4684 curl_shape(9,0) = -x;
4685 curl_shape(9,1) = y - 1.;
4686 curl_shape(9,2) = 0;
4688 curl_shape(10,0) = x;
4689 curl_shape(10,1) = -y;
4690 curl_shape(10,2) = 0.;
4692 curl_shape(11,0) = 1. - x;
4693 curl_shape(11,1) = y;
4694 curl_shape(11,2) = 0.;
4697 const double Nedelec1HexFiniteElement::tk[12][3] =
4698 {{1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
4699 {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
4700 {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}};
4706 #ifdef MFEM_THREAD_SAFE
4711 for (k = 0; k < 12; k++)
4714 for (j = 0; j < 12; j++)
4716 double d = (
vshape(j,0)*tk[k][0] +
vshape(j,1)*tk[k][1] +
4718 if (j == k) d -= 1.0;
4719 if (fabs(d) > 1.0e-12)
4721 cerr <<
"Nedelec1HexFiniteElement::GetLocalInterpolation (...)\n"
4722 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
4730 ip.
x = ip.
y = ip.
z = 0.0;
4737 for (k = 0; k < 12; k++)
4740 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
4743 vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
4744 vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
4745 vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
4746 for (j = 0; j < 12; j++)
4747 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
4748 vshape(j,2)*vk[2])) < 1.0e-12)
4760 for (
int k = 0; k < 12; k++)
4768 vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
4769 vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
4770 vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
4807 double x = ip.
x, y = ip.
y, z = ip.
z;
4809 shape(0,0) = 1. - y - z;
4814 shape(1,1) = 1. - x - z;
4819 shape(2,2) = 1. - x - y;
4838 curl_shape(0,0) = 0.;
4839 curl_shape(0,1) = -2.;
4840 curl_shape(0,2) = 2.;
4842 curl_shape(1,0) = 2.;
4843 curl_shape(1,1) = 0.;
4844 curl_shape(1,2) = -2.;
4846 curl_shape(2,0) = -2.;
4847 curl_shape(2,1) = 2.;
4848 curl_shape(2,2) = 0.;
4850 curl_shape(3,0) = 0.;
4851 curl_shape(3,1) = 0.;
4852 curl_shape(3,2) = 2.;
4854 curl_shape(4,0) = 0.;
4855 curl_shape(4,1) = -2.;
4856 curl_shape(4,2) = 0.;
4858 curl_shape(5,0) = 2.;
4859 curl_shape(5,1) = 0.;
4860 curl_shape(5,2) = 0.;
4863 const double Nedelec1TetFiniteElement::tk[6][3] =
4864 {{1,0,0}, {0,1,0}, {0,0,1}, {-1,1,0}, {-1,0,1}, {0,-1,1}};
4870 #ifdef MFEM_THREAD_SAFE
4875 for (k = 0; k < 6; k++)
4878 for (j = 0; j < 6; j++)
4880 double d = (
vshape(j,0)*tk[k][0] +
vshape(j,1)*tk[k][1] +
4882 if (j == k) d -= 1.0;
4883 if (fabs(d) > 1.0e-12)
4885 cerr <<
"Nedelec1TetFiniteElement::GetLocalInterpolation (...)\n"
4886 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
4894 ip.
x = ip.
y = ip.
z = 0.0;
4901 for (k = 0; k < 6; k++)
4904 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
4907 vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
4908 vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
4909 vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
4910 for (j = 0; j < 6; j++)
4911 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
4912 vshape(j,2)*vk[2])) < 1.0e-12)
4924 for (
int k = 0; k < 6; k++)
4932 vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
4933 vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
4934 vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
4971 double x = ip.
x, y = ip.
y, z = ip.
z;
4975 shape(0,2) = z - 1.;
4978 shape(1,1) = y - 1.;
4989 shape(4,0) = x - 1.;
5009 const double RT0HexFiniteElement::nk[6][3] =
5010 {{0,0,-1}, {0,-1,0}, {1,0,0}, {0,1,0}, {-1,0,0}, {0,0,1}};
5016 #ifdef MFEM_THREAD_SAFE
5022 for (k = 0; k < 6; k++)
5025 for (j = 0; j < 6; j++)
5027 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5029 if (j == k) d -= 1.0;
5030 if (fabs(d) > 1.0e-12)
5032 cerr <<
"RT0HexFiniteElement::GetLocalInterpolation (...)\n"
5033 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5041 ip.
x = ip.
y = ip.
z = 0.0;
5049 for (k = 0; k < 6; k++)
5052 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5055 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2];
5056 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2];
5057 vk[2] =
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2];
5058 for (j = 0; j < 6; j++)
5059 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5060 vshape(j,2)*vk[2])) < 1.0e-12)
5071 #ifdef MFEM_THREAD_SAFE
5075 for (
int k = 0; k < 6; k++)
5084 vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2] ) +
5085 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2] ) +
5086 vk[2] * (
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2] );
5215 double x = ip.
x, y = ip.
y, z = ip.
z;
5219 shape(2,2) = -(1. - 3.*z + 2.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5222 shape(3,2) = -(1. - 3.*z + 2.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5225 shape(0,2) = -(1. - 3.*z + 2.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5228 shape(1,2) = -(1. - 3.*z + 2.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5231 shape(4,1) = -(1. - 3.*y + 2.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5234 shape(5,1) = -(1. - 3.*y + 2.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5237 shape(6,1) = -(1. - 3.*y + 2.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5240 shape(7,1) = -(1. - 3.*y + 2.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5243 shape(8,0) = (-x + 2.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5246 shape(9,0) = (-x + 2.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5249 shape(10,0) = (-x + 2.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5252 shape(11,0) = (-x + 2.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5257 shape(13,1) = (-y + 2.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5260 shape(12,1) = (-y + 2.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5263 shape(15,1) = (-y + 2.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5266 shape(14,1) = (-y + 2.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5269 shape(17,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5272 shape(16,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5275 shape(19,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5278 shape(18,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5284 shape(20,2) = (-z + 2.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5287 shape(21,2) = (-z + 2.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5290 shape(22,2) = (-z + 2.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5293 shape(23,2) = (-z + 2.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5295 shape(24,0) = (4.*x - 4.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5298 shape(25,0) = (4.*x - 4.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5301 shape(26,0) = (4.*x - 4.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5304 shape(27,0) = (4.*x - 4.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5309 shape(28,1) = (4.*y - 4.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5312 shape(29,1) = (4.*y - 4.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5315 shape(30,1) = (4.*y - 4.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5318 shape(31,1) = (4.*y - 4.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5323 shape(32,2) = (4.*z - 4.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5326 shape(33,2) = (4.*z - 4.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5329 shape(34,2) = (4.*z - 4.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5332 shape(35,2) = (4.*z - 4.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5338 double x = ip.
x, y = ip.
y, z = ip.
z;
5340 divshape(2) = -(-3. + 4.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5341 divshape(3) = -(-3. + 4.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5342 divshape(0) = -(-3. + 4.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5343 divshape(1) = -(-3. + 4.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5345 divshape(4) = -(-3. + 4.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5346 divshape(5) = -(-3. + 4.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5347 divshape(6) = -(-3. + 4.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5348 divshape(7) = -(-3. + 4.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5350 divshape(8) = (-1. + 4.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5351 divshape(9) = (-1. + 4.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5352 divshape(10) = (-1. + 4.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5353 divshape(11) = (-1. + 4.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5355 divshape(13) = (-1. + 4.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5356 divshape(12) = (-1. + 4.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5357 divshape(15) = (-1. + 4.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5358 divshape(14) = (-1. + 4.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5360 divshape(17) = -(-3. + 4.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5361 divshape(16) = -(-3. + 4.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5362 divshape(19) = -(-3. + 4.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5363 divshape(18) = -(-3. + 4.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5365 divshape(20) = (-1. + 4.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5366 divshape(21) = (-1. + 4.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5367 divshape(22) = (-1. + 4.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5368 divshape(23) = (-1. + 4.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5370 divshape(24) = ( 4. - 8.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5371 divshape(25) = ( 4. - 8.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5372 divshape(26) = ( 4. - 8.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5373 divshape(27) = ( 4. - 8.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5375 divshape(28) = ( 4. - 8.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5376 divshape(29) = ( 4. - 8.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5377 divshape(30) = ( 4. - 8.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5378 divshape(31) = ( 4. - 8.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5380 divshape(32) = ( 4. - 8.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5381 divshape(33) = ( 4. - 8.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5382 divshape(34) = ( 4. - 8.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5383 divshape(35) = ( 4. - 8.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5386 const double RT1HexFiniteElement::nk[36][3] =
5388 {0, 0,-1}, {0, 0,-1}, {0, 0,-1}, {0, 0,-1},
5389 {0,-1, 0}, {0,-1, 0}, {0,-1, 0}, {0,-1, 0},
5390 {1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0},
5391 {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0},
5392 {-1,0, 0}, {-1,0, 0}, {-1,0, 0}, {-1,0, 0},
5393 {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1},
5394 {1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0},
5395 {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0},
5396 {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1}
5403 #ifdef MFEM_THREAD_SAFE
5409 for (k = 0; k < 36; k++)
5412 for (j = 0; j < 36; j++)
5414 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5416 if (j == k) d -= 1.0;
5417 if (fabs(d) > 1.0e-12)
5419 cerr <<
"RT0HexFiniteElement::GetLocalInterpolation (...)\n"
5420 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5428 ip.
x = ip.
y = ip.
z = 0.0;
5436 for (k = 0; k < 36; k++)
5439 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5442 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2];
5443 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2];
5444 vk[2] =
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2];
5445 for (j = 0; j < 36; j++)
5446 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5447 vshape(j,2)*vk[2])) < 1.0e-12)
5458 #ifdef MFEM_THREAD_SAFE
5462 for (
int k = 0; k < 36; k++)
5471 vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2] ) +
5472 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2] ) +
5473 vk[2] * (
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2] );
5501 double x2 = 2.0*ip.
x, y2 = 2.0*ip.
y, z2 = 2.0*ip.
z;
5507 shape(1,0) = x2 - 2.0;
5512 shape(2,1) = y2 - 2.0;
5517 shape(3,2) = z2 - 2.0;
5529 const double RT0TetFiniteElement::nk[4][3] =
5530 {{.5,.5,.5}, {-.5,0,0}, {0,-.5,0}, {0,0,-.5}};
5536 #ifdef MFEM_THREAD_SAFE
5542 for (k = 0; k < 4; k++)
5545 for (j = 0; j < 4; j++)
5547 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5549 if (j == k) d -= 1.0;
5550 if (fabs(d) > 1.0e-12)
5552 cerr <<
"RT0TetFiniteElement::GetLocalInterpolation (...)\n"
5553 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5561 ip.
x = ip.
y = ip.
z = 0.0;
5569 for (k = 0; k < 4; k++)
5572 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5575 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2];
5576 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2];
5577 vk[2] =
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2];
5578 for (j = 0; j < 4; j++)
5579 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5580 vshape(j,2)*vk[2])) < 1.0e-12)
5591 #ifdef MFEM_THREAD_SAFE
5595 for (
int k = 0; k < 4; k++)
5604 vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2] ) +
5605 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2] ) +
5606 vk[2] * (
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2] );
5641 double x = 2. * ip.
x - 1.;
5642 double y = 2. * ip.
y - 1.;
5643 double z = 2. * ip.
z - 1.;
5644 double f5 = x * x - y * y;
5645 double f6 = y * y - z * z;
5647 shape(0) = (1./6.) * (1. - 3. * z - f5 - 2. * f6);
5648 shape(1) = (1./6.) * (1. - 3. * y - f5 + f6);
5649 shape(2) = (1./6.) * (1. + 3. * x + 2. * f5 + f6);
5650 shape(3) = (1./6.) * (1. + 3. * y - f5 + f6);
5651 shape(4) = (1./6.) * (1. - 3. * x + 2. * f5 + f6);
5652 shape(5) = (1./6.) * (1. + 3. * z - f5 - 2. * f6);
5658 const double a = 2./3.;
5660 double xt = a * (1. - 2. * ip.
x);
5661 double yt = a * (1. - 2. * ip.
y);
5662 double zt = a * (1. - 2. * ip.
z);
5666 dshape(0,2) = -1. - 2. * zt;
5669 dshape(1,1) = -1. - 2. * yt;
5672 dshape(2,0) = 1. - 2. * xt;
5677 dshape(3,1) = 1. - 2. * yt;
5680 dshape(4,0) = -1. - 2. * xt;
5686 dshape(5,2) = 1. - 2. * zt;
5691 : x(p + 1), w(p + 1)
5697 for (
int i = 0; i <= p; i++)
5700 for (
int j = 0; j <= p; j++)
5711 for (
int i = 0; i <= p; i++)
5712 for (
int j = 0; j < i; j++)
5714 double xij = x(i) - x(j);
5718 for (
int i = 0; i <= p; i++)
5723 for (
int i = 0; i < p; i++)
5725 mfem_error(
"Poly_1D::Basis::Basis : nodes are not increasing!");
5739 int i, k, p = x.Size() - 1;
5749 for (k = 0; k < p; k++)
5750 if (y >= (x(k) + x(k+1))/2)
5754 for (i = k+1; i <= p; i++)
5758 l = lk * (y - x(k));
5760 for (i = 0; i < k; i++)
5761 u(i) = l * w(i) / (y - x(i));
5763 for (i++; i <= p; i++)
5764 u(i) = l * w(i) / (y - x(i));
5778 int i, k, p = x.Size() - 1;
5779 double l, lp, lk, sk, si;
5789 for (k = 0; k < p; k++)
5790 if (y >= (x(k) + x(k+1))/2)
5794 for (i = k+1; i <= p; i++)
5798 l = lk * (y - x(k));
5801 for (i = 0; i < k; i++)
5803 si = 1.0/(y - x(i));
5805 u(i) = l * si * w(i);
5808 for (i++; i <= p; i++)
5810 si = 1.0/(y - x(i));
5812 u(i) = l * si * w(i);
5816 for (i = 0; i < k; i++)
5817 d(i) = (lp * w(i) - u(i))/(y - x(i));
5819 for (i++; i <= p; i++)
5820 d(i) = (lp * w(i) - u(i))/(y - x(i));
5829 for (
int i = 0; i <= p; i++)
5831 binom(i,0) = binom(i,i) = 1;
5832 for (
int j = 1; j < i; j++)
5833 binom(i,j) = binom(i-1,j) + binom(i-1,j-1);
5847 for (
int i = 0; i <= p; i++)
5854 int m = (p+1)/2, odd_p = p%2;
5859 for (
int i = 0; i < m; i++)
5861 double z, y, p0, d0;
5862 z = cos(M_PI*(i + 0.75)/(p + 1.5));
5864 for (
int k = 0;
true; )
5874 for (
int n = 1;
true; n++)
5876 p0 = ((2*n+1)*z*p1 - n*p2)/(n + 1);
5886 if (fabs(p0/d0) < 2e-16)
break;
5890 std::cout <<
"Poly_1D::GaussPoints : No convergence!"
5891 " p = " << p <<
", i = " << i <<
", p0/d0 = " << p0/d0
5900 y = ((1 - z) + p0/d0)/2;
5921 int m = (p - 1)/2, odd_p = p%2;
5925 for (
int i = 0; i < m; )
5927 double y, z, d0, s0;
5928 z = cos(M_PI*(i + 1)/p);
5930 for (
int k = 0;
true; )
5937 double p0, p1, p2, d1;
5943 for (
int n = 1;
true; n++)
5945 p0 = ((2*n+1)*z*p1 - n*p2)/(n + 1);
5955 if (n == p - 1)
break;
5961 if (fabs(d0/s0) < 2e-16)
break;
5965 std::cout <<
"Poly_1D::GaussLobattoPoints : No convergence!"
5966 " p = " << p <<
", i = " << i <<
", d0/s0 = " << d0/s0
5974 y = ((1 - z) + d0/s0)/2;
5984 for (
int i = 0; i <= p; i++)
5987 double s = sin(M_PI_2*(i + 0.5)/(p + 1));
5992 void Poly_1D::CalcMono(
const int p,
const double x,
double *u)
5996 for (
int n = 1; n <= p; n++)
6000 void Poly_1D::CalcMono(
const int p,
const double x,
double *u,
double *d)
6005 for (
int n = 1; n <= p; n++)
6022 const int *b =
Binom(p);
6025 for (i = 1; i < p; i++)
6032 for (i--; i > 0; i--)
6042 double *u,
double *d)
6052 const int *b =
Binom(p);
6053 const double xpy = x + y, ptx = p*x;
6056 for (i = 1; i < p; i++)
6058 d[i] = b[i]*z*(i*xpy - ptx);
6065 for (i--; i > 0; i--)
6086 const int *b =
Binom(p);
6087 const double xpy = x + y, ptx = p*x;
6090 for (i = 1; i < p; i++)
6092 d[i] = b[i]*z*(i*xpy - ptx);
6097 for (i--; i > 0; i--)
6106 void Poly_1D::CalcLegendre(
const int p,
const double x,
double *u)
6113 u[1] = z = 2.*x - 1.;
6114 for (
int n = 1; n < p; n++)
6116 u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
6120 void Poly_1D::CalcLegendre(
const int p,
const double x,
double *u,
double *d)
6130 u[1] = z = 2.*x - 1.;
6132 for (
int n = 1; n < p; n++)
6134 u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
6135 d[n+1] = (4*n + 2)*u[n] + d[n-1];
6139 void Poly_1D::CalcChebyshev(
const int p,
const double x,
double *u)
6147 u[1] = z = 2.*x - 1.;
6148 for (
int n = 1; n < p; n++)
6150 u[n+1] = 2*z*u[n] - u[n-1];
6154 void Poly_1D::CalcChebyshev(
const int p,
const double x,
double *u,
double *d)
6168 u[1] = z = 2.*x - 1.;
6170 for (
int n = 1; n < p; n++)
6172 u[n+1] = 2*z*u[n] - u[n-1];
6173 d[n+1] = (n + 1)*(z*d[n]/n + 2*u[n]);
6181 if (open_pts.
Size() <= p)
6183 int i = open_pts.
Size();
6189 if ((op = open_pts[p]) != NULL)
6192 open_pts[p] = op =
new double[p + 1];
6202 if (closed_pts.
Size() <= p)
6204 int i = closed_pts.
Size();
6207 closed_pts[i] = NULL;
6210 if ((cp = closed_pts[p]) != NULL)
6213 closed_pts[p] = cp =
new double[p + 1];
6223 if (open_basis.Size() <= p)
6225 int i = open_basis.Size();
6226 open_basis.SetSize(p + 1);
6228 open_basis[i] = NULL;
6231 if ((ob = open_basis[p]) != NULL)
6242 if (closed_basis.Size() <= p)
6244 int i = closed_basis.Size();
6245 closed_basis.SetSize(p + 1);
6247 closed_basis[i] = NULL;
6250 if ((cb = closed_basis[p]) != NULL)
6259 for (
int i = 0; i < open_pts.
Size(); i++)
6260 delete [] open_pts[i];
6261 for (
int i = 0; i < closed_pts.
Size(); i++)
6262 delete [] closed_pts[i];
6263 for (
int i = 0; i < open_basis.Size(); i++)
6264 delete open_basis[i];
6265 for (
int i = 0; i < closed_basis.Size(); i++)
6266 delete closed_basis[i];
6279 #ifndef MFEM_THREAD_SAFE
6286 for (
int i = 1; i < p; i++)
6293 const int p =
Order;
6295 #ifdef MFEM_THREAD_SAFE
6299 basis1d.
Eval(ip.
x, shape_x);
6301 shape(0) = shape_x(0);
6302 shape(1) = shape_x(p);
6303 for (
int i = 1; i < p; i++)
6304 shape(i+1) = shape_x(i);
6310 const int p =
Order;
6312 #ifdef MFEM_THREAD_SAFE
6313 Vector shape_x(p+1), dshape_x(p+1);
6316 basis1d.
Eval(ip.
x, shape_x, dshape_x);
6318 dshape(0,0) = dshape_x(0);
6319 dshape(1,0) = dshape_x(p);
6320 for (
int i = 1; i < p; i++)
6321 dshape(i+1,0) = dshape_x(i);
6326 const int p =
Order;
6334 for (
int i = 1; i < p; i++)
6341 for (
int i = 1; i < p; i++)
6351 basis1d(
poly1d.ClosedBasis(p)), dof_map((p + 1)*(p + 1))
6355 const int p1 = p + 1;
6357 #ifndef MFEM_THREAD_SAFE
6365 dof_map[0 + 0*p1] = 0;
6366 dof_map[p + 0*p1] = 1;
6367 dof_map[p + p*p1] = 2;
6368 dof_map[0 + p*p1] = 3;
6372 for (
int i = 1; i < p; i++)
6373 dof_map[i + 0*p1] = o++;
6374 for (
int i = 1; i < p; i++)
6375 dof_map[p + i*p1] = o++;
6376 for (
int i = 1; i < p; i++)
6377 dof_map[(p-i) + p*p1] = o++;
6378 for (
int i = 1; i < p; i++)
6379 dof_map[0 + (p-i)*p1] = o++;
6382 for (
int j = 1; j < p; j++)
6383 for (
int i = 1; i < p; i++)
6384 dof_map[i + j*p1] = o++;
6387 for (
int j = 0; j <= p; j++)
6388 for (
int i = 0; i <= p; i++)
6395 const int p =
Order;
6397 #ifdef MFEM_THREAD_SAFE
6398 Vector shape_x(p+1), shape_y(p+1);
6401 basis1d.
Eval(ip.
x, shape_x);
6402 basis1d.
Eval(ip.
y, shape_y);
6404 for (
int o = 0, j = 0; j <= p; j++)
6405 for (
int i = 0; i <= p; i++)
6406 shape(dof_map[o++]) = shape_x(i)*shape_y(j);
6412 const int p =
Order;
6414 #ifdef MFEM_THREAD_SAFE
6415 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
6418 basis1d.
Eval(ip.
x, shape_x, dshape_x);
6419 basis1d.
Eval(ip.
y, shape_y, dshape_y);
6421 for (
int o = 0, j = 0; j <= p; j++)
6422 for (
int i = 0; i <= p; i++)
6424 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j);
6425 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j); o++;
6431 const int p =
Order;
6434 #ifdef MFEM_THREAD_SAFE
6435 Vector shape_x(p+1), shape_y(p+1);
6438 for (
int i = 0; i <= p; i++)
6447 for (
int o = 0, j = 0; j <= p; j++)
6448 for (
int i = 0; i <= p; i++)
6449 dofs(dof_map[o++]) = shape_x(i)*shape_x(j);
6452 for (
int o = 0, j = 0; j <= p; j++)
6453 for (
int i = 0; i <= p; i++)
6454 dofs(dof_map[o++]) = shape_y(i)*shape_x(j);
6457 for (
int o = 0, j = 0; j <= p; j++)
6458 for (
int i = 0; i <= p; i++)
6459 dofs(dof_map[o++]) = shape_y(i)*shape_y(j);
6462 for (
int o = 0, j = 0; j <= p; j++)
6463 for (
int i = 0; i <= p; i++)
6464 dofs(dof_map[o++]) = shape_x(i)*shape_y(j);
6473 basis1d(
poly1d.ClosedBasis(p)), dof_map((p + 1)*(p + 1)*(p + 1))
6477 const int p1 = p + 1;
6479 #ifndef MFEM_THREAD_SAFE
6489 dof_map[0 + (0 + 0*p1)*p1] = 0;
6490 dof_map[p + (0 + 0*p1)*p1] = 1;
6491 dof_map[p + (p + 0*p1)*p1] = 2;
6492 dof_map[0 + (p + 0*p1)*p1] = 3;
6493 dof_map[0 + (0 + p*p1)*p1] = 4;
6494 dof_map[p + (0 + p*p1)*p1] = 5;
6495 dof_map[p + (p + p*p1)*p1] = 6;
6496 dof_map[0 + (p + p*p1)*p1] = 7;
6500 for (
int i = 1; i < p; i++)
6501 dof_map[i + (0 + 0*p1)*p1] = o++;
6502 for (
int i = 1; i < p; i++)
6503 dof_map[p + (i + 0*p1)*p1] = o++;
6504 for (
int i = 1; i < p; i++)
6505 dof_map[i + (p + 0*p1)*p1] = o++;
6506 for (
int i = 1; i < p; i++)
6507 dof_map[0 + (i + 0*p1)*p1] = o++;
6508 for (
int i = 1; i < p; i++)
6509 dof_map[i + (0 + p*p1)*p1] = o++;
6510 for (
int i = 1; i < p; i++)
6511 dof_map[p + (i + p*p1)*p1] = o++;
6512 for (
int i = 1; i < p; i++)
6513 dof_map[i + (p + p*p1)*p1] = o++;
6514 for (
int i = 1; i < p; i++)
6515 dof_map[0 + (i + p*p1)*p1] = o++;
6516 for (
int i = 1; i < p; i++)
6517 dof_map[0 + (0 + i*p1)*p1] = o++;
6518 for (
int i = 1; i < p; i++)
6519 dof_map[p + (0 + i*p1)*p1] = o++;
6520 for (
int i = 1; i < p; i++)
6521 dof_map[p + (p + i*p1)*p1] = o++;
6522 for (
int i = 1; i < p; i++)
6523 dof_map[0 + (p + i*p1)*p1] = o++;
6526 for (
int j = 1; j < p; j++)
6527 for (
int i = 1; i < p; i++)
6528 dof_map[i + ((p-j) + 0*p1)*p1] = o++;
6529 for (
int j = 1; j < p; j++)
6530 for (
int i = 1; i < p; i++)
6531 dof_map[i + (0 + j*p1)*p1] = o++;
6532 for (
int j = 1; j < p; j++)
6533 for (
int i = 1; i < p; i++)
6534 dof_map[p + (i + j*p1)*p1] = o++;
6535 for (
int j = 1; j < p; j++)
6536 for (
int i = 1; i < p; i++)
6537 dof_map[(p-i) + (p + j*p1)*p1] = o++;
6538 for (
int j = 1; j < p; j++)
6539 for (
int i = 1; i < p; i++)
6540 dof_map[0 + ((p-i) + j*p1)*p1] = o++;
6541 for (
int j = 1; j < p; j++)
6542 for (
int i = 1; i < p; i++)
6543 dof_map[i + (j + p*p1)*p1] = o++;
6546 for (
int k = 1; k < p; k++)
6547 for (
int j = 1; j < p; j++)
6548 for (
int i = 1; i < p; i++)
6549 dof_map[i + (j + k*p1)*p1] = o++;
6552 for (
int k = 0; k <= p; k++)
6553 for (
int j = 0; j <= p; j++)
6554 for (
int i = 0; i <= p; i++)
6561 const int p =
Order;
6563 #ifdef MFEM_THREAD_SAFE
6564 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
6567 basis1d.
Eval(ip.
x, shape_x);
6568 basis1d.
Eval(ip.
y, shape_y);
6569 basis1d.
Eval(ip.
z, shape_z);
6571 for (
int o = 0, k = 0; k <= p; k++)
6572 for (
int j = 0; j <= p; j++)
6573 for (
int i = 0; i <= p; i++)
6574 shape(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_z(k);
6580 const int p =
Order;
6582 #ifdef MFEM_THREAD_SAFE
6583 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
6584 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
6587 basis1d.
Eval(ip.
x, shape_x, dshape_x);
6588 basis1d.
Eval(ip.
y, shape_y, dshape_y);
6589 basis1d.
Eval(ip.
z, shape_z, dshape_z);
6591 for (
int o = 0, k = 0; k <= p; k++)
6592 for (
int j = 0; j <= p; j++)
6593 for (
int i = 0; i <= p; i++)
6595 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j)* shape_z(k);
6596 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j)* shape_z(k);
6597 dshape(dof_map[o],2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
6603 const int p =
Order;
6606 #ifdef MFEM_THREAD_SAFE
6607 Vector shape_x(p+1), shape_y(p+1);
6610 for (
int i = 0; i <= p; i++)
6619 for (
int o = 0, k = 0; k <= p; k++)
6620 for (
int j = 0; j <= p; j++)
6621 for (
int i = 0; i <= p; i++)
6622 dofs(dof_map[o++]) = shape_x(i)*shape_x(j)*shape_x(k);
6625 for (
int o = 0, k = 0; k <= p; k++)
6626 for (
int j = 0; j <= p; j++)
6627 for (
int i = 0; i <= p; i++)
6628 dofs(dof_map[o++]) = shape_y(i)*shape_x(j)*shape_x(k);
6631 for (
int o = 0, k = 0; k <= p; k++)
6632 for (
int j = 0; j <= p; j++)
6633 for (
int i = 0; i <= p; i++)
6634 dofs(dof_map[o++]) = shape_y(i)*shape_y(j)*shape_x(k);
6637 for (
int o = 0, k = 0; k <= p; k++)
6638 for (
int j = 0; j <= p; j++)
6639 for (
int i = 0; i <= p; i++)
6640 dofs(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_x(k);
6643 for (
int o = 0, k = 0; k <= p; k++)
6644 for (
int j = 0; j <= p; j++)
6645 for (
int i = 0; i <= p; i++)
6646 dofs(dof_map[o++]) = shape_x(i)*shape_x(j)*shape_y(k);
6649 for (
int o = 0, k = 0; k <= p; k++)
6650 for (
int j = 0; j <= p; j++)
6651 for (
int i = 0; i <= p; i++)
6652 dofs(dof_map[o++]) = shape_y(i)*shape_x(j)*shape_y(k);
6655 for (
int o = 0, k = 0; k <= p; k++)
6656 for (
int j = 0; j <= p; j++)
6657 for (
int i = 0; i <= p; i++)
6658 dofs(dof_map[o++]) = shape_y(i)*shape_y(j)*shape_y(k);
6661 for (
int o = 0, k = 0; k <= p; k++)
6662 for (
int j = 0; j <= p; j++)
6663 for (
int i = 0; i <= p; i++)
6664 dofs(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_y(k);
6673 #ifndef MFEM_THREAD_SAFE
6682 for (
int i = 1; i < p; i++)
6689 const int p =
Order;
6691 #ifdef MFEM_THREAD_SAFE
6698 shape(0) = shape_x(0);
6699 shape(1) = shape_x(p);
6700 for (
int i = 1; i < p; i++)
6701 shape(i+1) = shape_x(i);
6707 const int p =
Order;
6709 #ifdef MFEM_THREAD_SAFE
6710 Vector shape_x(p+1), dshape_x(p+1);
6716 dshape(0,0) = dshape_x(0);
6717 dshape(1,0) = dshape_x(p);
6718 for (
int i = 1; i < p; i++)
6719 dshape(i+1,0) = dshape_x(i);
6732 dof_map((p + 1)*(p + 1))
6734 const int p1 = p + 1;
6736 #ifndef MFEM_THREAD_SAFE
6745 dof_map[0 + 0*p1] = 0;
6746 dof_map[p + 0*p1] = 1;
6747 dof_map[p + p*p1] = 2;
6748 dof_map[0 + p*p1] = 3;
6752 for (
int i = 1; i < p; i++)
6753 dof_map[i + 0*p1] = o++;
6754 for (
int i = 1; i < p; i++)
6755 dof_map[p + i*p1] = o++;
6756 for (
int i = 1; i < p; i++)
6757 dof_map[(p-i) + p*p1] = o++;
6758 for (
int i = 1; i < p; i++)
6759 dof_map[0 + (p-i)*p1] = o++;
6762 for (
int j = 1; j < p; j++)
6763 for (
int i = 1; i < p; i++)
6764 dof_map[i + j*p1] = o++;
6767 for (
int j = 0; j <= p; j++)
6768 for (
int i = 0; i <= p; i++)
6775 const int p =
Order;
6777 #ifdef MFEM_THREAD_SAFE
6778 Vector shape_x(p+1), shape_y(p+1);
6785 for (
int o = 0, j = 0; j <= p; j++)
6786 for (
int i = 0; i <= p; i++)
6787 shape(dof_map[o++]) = shape_x(i)*shape_y(j);
6793 const int p =
Order;
6795 #ifdef MFEM_THREAD_SAFE
6796 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
6803 for (
int o = 0, j = 0; j <= p; j++)
6804 for (
int i = 0; i <= p; i++)
6806 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j);
6807 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j); o++;
6821 dof_map((p + 1)*(p + 1)*(p + 1))
6823 const int p1 = p + 1;
6825 #ifndef MFEM_THREAD_SAFE
6836 dof_map[0 + (0 + 0*p1)*p1] = 0;
6837 dof_map[p + (0 + 0*p1)*p1] = 1;
6838 dof_map[p + (p + 0*p1)*p1] = 2;
6839 dof_map[0 + (p + 0*p1)*p1] = 3;
6840 dof_map[0 + (0 + p*p1)*p1] = 4;
6841 dof_map[p + (0 + p*p1)*p1] = 5;
6842 dof_map[p + (p + p*p1)*p1] = 6;
6843 dof_map[0 + (p + p*p1)*p1] = 7;
6847 for (
int i = 1; i < p; i++)
6848 dof_map[i + (0 + 0*p1)*p1] = o++;
6849 for (
int i = 1; i < p; i++)
6850 dof_map[p + (i + 0*p1)*p1] = o++;
6851 for (
int i = 1; i < p; i++)
6852 dof_map[i + (p + 0*p1)*p1] = o++;
6853 for (
int i = 1; i < p; i++)
6854 dof_map[0 + (i + 0*p1)*p1] = o++;
6855 for (
int i = 1; i < p; i++)
6856 dof_map[i + (0 + p*p1)*p1] = o++;
6857 for (
int i = 1; i < p; i++)
6858 dof_map[p + (i + p*p1)*p1] = o++;
6859 for (
int i = 1; i < p; i++)
6860 dof_map[i + (p + p*p1)*p1] = o++;
6861 for (
int i = 1; i < p; i++)
6862 dof_map[0 + (i + p*p1)*p1] = o++;
6863 for (
int i = 1; i < p; i++)
6864 dof_map[0 + (0 + i*p1)*p1] = o++;
6865 for (
int i = 1; i < p; i++)
6866 dof_map[p + (0 + i*p1)*p1] = o++;
6867 for (
int i = 1; i < p; i++)
6868 dof_map[p + (p + i*p1)*p1] = o++;
6869 for (
int i = 1; i < p; i++)
6870 dof_map[0 + (p + i*p1)*p1] = o++;
6873 for (
int j = 1; j < p; j++)
6874 for (
int i = 1; i < p; i++)
6875 dof_map[i + ((p-j) + 0*p1)*p1] = o++;
6876 for (
int j = 1; j < p; j++)
6877 for (
int i = 1; i < p; i++)
6878 dof_map[i + (0 + j*p1)*p1] = o++;
6879 for (
int j = 1; j < p; j++)
6880 for (
int i = 1; i < p; i++)
6881 dof_map[p + (i + j*p1)*p1] = o++;
6882 for (
int j = 1; j < p; j++)
6883 for (
int i = 1; i < p; i++)
6884 dof_map[(p-i) + (p + j*p1)*p1] = o++;
6885 for (
int j = 1; j < p; j++)
6886 for (
int i = 1; i < p; i++)
6887 dof_map[0 + ((p-i) + j*p1)*p1] = o++;
6888 for (
int j = 1; j < p; j++)
6889 for (
int i = 1; i < p; i++)
6890 dof_map[i + (j + p*p1)*p1] = o++;
6893 for (
int k = 1; k < p; k++)
6894 for (
int j = 1; j < p; j++)
6895 for (
int i = 1; i < p; i++)
6896 dof_map[i + (j + k*p1)*p1] = o++;
6899 for (
int k = 0; k <= p; k++)
6900 for (
int j = 0; j <= p; j++)
6901 for (
int i = 0; i <= p; i++)
6909 const int p =
Order;
6911 #ifdef MFEM_THREAD_SAFE
6912 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
6919 for (
int o = 0, k = 0; k <= p; k++)
6920 for (
int j = 0; j <= p; j++)
6921 for (
int i = 0; i <= p; i++)
6922 shape(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_z(k);
6928 const int p =
Order;
6930 #ifdef MFEM_THREAD_SAFE
6931 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
6932 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
6939 for (
int o = 0, k = 0; k <= p; k++)
6940 for (
int j = 0; j <= p; j++)
6941 for (
int i = 0; i <= p; i++)
6943 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j)* shape_z(k);
6944 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j)* shape_z(k);
6945 dshape(dof_map[o],2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
6962 #ifndef MFEM_THREAD_SAFE
6972 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
6982 for (
int i = 1; i < p; i++)
6984 for (
int i = 1; i < p; i++)
6986 for (
int i = 1; i < p; i++)
6990 for (
int j = 1; j < p; j++)
6991 for (
int i = 1; i + j < p; i++)
6993 const double w = cp[i] + cp[j] + cp[p-i-j];
6997 for (
int k = 0; k <
Dof; k++)
7005 for (
int j = 0; j <= p; j++)
7006 for (
int i = 0; i + j <= p; i++)
7007 T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
7017 const int p =
Order;
7019 #ifdef MFEM_THREAD_SAFE
7020 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(
Dof);
7027 for (
int o = 0, j = 0; j <= p; j++)
7028 for (
int i = 0; i + j <= p; i++)
7029 u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
7037 const int p =
Order;
7039 #ifdef MFEM_THREAD_SAFE
7040 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
7041 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
7049 for (
int o = 0, j = 0; j <= p; j++)
7050 for (
int i = 0; i + j <= p; i++)
7053 du(o,0) = ((dshape_x(i)* shape_l(k)) -
7054 ( shape_x(i)*dshape_l(k)))*shape_y(j);
7055 du(o,1) = ((dshape_y(j)* shape_l(k)) -
7056 ( shape_y(j)*dshape_l(k)))*shape_x(i);
7060 Mult(T, du, dshape);
7070 #ifndef MFEM_THREAD_SAFE
7082 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7093 for (
int i = 1; i < p; i++)
7095 for (
int i = 1; i < p; i++)
7097 for (
int i = 1; i < p; i++)
7099 for (
int i = 1; i < p; i++)
7101 for (
int i = 1; i < p; i++)
7103 for (
int i = 1; i < p; i++)
7107 for (
int j = 1; j < p; j++)
7108 for (
int i = 1; i + j < p; i++)
7110 double w = cp[i] + cp[j] + cp[p-i-j];
7113 for (
int j = 1; j < p; j++)
7114 for (
int i = 1; i + j < p; i++)
7116 double w = cp[i] + cp[j] + cp[p-i-j];
7119 for (
int j = 1; j < p; j++)
7120 for (
int i = 1; i + j < p; i++)
7122 double w = cp[i] + cp[j] + cp[p-i-j];
7125 for (
int j = 1; j < p; j++)
7126 for (
int i = 1; i + j < p; i++)
7128 double w = cp[i] + cp[j] + cp[p-i-j];
7133 for (
int k = 1; k < p; k++)
7134 for (
int j = 1; j + k < p; j++)
7135 for (
int i = 1; i + j + k < p; i++)
7137 double w = cp[i] + cp[j] + cp[k] + cp[p-i-j-k];
7141 for (
int m = 0; m <
Dof; m++)
7150 for (
int k = 0; k <= p; k++)
7151 for (
int j = 0; j + k <= p; j++)
7152 for (
int i = 0; i + j + k <= p; i++)
7153 T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
7163 const int p =
Order;
7165 #ifdef MFEM_THREAD_SAFE
7166 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7175 for (
int o = 0, k = 0; k <= p; k++)
7176 for (
int j = 0; j + k <= p; j++)
7177 for (
int i = 0; i + j + k <= p; i++)
7178 u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
7186 const int p =
Order;
7188 #ifdef MFEM_THREAD_SAFE
7189 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7190 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
7199 for (
int o = 0, k = 0; k <= p; k++)
7200 for (
int j = 0; j + k <= p; j++)
7201 for (
int i = 0; i + j + k <= p; i++)
7203 int l = p - i - j - k;
7204 du(o,0) = ((dshape_x(i)* shape_l(l)) -
7205 ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
7206 du(o,1) = ((dshape_y(j)* shape_l(l)) -
7207 ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
7208 du(o,2) = ((dshape_z(k)* shape_l(l)) -
7209 ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
7213 Mult(T, du, dshape);
7235 #ifndef MFEM_THREAD_SAFE
7240 for (
int i = 0; i <= p; i++)
7247 basis1d->
Eval(ip.
x, shape);
7253 #ifdef MFEM_THREAD_SAFE
7258 basis1d->
Eval(ip.
x, shape_x, dshape_x);
7264 const int p =
Order;
7277 for (
int i = 0; i <= p; i++)
7282 for (
int i = 0; i <= p; i++)
7292 #ifndef MFEM_THREAD_SAFE
7297 for (
int i = 0; i <= p; i++)
7310 #ifdef MFEM_THREAD_SAFE
7321 dofs[vertex*
Order] = 1.0;
7344 #ifndef MFEM_THREAD_SAFE
7351 for (
int o = 0, j = 0; j <= p; j++)
7352 for (
int i = 0; i <= p; i++)
7359 const int p =
Order;
7361 #ifdef MFEM_THREAD_SAFE
7362 Vector shape_x(p+1), shape_y(p+1);
7365 basis1d->
Eval(ip.
x, shape_x);
7366 basis1d->
Eval(ip.
y, shape_y);
7368 for (
int o = 0, j = 0; j <= p; j++)
7369 for (
int i = 0; i <= p; i++)
7370 shape(o++) = shape_x(i)*shape_y(j);
7376 const int p =
Order;
7378 #ifdef MFEM_THREAD_SAFE
7379 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
7382 basis1d->
Eval(ip.
x, shape_x, dshape_x);
7383 basis1d->
Eval(ip.
y, shape_y, dshape_y);
7385 for (
int o = 0, j = 0; j <= p; j++)
7386 for (
int i = 0; i <= p; i++)
7388 dshape(o,0) = dshape_x(i)* shape_y(j);
7389 dshape(o,1) = shape_x(i)*dshape_y(j); o++;
7395 const int p =
Order;
7405 #ifdef MFEM_THREAD_SAFE
7406 Vector shape_x(p+1), shape_y(p+1);
7409 for (
int i = 0; i <= p; i++)
7418 for (
int o = 0, j = 0; j <= p; j++)
7419 for (
int i = 0; i <= p; i++)
7420 dofs[o++] = shape_x(i)*shape_x(j);
7423 for (
int o = 0, j = 0; j <= p; j++)
7424 for (
int i = 0; i <= p; i++)
7425 dofs[o++] = shape_y(i)*shape_x(j);
7428 for (
int o = 0, j = 0; j <= p; j++)
7429 for (
int i = 0; i <= p; i++)
7430 dofs[o++] = shape_y(i)*shape_y(j);
7433 for (
int o = 0, j = 0; j <= p; j++)
7434 for (
int i = 0; i <= p; i++)
7435 dofs[o++] = shape_x(i)*shape_y(j);
7445 #ifndef MFEM_THREAD_SAFE
7452 for (
int o = 0, j = 0; j <= p; j++)
7453 for (
int i = 0; i <= p; i++)
7460 const int p =
Order;
7462 #ifdef MFEM_THREAD_SAFE
7463 Vector shape_x(p+1), shape_y(p+1);
7469 for (
int o = 0, j = 0; j <= p; j++)
7470 for (
int i = 0; i <= p; i++)
7471 shape(o++) = shape_x(i)*shape_y(j);
7477 const int p =
Order;
7479 #ifdef MFEM_THREAD_SAFE
7480 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
7486 for (
int o = 0, j = 0; j <= p; j++)
7487 for (
int i = 0; i <= p; i++)
7489 dshape(o,0) = dshape_x(i)* shape_y(j);
7490 dshape(o,1) = shape_x(i)*dshape_y(j); o++;
7496 const int p =
Order;
7501 case 0: dofs[0] = 1.0;
break;
7502 case 1: dofs[p] = 1.0;
break;
7503 case 2: dofs[p*(p + 2)] = 1.0;
break;
7504 case 3: dofs[p*(p + 1)] = 1.0;
break;
7528 #ifndef MFEM_THREAD_SAFE
7537 for (
int o = 0, k = 0; k <= p; k++)
7538 for (
int j = 0; j <= p; j++)
7539 for (
int i = 0; i <= p; i++)
7546 const int p =
Order;
7548 #ifdef MFEM_THREAD_SAFE
7549 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7552 basis1d->
Eval(ip.
x, shape_x);
7553 basis1d->
Eval(ip.
y, shape_y);
7554 basis1d->
Eval(ip.
z, shape_z);
7556 for (
int o = 0, k = 0; k <= p; k++)
7557 for (
int j = 0; j <= p; j++)
7558 for (
int i = 0; i <= p; i++)
7559 shape(o++) = shape_x(i)*shape_y(j)*shape_z(k);
7565 const int p =
Order;
7567 #ifdef MFEM_THREAD_SAFE
7568 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7569 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
7572 basis1d->
Eval(ip.
x, shape_x, dshape_x);
7573 basis1d->
Eval(ip.
y, shape_y, dshape_y);
7574 basis1d->
Eval(ip.
z, shape_z, dshape_z);
7576 for (
int o = 0, k = 0; k <= p; k++)
7577 for (
int j = 0; j <= p; j++)
7578 for (
int i = 0; i <= p; i++)
7580 dshape(o,0) = dshape_x(i)* shape_y(j)* shape_z(k);
7581 dshape(o,1) = shape_x(i)*dshape_y(j)* shape_z(k);
7582 dshape(o,2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
7588 const int p =
Order;
7598 #ifdef MFEM_THREAD_SAFE
7599 Vector shape_x(p+1), shape_y(p+1);
7602 for (
int i = 0; i <= p; i++)
7611 for (
int o = 0, k = 0; k <= p; k++)
7612 for (
int j = 0; j <= p; j++)
7613 for (
int i = 0; i <= p; i++)
7614 dofs[o++] = shape_x(i)*shape_x(j)*shape_x(k);
7617 for (
int o = 0, k = 0; k <= p; k++)
7618 for (
int j = 0; j <= p; j++)
7619 for (
int i = 0; i <= p; i++)
7620 dofs[o++] = shape_y(i)*shape_x(j)*shape_x(k);
7623 for (
int o = 0, k = 0; k <= p; k++)
7624 for (
int j = 0; j <= p; j++)
7625 for (
int i = 0; i <= p; i++)
7626 dofs[o++] = shape_y(i)*shape_y(j)*shape_x(k);
7629 for (
int o = 0, k = 0; k <= p; k++)
7630 for (
int j = 0; j <= p; j++)
7631 for (
int i = 0; i <= p; i++)
7632 dofs[o++] = shape_x(i)*shape_y(j)*shape_x(k);
7635 for (
int o = 0, k = 0; k <= p; k++)
7636 for (
int j = 0; j <= p; j++)
7637 for (
int i = 0; i <= p; i++)
7638 dofs[o++] = shape_x(i)*shape_x(j)*shape_y(k);
7641 for (
int o = 0, k = 0; k <= p; k++)
7642 for (
int j = 0; j <= p; j++)
7643 for (
int i = 0; i <= p; i++)
7644 dofs[o++] = shape_y(i)*shape_x(j)*shape_y(k);
7647 for (
int o = 0, k = 0; k <= p; k++)
7648 for (
int j = 0; j <= p; j++)
7649 for (
int i = 0; i <= p; i++)
7650 dofs[o++] = shape_y(i)*shape_y(j)*shape_y(k);
7653 for (
int o = 0, k = 0; k <= p; k++)
7654 for (
int j = 0; j <= p; j++)
7655 for (
int i = 0; i <= p; i++)
7656 dofs[o++] = shape_x(i)*shape_y(j)*shape_y(k);
7666 #ifndef MFEM_THREAD_SAFE
7675 for (
int o = 0, k = 0; k <= p; k++)
7676 for (
int j = 0; j <= p; j++)
7677 for (
int i = 0; i <= p; i++)
7684 const int p =
Order;
7686 #ifdef MFEM_THREAD_SAFE
7687 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7694 for (
int o = 0, k = 0; k <= p; k++)
7695 for (
int j = 0; j <= p; j++)
7696 for (
int i = 0; i <= p; i++)
7697 shape(o++) = shape_x(i)*shape_y(j)*shape_z(k);
7703 const int p =
Order;
7705 #ifdef MFEM_THREAD_SAFE
7706 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7707 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
7714 for (
int o = 0, k = 0; k <= p; k++)
7715 for (
int j = 0; j <= p; j++)
7716 for (
int i = 0; i <= p; i++)
7718 dshape(o,0) = dshape_x(i)* shape_y(j)* shape_z(k);
7719 dshape(o,1) = shape_x(i)*dshape_y(j)* shape_z(k);
7720 dshape(o,2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
7726 const int p =
Order;
7731 case 0: dofs[0] = 1.0;
break;
7732 case 1: dofs[p] = 1.0;
break;
7733 case 2: dofs[p*(p + 2)] = 1.0;
break;
7734 case 3: dofs[p*(p + 1)] = 1.0;
break;
7735 case 4: dofs[p*(p + 1)*(p + 1)] = 1.0;
break;
7736 case 5: dofs[p + p*(p + 1)*(p + 1)] = 1.0;
break;
7737 case 6: dofs[
Dof - 1] = 1.0;
break;
7738 case 7: dofs[
Dof - p - 1] = 1.0;
break;
7757 #ifndef MFEM_THREAD_SAFE
7767 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
7770 for (
int o = 0, j = 0; j <= p; j++)
7771 for (
int i = 0; i + j <= p; i++)
7773 double w = op[i] + op[j] + op[p-i-j];
7777 for (
int k = 0; k <
Dof; k++)
7784 for (
int o = 0, j = 0; j <= p; j++)
7785 for (
int i = 0; i + j <= p; i++)
7786 T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
7796 const int p =
Order;
7798 #ifdef MFEM_THREAD_SAFE
7799 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(
Dof);
7806 for (
int o = 0, j = 0; j <= p; j++)
7807 for (
int i = 0; i + j <= p; i++)
7808 u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
7816 const int p =
Order;
7818 #ifdef MFEM_THREAD_SAFE
7819 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
7820 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
7828 for (
int o = 0, j = 0; j <= p; j++)
7829 for (
int i = 0; i + j <= p; i++)
7832 du(o,0) = ((dshape_x(i)* shape_l(k)) -
7833 ( shape_x(i)*dshape_l(k)))*shape_y(j);
7834 du(o,1) = ((dshape_y(j)* shape_l(k)) -
7835 ( shape_y(j)*dshape_l(k)))*shape_x(i);
7839 Mult(T, du, dshape);
7847 for (
int i = 0; i <
Dof; i++)
7850 dofs[i] = pow(1.0 - ip.
x - ip.
y,
Order);
7854 for (
int i = 0; i <
Dof; i++)
7857 dofs[i] = pow(ip.
x,
Order);
7861 for (
int i = 0; i <
Dof; i++)
7864 dofs[i] = pow(ip.
y,
Order);
7875 #ifndef MFEM_THREAD_SAFE
7879 for (
int o = 0, j = 0; j <= p; j++)
7880 for (
int i = 0; i + j <= p; i++)
7889 const int p =
Order;
7890 const double l1 = ip.
x, l2 = ip.
y, l3 = 1. - l1 - l2;
7900 for (
int o = 0, j = 0; j <= p; j++)
7904 for (
int i = 0; i <= p - j; i++)
7913 const int p =
Order;
7914 const double l1 = ip.
x, l2 = ip.
y, l3 = 1. - l1 - l2;
7916 #ifdef MFEM_THREAD_SAFE
7922 for (
int o = 0, j = 0; j <= p; j++)
7926 for (
int i = 0; i <= p - j; i++)
7927 dshape(o++,0) = s*dshape_1d(i);
7931 for (
int i = 0; i <= p; i++)
7935 for (
int o = i, j = 0; j <= p - i; j++)
7937 dshape(o,1) = s*dshape_1d(j);
7949 case 0: dofs[0] = 1.0;
break;
7950 case 1: dofs[
Order] = 1.0;
break;
7951 case 2: dofs[
Dof-1] = 1.0;
break;
7970 #ifndef MFEM_THREAD_SAFE
7982 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7985 for (
int o = 0, k = 0; k <= p; k++)
7986 for (
int j = 0; j + k <= p; j++)
7987 for (
int i = 0; i + j + k <= p; i++)
7989 double w = op[i] + op[j] + op[k] + op[p-i-j-k];
7993 for (
int m = 0; m <
Dof; m++)
8001 for (
int o = 0, k = 0; k <= p; k++)
8002 for (
int j = 0; j + k <= p; j++)
8003 for (
int i = 0; i + j + k <= p; i++)
8004 T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
8014 const int p =
Order;
8016 #ifdef MFEM_THREAD_SAFE
8017 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8026 for (
int o = 0, k = 0; k <= p; k++)
8027 for (
int j = 0; j + k <= p; j++)
8028 for (
int i = 0; i + j + k <= p; i++)
8029 u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
8037 const int p =
Order;
8039 #ifdef MFEM_THREAD_SAFE
8040 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8041 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
8050 for (
int o = 0, k = 0; k <= p; k++)
8051 for (
int j = 0; j + k <= p; j++)
8052 for (
int i = 0; i + j + k <= p; i++)
8054 int l = p - i - j - k;
8055 du(o,0) = ((dshape_x(i)* shape_l(l)) -
8056 ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
8057 du(o,1) = ((dshape_y(j)* shape_l(l)) -
8058 ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
8059 du(o,2) = ((dshape_z(k)* shape_l(l)) -
8060 ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
8064 Mult(T, du, dshape);
8072 for (
int i = 0; i <
Dof; i++)
8075 dofs[i] = pow(1.0 - ip.
x - ip.
y - ip.
z,
Order);
8079 for (
int i = 0; i <
Dof; i++)
8082 dofs[i] = pow(ip.
x,
Order);
8086 for (
int i = 0; i <
Dof; i++)
8089 dofs[i] = pow(ip.
y,
Order);
8092 for (
int i = 0; i <
Dof; i++)
8095 dofs[i] = pow(ip.
z,
Order);
8106 #ifndef MFEM_THREAD_SAFE
8110 for (
int o = 0, k = 0; k <= p; k++)
8111 for (
int j = 0; j + k <= p; j++)
8112 for (
int i = 0; i + j + k <= p; i++)
8121 const int p =
Order;
8122 const double l1 = ip.
x, l2 = ip.
y, l3 = ip.
z, l4 = 1. - l1 - l2 - l3;
8131 for (
int o = 0, k = 0; k <= p; k++)
8134 const double ek = bp[k]*l3k;
8136 for (
int j = 0; j <= p - k; j++)
8139 double ekj = ek*bpk[j]*l2j;
8140 for (
int i = 0; i <= p - k - j; i++)
8151 const int p =
Order;
8152 const double l1 = ip.
x, l2 = ip.
y, l3 = ip.
z, l4 = 1. - l1 - l2 - l3;
8154 #ifdef MFEM_THREAD_SAFE
8164 for (
int o = 0, k = 0; k <= p; k++)
8167 const double ek = bp[k]*l3k;
8169 for (
int j = 0; j <= p - k; j++)
8172 double ekj = ek*bpk[j]*l2j;
8173 for (
int i = 0; i <= p - k - j; i++)
8174 dshape(o++,0) = dshape_1d(i)*ekj;
8184 for (
int ok = 0, k = 0; k <= p; k++)
8187 const double ek = bp[k]*l3k;
8189 for (
int i = 0; i <= p - k; i++)
8192 double eki = ek*bpk[i]*l1i;
8194 for (
int j = 0; j <= p - k - i; j++)
8196 dshape(o,1) = dshape_1d(j)*eki;
8202 ok += ((p - k + 2)*(p - k + 1))/2;
8209 for (
int j = 0; j <= p; j++)
8212 const double ej = bp[j]*l2j;
8214 for (
int i = 0; i <= p - j; i++)
8217 double eji = ej*bpj[i]*l1i;
8218 int m = ((p + 2)*(p + 1))/2;
8219 int n = ((p - j + 2)*(p - j + 1))/2;
8220 for (
int o = i, k = 0; k <= p - j - i; k++)
8225 dshape(o - n,2) = dshape_1d(k)*eji;
8240 case 0: dofs[0] = 1.0;
break;
8241 case 1: dofs[
Order] = 1.0;
break;
8242 case 2: dofs[(
Order*(
Order+3))/2] = 1.0;
break;
8243 case 3: dofs[
Dof-1] = 1.0;
break;
8248 const double RT_QuadrilateralElement::nk[8] =
8249 { 0., -1., 1., 0., 0., 1., -1., 0. };
8254 cbasis1d(
poly1d.ClosedBasis(p + 1)), obasis1d(
poly1d.OpenBasis(p)),
8255 dof_map(Dof), dof2nk(Dof)
8259 const int dof2 =
Dof/2;
8261 #ifndef MFEM_THREAD_SAFE
8272 for (
int i = 0; i <= p; i++)
8273 dof_map[1*dof2 + i + 0*(p + 1)] = o++;
8274 for (
int i = 0; i <= p; i++)
8275 dof_map[0*dof2 + (p + 1) + i*(p + 2)] = o++;
8276 for (
int i = 0; i <= p; i++)
8277 dof_map[1*dof2 + (p - i) + (p + 1)*(p + 1)] = o++;
8278 for (
int i = 0; i <= p; i++)
8279 dof_map[0*dof2 + 0 + (p - i)*(p + 2)] = o++;
8282 for (
int j = 0; j <= p; j++)
8283 for (
int i = 1; i <= p; i++)
8284 dof_map[0*dof2 + i + j*(p + 2)] = o++;
8285 for (
int j = 1; j <= p; j++)
8286 for (
int i = 0; i <= p; i++)
8287 dof_map[1*dof2 + i + j*(p + 1)] = o++;
8291 for (
int j = 0; j <= p; j++)
8292 for (
int i = 0; i <= p/2; i++)
8294 int idx = 0*dof2 + i + j*(p + 2);
8295 dof_map[idx] = -1 - dof_map[idx];
8298 for (
int j = p/2 + 1; j <= p; j++)
8300 int idx = 0*dof2 + (p/2 + 1) + j*(p + 2);
8301 dof_map[idx] = -1 - dof_map[idx];
8304 for (
int j = 0; j <= p/2; j++)
8305 for (
int i = 0; i <= p; i++)
8307 int idx = 1*dof2 + i + j*(p + 1);
8308 dof_map[idx] = -1 - dof_map[idx];
8311 for (
int i = 0; i <= p/2; i++)
8313 int idx = 1*dof2 + i + (p/2 + 1)*(p + 1);
8314 dof_map[idx] = -1 - dof_map[idx];
8318 for (
int j = 0; j <= p; j++)
8319 for (
int i = 0; i <= p + 1; i++)
8322 if ((idx = dof_map[o++]) < 0)
8333 for (
int j = 0; j <= p + 1; j++)
8334 for (
int i = 0; i <= p; i++)
8337 if ((idx = dof_map[o++]) < 0)
8353 const int pp1 =
Order;
8355 #ifdef MFEM_THREAD_SAFE
8356 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
8359 cbasis1d.
Eval(ip.
x, shape_cx);
8360 obasis1d.
Eval(ip.
x, shape_ox);
8361 cbasis1d.
Eval(ip.
y, shape_cy);
8362 obasis1d.
Eval(ip.
y, shape_oy);
8365 for (
int j = 0; j < pp1; j++)
8366 for (
int i = 0; i <= pp1; i++)
8369 if ((idx = dof_map[o++]) < 0)
8370 idx = -1 - idx, s = -1;
8373 shape(idx,0) = s*shape_cx(i)*shape_oy(j);
8376 for (
int j = 0; j <= pp1; j++)
8377 for (
int i = 0; i < pp1; i++)
8380 if ((idx = dof_map[o++]) < 0)
8381 idx = -1 - idx, s = -1;
8385 shape(idx,1) = s*shape_ox(i)*shape_cy(j);
8392 const int pp1 =
Order;
8394 #ifdef MFEM_THREAD_SAFE
8395 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
8396 Vector dshape_cx(pp1 + 1), dshape_cy(pp1 + 1);
8399 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
8400 obasis1d.
Eval(ip.
x, shape_ox);
8401 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
8402 obasis1d.
Eval(ip.
y, shape_oy);
8405 for (
int j = 0; j < pp1; j++)
8406 for (
int i = 0; i <= pp1; i++)
8409 if ((idx = dof_map[o++]) < 0)
8410 idx = -1 - idx, s = -1;
8413 divshape(idx) = s*dshape_cx(i)*shape_oy(j);
8415 for (
int j = 0; j <= pp1; j++)
8416 for (
int i = 0; i < pp1; i++)
8419 if ((idx = dof_map[o++]) < 0)
8420 idx = -1 - idx, s = -1;
8423 divshape(idx) = s*shape_ox(i)*dshape_cy(j);
8428 const double RT_HexahedronElement::nk[18] =
8429 { 0.,0.,-1., 0.,-1.,0., 1.,0.,0., 0.,1.,0., -1.,0.,0., 0.,0.,1. };
8434 cbasis1d(
poly1d.ClosedBasis(p + 1)), obasis1d(
poly1d.OpenBasis(p)),
8435 dof_map(Dof), dof2nk(Dof)
8439 const int dof3 =
Dof/3;
8441 #ifndef MFEM_THREAD_SAFE
8455 for (
int j = 0; j <= p; j++)
8456 for (
int i = 0; i <= p; i++)
8457 dof_map[2*dof3 + i + ((p - j) + 0*(p + 1))*(p + 1)] = o++;
8458 for (
int j = 0; j <= p; j++)
8459 for (
int i = 0; i <= p; i++)
8460 dof_map[1*dof3 + i + (0 + j*(p + 2))*(p + 1)] = o++;
8461 for (
int j = 0; j <= p; j++)
8462 for (
int i = 0; i <= p; i++)
8463 dof_map[0*dof3 + (p + 1) + (i + j*(p + 1))*(p + 2)] = o++;
8464 for (
int j = 0; j <= p; j++)
8465 for (
int i = 0; i <= p; i++)
8466 dof_map[1*dof3 + (p - i) + ((p + 1) + j*(p + 2))*(p + 1)] = o++;
8467 for (
int j = 0; j <= p; j++)
8468 for (
int i = 0; i <= p; i++)
8469 dof_map[0*dof3 + 0 + ((p - i) + j*(p + 1))*(p + 2)] = o++;
8470 for (
int j = 0; j <= p; j++)
8471 for (
int i = 0; i <= p; i++)
8472 dof_map[2*dof3 + i + (j + (p + 1)*(p + 1))*(p + 1)] = o++;
8476 for (
int k = 0; k <= p; k++)
8477 for (
int j = 0; j <= p; j++)
8478 for (
int i = 1; i <= p; i++)
8479 dof_map[0*dof3 + i + (j + k*(p + 1))*(p + 2)] = o++;
8481 for (
int k = 0; k <= p; k++)
8482 for (
int j = 1; j <= p; j++)
8483 for (
int i = 0; i <= p; i++)
8484 dof_map[1*dof3 + i + (j + k*(p + 2))*(p + 1)] = o++;
8486 for (
int k = 1; k <= p; k++)
8487 for (
int j = 0; j <= p; j++)
8488 for (
int i = 0; i <= p; i++)
8489 dof_map[2*dof3 + i + (j + k*(p + 1))*(p + 1)] = o++;
8496 for (
int k = 0; k <= p; k++)
8497 for (
int j = 0; j <= p; j++)
8498 for (
int i = 0; i <= p/2; i++)
8500 int idx = 0*dof3 + i + (j + k*(p + 1))*(p + 2);
8501 dof_map[idx] = -1 - dof_map[idx];
8504 for (
int k = 0; k <= p; k++)
8505 for (
int j = 0; j <= p/2; j++)
8506 for (
int i = 0; i <= p; i++)
8508 int idx = 1*dof3 + i + (j + k*(p + 2))*(p + 1);
8509 dof_map[idx] = -1 - dof_map[idx];
8512 for (
int k = 0; k <= p/2; k++)
8513 for (
int j = 0; j <= p; j++)
8514 for (
int i = 0; i <= p; i++)
8516 int idx = 2*dof3 + i + (j + k*(p + 1))*(p + 1);
8517 dof_map[idx] = -1 - dof_map[idx];
8522 for (
int k = 0; k <= p; k++)
8523 for (
int j = 0; j <= p; j++)
8524 for (
int i = 0; i <= p + 1; i++)
8527 if ((idx = dof_map[o++]) < 0)
8539 for (
int k = 0; k <= p; k++)
8540 for (
int j = 0; j <= p + 1; j++)
8541 for (
int i = 0; i <= p; i++)
8544 if ((idx = dof_map[o++]) < 0)
8556 for (
int k = 0; k <= p + 1; k++)
8557 for (
int j = 0; j <= p; j++)
8558 for (
int i = 0; i <= p; i++)
8561 if ((idx = dof_map[o++]) < 0)
8577 const int pp1 =
Order;
8579 #ifdef MFEM_THREAD_SAFE
8580 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
8581 Vector shape_cz(pp1 + 1), shape_oz(pp1);
8584 cbasis1d.
Eval(ip.
x, shape_cx);
8585 obasis1d.
Eval(ip.
x, shape_ox);
8586 cbasis1d.
Eval(ip.
y, shape_cy);
8587 obasis1d.
Eval(ip.
y, shape_oy);
8588 cbasis1d.
Eval(ip.
z, shape_cz);
8589 obasis1d.
Eval(ip.
z, shape_oz);
8593 for (
int k = 0; k < pp1; k++)
8594 for (
int j = 0; j < pp1; j++)
8595 for (
int i = 0; i <= pp1; i++)
8598 if ((idx = dof_map[o++]) < 0)
8599 idx = -1 - idx, s = -1;
8602 shape(idx,0) = s*shape_cx(i)*shape_oy(j)*shape_oz(k);
8607 for (
int k = 0; k < pp1; k++)
8608 for (
int j = 0; j <= pp1; j++)
8609 for (
int i = 0; i < pp1; i++)
8612 if ((idx = dof_map[o++]) < 0)
8613 idx = -1 - idx, s = -1;
8617 shape(idx,1) = s*shape_ox(i)*shape_cy(j)*shape_oz(k);
8621 for (
int k = 0; k <= pp1; k++)
8622 for (
int j = 0; j < pp1; j++)
8623 for (
int i = 0; i < pp1; i++)
8626 if ((idx = dof_map[o++]) < 0)
8627 idx = -1 - idx, s = -1;
8632 shape(idx,2) = s*shape_ox(i)*shape_oy(j)*shape_cz(k);
8639 const int pp1 =
Order;
8641 #ifdef MFEM_THREAD_SAFE
8642 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
8643 Vector shape_cz(pp1 + 1), shape_oz(pp1);
8644 Vector dshape_cx(pp1 + 1), dshape_cy(pp1 + 1), dshape_cz(pp1 + 1);
8647 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
8648 obasis1d.
Eval(ip.
x, shape_ox);
8649 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
8650 obasis1d.
Eval(ip.
y, shape_oy);
8651 cbasis1d.
Eval(ip.
z, shape_cz, dshape_cz);
8652 obasis1d.
Eval(ip.
z, shape_oz);
8656 for (
int k = 0; k < pp1; k++)
8657 for (
int j = 0; j < pp1; j++)
8658 for (
int i = 0; i <= pp1; i++)
8661 if ((idx = dof_map[o++]) < 0)
8662 idx = -1 - idx, s = -1;
8665 divshape(idx) = s*dshape_cx(i)*shape_oy(j)*shape_oz(k);
8668 for (
int k = 0; k < pp1; k++)
8669 for (
int j = 0; j <= pp1; j++)
8670 for (
int i = 0; i < pp1; i++)
8673 if ((idx = dof_map[o++]) < 0)
8674 idx = -1 - idx, s = -1;
8677 divshape(idx) = s*shape_ox(i)*dshape_cy(j)*shape_oz(k);
8680 for (
int k = 0; k <= pp1; k++)
8681 for (
int j = 0; j < pp1; j++)
8682 for (
int i = 0; i < pp1; i++)
8685 if ((idx = dof_map[o++]) < 0)
8686 idx = -1 - idx, s = -1;
8689 divshape(idx) = s*shape_ox(i)*shape_oy(j)*dshape_cz(k);
8694 const double RT_TriangleElement::nk[6] =
8695 { 0., -1., 1., 1., -1., 0. };
8697 const double RT_TriangleElement::c = 1./3.;
8706 #ifndef MFEM_THREAD_SAFE
8716 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
8721 for (
int i = 0; i <= p; i++)
8726 for (
int i = 0; i <= p; i++)
8731 for (
int i = 0; i <= p; i++)
8738 for (
int j = 0; j < p; j++)
8739 for (
int i = 0; i + j < p; i++)
8741 double w = iop[i] + iop[j] + iop[p-1-i-j];
8748 for (
int k = 0; k <
Dof; k++)
8754 const double *n_k = nk + 2*dof2nk[k];
8757 for (
int j = 0; j <= p; j++)
8758 for (
int i = 0; i + j <= p; i++)
8760 double s = shape_x(i)*shape_y(j)*shape_l(p-i-j);
8761 T(o++, k) = s*n_k[0];
8762 T(o++, k) = s*n_k[1];
8764 for (
int i = 0; i <= p; i++)
8766 double s = shape_x(i)*shape_y(p-i);
8767 T(o++, k) = s*((ip.
x - c)*n_k[0] + (ip.
y - c)*n_k[1]);
8778 const int p =
Order - 1;
8780 #ifdef MFEM_THREAD_SAFE
8781 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
8790 for (
int j = 0; j <= p; j++)
8791 for (
int i = 0; i + j <= p; i++)
8793 double s = shape_x(i)*shape_y(j)*shape_l(p-i-j);
8794 u(o,0) = s; u(o,1) = 0; o++;
8795 u(o,0) = 0; u(o,1) = s; o++;
8797 for (
int i = 0; i <= p; i++)
8799 double s = shape_x(i)*shape_y(p-i);
8800 u(o,0) = (ip.
x - c)*s; u(o,1) = (ip.
y - c)*s; o++;
8809 const int p =
Order - 1;
8811 #ifdef MFEM_THREAD_SAFE
8812 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
8813 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
8822 for (
int j = 0; j <= p; j++)
8823 for (
int i = 0; i + j <= p; i++)
8826 divu(o++) = (dshape_x(i)*shape_l(k) -
8827 shape_x(i)*dshape_l(k))*shape_y(j);
8828 divu(o++) = (dshape_y(j)*shape_l(k) -
8829 shape_y(j)*dshape_l(k))*shape_x(i);
8831 for (
int i = 0; i <= p; i++)
8834 divu(o++) = ((shape_x(i) + (ip.
x - c)*dshape_x(i))*shape_y(j) +
8835 (shape_y(j) + (ip.
y - c)*dshape_y(j))*shape_x(i));
8838 T.
Mult(divu, divshape);
8842 const double RT_TetrahedronElement::nk[12] =
8843 { 1,1,1, -1,0,0, 0,-1,0, 0,0,-1 };
8846 const double RT_TetrahedronElement::c = 1./4.;
8855 #ifndef MFEM_THREAD_SAFE
8867 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8873 for (
int j = 0; j <= p; j++)
8874 for (
int i = 0; i + j <= p; i++)
8876 double w = bop[i] + bop[j] + bop[p-i-j];
8880 for (
int j = 0; j <= p; j++)
8881 for (
int i = 0; i + j <= p; i++)
8883 double w = bop[i] + bop[j] + bop[p-i-j];
8887 for (
int j = 0; j <= p; j++)
8888 for (
int i = 0; i + j <= p; i++)
8890 double w = bop[i] + bop[j] + bop[p-i-j];
8894 for (
int j = 0; j <= p; j++)
8895 for (
int i = 0; i + j <= p; i++)
8897 double w = bop[i] + bop[j] + bop[p-i-j];
8903 for (
int k = 0; k < p; k++)
8904 for (
int j = 0; j + k < p; j++)
8905 for (
int i = 0; i + j + k < p; i++)
8907 double w = iop[i] + iop[j] + iop[k] + iop[p-1-i-j-k];
8916 for (
int m = 0; m <
Dof; m++)
8923 const double *nm = nk + 3*dof2nk[m];
8926 for (
int k = 0; k <= p; k++)
8927 for (
int j = 0; j + k <= p; j++)
8928 for (
int i = 0; i + j + k <= p; i++)
8930 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
8931 T(o++, m) = s * nm[0];
8932 T(o++, m) = s * nm[1];
8933 T(o++, m) = s * nm[2];
8935 for (
int j = 0; j <= p; j++)
8936 for (
int i = 0; i + j <= p; i++)
8938 double s = shape_x(i)*shape_y(j)*shape_z(p-i-j);
8939 T(o++, m) = s*((ip.
x - c)*nm[0] + (ip.
y - c)*nm[1] +
8951 const int p =
Order - 1;
8953 #ifdef MFEM_THREAD_SAFE
8954 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8964 for (
int k = 0; k <= p; k++)
8965 for (
int j = 0; j + k <= p; j++)
8966 for (
int i = 0; i + j + k <= p; i++)
8968 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
8969 u(o,0) = s; u(o,1) = 0; u(o,2) = 0; o++;
8970 u(o,0) = 0; u(o,1) = s; u(o,2) = 0; o++;
8971 u(o,0) = 0; u(o,1) = 0; u(o,2) = s; o++;
8973 for (
int j = 0; j <= p; j++)
8974 for (
int i = 0; i + j <= p; i++)
8976 double s = shape_x(i)*shape_y(j)*shape_z(p-i-j);
8977 u(o,0) = (ip.
x - c)*s; u(o,1) = (ip.
y - c)*s; u(o,2) = (ip.
z - c)*s;
8987 const int p =
Order - 1;
8989 #ifdef MFEM_THREAD_SAFE
8990 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8991 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
9001 for (
int k = 0; k <= p; k++)
9002 for (
int j = 0; j + k <= p; j++)
9003 for (
int i = 0; i + j + k <= p; i++)
9005 int l = p - i - j - k;
9006 divu(o++) = (dshape_x(i)*shape_l(l) -
9007 shape_x(i)*dshape_l(l))*shape_y(j)*shape_z(k);
9008 divu(o++) = (dshape_y(j)*shape_l(l) -
9009 shape_y(j)*dshape_l(l))*shape_x(i)*shape_z(k);
9010 divu(o++) = (dshape_z(k)*shape_l(l) -
9011 shape_z(k)*dshape_l(l))*shape_x(i)*shape_y(j);
9013 for (
int j = 0; j <= p; j++)
9014 for (
int i = 0; i + j <= p; i++)
9018 (shape_x(i) + (ip.
x - c)*dshape_x(i))*shape_y(j)*shape_z(k) +
9019 (shape_y(j) + (ip.
y - c)*dshape_y(j))*shape_x(i)*shape_z(k) +
9020 (shape_z(k) + (ip.
z - c)*dshape_z(k))*shape_x(i)*shape_y(j);
9023 T.
Mult(divu, divshape);
9027 const double ND_HexahedronElement::tk[18] =
9028 { 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,0.,0., 0.,-1.,0., 0.,0.,-1. };
9033 cbasis1d(
poly1d.ClosedBasis(p)), obasis1d(
poly1d.OpenBasis(p - 1)),
9034 dof_map(Dof), dof2tk(Dof)
9038 const int dof3 =
Dof/3;
9040 #ifndef MFEM_THREAD_SAFE
9054 for (
int i = 0; i < p; i++)
9055 dof_map[0*dof3 + i + (0 + 0*(p + 1))*p] = o++;
9056 for (
int i = 0; i < p; i++)
9057 dof_map[1*dof3 + p + (i + 0*p)*(p + 1)] = o++;
9058 for (
int i = 0; i < p; i++)
9059 dof_map[0*dof3 + i + (p + 0*(p + 1))*p] = o++;
9060 for (
int i = 0; i < p; i++)
9061 dof_map[1*dof3 + 0 + (i + 0*p)*(p + 1)] = o++;
9062 for (
int i = 0; i < p; i++)
9063 dof_map[0*dof3 + i + (0 + p*(p + 1))*p] = o++;
9064 for (
int i = 0; i < p; i++)
9065 dof_map[1*dof3 + p + (i + p*p)*(p + 1)] = o++;
9066 for (
int i = 0; i < p; i++)
9067 dof_map[0*dof3 + i + (p + p*(p + 1))*p] = o++;
9068 for (
int i = 0; i < p; i++)
9069 dof_map[1*dof3 + 0 + (i + p*p)*(p + 1)] = o++;
9070 for (
int i = 0; i < p; i++)
9071 dof_map[2*dof3 + 0 + (0 + i*(p + 1))*(p + 1)] = o++;
9072 for (
int i = 0; i < p; i++)
9073 dof_map[2*dof3 + p + (0 + i*(p + 1))*(p + 1)] = o++;
9074 for (
int i = 0; i < p; i++)
9075 dof_map[2*dof3 + p + (p + i*(p + 1))*(p + 1)] = o++;
9076 for (
int i = 0; i < p; i++)
9077 dof_map[2*dof3 + 0 + (p + i*(p + 1))*(p + 1)] = o++;
9081 for (
int j = 1; j < p; j++)
9082 for (
int i = 0; i < p; i++)
9083 dof_map[0*dof3 + i + ((p - j) + 0*(p + 1))*p] = o++;
9084 for (
int j = 0; j < p; j++)
9085 for (
int i = 1; i < p; i++)
9086 dof_map[1*dof3 + i + ((p - 1 - j) + 0*p)*(p + 1)] = -1 - (o++);
9088 for (
int k = 1; k < p; k++)
9089 for (
int i = 0; i < p; i++)
9090 dof_map[0*dof3 + i + (0 + k*(p + 1))*p] = o++;
9091 for (
int k = 0; k < p; k++)
9092 for (
int i = 1; i < p; i++ )
9093 dof_map[2*dof3 + i + (0 + k*(p + 1))*(p + 1)] = o++;
9095 for (
int k = 1; k < p; k++)
9096 for (
int j = 0; j < p; j++)
9097 dof_map[1*dof3 + p + (j + k*p)*(p + 1)] = o++;
9098 for (
int k = 0; k < p; k++)
9099 for (
int j = 1; j < p; j++)
9100 dof_map[2*dof3 + p + (j + k*(p + 1))*(p + 1)] = o++;
9102 for (
int k = 1; k < p; k++)
9103 for (
int i = 0; i < p; i++)
9104 dof_map[0*dof3 + (p - 1 - i) + (p + k*(p + 1))*p] = -1 - (o++);
9105 for (
int k = 0; k < p; k++)
9106 for (
int i = 1; i < p; i++)
9107 dof_map[2*dof3 + (p - i) + (p + k*(p + 1))*(p + 1)] = o++;
9109 for (
int k = 1; k < p; k++)
9110 for (
int j = 0; j < p; j++)
9111 dof_map[1*dof3 + 0 + ((p - 1 - j) + k*p)*(p + 1)] = -1 - (o++);
9112 for (
int k = 0; k < p; k++)
9113 for (
int j = 1; j < p; j++)
9114 dof_map[2*dof3 + 0 + ((p - j) + k*(p + 1))*(p + 1)] = o++;
9116 for (
int j = 1; j < p; j++)
9117 for (
int i = 0; i < p; i++)
9118 dof_map[0*dof3 + i + (j + p*(p + 1))*p] = o++;
9119 for (
int j = 0; j < p; j++)
9120 for (
int i = 1; i < p; i++)
9121 dof_map[1*dof3 + i + (j + p*p)*(p + 1)] = o++;
9125 for (
int k = 1; k < p; k++)
9126 for (
int j = 1; j < p; j++)
9127 for (
int i = 0; i < p; i++)
9128 dof_map[0*dof3 + i + (j + k*(p + 1))*p] = o++;
9130 for (
int k = 1; k < p; k++)
9131 for (
int j = 0; j < p; j++)
9132 for (
int i = 1; i < p; i++)
9133 dof_map[1*dof3 + i + (j + k*p)*(p + 1)] = o++;
9135 for (
int k = 0; k < p; k++)
9136 for (
int j = 1; j < p; j++)
9137 for (
int i = 1; i < p; i++)
9138 dof_map[2*dof3 + i + (j + k*(p + 1))*(p + 1)] = o++;
9143 for (
int k = 0; k <= p; k++)
9144 for (
int j = 0; j <= p; j++)
9145 for (
int i = 0; i < p; i++)
9148 if ((idx = dof_map[o++]) < 0)
9149 dof2tk[idx = -1 - idx] = 3;
9155 for (
int k = 0; k <= p; k++)
9156 for (
int j = 0; j < p; j++)
9157 for (
int i = 0; i <= p; i++)
9160 if ((idx = dof_map[o++]) < 0)
9161 dof2tk[idx = -1 - idx] = 4;
9167 for (
int k = 0; k < p; k++)
9168 for (
int j = 0; j <= p; j++)
9169 for (
int i = 0; i <= p; i++)
9172 if ((idx = dof_map[o++]) < 0)
9173 dof2tk[idx = -1 - idx] = 5;
9183 const int p =
Order;
9185 #ifdef MFEM_THREAD_SAFE
9186 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
9187 Vector shape_cz(p + 1), shape_oz(p);
9190 cbasis1d.
Eval(ip.
x, shape_cx);
9191 obasis1d.
Eval(ip.
x, shape_ox);
9192 cbasis1d.
Eval(ip.
y, shape_cy);
9193 obasis1d.
Eval(ip.
y, shape_oy);
9194 cbasis1d.
Eval(ip.
z, shape_cz);
9195 obasis1d.
Eval(ip.
z, shape_oz);
9199 for (
int k = 0; k <= p; k++)
9200 for (
int j = 0; j <= p; j++)
9201 for (
int i = 0; i < p; i++)
9204 if ((idx = dof_map[o++]) < 0)
9205 idx = -1 - idx, s = -1;
9208 shape(idx,0) = s*shape_ox(i)*shape_cy(j)*shape_cz(k);
9213 for (
int k = 0; k <= p; k++)
9214 for (
int j = 0; j < p; j++)
9215 for (
int i = 0; i <= p; i++)
9218 if ((idx = dof_map[o++]) < 0)
9219 idx = -1 - idx, s = -1;
9223 shape(idx,1) = s*shape_cx(i)*shape_oy(j)*shape_cz(k);
9227 for (
int k = 0; k < p; k++)
9228 for (
int j = 0; j <= p; j++)
9229 for (
int i = 0; i <= p; i++)
9232 if ((idx = dof_map[o++]) < 0)
9233 idx = -1 - idx, s = -1;
9238 shape(idx,2) = s*shape_cx(i)*shape_cy(j)*shape_oz(k);
9245 const int p =
Order;
9247 #ifdef MFEM_THREAD_SAFE
9248 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
9249 Vector shape_cz(p + 1), shape_oz(p);
9250 Vector dshape_cx(p + 1), dshape_cy(p + 1), dshape_cz(p + 1);
9253 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
9254 obasis1d.
Eval(ip.
x, shape_ox);
9255 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
9256 obasis1d.
Eval(ip.
y, shape_oy);
9257 cbasis1d.
Eval(ip.
z, shape_cz, dshape_cz);
9258 obasis1d.
Eval(ip.
z, shape_oz);
9262 for (
int k = 0; k <= p; k++)
9263 for (
int j = 0; j <= p; j++)
9264 for (
int i = 0; i < p; i++)
9267 if ((idx = dof_map[o++]) < 0)
9268 idx = -1 - idx, s = -1;
9271 curl_shape(idx,0) = 0.;
9272 curl_shape(idx,1) = s*shape_ox(i)* shape_cy(j)*dshape_cz(k);
9273 curl_shape(idx,2) = -s*shape_ox(i)*dshape_cy(j)* shape_cz(k);
9276 for (
int k = 0; k <= p; k++)
9277 for (
int j = 0; j < p; j++)
9278 for (
int i = 0; i <= p; i++)
9281 if ((idx = dof_map[o++]) < 0)
9282 idx = -1 - idx, s = -1;
9285 curl_shape(idx,0) = -s* shape_cx(i)*shape_oy(j)*dshape_cz(k);
9286 curl_shape(idx,1) = 0.;
9287 curl_shape(idx,2) = s*dshape_cx(i)*shape_oy(j)* shape_cz(k);
9290 for (
int k = 0; k < p; k++)
9291 for (
int j = 0; j <= p; j++)
9292 for (
int i = 0; i <= p; i++)
9295 if ((idx = dof_map[o++]) < 0)
9296 idx = -1 - idx, s = -1;
9299 curl_shape(idx,0) = s* shape_cx(i)*dshape_cy(j)*shape_oz(k);
9300 curl_shape(idx,1) = -s*dshape_cx(i)* shape_cy(j)*shape_oz(k);
9301 curl_shape(idx,2) = 0.;
9306 const double ND_QuadrilateralElement::tk[8] =
9307 { 1.,0., 0.,1., -1.,0., 0.,-1. };
9312 cbasis1d(
poly1d.ClosedBasis(p)), obasis1d(
poly1d.OpenBasis(p - 1)),
9313 dof_map(Dof), dof2tk(Dof)
9317 const int dof2 =
Dof/2;
9319 #ifndef MFEM_THREAD_SAFE
9330 for (
int i = 0; i < p; i++)
9331 dof_map[0*dof2 + i + 0*p] = o++;
9332 for (
int j = 0; j < p; j++)
9333 dof_map[1*dof2 + p + j*(p + 1)] = o++;
9334 for (
int i = 0; i < p; i++)
9335 dof_map[0*dof2 + (p - 1 - i) + p*p] = -1 - (o++);
9336 for (
int j = 0; j < p; j++)
9337 dof_map[1*dof2 + 0 + (p - 1 - j)*(p + 1)] = -1 - (o++);
9341 for (
int j = 1; j < p; j++)
9342 for (
int i = 0; i < p; i++)
9343 dof_map[0*dof2 + i + j*p] = o++;
9345 for (
int j = 0; j < p; j++)
9346 for (
int i = 1; i < p; i++)
9347 dof_map[1*dof2 + i + j*(p + 1)] = o++;
9352 for (
int j = 0; j <= p; j++)
9353 for (
int i = 0; i < p; i++)
9356 if ((idx = dof_map[o++]) < 0)
9357 dof2tk[idx = -1 - idx] = 2;
9363 for (
int j = 0; j < p; j++)
9364 for (
int i = 0; i <= p; i++)
9367 if ((idx = dof_map[o++]) < 0)
9368 dof2tk[idx = -1 - idx] = 3;
9378 const int p =
Order;
9380 #ifdef MFEM_THREAD_SAFE
9381 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
9384 cbasis1d.
Eval(ip.
x, shape_cx);
9385 obasis1d.
Eval(ip.
x, shape_ox);
9386 cbasis1d.
Eval(ip.
y, shape_cy);
9387 obasis1d.
Eval(ip.
y, shape_oy);
9391 for (
int j = 0; j <= p; j++)
9392 for (
int i = 0; i < p; i++)
9395 if ((idx = dof_map[o++]) < 0)
9396 idx = -1 - idx, s = -1;
9399 shape(idx,0) = s*shape_ox(i)*shape_cy(j);
9403 for (
int j = 0; j < p; j++)
9404 for (
int i = 0; i <= p; i++)
9407 if ((idx = dof_map[o++]) < 0)
9408 idx = -1 - idx, s = -1;
9412 shape(idx,1) = s*shape_cx(i)*shape_oy(j);
9419 mfem_error(
"ND_QuadrilateralElement::CalcCurlShape");
9423 const double ND_TetrahedronElement::tk[18] =
9424 { 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,1.,0., -1.,0.,1., 0.,-1.,1. };
9426 const double ND_TetrahedronElement::c = 1./4.;
9436 const int pm1 = p - 1, pm2 = p - 2, pm3 = p - 3;
9438 #ifndef MFEM_THREAD_SAFE
9449 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
9454 for (
int i = 0; i < p; i++)
9459 for (
int i = 0; i < p; i++)
9464 for (
int i = 0; i < p; i++)
9469 for (
int i = 0; i < p; i++)
9474 for (
int i = 0; i < p; i++)
9479 for (
int i = 0; i < p; i++)
9486 for (
int j = 0; j <= pm2; j++)
9487 for (
int i = 0; i + j <= pm2; i++)
9489 double w = fop[i] + fop[j] + fop[pm2-i-j];
9495 for (
int j = 0; j <= pm2; j++)
9496 for (
int i = 0; i + j <= pm2; i++)
9498 double w = fop[i] + fop[j] + fop[pm2-i-j];
9504 for (
int j = 0; j <= pm2; j++)
9505 for (
int i = 0; i + j <= pm2; i++)
9507 double w = fop[i] + fop[j] + fop[pm2-i-j];
9513 for (
int j = 0; j <= pm2; j++)
9514 for (
int i = 0; i + j <= pm2; i++)
9516 double w = fop[i] + fop[j] + fop[pm2-i-j];
9524 for (
int k = 0; k <= pm3; k++)
9525 for (
int j = 0; j + k <= pm3; j++)
9526 for (
int i = 0; i + j + k <= pm3; i++)
9528 double w = iop[i] + iop[j] + iop[k] + iop[pm3-i-j-k];
9537 for (
int m = 0; m <
Dof; m++)
9540 const double *tm = tk + 3*dof2tk[m];
9548 for (
int k = 0; k <= pm1; k++)
9549 for (
int j = 0; j + k <= pm1; j++)
9550 for (
int i = 0; i + j + k <= pm1; i++)
9552 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
9553 T(o++, m) = s * tm[0];
9554 T(o++, m) = s * tm[1];
9555 T(o++, m) = s * tm[2];
9557 for (
int k = 0; k <= pm1; k++)
9558 for (
int j = 0; j + k <= pm1; j++)
9560 double s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
9561 T(o++, m) = s*((ip.
y - c)*tm[0] - (ip.
x - c)*tm[1]);
9562 T(o++, m) = s*((ip.
z - c)*tm[0] - (ip.
x - c)*tm[2]);
9564 for (
int k = 0; k <= pm1; k++)
9567 shape_y(pm1-k)*shape_z(k)*((ip.
z - c)*tm[1] - (ip.
y - c)*tm[2]);
9578 const int pm1 =
Order - 1;
9580 #ifdef MFEM_THREAD_SAFE
9581 const int p =
Order;
9582 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
9592 for (
int k = 0; k <= pm1; k++)
9593 for (
int j = 0; j + k <= pm1; j++)
9594 for (
int i = 0; i + j + k <= pm1; i++)
9596 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
9597 u(n,0) = s; u(n,1) = 0.; u(n,2) = 0.; n++;
9598 u(n,0) = 0.; u(n,1) = s; u(n,2) = 0.; n++;
9599 u(n,0) = 0.; u(n,1) = 0.; u(n,2) = s; n++;
9601 for (
int k = 0; k <= pm1; k++)
9602 for (
int j = 0; j + k <= pm1; j++)
9604 double s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
9605 u(n,0) = s*(ip.
y - c); u(n,1) = -s*(ip.
x - c); u(n,2) = 0.; n++;
9606 u(n,0) = s*(ip.
z - c); u(n,1) = 0.; u(n,2) = -s*(ip.
x - c); n++;
9608 for (
int k = 0; k <= pm1; k++)
9610 double s = shape_y(pm1-k)*shape_z(k);
9611 u(n,0) = 0.; u(n,1) = s*(ip.
z - c); u(n,2) = -s*(ip.
y - c); n++;
9620 const int pm1 =
Order - 1;
9622 #ifdef MFEM_THREAD_SAFE
9623 const int p =
Order;
9624 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
9625 Vector dshape_x(p), dshape_y(p), dshape_z(p), dshape_l(p);
9635 for (
int k = 0; k <= pm1; k++)
9636 for (
int j = 0; j + k <= pm1; j++)
9637 for (
int i = 0; i + j + k <= pm1; i++)
9640 const double dx = (dshape_x(i)*shape_l(l) -
9641 shape_x(i)*dshape_l(l))*shape_y(j)*shape_z(k);
9642 const double dy = (dshape_y(j)*shape_l(l) -
9643 shape_y(j)*dshape_l(l))*shape_x(i)*shape_z(k);
9644 const double dz = (dshape_z(k)*shape_l(l) -
9645 shape_z(k)*dshape_l(l))*shape_x(i)*shape_y(j);
9647 u(n,0) = 0.; u(n,1) = dz; u(n,2) = -dy; n++;
9648 u(n,0) = -dz; u(n,1) = 0.; u(n,2) = dx; n++;
9649 u(n,0) = dy; u(n,1) = -dx; u(n,2) = 0.; n++;
9651 for (
int k = 0; k <= pm1; k++)
9652 for (
int j = 0; j + k <= pm1; j++)
9654 int i = pm1 - j - k;
9657 u(n,0) = shape_x(i)*(ip.
x - c)*shape_y(j)*dshape_z(k);
9658 u(n,1) = shape_x(i)*shape_y(j)*(ip.
y - c)*dshape_z(k);
9660 -((dshape_x(i)*(ip.
x - c) + shape_x(i))*shape_y(j)*shape_z(k) +
9661 (dshape_y(j)*(ip.
y - c) + shape_y(j))*shape_x(i)*shape_z(k));
9664 u(n,0) = -shape_x(i)*(ip.
x - c)*dshape_y(j)*shape_z(k);
9665 u(n,1) = (shape_x(i)*shape_y(j)*(dshape_z(k)*(ip.
z - c) + shape_z(k)) +
9666 (dshape_x(i)*(ip.
x - c) + shape_x(i))*shape_y(j)*shape_z(k));
9667 u(n,2) = -shape_x(i)*dshape_y(j)*shape_z(k)*(ip.
z - c);
9670 for (
int k = 0; k <= pm1; k++)
9674 u(n,0) = -((dshape_y(j)*(ip.
y - c) + shape_y(j))*shape_z(k) +
9675 shape_y(j)*(dshape_z(k)*(ip.
z - c) + shape_z(k)));
9680 Mult(T, u, curl_shape);
9684 const double ND_TriangleElement::tk[8] =
9685 { 1.,0., -1.,1., 0.,-1., 0.,1. };
9687 const double ND_TriangleElement::c = 1./3.;
9696 const int pm1 = p - 1, pm2 = p - 2;
9698 #ifndef MFEM_THREAD_SAFE
9707 Vector shape_x(p), shape_y(p), shape_l(p);
9712 for (
int i = 0; i < p; i++)
9717 for (
int i = 0; i < p; i++)
9722 for (
int i = 0; i < p; i++)
9729 for (
int j = 0; j <= pm2; j++)
9730 for (
int i = 0; i + j <= pm2; i++)
9732 double w = iop[i] + iop[j] + iop[pm2-i-j];
9739 for (
int m = 0; m <
Dof; m++)
9742 const double *tm = tk + 2*dof2tk[m];
9749 for (
int j = 0; j <= pm1; j++)
9750 for (
int i = 0; i + j <= pm1; i++)
9752 double s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
9753 T(n++, m) = s * tm[0];
9754 T(n++, m) = s * tm[1];
9756 for (
int j = 0; j <= pm1; j++)
9759 shape_x(pm1-j)*shape_y(j)*((ip.
y - c)*tm[0] - (ip.
x - c)*tm[1]);
9770 const int pm1 =
Order - 1;
9772 #ifdef MFEM_THREAD_SAFE
9773 const int p =
Order;
9774 Vector shape_x(p), shape_y(p), shape_l(p);
9783 for (
int j = 0; j <= pm1; j++)
9784 for (
int i = 0; i + j <= pm1; i++)
9786 double s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
9787 u(n,0) = s; u(n,1) = 0.; n++;
9788 u(n,0) = 0.; u(n,1) = s; n++;
9790 for (
int j = 0; j <= pm1; j++)
9792 double s = shape_x(pm1-j)*shape_y(j);
9793 u(n,0) = s*(ip.
y - c);
9794 u(n,1) = -s*(ip.
x - c);
9804 mfem_error(
"ND_TriangleElement::CalcCurlShape");
9811 kv[0]->CalcShape(shape,
ijk[0], ip.
x);
9814 for (
int i = 0; i <=
Order; i++)
9815 sum += (shape(i) *=
weights(i));
9826 kv[0]->CalcDShape(grad,
ijk[0], ip.
x);
9828 double sum = 0.0, dsum = 0.0;
9829 for (
int i = 0; i <=
Order; i++)
9832 dsum += ( grad(i) *=
weights(i));
9836 add(sum, grad, -dsum*sum*sum,
shape_x, grad);
9846 for (
int o = 0, j = 0; j <=
Order; j++)
9849 for (
int i = 0; i <=
Order; i++, o++)
9861 double sum, dsum[2];
9869 sum = dsum[0] = dsum[1] = 0.0;
9870 for (
int o = 0, j = 0; j <=
Order; j++)
9873 for (
int i = 0; i <=
Order; i++, o++)
9886 for (
int o = 0; o <
Dof; o++)
9888 dshape(o,0) = dshape(o,0)*sum -
u(o)*dsum[0];
9889 dshape(o,1) = dshape(o,1)*sum -
u(o)*dsum[1];
9901 for (
int o = 0, k = 0; k <=
Order; k++)
9904 for (
int j = 0; j <=
Order; j++)
9906 const double sy_sz =
shape_y(j)*sz;
9907 for (
int i = 0; i <=
Order; i++, o++)
9920 double sum, dsum[3];
9930 sum = dsum[0] = dsum[1] = dsum[2] = 0.0;
9931 for (
int o = 0, k = 0; k <=
Order; k++)
9934 for (
int j = 0; j <=
Order; j++)
9936 const double sy_sz =
shape_y(j)* sz;
9937 const double dsy_sz =
dshape_y(j)* sz;
9938 const double sy_dsz =
shape_y(j)*dsz;
9939 for (
int i = 0; i <=
Order; i++, o++)
9955 for (
int o = 0; o <
Dof; o++)
9957 dshape(o,0) = dshape(o,0)*sum -
u(o)*dsum[0];
9958 dshape(o,1) = dshape(o,1)*sum -
u(o)*dsum[1];
9959 dshape(o,2) = dshape(o,2)*sum -
u(o)*dsum[2];
Abstract class for Finite Elements.
RefinedLinear3DFiniteElement()
Construct a quadratic FE on tetrahedron.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
void Set(const double *p, const int dim)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
Multiply a matrix A with the transpose of a matrix B: A*Bt.
RT0TriangleFiniteElement()
int Size() const
Logical size of the array.
RT_QuadrilateralElement(const int p)
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
RT2TriangleFiniteElement()
int GetDim() const
Returns the space dimension for the finite element.
GaussQuad2DFiniteElement()
virtual void ProjectDelta(int vertex, Vector &dofs) const
L2_TetrahedronElement(const int p, const int _type=0)
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
CrouzeixRaviartQuadFiniteElement()
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Linear3DFiniteElement()
Construct a linear FE on tetrahedron.
static double CalcDelta(const int p, const double x)
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
RT1TriangleFiniteElement()
H1_HexahedronElement(const int p)
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
RefinedBiLinear2DFiniteElement()
Construct a biquadratic FE on quadrilateral.
ND_QuadrilateralElement(const int p)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
L2Pos_TriangleElement(const int p)
H1Pos_SegmentElement(const int p)
virtual void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
L2Pos_TetrahedronElement(const int p)
void CalcVShape_RT(ElementTransformation &Trans, DenseMatrix &shape) const
void SetSize(int s)
Resizes the vector if the new size is different.
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
void LocalInterpolation_ND(const double *tk, const Array< int > &d2t, ElementTransformation &Trans, DenseMatrix &I) const
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &h) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
static void UniformPoints(const int p, double *x)
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
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
double InnerProduct(const double *x, const double *y) const
Compute y^t A x.
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
LagrangeHexFiniteElement(int degree)
BiQuadPos2DFiniteElement()
int GetOrder() const
Returns the order of the finite element.
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
GaussBiQuad2DFiniteElement()
virtual void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Basis & OpenBasis(const int p)
const double * ClosedPoints(const int p)
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
RefinedTriLinear3DFiniteElement()
Construct a biquadratic FE on quadrilateral.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Quadratic3DFiniteElement()
Construct a quadratic FE on tetrahedron.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
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
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
const IntegrationPoint & GetCenter(int GeomType)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
void ProjectCurl_RT(const double *nk, const Array< int > &d2n, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
H1Pos_HexahedronElement(const int p)
virtual void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
RT_HexahedronElement(const int p)
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &h) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Linear2DFiniteElement()
Construct a linear FE on triangle.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
void add(const Vector &v1, const Vector &v2, Vector &v)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
L2_SegmentElement(const int p, const int _type=0)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
L2Pos_HexahedronElement(const int p)
ND_TriangleElement(const int p)
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
static void CalcBasis(const int p, const double x, double *u)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
evaluate shape function - constant 1
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &h) const
static void CalcBinomTerms(const int p, const double x, const double y, double *u)
Compute the terms in the expansion of the binomial (x + y)^p.
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &h) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Cubic3DFiniteElement()
Construct a cubic FE on tetrahedron.
Linear1DFiniteElement()
Construct a linear FE on interval.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
virtual void ProjectCurl(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
CrouzeixRaviartFiniteElement()
friend void Mult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A = B * C.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
void SetSize(int m, int n)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat)
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
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.
H1_TriangleElement(const int p)
RefinedLinear1DFiniteElement()
Construct a quadratic FE on interval.
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
void LocalInterpolation_RT(const double *nk, const Array< int > &d2n, ElementTransformation &Trans, DenseMatrix &I) const
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
static void ChebyshevPoints(const int p, double *x)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
ND_TetrahedronElement(const int p)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Nedelec1TetFiniteElement()
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
GaussLinear2DFiniteElement()
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
void Invert()
Replaces the current matrix with its inverse.
const IntegrationRule & GetNodes() const
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
L2Pos_QuadrilateralElement(const int p)
P1TetNonConfFiniteElement()
P0SegmentFiniteElement(int Ord=0)
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Nedelec1HexFiniteElement()
void Set2(const double x1, const double x2)
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
int GetVDim()
Returns dimension of the vector.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &h) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
void ProjectGrad_ND(const double *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
void Set3(const double x1, const double x2, const double x3)
ND_HexahedronElement(const int p)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
double * Data() const
Returns vector of the elements.
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
H1_QuadrilateralElement(const int p)
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
~LagrangeHexFiniteElement()
void CalcVShape_ND(ElementTransformation &Trans, DenseMatrix &shape) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
int GetDof() const
Returns the degrees of freedom in the FE space.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Base class Coefficient that may optionally depend on time.
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
void mfem_error(const char *msg)
void SetSize(int nsize)
Change 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
static const int * Binom(const int p)
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
static void GaussLobattoPoints(const int p, double *x)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
void SetDataAndSize(double *d, int s)
FiniteElement(int D, int G, int Do, int O, int F=FunctionSpace::Pk)
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Quad1DFiniteElement()
Construct a quadratic FE on interval.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
H1Pos_QuadrilateralElement(const int p)
Class for integration point with weight.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
void Project_RT(const double *nk, const Array< int > &d2n, VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
virtual void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &elmat)
virtual void ProjectDelta(int vertex, Vector &dofs) const
RefinedLinear2DFiniteElement()
Construct a quadratic FE on triangle.
void NodalLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const NodalFiniteElement &fine_fe) const
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
RT_TetrahedronElement(const int p)
L2_TriangleElement(const int p, const int _type=0)
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
H1_TetrahedronElement(const int p)
BiQuad2DFiniteElement()
Construct a biquadratic FE on quadrilateral.
virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const
Basis & ClosedBasis(const int p)
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
P0TriangleFiniteElement()
Construct P0 triangle finite element.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
RT_TriangleElement(const int p)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
void Eval(const double x, Vector &u) const
L2Pos_SegmentElement(const int p)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
static void CalcBernstein(const int p, const double x, double *u)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Quad2DFiniteElement()
Construct a quadratic FE on triangle.
Describes the space on each element.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
const double * OpenPoints(const int p)
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
GaussBiLinear2DFiniteElement()
static void GaussPoints(const int p, double *x)
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
H1_SegmentElement(const int p)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
void SetSize(int s)
If the matrix is not a square matrix of size s then recreate it.
virtual void ProjectDelta(int vertex, Vector &dofs) const
L2_HexahedronElement(const int p, const int _type=0)
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
RotTriLinearHexFiniteElement()
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
void Project_ND(const double *tk, const Array< int > &d2t, VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
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
Basis(const int p, const double *nodes, const int _mode=1)
static void CalcDBinomTerms(const int p, const double x, const double y, double *d)
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
TriLinear3DFiniteElement()
Construct a tri-linear FE on cube.
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
BiLinear2DFiniteElement()
Construct a bilinear FE on quadrilateral.
L2_QuadrilateralElement(const int p, const int _type=0)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Lagrange1DFiniteElement(int degree)
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
void ProjectGrad_RT(const double *nk, const Array< int > &d2n, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const