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)
167 for (
int i = 0; i <
Dof; i++)
170 for (
int j = 0; j < fe.
GetDof(); j++)
172 curl(i,j) = curl_shape(j,0);
180 for (
int i = 0; i <
Dof; i++)
186 dofs(i) = coeff.
Eval (Trans, ip);
189 dofs(i) *= Trans.
Weight();
197 MFEM_ASSERT(vc.
GetVDim() <= 3,
"");
202 for (
int i = 0; i <
Dof; i++)
206 vc.
Eval (x, Trans, ip);
211 for (
int j = 0; j < x.Size(); j++)
213 dofs(Dof*j+i) = v[j];
228 for (
int k = 0; k <
Dof; k++)
231 for (
int j = 0; j < shape.Size(); j++)
233 I(k,j) = (fabs(shape(j)) < 1e-12) ? 0.0 : shape(j);
242 for (
int k = 0; k <
Dof; k++)
250 for (
int j = 0; j < vshape.Height(); j++)
251 for (
int d = 0; d < vshape.Width(); d++)
253 I(k+d*Dof,j) = vshape(j,d);
269 for (
int k = 0; k <
Dof; k++)
275 Mult(dshape, Jinv, grad_k);
280 for (
int j = 0; j < grad_k.Height(); j++)
281 for (
int d = 0; d <
Dim; d++)
283 grad(k+d*Dof,j) = grad_k(j,d);
296 for (
int k = 0; k <
Dof; k++)
304 for (
int j = 0; j < div_shape.Size(); j++)
306 div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j)/detJ;
311 for (
int j = 0; j < div_shape.Size(); j++)
313 div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j);
323 for (
int i = 0; i <
Dof; i++)
327 dofs(i) = coeff.
Eval(Trans, ip);
353 pos_mass_inv.
Mult(mixed_mass, I);
358 void VectorFiniteElement::CalcShape (
361 mfem_error (
"Error: Cannot use scalar CalcShape(...) function with\n"
362 " VectorFiniteElements!");
365 void VectorFiniteElement::CalcDShape (
366 const IntegrationPoint &ip, DenseMatrix &dshape )
const
368 mfem_error (
"Error: Cannot use scalar CalcDShape(...) function with\n"
369 " VectorFiniteElements!");
375 #ifdef MFEM_THREAD_SAFE
382 shape *= (1.0 / Trans.
Weight());
390 #ifdef MFEM_THREAD_SAFE
410 #ifdef MFEM_THREAD_SAFE
417 for (
int k = 0; k <
Dof; k++)
427 if (!square_J) { dofs(k) /= Trans.
Weight(); }
440 #ifdef MFEM_THREAD_SAFE
445 for (
int k = 0; k <
Dof; k++)
455 double w = 1.0/Trans.
Weight();
456 for (
int d = 0; d <
Dim; d++)
462 for (
int j = 0; j < shape.Size(); j++)
469 for (
int d = 0; d < sdim; d++)
471 I(k,j+d*shape.Size()) = s*vk[d];
478 mfem_error(
"VectorFiniteElement::Project_RT (fe version)");
488 mfem_error(
"VectorFiniteElement::ProjectGrad_RT works only in 2D!");
496 for (
int k = 0; k <
Dof; k++)
499 tk[0] = nk[d2n[k]*
Dim+1];
500 tk[1] = -nk[d2n[k]*
Dim];
501 dshape.Mult(tk, grad_k);
502 for (
int j = 0; j < grad_k.Size(); j++)
504 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
513 #ifdef MFEM_THREAD_SAFE
526 for (
int k = 0; k <
Dof; k++)
533 J *= 1.0 / Trans.
Weight();
540 for (
int j = 0; j < curl_k.Size(); j++)
542 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
555 for (
int k = 0; k <
Dof; k++)
558 curl_shape.Mult(nk + d2n[k]*
Dim, curl_k);
559 for (
int j = 0; j < curl_k.Size(); j++)
561 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
573 for (
int k = 0; k <
Dof; k++)
594 for (
int k = 0; k <
Dof; k++)
603 double w = 1.0/Trans.
Weight();
604 for (
int d = 0; d < sdim; d++)
610 for (
int j = 0; j < shape.Size(); j++)
617 for (
int d = 0; d < sdim; d++)
619 I(k, j + d*shape.Size()) = s*vk[d];
626 mfem_error(
"VectorFiniteElement::Project_ND (fe version)");
640 for (
int k = 0; k <
Dof; k++)
643 dshape.Mult(tk + d2t[k]*
Dim, grad_k);
644 for (
int j = 0; j < grad_k.Size(); j++)
646 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
658 #ifdef MFEM_THREAD_SAFE
667 for (
int k = 0; k <
Dof; k++)
673 Jinv.
Mult(nk + d2n[k]*
Dim, vk);
675 for (
int j = 0; j <
Dof; j++)
678 for (
int i = 0; i <
Dim; i++)
680 Ikj +=
vshape(j, i) * vk[i];
682 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
694 #ifdef MFEM_THREAD_SAFE
701 for (
int k = 0; k <
Dof; k++)
709 for (
int j = 0; j <
Dof; j++)
712 for (
int i = 0; i <
Dim; i++)
714 Ikj +=
vshape(j, i) * vk[i];
716 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
750 shape(0) = 1. - ip.
x;
775 shape(0) = 1. - ip.
x - ip.
y;
783 dshape(0,0) = -1.; dshape(0,1) = -1.;
784 dshape(1,0) = 1.; dshape(1,1) = 0.;
785 dshape(2,0) = 0.; dshape(2,1) = 1.;
804 shape(0) = (1. - ip.
x) * (1. - ip.
y) ;
805 shape(1) = ip.
x * (1. - ip.
y) ;
806 shape(2) = ip.
x * ip.
y ;
807 shape(3) = (1. - ip.
x) * ip.
y ;
813 dshape(0,0) = -1. + ip.
y; dshape(0,1) = -1. + ip.
x ;
814 dshape(1,0) = 1. - ip.
y; dshape(1,1) = -ip.
x ;
815 dshape(2,0) = ip.
y ; dshape(2,1) = ip.
x ;
816 dshape(3,0) = -ip.
y ; dshape(3,1) = 1. - ip.
x ;
822 h( 0,0) = 0.; h( 0,1) = 1.; h( 0,2) = 0.;
823 h( 1,0) = 0.; h( 1,1) = -1.; h( 1,2) = 0.;
824 h( 2,0) = 0.; h( 2,1) = 1.; h( 2,2) = 0.;
825 h( 3,0) = 0.; h( 3,1) = -1.; h( 3,2) = 0.;
843 const double x = ip.
x, y = ip.
y;
845 shape(0) = 5./3. - 2. * (x + y);
846 shape(1) = 2. * (x - 1./6.);
847 shape(2) = 2. * (y - 1./6.);
853 dshape(0,0) = -2.; dshape(0,1) = -2.;
854 dshape(1,0) = 2.; dshape(1,1) = 0.;
855 dshape(2,0) = 0.; dshape(2,1) = 2.;
860 dofs(vertex) = 2./3.;
861 dofs((vertex+1)%3) = 1./6.;
862 dofs((vertex+2)%3) = 1./6.;
867 const double GaussBiLinear2DFiniteElement::p[] =
868 { 0.2113248654051871177454256, 0.7886751345948128822545744 };
886 const double x = ip.
x, y = ip.
y;
888 shape(0) = 3. * (p[1] - x) * (p[1] - y);
889 shape(1) = 3. * (x - p[0]) * (p[1] - y);
890 shape(2) = 3. * (x - p[0]) * (y - p[0]);
891 shape(3) = 3. * (p[1] - x) * (y - p[0]);
897 const double x = ip.
x, y = ip.
y;
899 dshape(0,0) = 3. * (y - p[1]); dshape(0,1) = 3. * (x - p[1]);
900 dshape(1,0) = 3. * (p[1] - y); dshape(1,1) = 3. * (p[0] - x);
901 dshape(2,0) = 3. * (y - p[0]); dshape(2,1) = 3. * (x - p[0]);
902 dshape(3,0) = 3. * (p[0] - y); dshape(3,1) = 3. * (p[1] - x);
908 dofs(vertex) = p[1]*p[1];
909 dofs((vertex+1)%4) = p[0]*p[1];
910 dofs((vertex+2)%4) = p[0]*p[0];
911 dofs((vertex+3)%4) = p[0]*p[1];
932 shape(0) = 1. - ip.
x - ip.
y;
940 dshape(0,0) = -1.; dshape(0,1) = -1.;
941 dshape(1,0) = 1.; dshape(1,1) = 0.;
942 dshape(2,0) = 0.; dshape(2,1) = 1.;
958 double l1 = 1.0 - x, l2 = x, l3 = 2. * x - 1.;
960 shape(0) = l1 * (-l3);
962 shape(2) = 4. * l1 * l2;
970 dshape(0,0) = 4. * x - 3.;
971 dshape(1,0) = 4. * x - 1.;
972 dshape(2,0) = 4. - 8. * x;
987 const double x = ip.
x, x1 = 1. - x;
991 shape(2) = 2. * x * x1;
997 const double x = ip.
x;
999 dshape(0,0) = 2. * x - 2.;
1000 dshape(1,0) = 2. * x;
1001 dshape(2,0) = 2. - 4. * x;
1024 double x = ip.
x, y = ip.
y;
1025 double l1 = 1.-x-y, l2 = x, l3 = y;
1027 shape(0) = l1 * (2. * l1 - 1.);
1028 shape(1) = l2 * (2. * l2 - 1.);
1029 shape(2) = l3 * (2. * l3 - 1.);
1030 shape(3) = 4. * l1 * l2;
1031 shape(4) = 4. * l2 * l3;
1032 shape(5) = 4. * l3 * l1;
1038 double x = ip.
x, y = ip.
y;
1041 dshape(0,1) = 4. * (x + y) - 3.;
1043 dshape(1,0) = 4. * x - 1.;
1047 dshape(2,1) = 4. * y - 1.;
1049 dshape(3,0) = -4. * (2. * x + y - 1.);
1050 dshape(3,1) = -4. * x;
1052 dshape(4,0) = 4. * y;
1053 dshape(4,1) = 4. * x;
1055 dshape(5,0) = -4. * y;
1056 dshape(5,1) = -4. * (x + 2. * y - 1.);
1096 case 0: dofs(3) = 0.25; dofs(5) = 0.25;
break;
1097 case 1: dofs(3) = 0.25; dofs(4) = 0.25;
break;
1098 case 2: dofs(4) = 0.25; dofs(5) = 0.25;
break;
1104 const double GaussQuad2DFiniteElement::p[] =
1105 { 0.0915762135097707434595714634022015, 0.445948490915964886318329253883051 };
1123 for (
int i = 0; i < 6; i++)
1140 const double x = ip.
x, y = ip.
y;
1154 const double x = ip.
x, y = ip.
y;
1155 D(0,0) = 0.; D(0,1) = 0.;
1156 D(1,0) = 1.; D(1,1) = 0.;
1157 D(2,0) = 0.; D(2,1) = 1.;
1158 D(3,0) = 2. * x; D(3,1) = 0.;
1159 D(4,0) = y; D(4,1) = x;
1160 D(5,0) = 0.; D(5,1) = 2. * y;
1192 double x = ip.
x, y = ip.
y;
1193 double l1x, l2x, l3x, l1y, l2y, l3y;
1195 l1x = (x - 1.) * (2. * x - 1);
1196 l2x = 4. * x * (1. - x);
1197 l3x = x * (2. * x - 1.);
1198 l1y = (y - 1.) * (2. * y - 1);
1199 l2y = 4. * y * (1. - y);
1200 l3y = y * (2. * y - 1.);
1202 shape(0) = l1x * l1y;
1203 shape(4) = l2x * l1y;
1204 shape(1) = l3x * l1y;
1205 shape(7) = l1x * l2y;
1206 shape(8) = l2x * l2y;
1207 shape(5) = l3x * l2y;
1208 shape(3) = l1x * l3y;
1209 shape(6) = l2x * l3y;
1210 shape(2) = l3x * l3y;
1216 double x = ip.
x, y = ip.
y;
1217 double l1x, l2x, l3x, l1y, l2y, l3y;
1218 double d1x, d2x, d3x, d1y, d2y, d3y;
1220 l1x = (x - 1.) * (2. * x - 1);
1221 l2x = 4. * x * (1. - x);
1222 l3x = x * (2. * x - 1.);
1223 l1y = (y - 1.) * (2. * y - 1);
1224 l2y = 4. * y * (1. - y);
1225 l3y = y * (2. * y - 1.);
1234 dshape(0,0) = d1x * l1y;
1235 dshape(0,1) = l1x * d1y;
1237 dshape(4,0) = d2x * l1y;
1238 dshape(4,1) = l2x * d1y;
1240 dshape(1,0) = d3x * l1y;
1241 dshape(1,1) = l3x * d1y;
1243 dshape(7,0) = d1x * l2y;
1244 dshape(7,1) = l1x * d2y;
1246 dshape(8,0) = d2x * l2y;
1247 dshape(8,1) = l2x * d2y;
1249 dshape(5,0) = d3x * l2y;
1250 dshape(5,1) = l3x * d2y;
1252 dshape(3,0) = d1x * l3y;
1253 dshape(3,1) = l1x * d3y;
1255 dshape(6,0) = d2x * l3y;
1256 dshape(6,1) = l2x * d3y;
1258 dshape(2,0) = d3x * l3y;
1259 dshape(2,1) = l3x * d3y;
1271 case 0: dofs(4) = 0.25; dofs(7) = 0.25;
break;
1272 case 1: dofs(4) = 0.25; dofs(5) = 0.25;
break;
1273 case 2: dofs(5) = 0.25; dofs(6) = 0.25;
break;
1274 case 3: dofs(6) = 0.25; dofs(7) = 0.25;
break;
1306 double x = ip.
x, y = ip.
y;
1307 double l1x, l2x, l3x, l1y, l2y, l3y;
1309 l1x = (1. - x) * (1. - x);
1310 l2x = 2. * x * (1. - x);
1312 l1y = (1. - y) * (1. - y);
1313 l2y = 2. * y * (1. - y);
1316 shape(0) = l1x * l1y;
1317 shape(4) = l2x * l1y;
1318 shape(1) = l3x * l1y;
1319 shape(7) = l1x * l2y;
1320 shape(8) = l2x * l2y;
1321 shape(5) = l3x * l2y;
1322 shape(3) = l1x * l3y;
1323 shape(6) = l2x * l3y;
1324 shape(2) = l3x * l3y;
1330 double x = ip.
x, y = ip.
y;
1331 double l1x, l2x, l3x, l1y, l2y, l3y;
1332 double d1x, d2x, d3x, d1y, d2y, d3y;
1334 l1x = (1. - x) * (1. - x);
1335 l2x = 2. * x * (1. - x);
1337 l1y = (1. - y) * (1. - y);
1338 l2y = 2. * y * (1. - y);
1348 dshape(0,0) = d1x * l1y;
1349 dshape(0,1) = l1x * d1y;
1351 dshape(4,0) = d2x * l1y;
1352 dshape(4,1) = l2x * d1y;
1354 dshape(1,0) = d3x * l1y;
1355 dshape(1,1) = l3x * d1y;
1357 dshape(7,0) = d1x * l2y;
1358 dshape(7,1) = l1x * d2y;
1360 dshape(8,0) = d2x * l2y;
1361 dshape(8,1) = l2x * d2y;
1363 dshape(5,0) = d3x * l2y;
1364 dshape(5,1) = l3x * d2y;
1366 dshape(3,0) = d1x * l3y;
1367 dshape(3,1) = l1x * d3y;
1369 dshape(6,0) = d2x * l3y;
1370 dshape(6,1) = l2x * d3y;
1372 dshape(2,0) = d3x * l3y;
1373 dshape(2,1) = l3x * d3y;
1381 Vector xx(&tr_ip.
x, 2), shape(s, 9);
1383 for (
int i = 0; i < 9; i++)
1387 for (
int j = 0; j < 9; j++)
1388 if (fabs(I(i,j) = s[j]) < 1.0e-12)
1393 for (
int i = 0; i < 9; i++)
1395 double *d = &I(0,i);
1396 d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
1397 d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
1398 d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
1399 d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
1400 d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
1401 0.25 * (d[0] + d[1] + d[2] + d[3]);
1410 for (
int i = 0; i < 9; i++)
1414 d[i] = coeff.
Eval(Trans, ip);
1416 d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
1417 d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
1418 d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
1419 d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
1420 d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
1421 0.25 * (d[0] + d[1] + d[2] + d[3]);
1431 for (
int i = 0; i < 9; i++)
1435 vc.
Eval (x, Trans, ip);
1436 for (
int j = 0; j < x.Size(); j++)
1441 for (
int j = 0; j < x.Size(); j++)
1443 double *d = &dofs(9*j);
1445 d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
1446 d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
1447 d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
1448 d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
1449 d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
1450 0.25 * (d[0] + d[1] + d[2] + d[3]);
1458 const double p1 = 0.5*(1.-sqrt(3./5.));
1483 const double a = sqrt(5./3.);
1484 const double p1 = 0.5*(1.-sqrt(3./5.));
1486 double x = a*(ip.
x-p1), y = a*(ip.
y-p1);
1487 double l1x, l2x, l3x, l1y, l2y, l3y;
1489 l1x = (x - 1.) * (2. * x - 1);
1490 l2x = 4. * x * (1. - x);
1491 l3x = x * (2. * x - 1.);
1492 l1y = (y - 1.) * (2. * y - 1);
1493 l2y = 4. * y * (1. - y);
1494 l3y = y * (2. * y - 1.);
1496 shape(0) = l1x * l1y;
1497 shape(4) = l2x * l1y;
1498 shape(1) = l3x * l1y;
1499 shape(7) = l1x * l2y;
1500 shape(8) = l2x * l2y;
1501 shape(5) = l3x * l2y;
1502 shape(3) = l1x * l3y;
1503 shape(6) = l2x * l3y;
1504 shape(2) = l3x * l3y;
1510 const double a = sqrt(5./3.);
1511 const double p1 = 0.5*(1.-sqrt(3./5.));
1513 double x = a*(ip.
x-p1), y = a*(ip.
y-p1);
1514 double l1x, l2x, l3x, l1y, l2y, l3y;
1515 double d1x, d2x, d3x, d1y, d2y, d3y;
1517 l1x = (x - 1.) * (2. * x - 1);
1518 l2x = 4. * x * (1. - x);
1519 l3x = x * (2. * x - 1.);
1520 l1y = (y - 1.) * (2. * y - 1);
1521 l2y = 4. * y * (1. - y);
1522 l3y = y * (2. * y - 1.);
1524 d1x = a * (4. * x - 3.);
1525 d2x = a * (4. - 8. * x);
1526 d3x = a * (4. * x - 1.);
1527 d1y = a * (4. * y - 3.);
1528 d2y = a * (4. - 8. * y);
1529 d3y = a * (4. * y - 1.);
1531 dshape(0,0) = d1x * l1y;
1532 dshape(0,1) = l1x * d1y;
1534 dshape(4,0) = d2x * l1y;
1535 dshape(4,1) = l2x * d1y;
1537 dshape(1,0) = d3x * l1y;
1538 dshape(1,1) = l3x * d1y;
1540 dshape(7,0) = d1x * l2y;
1541 dshape(7,1) = l1x * d2y;
1543 dshape(8,0) = d2x * l2y;
1544 dshape(8,1) = l2x * d2y;
1546 dshape(5,0) = d3x * l2y;
1547 dshape(5,1) = l3x * d2y;
1549 dshape(3,0) = d1x * l3y;
1550 dshape(3,1) = l1x * d3y;
1552 dshape(6,0) = d2x * l3y;
1553 dshape(6,1) = l2x * d3y;
1555 dshape(2,0) = d3x * l3y;
1556 dshape(2,1) = l3x * d3y;
1599 double x = ip.
x, y = ip.
y;
1601 double w1x, w2x, w3x, w1y, w2y, w3y;
1602 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
1604 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
1605 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
1607 l0x = (- 4.5) * w1x * w2x * w3x;
1608 l1x = ( 13.5) * x * w2x * w3x;
1609 l2x = (-13.5) * x * w1x * w3x;
1610 l3x = ( 4.5) * x * w1x * w2x;
1612 l0y = (- 4.5) * w1y * w2y * w3y;
1613 l1y = ( 13.5) * y * w2y * w3y;
1614 l2y = (-13.5) * y * w1y * w3y;
1615 l3y = ( 4.5) * y * w1y * w2y;
1617 shape(0) = l0x * l0y;
1618 shape(1) = l3x * l0y;
1619 shape(2) = l3x * l3y;
1620 shape(3) = l0x * l3y;
1621 shape(4) = l1x * l0y;
1622 shape(5) = l2x * l0y;
1623 shape(6) = l3x * l1y;
1624 shape(7) = l3x * l2y;
1625 shape(8) = l2x * l3y;
1626 shape(9) = l1x * l3y;
1627 shape(10) = l0x * l2y;
1628 shape(11) = l0x * l1y;
1629 shape(12) = l1x * l1y;
1630 shape(13) = l2x * l1y;
1631 shape(14) = l1x * l2y;
1632 shape(15) = l2x * l2y;
1638 double x = ip.
x, y = ip.
y;
1640 double w1x, w2x, w3x, w1y, w2y, w3y;
1641 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
1642 double d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
1644 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
1645 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
1647 l0x = (- 4.5) * w1x * w2x * w3x;
1648 l1x = ( 13.5) * x * w2x * w3x;
1649 l2x = (-13.5) * x * w1x * w3x;
1650 l3x = ( 4.5) * x * w1x * w2x;
1652 l0y = (- 4.5) * w1y * w2y * w3y;
1653 l1y = ( 13.5) * y * w2y * w3y;
1654 l2y = (-13.5) * y * w1y * w3y;
1655 l3y = ( 4.5) * y * w1y * w2y;
1657 d0x = -5.5 + ( 18. - 13.5 * x) * x;
1658 d1x = 9. + (-45. + 40.5 * x) * x;
1659 d2x = -4.5 + ( 36. - 40.5 * x) * x;
1660 d3x = 1. + (- 9. + 13.5 * x) * x;
1662 d0y = -5.5 + ( 18. - 13.5 * y) * y;
1663 d1y = 9. + (-45. + 40.5 * y) * y;
1664 d2y = -4.5 + ( 36. - 40.5 * y) * y;
1665 d3y = 1. + (- 9. + 13.5 * y) * y;
1667 dshape( 0,0) = d0x * l0y; dshape( 0,1) = l0x * d0y;
1668 dshape( 1,0) = d3x * l0y; dshape( 1,1) = l3x * d0y;
1669 dshape( 2,0) = d3x * l3y; dshape( 2,1) = l3x * d3y;
1670 dshape( 3,0) = d0x * l3y; dshape( 3,1) = l0x * d3y;
1671 dshape( 4,0) = d1x * l0y; dshape( 4,1) = l1x * d0y;
1672 dshape( 5,0) = d2x * l0y; dshape( 5,1) = l2x * d0y;
1673 dshape( 6,0) = d3x * l1y; dshape( 6,1) = l3x * d1y;
1674 dshape( 7,0) = d3x * l2y; dshape( 7,1) = l3x * d2y;
1675 dshape( 8,0) = d2x * l3y; dshape( 8,1) = l2x * d3y;
1676 dshape( 9,0) = d1x * l3y; dshape( 9,1) = l1x * d3y;
1677 dshape(10,0) = d0x * l2y; dshape(10,1) = l0x * d2y;
1678 dshape(11,0) = d0x * l1y; dshape(11,1) = l0x * d1y;
1679 dshape(12,0) = d1x * l1y; dshape(12,1) = l1x * d1y;
1680 dshape(13,0) = d2x * l1y; dshape(13,1) = l2x * d1y;
1681 dshape(14,0) = d1x * l2y; dshape(14,1) = l1x * d2y;
1682 dshape(15,0) = d2x * l2y; dshape(15,1) = l2x * d2y;
1688 double x = ip.
x, y = ip.
y;
1690 double w1x, w2x, w3x, w1y, w2y, w3y;
1691 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
1692 double d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
1693 double h0x, h1x, h2x, h3x, h0y, h1y, h2y, h3y;
1695 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
1696 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
1698 l0x = (- 4.5) * w1x * w2x * w3x;
1699 l1x = ( 13.5) * x * w2x * w3x;
1700 l2x = (-13.5) * x * w1x * w3x;
1701 l3x = ( 4.5) * x * w1x * w2x;
1703 l0y = (- 4.5) * w1y * w2y * w3y;
1704 l1y = ( 13.5) * y * w2y * w3y;
1705 l2y = (-13.5) * y * w1y * w3y;
1706 l3y = ( 4.5) * y * w1y * w2y;
1708 d0x = -5.5 + ( 18. - 13.5 * x) * x;
1709 d1x = 9. + (-45. + 40.5 * x) * x;
1710 d2x = -4.5 + ( 36. - 40.5 * x) * x;
1711 d3x = 1. + (- 9. + 13.5 * x) * x;
1713 d0y = -5.5 + ( 18. - 13.5 * y) * y;
1714 d1y = 9. + (-45. + 40.5 * y) * y;
1715 d2y = -4.5 + ( 36. - 40.5 * y) * y;
1716 d3y = 1. + (- 9. + 13.5 * y) * y;
1718 h0x = -27. * x + 18.;
1719 h1x = 81. * x - 45.;
1720 h2x = -81. * x + 36.;
1723 h0y = -27. * y + 18.;
1724 h1y = 81. * y - 45.;
1725 h2y = -81. * y + 36.;
1728 h( 0,0) = h0x * l0y; h( 0,1) = d0x * d0y; h( 0,2) = l0x * h0y;
1729 h( 1,0) = h3x * l0y; h( 1,1) = d3x * d0y; h( 1,2) = l3x * h0y;
1730 h( 2,0) = h3x * l3y; h( 2,1) = d3x * d3y; h( 2,2) = l3x * h3y;
1731 h( 3,0) = h0x * l3y; h( 3,1) = d0x * d3y; h( 3,2) = l0x * h3y;
1732 h( 4,0) = h1x * l0y; h( 4,1) = d1x * d0y; h( 4,2) = l1x * h0y;
1733 h( 5,0) = h2x * l0y; h( 5,1) = d2x * d0y; h( 5,2) = l2x * h0y;
1734 h( 6,0) = h3x * l1y; h( 6,1) = d3x * d1y; h( 6,2) = l3x * h1y;
1735 h( 7,0) = h3x * l2y; h( 7,1) = d3x * d2y; h( 7,2) = l3x * h2y;
1736 h( 8,0) = h2x * l3y; h( 8,1) = d2x * d3y; h( 8,2) = l2x * h3y;
1737 h( 9,0) = h1x * l3y; h( 9,1) = d1x * d3y; h( 9,2) = l1x * h3y;
1738 h(10,0) = h0x * l2y; h(10,1) = d0x * d2y; h(10,2) = l0x * h2y;
1739 h(11,0) = h0x * l1y; h(11,1) = d0x * d1y; h(11,2) = l0x * h1y;
1740 h(12,0) = h1x * l1y; h(12,1) = d1x * d1y; h(12,2) = l1x * h1y;
1741 h(13,0) = h2x * l1y; h(13,1) = d2x * d1y; h(13,2) = l2x * h1y;
1742 h(14,0) = h1x * l2y; h(14,1) = d1x * d2y; h(14,2) = l1x * h2y;
1743 h(15,0) = h2x * l2y; h(15,1) = d2x * d2y; h(15,2) = l2x * h2y;
1762 l3 = (0.33333333333333333333-x),
1763 l4 = (0.66666666666666666667-x);
1765 shape(0) = 4.5 * l2 * l3 * l4;
1766 shape(1) = 4.5 * l1 * l3 * l4;
1767 shape(2) = 13.5 * l1 * l2 * l4;
1768 shape(3) = -13.5 * l1 * l2 * l3;
1776 dshape(0,0) = -5.5 + x * (18. - 13.5 * x);
1777 dshape(1,0) = 1. - x * (9. - 13.5 * x);
1778 dshape(2,0) = 9. - x * (45. - 40.5 * x);
1779 dshape(3,0) = -4.5 + x * (36. - 40.5 * x);
1811 double x = ip.
x, y = ip.
y;
1812 double l1 = (-1. + x + y),
1816 shape(0) = -0.5*l1*(3.*l1 + 1.)*(3.*l1 + 2.);
1817 shape(1) = 0.5*x*(lx - 1.)*lx;
1818 shape(2) = 0.5*y*(-1. + ly)*ly;
1819 shape(3) = 4.5*x*l1*(3.*l1 + 1.);
1820 shape(4) = -4.5*x*lx*l1;
1821 shape(5) = 4.5*x*lx*y;
1822 shape(6) = 4.5*x*y*ly;
1823 shape(7) = -4.5*y*l1*ly;
1824 shape(8) = 4.5*y*l1*(1. + 3.*l1);
1825 shape(9) = -27.*x*y*l1;
1831 double x = ip.
x, y = ip.
y;
1833 dshape(0,0) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
1834 dshape(1,0) = 1. + 4.5*x*(-2. + 3.*x);
1836 dshape(3,0) = 4.5*(2. + 9.*x*x - 5.*y + 3.*y*y + 2.*x*(-5. + 6.*y));
1837 dshape(4,0) = -4.5*(1. - 1.*y + x*(-8. + 9.*x + 6.*y));
1838 dshape(5,0) = 4.5*(-1. + 6.*x)*y;
1839 dshape(6,0) = 4.5*y*(-1. + 3.*y);
1840 dshape(7,0) = 4.5*(1. - 3.*y)*y;
1841 dshape(8,0) = 4.5*y*(-5. + 6.*x + 6.*y);
1842 dshape(9,0) = -27.*y*(-1. + 2.*x + y);
1844 dshape(0,1) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
1846 dshape(2,1) = 1. + 4.5*y*(-2. + 3.*y);
1847 dshape(3,1) = 4.5*x*(-5. + 6.*x + 6.*y);
1848 dshape(4,1) = 4.5*(1. - 3.*x)*x;
1849 dshape(5,1) = 4.5*x*(-1. + 3.*x);
1850 dshape(6,1) = 4.5*x*(-1. + 6.*y);
1851 dshape(7,1) = -4.5*(1. + x*(-1. + 6.*y) + y*(-8. + 9.*y));
1852 dshape(8,1) = 4.5*(2. + 3.*x*x + y*(-10. + 9.*y) + x*(-5. + 12.*y));
1853 dshape(9,1) = -27.*x*(-1. + x + 2.*y);
1859 double x = ip.
x, y = ip.
y;
1861 h(0,0) = 18.-27.*(x+y);
1862 h(0,1) = 18.-27.*(x+y);
1863 h(0,2) = 18.-27.*(x+y);
1873 h(3,0) = -45.+81.*x+54.*y;
1874 h(3,1) = -22.5+54.*x+27.*y;
1877 h(4,0) = 36.-81.*x-27.*y;
1882 h(5,1) = -4.5+27.*x;
1886 h(6,1) = -4.5+27.*y;
1891 h(7,2) = 36.-27.*x-81.*y;
1894 h(8,1) = -22.5+27.*x+54.*y;
1895 h(8,2) = -45.+54.*x+81.*y;
1898 h(9,1) = 27.-54.*(x+y);
1971 double x = ip.
x, y = ip.
y, z = ip.
z;
1973 shape(0) = -((-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z)*
1974 (-1 + 3*x + 3*y + 3*z))/2.;
1975 shape(4) = (9*x*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
1976 shape(5) = (-9*x*(-1 + 3*x)*(-1 + x + y + z))/2.;
1977 shape(1) = (x*(2 + 9*(-1 + x)*x))/2.;
1978 shape(6) = (9*y*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
1979 shape(19) = -27*x*y*(-1 + x + y + z);
1980 shape(10) = (9*x*(-1 + 3*x)*y)/2.;
1981 shape(7) = (-9*y*(-1 + 3*y)*(-1 + x + y + z))/2.;
1982 shape(11) = (9*x*y*(-1 + 3*y))/2.;
1983 shape(2) = (y*(2 + 9*(-1 + y)*y))/2.;
1984 shape(8) = (9*z*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
1985 shape(18) = -27*x*z*(-1 + x + y + z);
1986 shape(12) = (9*x*(-1 + 3*x)*z)/2.;
1987 shape(17) = -27*y*z*(-1 + x + y + z);
1988 shape(16) = 27*x*y*z;
1989 shape(14) = (9*y*(-1 + 3*y)*z)/2.;
1990 shape(9) = (-9*z*(-1 + x + y + z)*(-1 + 3*z))/2.;
1991 shape(13) = (9*x*z*(-1 + 3*z))/2.;
1992 shape(15) = (9*y*z*(-1 + 3*z))/2.;
1993 shape(3) = (z*(2 + 9*(-1 + z)*z))/2.;
1999 double x = ip.
x, y = ip.
y, z = ip.
z;
2001 dshape(0,0) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
2002 x*(-4 + 6*y + 6*z)))/2.;
2003 dshape(0,1) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
2004 x*(-4 + 6*y + 6*z)))/2.;
2005 dshape(0,2) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
2006 x*(-4 + 6*y + 6*z)))/2.;
2007 dshape(4,0) = (9*(9*pow(x,2) + (-1 + y + z)*(-2 + 3*y + 3*z) +
2008 2*x*(-5 + 6*y + 6*z)))/2.;
2009 dshape(4,1) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
2010 dshape(4,2) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
2011 dshape(5,0) = (-9*(1 - y - z + x*(-8 + 9*x + 6*y + 6*z)))/2.;
2012 dshape(5,1) = (9*(1 - 3*x)*x)/2.;
2013 dshape(5,2) = (9*(1 - 3*x)*x)/2.;
2014 dshape(1,0) = 1 + (9*x*(-2 + 3*x))/2.;
2017 dshape(6,0) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
2018 dshape(6,1) = (9*(2 + 3*pow(x,2) - 10*y - 5*z + 3*(y + z)*(3*y + z) +
2019 x*(-5 + 12*y + 6*z)))/2.;
2020 dshape(6,2) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
2021 dshape(19,0) = -27*y*(-1 + 2*x + y + z);
2022 dshape(19,1) = -27*x*(-1 + x + 2*y + z);
2023 dshape(19,2) = -27*x*y;
2024 dshape(10,0) = (9*(-1 + 6*x)*y)/2.;
2025 dshape(10,1) = (9*x*(-1 + 3*x))/2.;
2027 dshape(7,0) = (9*(1 - 3*y)*y)/2.;
2028 dshape(7,1) = (-9*(1 + x*(-1 + 6*y) - z + y*(-8 + 9*y + 6*z)))/2.;
2029 dshape(7,2) = (9*(1 - 3*y)*y)/2.;
2030 dshape(11,0) = (9*y*(-1 + 3*y))/2.;
2031 dshape(11,1) = (9*x*(-1 + 6*y))/2.;
2034 dshape(2,1) = 1 + (9*y*(-2 + 3*y))/2.;
2036 dshape(8,0) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
2037 dshape(8,1) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
2038 dshape(8,2) = (9*(2 + 3*pow(x,2) - 5*y - 10*z + 3*(y + z)*(y + 3*z) +
2039 x*(-5 + 6*y + 12*z)))/2.;
2040 dshape(18,0) = -27*z*(-1 + 2*x + y + z);
2041 dshape(18,1) = -27*x*z;
2042 dshape(18,2) = -27*x*(-1 + x + y + 2*z);
2043 dshape(12,0) = (9*(-1 + 6*x)*z)/2.;
2045 dshape(12,2) = (9*x*(-1 + 3*x))/2.;
2046 dshape(17,0) = -27*y*z;
2047 dshape(17,1) = -27*z*(-1 + x + 2*y + z);
2048 dshape(17,2) = -27*y*(-1 + x + y + 2*z);
2049 dshape(16,0) = 27*y*z;
2050 dshape(16,1) = 27*x*z;
2051 dshape(16,2) = 27*x*y;
2053 dshape(14,1) = (9*(-1 + 6*y)*z)/2.;
2054 dshape(14,2) = (9*y*(-1 + 3*y))/2.;
2055 dshape(9,0) = (9*(1 - 3*z)*z)/2.;
2056 dshape(9,1) = (9*(1 - 3*z)*z)/2.;
2057 dshape(9,2) = (9*(-1 + x + y + 8*z - 6*(x + y)*z - 9*pow(z,2)))/2.;
2058 dshape(13,0) = (9*z*(-1 + 3*z))/2.;
2060 dshape(13,2) = (9*x*(-1 + 6*z))/2.;
2062 dshape(15,1) = (9*z*(-1 + 3*z))/2.;
2063 dshape(15,2) = (9*y*(-1 + 6*z))/2.;
2066 dshape(3,2) = 1 + (9*z*(-2 + 3*z))/2.;
2132 shape(0) = 1. - ip.
x - ip.
y - ip.
z;
2141 if (dshape.
Height() == 4)
2143 double *A = &dshape(0,0);
2144 A[0] = -1.; A[4] = -1.; A[8] = -1.;
2145 A[1] = 1.; A[5] = 0.; A[9] = 0.;
2146 A[2] = 0.; A[6] = 1.; A[10] = 0.;
2147 A[3] = 0.; A[7] = 0.; A[11] = 1.;
2151 dshape(0,0) = -1.; dshape(0,1) = -1.; dshape(0,2) = -1.;
2152 dshape(1,0) = 1.; dshape(1,1) = 0.; dshape(1,2) = 0.;
2153 dshape(2,0) = 0.; dshape(2,1) = 1.; dshape(2,2) = 0.;
2154 dshape(3,0) = 0.; dshape(3,1) = 0.; dshape(3,2) = 1.;
2161 static int face_dofs[4][3] = {{1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {0, 1, 2}};
2164 *dofs = face_dofs[face];
2206 double L0, L1, L2, L3;
2208 L0 = 1. - ip.
x - ip.
y - ip.
z;
2213 shape(0) = L0 * ( 2.0 * L0 - 1.0 );
2214 shape(1) = L1 * ( 2.0 * L1 - 1.0 );
2215 shape(2) = L2 * ( 2.0 * L2 - 1.0 );
2216 shape(3) = L3 * ( 2.0 * L3 - 1.0 );
2217 shape(4) = 4.0 * L0 * L1;
2218 shape(5) = 4.0 * L0 * L2;
2219 shape(6) = 4.0 * L0 * L3;
2220 shape(7) = 4.0 * L1 * L2;
2221 shape(8) = 4.0 * L1 * L3;
2222 shape(9) = 4.0 * L2 * L3;
2233 L0 = 1.0 - x - y - z;
2235 dshape(0,0) = dshape(0,1) = dshape(0,2) = 1.0 - 4.0 * L0;
2236 dshape(1,0) = -1.0 + 4.0 * x; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
2237 dshape(2,0) = 0.0; dshape(2,1) = -1.0 + 4.0 * y; dshape(2,2) = 0.0;
2238 dshape(3,0) = dshape(3,1) = 0.0; dshape(3,2) = -1.0 + 4.0 * z;
2239 dshape(4,0) = 4.0 * (L0 - x); dshape(4,1) = dshape(4,2) = -4.0 * x;
2240 dshape(5,0) = dshape(5,2) = -4.0 * y; dshape(5,1) = 4.0 * (L0 - y);
2241 dshape(6,0) = dshape(6,1) = -4.0 * z; dshape(6,2) = 4.0 * (L0 - z);
2242 dshape(7,0) = 4.0 * y; dshape(7,1) = 4.0 * x; dshape(7,2) = 0.0;
2243 dshape(8,0) = 4.0 * z; dshape(8,1) = 0.0; dshape(8,2) = 4.0 * x;
2244 dshape(9,0) = 0.0; dshape(9,1) = 4.0 * z; dshape(9,2) = 4.0 * y;
2286 double x = ip.
x, y = ip.
y, z = ip.
z;
2287 double ox = 1.-x, oy = 1.-y, oz = 1.-z;
2289 shape(0) = ox * oy * oz;
2290 shape(1) = x * oy * oz;
2291 shape(2) = x * y * oz;
2292 shape(3) = ox * y * oz;
2293 shape(4) = ox * oy * z;
2294 shape(5) = x * oy * z;
2295 shape(6) = x * y * z;
2296 shape(7) = ox * y * z;
2302 double x = ip.
x, y = ip.
y, z = ip.
z;
2303 double ox = 1.-x, oy = 1.-y, oz = 1.-z;
2305 dshape(0,0) = - oy * oz;
2306 dshape(0,1) = - ox * oz;
2307 dshape(0,2) = - ox * oy;
2309 dshape(1,0) = oy * oz;
2310 dshape(1,1) = - x * oz;
2311 dshape(1,2) = - x * oy;
2313 dshape(2,0) = y * oz;
2314 dshape(2,1) = x * oz;
2315 dshape(2,2) = - x * y;
2317 dshape(3,0) = - y * oz;
2318 dshape(3,1) = ox * oz;
2319 dshape(3,2) = - ox * y;
2321 dshape(4,0) = - oy * z;
2322 dshape(4,1) = - ox * z;
2323 dshape(4,2) = ox * oy;
2325 dshape(5,0) = oy * z;
2326 dshape(5,1) = - x * z;
2327 dshape(5,2) = x * oy;
2329 dshape(6,0) = y * z;
2330 dshape(6,1) = x * z;
2331 dshape(6,2) = x * y;
2333 dshape(7,0) = - y * z;
2334 dshape(7,1) = ox * z;
2335 dshape(7,2) = ox * y;
2370 shape(0) = 1.0 - 2.0 * ip.
y;
2371 shape(1) = -1.0 + 2.0 * ( ip.
x + ip.
y );
2372 shape(2) = 1.0 - 2.0 * ip.
x;
2378 dshape(0,0) = 0.0; dshape(0,1) = -2.0;
2379 dshape(1,0) = 2.0; dshape(1,1) = 2.0;
2380 dshape(2,0) = -2.0; dshape(2,1) = 0.0;
2402 const double l1 = ip.
x+ip.
y-0.5, l2 = 1.-l1, l3 = ip.
x-ip.
y+0.5, l4 = 1.-l3;
2413 const double x2 = 2.*ip.
x, y2 = 2.*ip.
y;
2415 dshape(0,0) = 1. - x2; dshape(0,1) = -2. + y2;
2416 dshape(1,0) = x2; dshape(1,1) = 1. - y2;
2417 dshape(2,0) = 1. - x2; dshape(2,1) = y2;
2418 dshape(3,0) = -2. + x2; dshape(3,1) = 1. - y2;
2436 double x = ip.
x, y = ip.
y;
2439 shape(0,1) = y - 1.;
2442 shape(2,0) = x - 1.;
2454 const double RT0TriangleFiniteElement::nk[3][2] =
2455 { {0, -1}, {1, 1}, {-1, 0} };
2461 #ifdef MFEM_THREAD_SAFE
2467 for (k = 0; k < 3; k++)
2470 for (j = 0; j < 3; j++)
2473 if (j == k) { d -= 1.0; }
2474 if (fabs(d) > 1.0e-12)
2476 cerr <<
"RT0TriangleFiniteElement::GetLocalInterpolation (...)\n"
2477 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2493 for (k = 0; k < 3; k++)
2496 ip.
x = vk[0]; ip.
y = vk[1];
2499 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2500 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2501 for (j = 0; j < 3; j++)
2502 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2515 #ifdef MFEM_THREAD_SAFE
2519 for (
int k = 0; k < 3; k++)
2527 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2528 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
2548 double x = ip.
x, y = ip.
y;
2551 shape(0,1) = y - 1.;
2556 shape(3,0) = x - 1.;
2569 const double RT0QuadFiniteElement::nk[4][2] =
2570 { {0, -1}, {1, 0}, {0, 1}, {-1, 0} };
2576 #ifdef MFEM_THREAD_SAFE
2582 for (k = 0; k < 4; k++)
2585 for (j = 0; j < 4; j++)
2588 if (j == k) { d -= 1.0; }
2589 if (fabs(d) > 1.0e-12)
2591 cerr <<
"RT0QuadFiniteElement::GetLocalInterpolation (...)\n"
2592 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2608 for (k = 0; k < 4; k++)
2611 ip.
x = vk[0]; ip.
y = vk[1];
2614 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2615 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2616 for (j = 0; j < 4; j++)
2617 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2630 #ifdef MFEM_THREAD_SAFE
2634 for (
int k = 0; k < 4; k++)
2642 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2643 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
2671 double x = ip.
x, y = ip.
y;
2673 shape(0,0) = -2 * x * (-1 + x + 2 * y);
2674 shape(0,1) = -2 * (-1 + y) * (-1 + x + 2 * y);
2675 shape(1,0) = 2 * x * (x - y);
2676 shape(1,1) = 2 * (x - y) * (-1 + y);
2677 shape(2,0) = 2 * x * (-1 + 2 * x + y);
2678 shape(2,1) = 2 * y * (-1 + 2 * x + y);
2679 shape(3,0) = 2 * x * (-1 + x + 2 * y);
2680 shape(3,1) = 2 * y * (-1 + x + 2 * y);
2681 shape(4,0) = -2 * (-1 + x) * (x - y);
2682 shape(4,1) = 2 * y * (-x + y);
2683 shape(5,0) = -2 * (-1 + x) * (-1 + 2 * x + y);
2684 shape(5,1) = -2 * y * (-1 + 2 * x + y);
2685 shape(6,0) = -3 * x * (-2 + 2 * x + y);
2686 shape(6,1) = -3 * y * (-1 + 2 * x + y);
2687 shape(7,0) = -3 * x * (-1 + x + 2 * y);
2688 shape(7,1) = -3 * y * (-2 + x + 2 * y);
2694 double x = ip.
x, y = ip.
y;
2696 divshape(0) = -2 * (-4 + 3 * x + 6 * y);
2697 divshape(1) = 2 + 6 * x - 6 * y;
2698 divshape(2) = -4 + 12 * x + 6 * y;
2699 divshape(3) = -4 + 6 * x + 12 * y;
2700 divshape(4) = 2 - 6 * x + 6 * y;
2701 divshape(5) = -2 * (-4 + 6 * x + 3 * y);
2702 divshape(6) = -9 * (-1 + 2 * x + y);
2703 divshape(7) = -9 * (-1 + x + 2 * y);
2706 const double RT1TriangleFiniteElement::nk[8][2] =
2718 #ifdef MFEM_THREAD_SAFE
2724 for (k = 0; k < 8; k++)
2727 for (j = 0; j < 8; j++)
2730 if (j == k) { d -= 1.0; }
2731 if (fabs(d) > 1.0e-12)
2733 cerr <<
"RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
2734 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2750 for (k = 0; k < 8; k++)
2753 ip.
x = vk[0]; ip.
y = vk[1];
2756 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2757 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2758 for (j = 0; j < 8; j++)
2759 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2771 #ifdef MFEM_THREAD_SAFE
2775 for (
int k = 0; k < 8; k++)
2783 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2784 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
2827 double x = ip.
x, y = ip.
y;
2831 shape(0,1) = -( 1. - 3.*y + 2.*y*y)*( 2. - 3.*x);
2833 shape(1,1) = -( 1. - 3.*y + 2.*y*y)*(-1. + 3.*x);
2835 shape(2,0) = (-x + 2.*x*x)*( 2. - 3.*y);
2837 shape(3,0) = (-x + 2.*x*x)*(-1. + 3.*y);
2841 shape(4,1) = (-y + 2.*y*y)*(-1. + 3.*x);
2843 shape(5,1) = (-y + 2.*y*y)*( 2. - 3.*x);
2845 shape(6,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y);
2847 shape(7,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y);
2850 shape(8,0) = (4.*x - 4.*x*x)*( 2. - 3.*y);
2852 shape(9,0) = (4.*x - 4.*x*x)*(-1. + 3.*y);
2856 shape(10,1) = (4.*y - 4.*y*y)*( 2. - 3.*x);
2858 shape(11,1) = (4.*y - 4.*y*y)*(-1. + 3.*x);
2864 double x = ip.
x, y = ip.
y;
2866 divshape(0) = -(-3. + 4.*y)*( 2. - 3.*x);
2867 divshape(1) = -(-3. + 4.*y)*(-1. + 3.*x);
2868 divshape(2) = (-1. + 4.*x)*( 2. - 3.*y);
2869 divshape(3) = (-1. + 4.*x)*(-1. + 3.*y);
2870 divshape(4) = (-1. + 4.*y)*(-1. + 3.*x);
2871 divshape(5) = (-1. + 4.*y)*( 2. - 3.*x);
2872 divshape(6) = -(-3. + 4.*x)*(-1. + 3.*y);
2873 divshape(7) = -(-3. + 4.*x)*( 2. - 3.*y);
2874 divshape(8) = ( 4. - 8.*x)*( 2. - 3.*y);
2875 divshape(9) = ( 4. - 8.*x)*(-1. + 3.*y);
2876 divshape(10) = ( 4. - 8.*y)*( 2. - 3.*x);
2877 divshape(11) = ( 4. - 8.*y)*(-1. + 3.*x);
2880 const double RT1QuadFiniteElement::nk[12][2] =
2900 #ifdef MFEM_THREAD_SAFE
2906 for (k = 0; k < 12; k++)
2909 for (j = 0; j < 12; j++)
2912 if (j == k) { d -= 1.0; }
2913 if (fabs(d) > 1.0e-12)
2915 cerr <<
"RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
2916 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2932 for (k = 0; k < 12; k++)
2935 ip.
x = vk[0]; ip.
y = vk[1];
2938 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2939 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2940 for (j = 0; j < 12; j++)
2941 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2953 #ifdef MFEM_THREAD_SAFE
2957 for (
int k = 0; k < 12; k++)
2965 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2966 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
2970 const double RT2TriangleFiniteElement::M[15][15] =
2973 0, -5.3237900077244501311, 5.3237900077244501311, 16.647580015448900262,
2974 0, 24.442740046346700787, -16.647580015448900262, -12.,
2975 -19.118950038622250656, -47.237900077244501311, 0, -34.414110069520051180,
2976 12., 30.590320061795601049, 15.295160030897800524
2979 0, 1.5, -1.5, -15., 0, 2.625, 15., 15., -4.125, 30., 0, -14.625, -15.,
2983 0, -0.67620999227554986889, 0.67620999227554986889, 7.3524199845510997378,
2984 0, -3.4427400463467007866, -7.3524199845510997378, -12.,
2985 4.1189500386222506555, -0.76209992275549868892, 0, 7.4141100695200511800,
2986 12., -6.5903200617956010489, -3.2951600308978005244
2989 0, 0, 1.5, 0, 0, 1.5, -11.471370023173350393, 0, 2.4713700231733503933,
2990 -11.471370023173350393, 0, 2.4713700231733503933, 15.295160030897800524,
2991 0, -3.2951600308978005244
2994 0, 0, 4.875, 0, 0, 4.875, -16.875, 0, -16.875, -16.875, 0, -16.875, 10.5,
2998 0, 0, 1.5, 0, 0, 1.5, 2.4713700231733503933, 0, -11.471370023173350393,
2999 2.4713700231733503933, 0, -11.471370023173350393, -3.2951600308978005244,
3000 0, 15.295160030897800524
3003 -0.67620999227554986889, 0, -3.4427400463467007866, 0,
3004 7.3524199845510997378, 0.67620999227554986889, 7.4141100695200511800, 0,
3005 -0.76209992275549868892, 4.1189500386222506555, -12.,
3006 -7.3524199845510997378, -3.2951600308978005244, -6.5903200617956010489,
3010 1.5, 0, 2.625, 0, -15., -1.5, -14.625, 0, 30., -4.125, 15., 15., 10.5,
3014 -5.3237900077244501311, 0, 24.442740046346700787, 0, 16.647580015448900262,
3015 5.3237900077244501311, -34.414110069520051180, 0, -47.237900077244501311,
3016 -19.118950038622250656, -12., -16.647580015448900262, 15.295160030897800524,
3017 30.590320061795601049, 12.
3019 { 0, 0, 18., 0, 0, 6., -42., 0, -30., -26., 0, -14., 24., 32., 8.},
3020 { 0, 0, 6., 0, 0, 18., -14., 0, -26., -30., 0, -42., 8., 32., 24.},
3021 { 0, 0, -6., 0, 0, -4., 30., 0, 4., 22., 0, 4., -24., -16., 0},
3022 { 0, 0, -4., 0, 0, -8., 20., 0, 8., 36., 0, 8., -16., -32., 0},
3023 { 0, 0, -8., 0, 0, -4., 8., 0, 36., 8., 0, 20., 0, -32., -16.},
3024 { 0, 0, -4., 0, 0, -6., 4., 0, 22., 4., 0, 30., 0, -16., -24.}
3030 const double p = 0.11270166537925831148;
3067 double x = ip.
x, y = ip.
y;
3069 double Bx[15] = {1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y, 0., x*x*x,
3072 double By[15] = {0., 1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y,
3076 for (
int i = 0; i < 15; i++)
3078 double cx = 0.0, cy = 0.0;
3079 for (
int j = 0; j < 15; j++)
3081 cx += M[i][j] * Bx[j];
3082 cy += M[i][j] * By[j];
3092 double x = ip.
x, y = ip.
y;
3094 double DivB[15] = {0., 0., 1., 0., 0., 1., 2.*x, 0., y, x, 0., 2.*y,
3095 4.*x*x, 4.*x*y, 4.*y*y
3098 for (
int i = 0; i < 15; i++)
3101 for (
int j = 0; j < 15; j++)
3103 div += M[i][j] * DivB[j];
3109 const double RT2QuadFiniteElement::pt[4] = {0.,1./3.,2./3.,1.};
3111 const double RT2QuadFiniteElement::dpt[3] = {0.25,0.5,0.75};
3153 double x = ip.
x, y = ip.
y;
3155 double ax0 = pt[0] - x;
3156 double ax1 = pt[1] - x;
3157 double ax2 = pt[2] - x;
3158 double ax3 = pt[3] - x;
3160 double by0 = dpt[0] - y;
3161 double by1 = dpt[1] - y;
3162 double by2 = dpt[2] - y;
3164 double ay0 = pt[0] - y;
3165 double ay1 = pt[1] - y;
3166 double ay2 = pt[2] - y;
3167 double ay3 = pt[3] - y;
3169 double bx0 = dpt[0] - x;
3170 double bx1 = dpt[1] - x;
3171 double bx2 = dpt[2] - x;
3173 double A01 = pt[0] - pt[1];
3174 double A02 = pt[0] - pt[2];
3175 double A12 = pt[1] - pt[2];
3176 double A03 = pt[0] - pt[3];
3177 double A13 = pt[1] - pt[3];
3178 double A23 = pt[2] - pt[3];
3180 double B01 = dpt[0] - dpt[1];
3181 double B02 = dpt[0] - dpt[2];
3182 double B12 = dpt[1] - dpt[2];
3184 double tx0 = (bx1*bx2)/(B01*B02);
3185 double tx1 = -(bx0*bx2)/(B01*B12);
3186 double tx2 = (bx0*bx1)/(B02*B12);
3188 double ty0 = (by1*by2)/(B01*B02);
3189 double ty1 = -(by0*by2)/(B01*B12);
3190 double ty2 = (by0*by1)/(B02*B12);
3194 shape(0, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx0;
3196 shape(1, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx1;
3198 shape(2, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx2;
3200 shape(3, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty0;
3202 shape(4, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty1;
3204 shape(5, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty2;
3208 shape(6, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx2;
3210 shape(7, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx1;
3212 shape(8, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx0;
3214 shape(9, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty2;
3216 shape(10, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty1;
3218 shape(11, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty0;
3221 shape(12, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty0;
3223 shape(13, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty1;
3225 shape(14, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty2;
3228 shape(15, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty0;
3230 shape(16, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty1;
3232 shape(17, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty2;
3236 shape(18, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx0;
3238 shape(19, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx1;
3240 shape(20, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx2;
3243 shape(21, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx0;
3245 shape(22, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx1;
3247 shape(23, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx2;
3253 double x = ip.
x, y = ip.
y;
3255 double a01 = pt[0]*pt[1];
3256 double a02 = pt[0]*pt[2];
3257 double a12 = pt[1]*pt[2];
3258 double a03 = pt[0]*pt[3];
3259 double a13 = pt[1]*pt[3];
3260 double a23 = pt[2]*pt[3];
3262 double bx0 = dpt[0] - x;
3263 double bx1 = dpt[1] - x;
3264 double bx2 = dpt[2] - x;
3266 double by0 = dpt[0] - y;
3267 double by1 = dpt[1] - y;
3268 double by2 = dpt[2] - y;
3270 double A01 = pt[0] - pt[1];
3271 double A02 = pt[0] - pt[2];
3272 double A12 = pt[1] - pt[2];
3273 double A03 = pt[0] - pt[3];
3274 double A13 = pt[1] - pt[3];
3275 double A23 = pt[2] - pt[3];
3277 double A012 = pt[0] + pt[1] + pt[2];
3278 double A013 = pt[0] + pt[1] + pt[3];
3279 double A023 = pt[0] + pt[2] + pt[3];
3280 double A123 = pt[1] + pt[2] + pt[3];
3282 double B01 = dpt[0] - dpt[1];
3283 double B02 = dpt[0] - dpt[2];
3284 double B12 = dpt[1] - dpt[2];
3286 double tx0 = (bx1*bx2)/(B01*B02);
3287 double tx1 = -(bx0*bx2)/(B01*B12);
3288 double tx2 = (bx0*bx1)/(B02*B12);
3290 double ty0 = (by1*by2)/(B01*B02);
3291 double ty1 = -(by0*by2)/(B01*B12);
3292 double ty2 = (by0*by1)/(B02*B12);
3295 divshape(0) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx0;
3296 divshape(1) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx1;
3297 divshape(2) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx2;
3299 divshape(3) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty0;
3300 divshape(4) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty1;
3301 divshape(5) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty2;
3303 divshape(6) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx2;
3304 divshape(7) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx1;
3305 divshape(8) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx0;
3307 divshape(9) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty2;
3308 divshape(10) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty1;
3309 divshape(11) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty0;
3311 divshape(12) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty0;
3312 divshape(13) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty1;
3313 divshape(14) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty2;
3315 divshape(15) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty0;
3316 divshape(16) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty1;
3317 divshape(17) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty2;
3319 divshape(18) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx0;
3320 divshape(19) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx1;
3321 divshape(20) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx2;
3323 divshape(21) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx0;
3324 divshape(22) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx1;
3325 divshape(23) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx2;
3328 const double RT2QuadFiniteElement::nk[24][2] =
3331 {0,-1}, {0,-1}, {0,-1},
3333 {1, 0}, {1, 0}, {1, 0},
3335 {0, 1}, {0, 1}, {0, 1},
3337 {-1,0}, {-1,0}, {-1,0},
3339 {1, 0}, {1, 0}, {1, 0},
3341 {1, 0}, {1, 0}, {1, 0},
3343 {0, 1}, {0, 1}, {0, 1},
3345 {0, 1}, {0, 1}, {0, 1}
3352 #ifdef MFEM_THREAD_SAFE
3358 for (k = 0; k < 24; k++)
3361 for (j = 0; j < 24; j++)
3364 if (j == k) { d -= 1.0; }
3365 if (fabs(d) > 1.0e-12)
3367 cerr <<
"RT2QuadFiniteElement::GetLocalInterpolation (...)\n"
3368 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
3384 for (k = 0; k < 24; k++)
3387 ip.
x = vk[0]; ip.
y = vk[1];
3390 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
3391 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
3392 for (j = 0; j < 24; j++)
3393 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
3405 #ifdef MFEM_THREAD_SAFE
3409 for (
int k = 0; k < 24; k++)
3417 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
3418 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
3434 shape(0) = 2. - 3. * x;
3435 shape(1) = 3. * x - 1.;
3449 const double p = 0.11270166537925831148;
3459 const double p = 0.11270166537925831148;
3460 const double w = 1./((1-2*p)*(1-2*p));
3463 shape(0) = (2*x-1)*(x-1+p)*w;
3464 shape(1) = 4*(x-1+p)*(p-x)*w;
3465 shape(2) = (2*x-1)*(x-p)*w;
3471 const double p = 0.11270166537925831148;
3472 const double w = 1./((1-2*p)*(1-2*p));
3475 dshape(0,0) = (-3+4*x+2*p)*w;
3476 dshape(1,0) = (4-8*x)*w;
3477 dshape(2,0) = (-1+4*x-2*p)*w;
3488 for (i = 1; i < m; i++)
3494 #ifndef MFEM_THREAD_SAFE
3499 for (i = 1; i <= m; i++)
3501 rwk(i) = rwk(i-1) * ( (double)(m) / (double)(i) );
3503 for (i = 0; i < m/2+1; i++)
3505 rwk(m-i) = ( rwk(i) *= rwk(m-i) );
3507 for (i = m-1; i >= 0; i -= 2)
3516 double w, wk, x = ip.
x;
3519 #ifdef MFEM_THREAD_SAFE
3523 k = (int) floor ( m * x + 0.5 );
3525 for (i = 0; i <= m; i++)
3528 wk *= ( rxxk(i) = x - (double)(i) / m );
3530 w = wk * ( rxxk(k) = x - (double)(k) / m );
3534 shape(0) = w * rwk(0) / rxxk(0);
3538 shape(0) = wk * rwk(0);
3542 shape(1) = w * rwk(m) / rxxk(m);
3546 shape(1) = wk * rwk(k);
3548 for (i = 1; i < m; i++)
3551 shape(i+1) = w * rwk(i) / rxxk(i);
3555 shape(k+1) = wk * rwk(k);
3562 double s, srx, w, wk, x = ip.
x;
3565 #ifdef MFEM_THREAD_SAFE
3569 k = (int) floor ( m * x + 0.5 );
3571 for (i = 0; i <= m; i++)
3574 wk *= ( rxxk(i) = x - (double)(i) / m );
3576 w = wk * ( rxxk(k) = x - (double)(k) / m );
3578 for (i = 0; i <= m; i++)
3580 rxxk(i) = 1.0 / rxxk(i);
3583 for (i = 0; i <= m; i++)
3592 dshape(0,0) = (s - w * rxxk(0)) * rwk(0) * rxxk(0);
3596 dshape(0,0) = wk * srx * rwk(0);
3600 dshape(1,0) = (s - w * rxxk(m)) * rwk(m) * rxxk(m);
3604 dshape(1,0) = wk * srx * rwk(k);
3606 for (i = 1; i < m; i++)
3609 dshape(i+1,0) = (s - w * rxxk(i)) * rwk(i) * rxxk(i);
3613 dshape(k+1,0) = wk * srx * rwk(k);
3642 double L0, L1, L2, L3;
3644 L1 = ip.
x; L2 = ip.
y; L3 = ip.
z; L0 = 1.0 - L1 - L2 - L3;
3645 shape(0) = 1.0 - 3.0 * L0;
3646 shape(1) = 1.0 - 3.0 * L1;
3647 shape(2) = 1.0 - 3.0 * L2;
3648 shape(3) = 1.0 - 3.0 * L3;
3654 dshape(0,0) = 3.0; dshape(0,1) = 3.0; dshape(0,2) = 3.0;
3655 dshape(1,0) = -3.0; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
3656 dshape(2,0) = 0.0; dshape(2,1) = -3.0; dshape(2,2) = 0.0;
3657 dshape(3,0) = 0.0; dshape(3,1) = 0.0; dshape(3,2) = -3.0;
3678 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3699 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3713 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3714 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3715 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3716 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3717 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3718 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3719 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3720 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3722 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3723 I[ 9] = 1; J[ 9] = 2; K[ 9] = 0;
3724 I[10] = 2; J[10] = 1; K[10] = 0;
3725 I[11] = 0; J[11] = 2; K[11] = 0;
3726 I[12] = 2; J[12] = 0; K[12] = 1;
3727 I[13] = 1; J[13] = 2; K[13] = 1;
3728 I[14] = 2; J[14] = 1; K[14] = 1;
3729 I[15] = 0; J[15] = 2; K[15] = 1;
3730 I[16] = 0; J[16] = 0; K[16] = 2;
3731 I[17] = 1; J[17] = 0; K[17] = 2;
3732 I[18] = 1; J[18] = 1; K[18] = 2;
3733 I[19] = 0; J[19] = 1; K[19] = 2;
3735 I[20] = 2; J[20] = 2; K[20] = 0;
3736 I[21] = 2; J[21] = 0; K[21] = 2;
3737 I[22] = 1; J[22] = 2; K[22] = 2;
3738 I[23] = 2; J[23] = 1; K[23] = 2;
3739 I[24] = 0; J[24] = 2; K[24] = 2;
3740 I[25] = 2; J[25] = 2; K[25] = 1;
3742 I[26] = 2; J[26] = 2; K[26] = 2;
3744 else if (degree == 3)
3750 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3751 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3752 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3753 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3754 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3755 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3756 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3757 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3759 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3760 I[ 9] = 3; J[ 9] = 0; K[ 9] = 0;
3761 I[10] = 1; J[10] = 2; K[10] = 0;
3762 I[11] = 1; J[11] = 3; K[11] = 0;
3763 I[12] = 2; J[12] = 1; K[12] = 0;
3764 I[13] = 3; J[13] = 1; K[13] = 0;
3765 I[14] = 0; J[14] = 2; K[14] = 0;
3766 I[15] = 0; J[15] = 3; K[15] = 0;
3767 I[16] = 2; J[16] = 0; K[16] = 1;
3768 I[17] = 3; J[17] = 0; K[17] = 1;
3769 I[18] = 1; J[18] = 2; K[18] = 1;
3770 I[19] = 1; J[19] = 3; K[19] = 1;
3771 I[20] = 2; J[20] = 1; K[20] = 1;
3772 I[21] = 3; J[21] = 1; K[21] = 1;
3773 I[22] = 0; J[22] = 2; K[22] = 1;
3774 I[23] = 0; J[23] = 3; K[23] = 1;
3775 I[24] = 0; J[24] = 0; K[24] = 2;
3776 I[25] = 0; J[25] = 0; K[25] = 3;
3777 I[26] = 1; J[26] = 0; K[26] = 2;
3778 I[27] = 1; J[27] = 0; K[27] = 3;
3779 I[28] = 1; J[28] = 1; K[28] = 2;
3780 I[29] = 1; J[29] = 1; K[29] = 3;
3781 I[30] = 0; J[30] = 1; K[30] = 2;
3782 I[31] = 0; J[31] = 1; K[31] = 3;
3784 I[32] = 2; J[32] = 3; K[32] = 0;
3785 I[33] = 3; J[33] = 3; K[33] = 0;
3786 I[34] = 2; J[34] = 2; K[34] = 0;
3787 I[35] = 3; J[35] = 2; K[35] = 0;
3788 I[36] = 2; J[36] = 0; K[36] = 2;
3789 I[37] = 3; J[37] = 0; K[37] = 2;
3790 I[38] = 2; J[38] = 0; K[38] = 3;
3791 I[39] = 3; J[39] = 0; K[39] = 3;
3792 I[40] = 1; J[40] = 2; K[40] = 2;
3793 I[41] = 1; J[41] = 3; K[41] = 2;
3794 I[42] = 1; J[42] = 2; K[42] = 3;
3795 I[43] = 1; J[43] = 3; K[43] = 3;
3796 I[44] = 3; J[44] = 1; K[44] = 2;
3797 I[45] = 2; J[45] = 1; K[45] = 2;
3798 I[46] = 3; J[46] = 1; K[46] = 3;
3799 I[47] = 2; J[47] = 1; K[47] = 3;
3800 I[48] = 0; J[48] = 3; K[48] = 2;
3801 I[49] = 0; J[49] = 2; K[49] = 2;
3802 I[50] = 0; J[50] = 3; K[50] = 3;
3803 I[51] = 0; J[51] = 2; K[51] = 3;
3804 I[52] = 2; J[52] = 2; K[52] = 1;
3805 I[53] = 3; J[53] = 2; K[53] = 1;
3806 I[54] = 2; J[54] = 3; K[54] = 1;
3807 I[55] = 3; J[55] = 3; K[55] = 1;
3809 I[56] = 2; J[56] = 2; K[56] = 2;
3810 I[57] = 3; J[57] = 2; K[57] = 2;
3811 I[58] = 3; J[58] = 3; K[58] = 2;
3812 I[59] = 2; J[59] = 3; K[59] = 2;
3813 I[60] = 2; J[60] = 2; K[60] = 3;
3814 I[61] = 3; J[61] = 2; K[61] = 3;
3815 I[62] = 3; J[62] = 3; K[62] = 3;
3816 I[63] = 2; J[63] = 3; K[63] = 3;
3820 mfem_error (
"LagrangeHexFiniteElement::LagrangeHexFiniteElement");
3824 dof1d = fe1d ->
GetDof();
3826 #ifndef MFEM_THREAD_SAFE
3836 for (
int n = 0; n <
Dof; n++)
3851 #ifdef MFEM_THREAD_SAFE
3852 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
3859 for (
int n = 0; n <
Dof; n++)
3861 shape(n) = shape1dx(I[n]) * shape1dy(J[n]) * shape1dz(K[n]);
3872 #ifdef MFEM_THREAD_SAFE
3873 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
3874 DenseMatrix dshape1dx(dof1d,1), dshape1dy(dof1d,1), dshape1dz(dof1d,1);
3885 for (
int n = 0; n <
Dof; n++)
3887 dshape(n,0) = dshape1dx(I[n],0) * shape1dy(J[n]) * shape1dz(K[n]);
3888 dshape(n,1) = shape1dx(I[n]) * dshape1dy(J[n],0) * shape1dz(K[n]);
3889 dshape(n,2) = shape1dx(I[n]) * shape1dy(J[n]) * dshape1dz(K[n],0);
3918 shape(0) = 1.0 - 2.0 * x;
3925 shape(1) = 2.0 * x - 1.0;
3926 shape(2) = 2.0 - 2.0 * x;
3937 dshape(0,0) = - 2.0;
3945 dshape(2,0) = - 2.0;
3972 L0 = 2.0 * ( 1. - ip.
x - ip.
y );
3973 L1 = 2.0 * ( ip.
x );
3974 L2 = 2.0 * ( ip.
y );
3983 for (i = 0; i < 6; i++)
3990 shape(0) = L0 - 1.0;
3997 shape(1) = L1 - 1.0;
4004 shape(2) = L2 - 1.0;
4008 shape(3) = 1.0 - L2;
4009 shape(4) = 1.0 - L0;
4010 shape(5) = 1.0 - L1;
4020 L0 = 2.0 * ( 1. - ip.
x - ip.
y );
4021 L1 = 2.0 * ( ip.
x );
4022 L2 = 2.0 * ( ip.
y );
4024 double DL0[2], DL1[2], DL2[2];
4025 DL0[0] = -2.0; DL0[1] = -2.0;
4026 DL1[0] = 2.0; DL1[1] = 0.0;
4027 DL2[0] = 0.0; DL2[1] = 2.0;
4029 for (i = 0; i < 6; i++)
4030 for (j = 0; j < 2; j++)
4037 for (j = 0; j < 2; j++)
4039 dshape(0,j) = DL0[j];
4040 dshape(3,j) = DL1[j];
4041 dshape(5,j) = DL2[j];
4046 for (j = 0; j < 2; j++)
4048 dshape(3,j) = DL0[j];
4049 dshape(1,j) = DL1[j];
4050 dshape(4,j) = DL2[j];
4055 for (j = 0; j < 2; j++)
4057 dshape(5,j) = DL0[j];
4058 dshape(4,j) = DL1[j];
4059 dshape(2,j) = DL2[j];
4064 for (j = 0; j < 2; j++)
4066 dshape(3,j) = - DL2[j];
4067 dshape(4,j) = - DL0[j];
4068 dshape(5,j) = - DL1[j];
4113 double L0, L1, L2, L3, L4, L5;
4114 L0 = 2.0 * ( 1. - ip.
x - ip.
y - ip.
z );
4115 L1 = 2.0 * ( ip.
x );
4116 L2 = 2.0 * ( ip.
y );
4117 L3 = 2.0 * ( ip.
z );
4118 L4 = 2.0 * ( ip.
x + ip.
y );
4119 L5 = 2.0 * ( ip.
y + ip.
z );
4132 for (i = 0; i < 10; i++)
4139 shape(0) = L0 - 1.0;
4147 shape(1) = L1 - 1.0;
4155 shape(2) = L2 - 1.0;
4163 shape(3) = L3 - 1.0;
4165 else if ((L4 <= 1.0) && (L5 <= 1.0))
4167 shape(4) = 1.0 - L5;
4169 shape(6) = 1.0 - L4;
4170 shape(8) = 1.0 - L0;
4172 else if ((L4 >= 1.0) && (L5 <= 1.0))
4174 shape(4) = 1.0 - L5;
4175 shape(5) = 1.0 - L1;
4176 shape(7) = L4 - 1.0;
4179 else if ((L4 <= 1.0) && (L5 >= 1.0))
4181 shape(5) = 1.0 - L3;
4182 shape(6) = 1.0 - L4;
4184 shape(9) = L5 - 1.0;
4186 else if ((L4 >= 1.0) && (L5 >= 1.0))
4189 shape(7) = L4 - 1.0;
4190 shape(8) = 1.0 - L2;
4191 shape(9) = L5 - 1.0;
4200 double L0, L1, L2, L3, L4, L5;
4201 L0 = 2.0 * ( 1. - ip.
x - ip.
y - ip.
z );
4202 L1 = 2.0 * ( ip.
x );
4203 L2 = 2.0 * ( ip.
y );
4204 L3 = 2.0 * ( ip.
z );
4205 L4 = 2.0 * ( ip.
x + ip.
y );
4206 L5 = 2.0 * ( ip.
y + ip.
z );
4208 double DL0[3], DL1[3], DL2[3], DL3[3], DL4[3], DL5[3];
4209 DL0[0] = -2.0; DL0[1] = -2.0; DL0[2] = -2.0;
4210 DL1[0] = 2.0; DL1[1] = 0.0; DL1[2] = 0.0;
4211 DL2[0] = 0.0; DL2[1] = 2.0; DL2[2] = 0.0;
4212 DL3[0] = 0.0; DL3[1] = 0.0; DL3[2] = 2.0;
4213 DL4[0] = 2.0; DL4[1] = 2.0; DL4[2] = 0.0;
4214 DL5[0] = 0.0; DL5[1] = 2.0; DL5[2] = 2.0;
4216 for (i = 0; i < 10; i++)
4217 for (j = 0; j < 3; j++)
4224 for (j = 0; j < 3; j++)
4226 dshape(0,j) = DL0[j];
4227 dshape(4,j) = DL1[j];
4228 dshape(5,j) = DL2[j];
4229 dshape(6,j) = DL3[j];
4234 for (j = 0; j < 3; j++)
4236 dshape(4,j) = DL0[j];
4237 dshape(1,j) = DL1[j];
4238 dshape(7,j) = DL2[j];
4239 dshape(8,j) = DL3[j];
4244 for (j = 0; j < 3; j++)
4246 dshape(5,j) = DL0[j];
4247 dshape(7,j) = DL1[j];
4248 dshape(2,j) = DL2[j];
4249 dshape(9,j) = DL3[j];
4254 for (j = 0; j < 3; j++)
4256 dshape(6,j) = DL0[j];
4257 dshape(8,j) = DL1[j];
4258 dshape(9,j) = DL2[j];
4259 dshape(3,j) = DL3[j];
4262 else if ((L4 <= 1.0) && (L5 <= 1.0))
4264 for (j = 0; j < 3; j++)
4266 dshape(4,j) = - DL5[j];
4267 dshape(5,j) = DL2[j];
4268 dshape(6,j) = - DL4[j];
4269 dshape(8,j) = - DL0[j];
4272 else if ((L4 >= 1.0) && (L5 <= 1.0))
4274 for (j = 0; j < 3; j++)
4276 dshape(4,j) = - DL5[j];
4277 dshape(5,j) = - DL1[j];
4278 dshape(7,j) = DL4[j];
4279 dshape(8,j) = DL3[j];
4282 else if ((L4 <= 1.0) && (L5 >= 1.0))
4284 for (j = 0; j < 3; j++)
4286 dshape(5,j) = - DL3[j];
4287 dshape(6,j) = - DL4[j];
4288 dshape(8,j) = DL1[j];
4289 dshape(9,j) = DL5[j];
4292 else if ((L4 >= 1.0) && (L5 >= 1.0))
4294 for (j = 0; j < 3; j++)
4296 dshape(5,j) = DL0[j];
4297 dshape(7,j) = DL4[j];
4298 dshape(8,j) = - DL2[j];
4299 dshape(9,j) = DL5[j];
4332 double x = ip.
x, y = ip.
y;
4334 Lx = 2.0 * ( 1. - x );
4335 Ly = 2.0 * ( 1. - y );
4344 for (i = 0; i < 9; i++)
4349 if ((x <= 0.5) && (y <= 0.5))
4351 shape(0) = (Lx - 1.0) * (Ly - 1.0);
4352 shape(4) = (2.0 - Lx) * (Ly - 1.0);
4353 shape(8) = (2.0 - Lx) * (2.0 - Ly);
4354 shape(7) = (Lx - 1.0) * (2.0 - Ly);
4356 else if ((x >= 0.5) && (y <= 0.5))
4358 shape(4) = Lx * (Ly - 1.0);
4359 shape(1) = (1.0 - Lx) * (Ly - 1.0);
4360 shape(5) = (1.0 - Lx) * (2.0 - Ly);
4361 shape(8) = Lx * (2.0 - Ly);
4363 else if ((x >= 0.5) && (y >= 0.5))
4365 shape(8) = Lx * Ly ;
4366 shape(5) = (1.0 - Lx) * Ly ;
4367 shape(2) = (1.0 - Lx) * (1.0 - Ly);
4368 shape(6) = Lx * (1.0 - Ly);
4370 else if ((x <= 0.5) && (y >= 0.5))
4372 shape(7) = (Lx - 1.0) * Ly ;
4373 shape(8) = (2.0 - Lx) * Ly ;
4374 shape(6) = (2.0 - Lx) * (1.0 - Ly);
4375 shape(3) = (Lx - 1.0) * (1.0 - Ly);
4383 double x = ip.
x, y = ip.
y;
4385 Lx = 2.0 * ( 1. - x );
4386 Ly = 2.0 * ( 1. - y );
4388 for (i = 0; i < 9; i++)
4389 for (j = 0; j < 2; j++)
4394 if ((x <= 0.5) && (y <= 0.5))
4396 dshape(0,0) = 2.0 * (1.0 - Ly);
4397 dshape(0,1) = 2.0 * (1.0 - Lx);
4399 dshape(4,0) = 2.0 * (Ly - 1.0);
4400 dshape(4,1) = -2.0 * (2.0 - Lx);
4402 dshape(8,0) = 2.0 * (2.0 - Ly);
4403 dshape(8,1) = 2.0 * (2.0 - Lx);
4405 dshape(7,0) = -2.0 * (2.0 - Ly);
4406 dshape(7,0) = 2.0 * (Lx - 1.0);
4408 else if ((x >= 0.5) && (y <= 0.5))
4410 dshape(4,0) = -2.0 * (Ly - 1.0);
4411 dshape(4,1) = -2.0 * Lx;
4413 dshape(1,0) = 2.0 * (Ly - 1.0);
4414 dshape(1,1) = -2.0 * (1.0 - Lx);
4416 dshape(5,0) = 2.0 * (2.0 - Ly);
4417 dshape(5,1) = 2.0 * (1.0 - Lx);
4419 dshape(8,0) = -2.0 * (2.0 - Ly);
4420 dshape(8,1) = 2.0 * Lx;
4422 else if ((x >= 0.5) && (y >= 0.5))
4424 dshape(8,0) = -2.0 * Ly;
4425 dshape(8,1) = -2.0 * Lx;
4427 dshape(5,0) = 2.0 * Ly;
4428 dshape(5,1) = -2.0 * (1.0 - Lx);
4430 dshape(2,0) = 2.0 * (1.0 - Ly);
4431 dshape(2,1) = 2.0 * (1.0 - Lx);
4433 dshape(6,0) = -2.0 * (1.0 - Ly);
4434 dshape(6,1) = 2.0 * Lx;
4436 else if ((x <= 0.5) && (y >= 0.5))
4438 dshape(7,0) = -2.0 * Ly;
4439 dshape(7,1) = -2.0 * (Lx - 1.0);
4441 dshape(8,0) = 2.0 * Ly ;
4442 dshape(8,1) = -2.0 * (2.0 - Lx);
4444 dshape(6,0) = 2.0 * (1.0 - Ly);
4445 dshape(6,1) = 2.0 * (2.0 - Lx);
4447 dshape(3,0) = -2.0 * (1.0 - Ly);
4448 dshape(3,1) = 2.0 * (Lx - 1.0);
4459 I[ 0] = 0.0; J[ 0] = 0.0; K[ 0] = 0.0;
4460 I[ 1] = 1.0; J[ 1] = 0.0; K[ 1] = 0.0;
4461 I[ 2] = 1.0; J[ 2] = 1.0; K[ 2] = 0.0;
4462 I[ 3] = 0.0; J[ 3] = 1.0; K[ 3] = 0.0;
4463 I[ 4] = 0.0; J[ 4] = 0.0; K[ 4] = 1.0;
4464 I[ 5] = 1.0; J[ 5] = 0.0; K[ 5] = 1.0;
4465 I[ 6] = 1.0; J[ 6] = 1.0; K[ 6] = 1.0;
4466 I[ 7] = 0.0; J[ 7] = 1.0; K[ 7] = 1.0;
4468 I[ 8] = 0.5; J[ 8] = 0.0; K[ 8] = 0.0;
4469 I[ 9] = 1.0; J[ 9] = 0.5; K[ 9] = 0.0;
4470 I[10] = 0.5; J[10] = 1.0; K[10] = 0.0;
4471 I[11] = 0.0; J[11] = 0.5; K[11] = 0.0;
4472 I[12] = 0.5; J[12] = 0.0; K[12] = 1.0;
4473 I[13] = 1.0; J[13] = 0.5; K[13] = 1.0;
4474 I[14] = 0.5; J[14] = 1.0; K[14] = 1.0;
4475 I[15] = 0.0; J[15] = 0.5; K[15] = 1.0;
4476 I[16] = 0.0; J[16] = 0.0; K[16] = 0.5;
4477 I[17] = 1.0; J[17] = 0.0; K[17] = 0.5;
4478 I[18] = 1.0; J[18] = 1.0; K[18] = 0.5;
4479 I[19] = 0.0; J[19] = 1.0; K[19] = 0.5;
4481 I[20] = 0.5; J[20] = 0.5; K[20] = 0.0;
4482 I[21] = 0.5; J[21] = 0.0; K[21] = 0.5;
4483 I[22] = 1.0; J[22] = 0.5; K[22] = 0.5;
4484 I[23] = 0.5; J[23] = 1.0; K[23] = 0.5;
4485 I[24] = 0.0; J[24] = 0.5; K[24] = 0.5;
4486 I[25] = 0.5; J[25] = 0.5; K[25] = 1.0;
4488 I[26] = 0.5; J[26] = 0.5; K[26] = 0.5;
4490 for (
int n = 0; n < 27; n++)
4503 double x = ip.
x, y = ip.
y, z = ip.
z;
4505 for (i = 0; i < 27; i++)
4510 if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5))
4525 else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5))
4540 else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5))
4555 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5))
4570 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5))
4585 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5))
4600 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5))
4631 shape(N[0]) = Lx * Ly * Lz;
4632 shape(N[1]) = (1 - Lx) * Ly * Lz;
4633 shape(N[2]) = (1 - Lx) * (1 - Ly) * Lz;
4634 shape(N[3]) = Lx * (1 - Ly) * Lz;
4635 shape(N[4]) = Lx * Ly * (1 - Lz);
4636 shape(N[5]) = (1 - Lx) * Ly * (1 - Lz);
4637 shape(N[6]) = (1 - Lx) * (1 - Ly) * (1 - Lz);
4638 shape(N[7]) = Lx * (1 - Ly) * (1 - Lz);
4646 double x = ip.
x, y = ip.
y, z = ip.
z;
4648 for (i = 0; i < 27; i++)
4649 for (j = 0; j < 3; j++)
4654 if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5))
4669 else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5))
4684 else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5))
4699 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5))
4714 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5))
4729 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5))
4744 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5))
4775 dshape(N[0],0) = -2.0 * Ly * Lz ;
4776 dshape(N[0],1) = -2.0 * Lx * Lz ;
4777 dshape(N[0],2) = -2.0 * Lx * Ly ;
4779 dshape(N[1],0) = 2.0 * Ly * Lz ;
4780 dshape(N[1],1) = -2.0 * (1 - Lx) * Lz ;
4781 dshape(N[1],2) = -2.0 * (1 - Lx) * Ly ;
4783 dshape(N[2],0) = 2.0 * (1 - Ly) * Lz ;
4784 dshape(N[2],1) = 2.0 * (1 - Lx) * Lz ;
4785 dshape(N[2],2) = -2.0 * (1 - Lx) * (1 - Ly);
4787 dshape(N[3],0) = -2.0 * (1 - Ly) * Lz ;
4788 dshape(N[3],1) = 2.0 * Lx * Lz ;
4789 dshape(N[3],2) = -2.0 * Lx * (1 - Ly);
4791 dshape(N[4],0) = -2.0 * Ly * (1 - Lz);
4792 dshape(N[4],1) = -2.0 * Lx * (1 - Lz);
4793 dshape(N[4],2) = 2.0 * Lx * Ly ;
4795 dshape(N[5],0) = 2.0 * Ly * (1 - Lz);
4796 dshape(N[5],1) = -2.0 * (1 - Lx) * (1 - Lz);
4797 dshape(N[5],2) = 2.0 * (1 - Lx) * Ly ;
4799 dshape(N[6],0) = 2.0 * (1 - Ly) * (1 - Lz);
4800 dshape(N[6],1) = 2.0 * (1 - Lx) * (1 - Lz);
4801 dshape(N[6],2) = 2.0 * (1 - Lx) * (1 - Ly);
4803 dshape(N[7],0) = -2.0 * (1 - Ly) * (1 - Lz);
4804 dshape(N[7],1) = 2.0 * Lx * (1 - Lz);
4805 dshape(N[7],2) = 2.0 * Lx * (1 - Ly);
4865 double x = ip.
x, y = ip.
y, z = ip.
z;
4867 shape(0,0) = (1. - y) * (1. - z);
4871 shape(2,0) = y * (1. - z);
4875 shape(4,0) = z * (1. - y);
4884 shape(1,1) = x * (1. - z);
4888 shape(3,1) = (1. - x) * (1. - z);
4896 shape(7,1) = (1. - x) * z;
4901 shape(8,2) = (1. - x) * (1. - y);
4905 shape(9,2) = x * (1. - y);
4909 shape(10,2) = x * y;
4913 shape(11,2) = y * (1. - x);
4921 double x = ip.
x, y = ip.
y, z = ip.
z;
4923 curl_shape(0,0) = 0.;
4924 curl_shape(0,1) = y - 1.;
4925 curl_shape(0,2) = 1. - z;
4927 curl_shape(2,0) = 0.;
4928 curl_shape(2,1) = -y;
4929 curl_shape(2,2) = z - 1.;
4931 curl_shape(4,0) = 0;
4932 curl_shape(4,1) = 1. - y;
4933 curl_shape(4,2) = z;
4935 curl_shape(6,0) = 0.;
4936 curl_shape(6,1) = y;
4937 curl_shape(6,2) = -z;
4939 curl_shape(1,0) = x;
4940 curl_shape(1,1) = 0.;
4941 curl_shape(1,2) = 1. - z;
4943 curl_shape(3,0) = 1. - x;
4944 curl_shape(3,1) = 0.;
4945 curl_shape(3,2) = z - 1.;
4947 curl_shape(5,0) = -x;
4948 curl_shape(5,1) = 0.;
4949 curl_shape(5,2) = z;
4951 curl_shape(7,0) = x - 1.;
4952 curl_shape(7,1) = 0.;
4953 curl_shape(7,2) = -z;
4955 curl_shape(8,0) = x - 1.;
4956 curl_shape(8,1) = 1. - y;
4957 curl_shape(8,2) = 0.;
4959 curl_shape(9,0) = -x;
4960 curl_shape(9,1) = y - 1.;
4961 curl_shape(9,2) = 0;
4963 curl_shape(10,0) = x;
4964 curl_shape(10,1) = -y;
4965 curl_shape(10,2) = 0.;
4967 curl_shape(11,0) = 1. - x;
4968 curl_shape(11,1) = y;
4969 curl_shape(11,2) = 0.;
4972 const double Nedelec1HexFiniteElement::tk[12][3] =
4974 {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
4975 {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
4976 {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}
4983 #ifdef MFEM_THREAD_SAFE
4988 for (k = 0; k < 12; k++)
4991 for (j = 0; j < 12; j++)
4993 double d = (
vshape(j,0)*tk[k][0] +
vshape(j,1)*tk[k][1] +
4995 if (j == k) { d -= 1.0; }
4996 if (fabs(d) > 1.0e-12)
4998 cerr <<
"Nedelec1HexFiniteElement::GetLocalInterpolation (...)\n"
4999 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5007 ip.
x = ip.
y = ip.
z = 0.0;
5014 for (k = 0; k < 12; k++)
5017 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5020 vk[0] =
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2];
5021 vk[1] =
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2];
5022 vk[2] =
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2];
5023 for (j = 0; j < 12; j++)
5024 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5025 vshape(j,2)*vk[2])) < 1.0e-12)
5039 for (
int k = 0; k < 12; k++)
5047 vk[0] * (
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2] ) +
5048 vk[1] * (
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2] ) +
5049 vk[2] * (
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2] );
5086 double x = ip.
x, y = ip.
y, z = ip.
z;
5088 shape(0,0) = 1. - y - z;
5093 shape(1,1) = 1. - x - z;
5098 shape(2,2) = 1. - x - y;
5117 curl_shape(0,0) = 0.;
5118 curl_shape(0,1) = -2.;
5119 curl_shape(0,2) = 2.;
5121 curl_shape(1,0) = 2.;
5122 curl_shape(1,1) = 0.;
5123 curl_shape(1,2) = -2.;
5125 curl_shape(2,0) = -2.;
5126 curl_shape(2,1) = 2.;
5127 curl_shape(2,2) = 0.;
5129 curl_shape(3,0) = 0.;
5130 curl_shape(3,1) = 0.;
5131 curl_shape(3,2) = 2.;
5133 curl_shape(4,0) = 0.;
5134 curl_shape(4,1) = -2.;
5135 curl_shape(4,2) = 0.;
5137 curl_shape(5,0) = 2.;
5138 curl_shape(5,1) = 0.;
5139 curl_shape(5,2) = 0.;
5142 const double Nedelec1TetFiniteElement::tk[6][3] =
5143 {{1,0,0}, {0,1,0}, {0,0,1}, {-1,1,0}, {-1,0,1}, {0,-1,1}};
5149 #ifdef MFEM_THREAD_SAFE
5154 for (k = 0; k < 6; k++)
5157 for (j = 0; j < 6; j++)
5159 double d = (
vshape(j,0)*tk[k][0] +
vshape(j,1)*tk[k][1] +
5161 if (j == k) { d -= 1.0; }
5162 if (fabs(d) > 1.0e-12)
5164 cerr <<
"Nedelec1TetFiniteElement::GetLocalInterpolation (...)\n"
5165 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5173 ip.
x = ip.
y = ip.
z = 0.0;
5180 for (k = 0; k < 6; k++)
5183 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5186 vk[0] =
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2];
5187 vk[1] =
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2];
5188 vk[2] =
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2];
5189 for (j = 0; j < 6; j++)
5190 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5191 vshape(j,2)*vk[2])) < 1.0e-12)
5205 for (
int k = 0; k < 6; k++)
5213 vk[0] * (
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2] ) +
5214 vk[1] * (
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2] ) +
5215 vk[2] * (
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2] );
5252 double x = ip.
x, y = ip.
y, z = ip.
z;
5256 shape(0,2) = z - 1.;
5259 shape(1,1) = y - 1.;
5270 shape(4,0) = x - 1.;
5290 const double RT0HexFiniteElement::nk[6][3] =
5291 {{0,0,-1}, {0,-1,0}, {1,0,0}, {0,1,0}, {-1,0,0}, {0,0,1}};
5297 #ifdef MFEM_THREAD_SAFE
5303 for (k = 0; k < 6; k++)
5306 for (j = 0; j < 6; j++)
5308 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5310 if (j == k) { d -= 1.0; }
5311 if (fabs(d) > 1.0e-12)
5313 cerr <<
"RT0HexFiniteElement::GetLocalInterpolation (...)\n"
5314 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5322 ip.
x = ip.
y = ip.
z = 0.0;
5330 for (k = 0; k < 6; k++)
5333 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5336 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2];
5337 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2];
5338 vk[2] =
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2];
5339 for (j = 0; j < 6; j++)
5340 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5341 vshape(j,2)*vk[2])) < 1.0e-12)
5354 #ifdef MFEM_THREAD_SAFE
5358 for (
int k = 0; k < 6; k++)
5367 vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2] ) +
5368 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2] ) +
5369 vk[2] * (
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2] );
5498 double x = ip.
x, y = ip.
y, z = ip.
z;
5502 shape(2,2) = -(1. - 3.*z + 2.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5505 shape(3,2) = -(1. - 3.*z + 2.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5508 shape(0,2) = -(1. - 3.*z + 2.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5511 shape(1,2) = -(1. - 3.*z + 2.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5514 shape(4,1) = -(1. - 3.*y + 2.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5517 shape(5,1) = -(1. - 3.*y + 2.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5520 shape(6,1) = -(1. - 3.*y + 2.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5523 shape(7,1) = -(1. - 3.*y + 2.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5526 shape(8,0) = (-x + 2.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5529 shape(9,0) = (-x + 2.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5532 shape(10,0) = (-x + 2.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5535 shape(11,0) = (-x + 2.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5540 shape(13,1) = (-y + 2.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5543 shape(12,1) = (-y + 2.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5546 shape(15,1) = (-y + 2.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5549 shape(14,1) = (-y + 2.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5552 shape(17,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5555 shape(16,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5558 shape(19,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5561 shape(18,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5567 shape(20,2) = (-z + 2.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5570 shape(21,2) = (-z + 2.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5573 shape(22,2) = (-z + 2.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5576 shape(23,2) = (-z + 2.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5578 shape(24,0) = (4.*x - 4.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5581 shape(25,0) = (4.*x - 4.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5584 shape(26,0) = (4.*x - 4.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5587 shape(27,0) = (4.*x - 4.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5592 shape(28,1) = (4.*y - 4.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5595 shape(29,1) = (4.*y - 4.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5598 shape(30,1) = (4.*y - 4.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5601 shape(31,1) = (4.*y - 4.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5606 shape(32,2) = (4.*z - 4.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5609 shape(33,2) = (4.*z - 4.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5612 shape(34,2) = (4.*z - 4.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5615 shape(35,2) = (4.*z - 4.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5621 double x = ip.
x, y = ip.
y, z = ip.
z;
5623 divshape(2) = -(-3. + 4.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5624 divshape(3) = -(-3. + 4.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5625 divshape(0) = -(-3. + 4.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5626 divshape(1) = -(-3. + 4.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5628 divshape(4) = -(-3. + 4.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5629 divshape(5) = -(-3. + 4.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5630 divshape(6) = -(-3. + 4.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5631 divshape(7) = -(-3. + 4.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5633 divshape(8) = (-1. + 4.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5634 divshape(9) = (-1. + 4.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5635 divshape(10) = (-1. + 4.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5636 divshape(11) = (-1. + 4.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5638 divshape(13) = (-1. + 4.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5639 divshape(12) = (-1. + 4.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5640 divshape(15) = (-1. + 4.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5641 divshape(14) = (-1. + 4.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5643 divshape(17) = -(-3. + 4.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5644 divshape(16) = -(-3. + 4.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5645 divshape(19) = -(-3. + 4.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5646 divshape(18) = -(-3. + 4.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5648 divshape(20) = (-1. + 4.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5649 divshape(21) = (-1. + 4.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5650 divshape(22) = (-1. + 4.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5651 divshape(23) = (-1. + 4.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5653 divshape(24) = ( 4. - 8.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5654 divshape(25) = ( 4. - 8.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5655 divshape(26) = ( 4. - 8.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5656 divshape(27) = ( 4. - 8.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5658 divshape(28) = ( 4. - 8.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5659 divshape(29) = ( 4. - 8.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5660 divshape(30) = ( 4. - 8.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5661 divshape(31) = ( 4. - 8.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5663 divshape(32) = ( 4. - 8.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5664 divshape(33) = ( 4. - 8.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5665 divshape(34) = ( 4. - 8.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5666 divshape(35) = ( 4. - 8.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5669 const double RT1HexFiniteElement::nk[36][3] =
5671 {0, 0,-1}, {0, 0,-1}, {0, 0,-1}, {0, 0,-1},
5672 {0,-1, 0}, {0,-1, 0}, {0,-1, 0}, {0,-1, 0},
5673 {1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0},
5674 {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0},
5675 {-1,0, 0}, {-1,0, 0}, {-1,0, 0}, {-1,0, 0},
5676 {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1},
5677 {1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0},
5678 {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0},
5679 {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1}
5686 #ifdef MFEM_THREAD_SAFE
5692 for (k = 0; k < 36; k++)
5695 for (j = 0; j < 36; j++)
5697 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5699 if (j == k) { d -= 1.0; }
5700 if (fabs(d) > 1.0e-12)
5702 cerr <<
"RT0HexFiniteElement::GetLocalInterpolation (...)\n"
5703 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5711 ip.
x = ip.
y = ip.
z = 0.0;
5719 for (k = 0; k < 36; k++)
5722 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5725 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2];
5726 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2];
5727 vk[2] =
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2];
5728 for (j = 0; j < 36; j++)
5729 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5730 vshape(j,2)*vk[2])) < 1.0e-12)
5743 #ifdef MFEM_THREAD_SAFE
5747 for (
int k = 0; k < 36; k++)
5756 vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2] ) +
5757 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2] ) +
5758 vk[2] * (
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2] );
5786 double x2 = 2.0*ip.
x, y2 = 2.0*ip.
y, z2 = 2.0*ip.
z;
5792 shape(1,0) = x2 - 2.0;
5797 shape(2,1) = y2 - 2.0;
5802 shape(3,2) = z2 - 2.0;
5814 const double RT0TetFiniteElement::nk[4][3] =
5815 {{.5,.5,.5}, {-.5,0,0}, {0,-.5,0}, {0,0,-.5}};
5821 #ifdef MFEM_THREAD_SAFE
5827 for (k = 0; k < 4; k++)
5830 for (j = 0; j < 4; j++)
5832 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5834 if (j == k) { d -= 1.0; }
5835 if (fabs(d) > 1.0e-12)
5837 cerr <<
"RT0TetFiniteElement::GetLocalInterpolation (...)\n"
5838 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5846 ip.
x = ip.
y = ip.
z = 0.0;
5854 for (k = 0; k < 4; k++)
5857 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5860 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2];
5861 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2];
5862 vk[2] =
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2];
5863 for (j = 0; j < 4; j++)
5864 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5865 vshape(j,2)*vk[2])) < 1.0e-12)
5878 #ifdef MFEM_THREAD_SAFE
5882 for (
int k = 0; k < 4; k++)
5891 vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2] ) +
5892 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2] ) +
5893 vk[2] * (
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2] );
5928 double x = 2. * ip.
x - 1.;
5929 double y = 2. * ip.
y - 1.;
5930 double z = 2. * ip.
z - 1.;
5931 double f5 = x * x - y * y;
5932 double f6 = y * y - z * z;
5934 shape(0) = (1./6.) * (1. - 3. * z - f5 - 2. * f6);
5935 shape(1) = (1./6.) * (1. - 3. * y - f5 + f6);
5936 shape(2) = (1./6.) * (1. + 3. * x + 2. * f5 + f6);
5937 shape(3) = (1./6.) * (1. + 3. * y - f5 + f6);
5938 shape(4) = (1./6.) * (1. - 3. * x + 2. * f5 + f6);
5939 shape(5) = (1./6.) * (1. + 3. * z - f5 - 2. * f6);
5945 const double a = 2./3.;
5947 double xt = a * (1. - 2. * ip.
x);
5948 double yt = a * (1. - 2. * ip.
y);
5949 double zt = a * (1. - 2. * ip.
z);
5953 dshape(0,2) = -1. - 2. * zt;
5956 dshape(1,1) = -1. - 2. * yt;
5959 dshape(2,0) = 1. - 2. * xt;
5964 dshape(3,1) = 1. - 2. * yt;
5967 dshape(4,0) = -1. - 2. * xt;
5973 dshape(5,2) = 1. - 2. * zt;
5978 : x(p + 1), w(p + 1)
5984 for (
int i = 0; i <= p; i++)
5987 for (
int j = 0; j <= p; j++)
6000 for (
int i = 0; i <= p; i++)
6001 for (
int j = 0; j < i; j++)
6003 double xij = x(i) - x(j);
6007 for (
int i = 0; i <= p; i++)
6014 for (
int i = 0; i < p; i++)
6017 mfem_error(
"Poly_1D::Basis::Basis : nodes are not increasing!");
6032 int i, k, p = x.Size() - 1;
6042 for (k = 0; k < p; k++)
6043 if (y >= (x(k) + x(k+1))/2)
6049 for (i = k+1; i <= p; i++)
6055 l = lk * (y - x(k));
6057 for (i = 0; i < k; i++)
6059 u(i) = l * w(i) / (y - x(i));
6062 for (i++; i <= p; i++)
6064 u(i) = l * w(i) / (y - x(i));
6079 int i, k, p = x.Size() - 1;
6080 double l, lp, lk, sk, si;
6090 for (k = 0; k < p; k++)
6091 if (y >= (x(k) + x(k+1))/2)
6097 for (i = k+1; i <= p; i++)
6103 l = lk * (y - x(k));
6106 for (i = 0; i < k; i++)
6108 si = 1.0/(y - x(i));
6110 u(i) = l * si * w(i);
6113 for (i++; i <= p; i++)
6115 si = 1.0/(y - x(i));
6117 u(i) = l * si * w(i);
6121 for (i = 0; i < k; i++)
6123 d(i) = (lp * w(i) - u(i))/(y - x(i));
6126 for (i++; i <= p; i++)
6128 d(i) = (lp * w(i) - u(i))/(y - x(i));
6138 for (
int i = 0; i <= p; i++)
6140 binom(i,0) = binom(i,i) = 1;
6141 for (
int j = 1; j < i; j++)
6143 binom(i,j) = binom(i-1,j) + binom(i-1,j-1);
6158 for (
int i = 0; i <= p; i++)
6167 int m = (p+1)/2, odd_p = p%2;
6174 for (
int i = 0; i < m; i++)
6176 double z, y, p0, d0;
6177 z = cos(M_PI*(i + 0.75)/(p + 1.5));
6179 for (
int k = 0;
true; )
6189 for (
int n = 1;
true; n++)
6191 p0 = ((2*n+1)*z*p1 - n*p2)/(n + 1);
6196 if (n == p) {
break; }
6203 if (fabs(p0/d0) < 2e-16) {
break; }
6207 std::cout <<
"Poly_1D::GaussPoints : No convergence!"
6208 " p = " << p <<
", i = " << i <<
", p0/d0 = " << p0/d0
6217 y = ((1 - z) + p0/d0)/2;
6235 if (p == 1) {
return; }
6238 int m = (p - 1)/2, odd_p = p%2;
6244 for (
int i = 0; i < m; )
6246 double y, z, d0, s0;
6247 z = cos(M_PI*(i + 1)/p);
6249 for (
int k = 0;
true; )
6256 double p0, p1, p2, d1;
6262 for (
int n = 1;
true; n++)
6264 p0 = ((2*n+1)*z*p1 - n*p2)/(n + 1);
6274 if (n == p - 1) {
break; }
6280 if (fabs(d0/s0) < 2e-16) {
break; }
6284 std::cout <<
"Poly_1D::GaussLobattoPoints : No convergence!"
6285 " p = " << p <<
", i = " << i <<
", d0/s0 = " << d0/s0
6293 y = ((1 - z) + d0/s0)/2;
6303 for (
int i = 0; i <= p; i++)
6306 double s = sin(M_PI_2*(i + 0.5)/(p + 1));
6311 void Poly_1D::CalcMono(
const int p,
const double x,
double *u)
6315 for (
int n = 1; n <= p; n++)
6321 void Poly_1D::CalcMono(
const int p,
const double x,
double *u,
double *d)
6326 for (
int n = 1; n <= p; n++)
6343 const int *b =
Binom(p);
6346 for (i = 1; i < p; i++)
6353 for (i--; i > 0; i--)
6363 double *u,
double *d)
6373 const int *b =
Binom(p);
6374 const double xpy = x + y, ptx = p*x;
6377 for (i = 1; i < p; i++)
6379 d[i] = b[i]*z*(i*xpy - ptx);
6386 for (i--; i > 0; i--)
6407 const int *b =
Binom(p);
6408 const double xpy = x + y, ptx = p*x;
6411 for (i = 1; i < p; i++)
6413 d[i] = b[i]*z*(i*xpy - ptx);
6418 for (i--; i > 0; i--)
6427 void Poly_1D::CalcLegendre(
const int p,
const double x,
double *u)
6433 if (p == 0) {
return; }
6434 u[1] = z = 2.*x - 1.;
6435 for (
int n = 1; n < p; n++)
6437 u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
6441 void Poly_1D::CalcLegendre(
const int p,
const double x,
double *u,
double *d)
6450 if (p == 0) {
return; }
6451 u[1] = z = 2.*x - 1.;
6453 for (
int n = 1; n < p; n++)
6455 u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
6456 d[n+1] = (4*n + 2)*u[n] + d[n-1];
6460 void Poly_1D::CalcChebyshev(
const int p,
const double x,
double *u)
6467 if (p == 0) {
return; }
6468 u[1] = z = 2.*x - 1.;
6469 for (
int n = 1; n < p; n++)
6471 u[n+1] = 2*z*u[n] - u[n-1];
6475 void Poly_1D::CalcChebyshev(
const int p,
const double x,
double *u,
double *d)
6488 if (p == 0) {
return; }
6489 u[1] = z = 2.*x - 1.;
6491 for (
int n = 1; n < p; n++)
6493 u[n+1] = 2*z*u[n] - u[n-1];
6494 d[n+1] = (n + 1)*(z*d[n]/n + 2*u[n]);
6502 if (open_pts.
Size() <= p)
6504 int i = open_pts.
Size();
6512 if ((op = open_pts[p]) != NULL)
6517 open_pts[p] = op =
new double[p + 1];
6527 if (closed_pts.
Size() <= p)
6529 int i = closed_pts.
Size();
6533 closed_pts[i] = NULL;
6537 if ((cp = closed_pts[p]) != NULL)
6542 closed_pts[p] = cp =
new double[p + 1];
6552 if (open_basis.Size() <= p)
6554 int i = open_basis.Size();
6555 open_basis.SetSize(p + 1);
6558 open_basis[i] = NULL;
6562 if ((ob = open_basis[p]) != NULL)
6575 if (closed_basis.Size() <= p)
6577 int i = closed_basis.Size();
6578 closed_basis.SetSize(p + 1);
6581 closed_basis[i] = NULL;
6585 if ((cb = closed_basis[p]) != NULL)
6596 for (
int i = 0; i < open_pts.
Size(); i++)
6598 delete [] open_pts[i];
6600 for (
int i = 0; i < closed_pts.
Size(); i++)
6602 delete [] closed_pts[i];
6604 for (
int i = 0; i < open_basis.Size(); i++)
6606 delete open_basis[i];
6608 for (
int i = 0; i < closed_basis.Size(); i++)
6610 delete closed_basis[i];
6625 #ifndef MFEM_THREAD_SAFE
6634 for (
int i = 1; i < p; i++)
6644 const int p =
Order;
6646 #ifdef MFEM_THREAD_SAFE
6650 basis1d.
Eval(ip.
x, shape_x);
6652 shape(0) = shape_x(0);
6653 shape(1) = shape_x(p);
6654 for (
int i = 1; i < p; i++)
6656 shape(i+1) = shape_x(i);
6663 const int p =
Order;
6665 #ifdef MFEM_THREAD_SAFE
6666 Vector shape_x(p+1), dshape_x(p+1);
6669 basis1d.
Eval(ip.
x, shape_x, dshape_x);
6671 dshape(0,0) = dshape_x(0);
6672 dshape(1,0) = dshape_x(p);
6673 for (
int i = 1; i < p; i++)
6675 dshape(i+1,0) = dshape_x(i);
6681 const int p =
Order;
6689 for (
int i = 1; i < p; i++)
6698 for (
int i = 1; i < p; i++)
6710 basis1d(
poly1d.ClosedBasis(p)), dof_map((p + 1)*(p + 1))
6714 const int p1 = p + 1;
6716 #ifndef MFEM_THREAD_SAFE
6724 dof_map[0 + 0*p1] = 0;
6725 dof_map[p + 0*p1] = 1;
6726 dof_map[p + p*p1] = 2;
6727 dof_map[0 + p*p1] = 3;
6731 for (
int i = 1; i < p; i++)
6733 dof_map[i + 0*p1] = o++;
6735 for (
int i = 1; i < p; i++)
6737 dof_map[p + i*p1] = o++;
6739 for (
int i = 1; i < p; i++)
6741 dof_map[(p-i) + p*p1] = o++;
6743 for (
int i = 1; i < p; i++)
6745 dof_map[0 + (p-i)*p1] = o++;
6749 for (
int j = 1; j < p; j++)
6750 for (
int i = 1; i < p; i++)
6752 dof_map[i + j*p1] = o++;
6756 for (
int j = 0; j <= p; j++)
6757 for (
int i = 0; i <= p; i++)
6766 const int p =
Order;
6768 #ifdef MFEM_THREAD_SAFE
6769 Vector shape_x(p+1), shape_y(p+1);
6772 basis1d.
Eval(ip.
x, shape_x);
6773 basis1d.
Eval(ip.
y, shape_y);
6775 for (
int o = 0, j = 0; j <= p; j++)
6776 for (
int i = 0; i <= p; i++)
6778 shape(dof_map[o++]) = shape_x(i)*shape_y(j);
6785 const int p =
Order;
6787 #ifdef MFEM_THREAD_SAFE
6788 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
6791 basis1d.
Eval(ip.
x, shape_x, dshape_x);
6792 basis1d.
Eval(ip.
y, shape_y, dshape_y);
6794 for (
int o = 0, j = 0; j <= p; j++)
6795 for (
int i = 0; i <= p; i++)
6797 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j);
6798 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j); o++;
6804 const int p =
Order;
6807 #ifdef MFEM_THREAD_SAFE
6808 Vector shape_x(p+1), shape_y(p+1);
6811 for (
int i = 0; i <= p; i++)
6820 for (
int o = 0, j = 0; j <= p; j++)
6821 for (
int i = 0; i <= p; i++)
6823 dofs(dof_map[o++]) = shape_x(i)*shape_x(j);
6827 for (
int o = 0, j = 0; j <= p; j++)
6828 for (
int i = 0; i <= p; i++)
6830 dofs(dof_map[o++]) = shape_y(i)*shape_x(j);
6834 for (
int o = 0, j = 0; j <= p; j++)
6835 for (
int i = 0; i <= p; i++)
6837 dofs(dof_map[o++]) = shape_y(i)*shape_y(j);
6841 for (
int o = 0, j = 0; j <= p; j++)
6842 for (
int i = 0; i <= p; i++)
6844 dofs(dof_map[o++]) = shape_x(i)*shape_y(j);
6854 basis1d(
poly1d.ClosedBasis(p)), dof_map((p + 1)*(p + 1)*(p + 1))
6858 const int p1 = p + 1;
6860 #ifndef MFEM_THREAD_SAFE
6870 dof_map[0 + (0 + 0*p1)*p1] = 0;
6871 dof_map[p + (0 + 0*p1)*p1] = 1;
6872 dof_map[p + (p + 0*p1)*p1] = 2;
6873 dof_map[0 + (p + 0*p1)*p1] = 3;
6874 dof_map[0 + (0 + p*p1)*p1] = 4;
6875 dof_map[p + (0 + p*p1)*p1] = 5;
6876 dof_map[p + (p + p*p1)*p1] = 6;
6877 dof_map[0 + (p + p*p1)*p1] = 7;
6881 for (
int i = 1; i < p; i++)
6883 dof_map[i + (0 + 0*p1)*p1] = o++;
6885 for (
int i = 1; i < p; i++)
6887 dof_map[p + (i + 0*p1)*p1] = o++;
6889 for (
int i = 1; i < p; i++)
6891 dof_map[i + (p + 0*p1)*p1] = o++;
6893 for (
int i = 1; i < p; i++)
6895 dof_map[0 + (i + 0*p1)*p1] = o++;
6897 for (
int i = 1; i < p; i++)
6899 dof_map[i + (0 + p*p1)*p1] = o++;
6901 for (
int i = 1; i < p; i++)
6903 dof_map[p + (i + p*p1)*p1] = o++;
6905 for (
int i = 1; i < p; i++)
6907 dof_map[i + (p + p*p1)*p1] = o++;
6909 for (
int i = 1; i < p; i++)
6911 dof_map[0 + (i + p*p1)*p1] = o++;
6913 for (
int i = 1; i < p; i++)
6915 dof_map[0 + (0 + i*p1)*p1] = o++;
6917 for (
int i = 1; i < p; i++)
6919 dof_map[p + (0 + i*p1)*p1] = o++;
6921 for (
int i = 1; i < p; i++)
6923 dof_map[p + (p + i*p1)*p1] = o++;
6925 for (
int i = 1; i < p; i++)
6927 dof_map[0 + (p + i*p1)*p1] = o++;
6931 for (
int j = 1; j < p; j++)
6932 for (
int i = 1; i < p; i++)
6934 dof_map[i + ((p-j) + 0*p1)*p1] = o++;
6936 for (
int j = 1; j < p; j++)
6937 for (
int i = 1; i < p; i++)
6939 dof_map[i + (0 + j*p1)*p1] = o++;
6941 for (
int j = 1; j < p; j++)
6942 for (
int i = 1; i < p; i++)
6944 dof_map[p + (i + j*p1)*p1] = o++;
6946 for (
int j = 1; j < p; j++)
6947 for (
int i = 1; i < p; i++)
6949 dof_map[(p-i) + (p + j*p1)*p1] = o++;
6951 for (
int j = 1; j < p; j++)
6952 for (
int i = 1; i < p; i++)
6954 dof_map[0 + ((p-i) + j*p1)*p1] = o++;
6956 for (
int j = 1; j < p; j++)
6957 for (
int i = 1; i < p; i++)
6959 dof_map[i + (j + p*p1)*p1] = o++;
6963 for (
int k = 1; k < p; k++)
6964 for (
int j = 1; j < p; j++)
6965 for (
int i = 1; i < p; i++)
6967 dof_map[i + (j + k*p1)*p1] = o++;
6971 for (
int k = 0; k <= p; k++)
6972 for (
int j = 0; j <= p; j++)
6973 for (
int i = 0; i <= p; i++)
6982 const int p =
Order;
6984 #ifdef MFEM_THREAD_SAFE
6985 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
6988 basis1d.
Eval(ip.
x, shape_x);
6989 basis1d.
Eval(ip.
y, shape_y);
6990 basis1d.
Eval(ip.
z, shape_z);
6992 for (
int o = 0, k = 0; k <= p; k++)
6993 for (
int j = 0; j <= p; j++)
6994 for (
int i = 0; i <= p; i++)
6996 shape(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_z(k);
7003 const int p =
Order;
7005 #ifdef MFEM_THREAD_SAFE
7006 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7007 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
7010 basis1d.
Eval(ip.
x, shape_x, dshape_x);
7011 basis1d.
Eval(ip.
y, shape_y, dshape_y);
7012 basis1d.
Eval(ip.
z, shape_z, dshape_z);
7014 for (
int o = 0, k = 0; k <= p; k++)
7015 for (
int j = 0; j <= p; j++)
7016 for (
int i = 0; i <= p; i++)
7018 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j)* shape_z(k);
7019 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j)* shape_z(k);
7020 dshape(dof_map[o],2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
7026 const int p =
Order;
7029 #ifdef MFEM_THREAD_SAFE
7030 Vector shape_x(p+1), shape_y(p+1);
7033 for (
int i = 0; i <= p; i++)
7042 for (
int o = 0, k = 0; k <= p; k++)
7043 for (
int j = 0; j <= p; j++)
7044 for (
int i = 0; i <= p; i++)
7046 dofs(dof_map[o++]) = shape_x(i)*shape_x(j)*shape_x(k);
7050 for (
int o = 0, k = 0; k <= p; k++)
7051 for (
int j = 0; j <= p; j++)
7052 for (
int i = 0; i <= p; i++)
7054 dofs(dof_map[o++]) = shape_y(i)*shape_x(j)*shape_x(k);
7058 for (
int o = 0, k = 0; k <= p; k++)
7059 for (
int j = 0; j <= p; j++)
7060 for (
int i = 0; i <= p; i++)
7062 dofs(dof_map[o++]) = shape_y(i)*shape_y(j)*shape_x(k);
7066 for (
int o = 0, k = 0; k <= p; k++)
7067 for (
int j = 0; j <= p; j++)
7068 for (
int i = 0; i <= p; i++)
7070 dofs(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_x(k);
7074 for (
int o = 0, k = 0; k <= p; k++)
7075 for (
int j = 0; j <= p; j++)
7076 for (
int i = 0; i <= p; i++)
7078 dofs(dof_map[o++]) = shape_x(i)*shape_x(j)*shape_y(k);
7082 for (
int o = 0, k = 0; k <= p; k++)
7083 for (
int j = 0; j <= p; j++)
7084 for (
int i = 0; i <= p; i++)
7086 dofs(dof_map[o++]) = shape_y(i)*shape_x(j)*shape_y(k);
7090 for (
int o = 0, k = 0; k <= p; k++)
7091 for (
int j = 0; j <= p; j++)
7092 for (
int i = 0; i <= p; i++)
7094 dofs(dof_map[o++]) = shape_y(i)*shape_y(j)*shape_y(k);
7098 for (
int o = 0, k = 0; k <= p; k++)
7099 for (
int j = 0; j <= p; j++)
7100 for (
int i = 0; i <= p; i++)
7102 dofs(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_y(k);
7113 #ifndef MFEM_THREAD_SAFE
7124 for (
int i = 1; i < p; i++)
7134 const int p =
Order;
7136 #ifdef MFEM_THREAD_SAFE
7143 shape(0) = shape_x(0);
7144 shape(1) = shape_x(p);
7145 for (
int i = 1; i < p; i++)
7147 shape(i+1) = shape_x(i);
7154 const int p =
Order;
7156 #ifdef MFEM_THREAD_SAFE
7157 Vector shape_x(p+1), dshape_x(p+1);
7163 dshape(0,0) = dshape_x(0);
7164 dshape(1,0) = dshape_x(p);
7165 for (
int i = 1; i < p; i++)
7167 dshape(i+1,0) = dshape_x(i);
7181 dof_map((p + 1)*(p + 1))
7183 const int p1 = p + 1;
7185 #ifndef MFEM_THREAD_SAFE
7194 dof_map[0 + 0*p1] = 0;
7195 dof_map[p + 0*p1] = 1;
7196 dof_map[p + p*p1] = 2;
7197 dof_map[0 + p*p1] = 3;
7201 for (
int i = 1; i < p; i++)
7203 dof_map[i + 0*p1] = o++;
7205 for (
int i = 1; i < p; i++)
7207 dof_map[p + i*p1] = o++;
7209 for (
int i = 1; i < p; i++)
7211 dof_map[(p-i) + p*p1] = o++;
7213 for (
int i = 1; i < p; i++)
7215 dof_map[0 + (p-i)*p1] = o++;
7219 for (
int j = 1; j < p; j++)
7220 for (
int i = 1; i < p; i++)
7222 dof_map[i + j*p1] = o++;
7226 for (
int j = 0; j <= p; j++)
7227 for (
int i = 0; i <= p; i++)
7236 const int p =
Order;
7238 #ifdef MFEM_THREAD_SAFE
7239 Vector shape_x(p+1), shape_y(p+1);
7246 for (
int o = 0, j = 0; j <= p; j++)
7247 for (
int i = 0; i <= p; i++)
7249 shape(dof_map[o++]) = shape_x(i)*shape_y(j);
7256 const int p =
Order;
7258 #ifdef MFEM_THREAD_SAFE
7259 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
7266 for (
int o = 0, j = 0; j <= p; j++)
7267 for (
int i = 0; i <= p; i++)
7269 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j);
7270 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j); o++;
7284 dof_map((p + 1)*(p + 1)*(p + 1))
7286 const int p1 = p + 1;
7288 #ifndef MFEM_THREAD_SAFE
7299 dof_map[0 + (0 + 0*p1)*p1] = 0;
7300 dof_map[p + (0 + 0*p1)*p1] = 1;
7301 dof_map[p + (p + 0*p1)*p1] = 2;
7302 dof_map[0 + (p + 0*p1)*p1] = 3;
7303 dof_map[0 + (0 + p*p1)*p1] = 4;
7304 dof_map[p + (0 + p*p1)*p1] = 5;
7305 dof_map[p + (p + p*p1)*p1] = 6;
7306 dof_map[0 + (p + p*p1)*p1] = 7;
7310 for (
int i = 1; i < p; i++)
7312 dof_map[i + (0 + 0*p1)*p1] = o++;
7314 for (
int i = 1; i < p; i++)
7316 dof_map[p + (i + 0*p1)*p1] = o++;
7318 for (
int i = 1; i < p; i++)
7320 dof_map[i + (p + 0*p1)*p1] = o++;
7322 for (
int i = 1; i < p; i++)
7324 dof_map[0 + (i + 0*p1)*p1] = o++;
7326 for (
int i = 1; i < p; i++)
7328 dof_map[i + (0 + p*p1)*p1] = o++;
7330 for (
int i = 1; i < p; i++)
7332 dof_map[p + (i + p*p1)*p1] = o++;
7334 for (
int i = 1; i < p; i++)
7336 dof_map[i + (p + p*p1)*p1] = o++;
7338 for (
int i = 1; i < p; i++)
7340 dof_map[0 + (i + p*p1)*p1] = o++;
7342 for (
int i = 1; i < p; i++)
7344 dof_map[0 + (0 + i*p1)*p1] = o++;
7346 for (
int i = 1; i < p; i++)
7348 dof_map[p + (0 + i*p1)*p1] = o++;
7350 for (
int i = 1; i < p; i++)
7352 dof_map[p + (p + i*p1)*p1] = o++;
7354 for (
int i = 1; i < p; i++)
7356 dof_map[0 + (p + i*p1)*p1] = o++;
7360 for (
int j = 1; j < p; j++)
7361 for (
int i = 1; i < p; i++)
7363 dof_map[i + ((p-j) + 0*p1)*p1] = o++;
7365 for (
int j = 1; j < p; j++)
7366 for (
int i = 1; i < p; i++)
7368 dof_map[i + (0 + j*p1)*p1] = o++;
7370 for (
int j = 1; j < p; j++)
7371 for (
int i = 1; i < p; i++)
7373 dof_map[p + (i + j*p1)*p1] = o++;
7375 for (
int j = 1; j < p; j++)
7376 for (
int i = 1; i < p; i++)
7378 dof_map[(p-i) + (p + j*p1)*p1] = o++;
7380 for (
int j = 1; j < p; j++)
7381 for (
int i = 1; i < p; i++)
7383 dof_map[0 + ((p-i) + j*p1)*p1] = o++;
7385 for (
int j = 1; j < p; j++)
7386 for (
int i = 1; i < p; i++)
7388 dof_map[i + (j + p*p1)*p1] = o++;
7392 for (
int k = 1; k < p; k++)
7393 for (
int j = 1; j < p; j++)
7394 for (
int i = 1; i < p; i++)
7396 dof_map[i + (j + k*p1)*p1] = o++;
7400 for (
int k = 0; k <= p; k++)
7401 for (
int j = 0; j <= p; j++)
7402 for (
int i = 0; i <= p; i++)
7410 const int p =
Order;
7412 #ifdef MFEM_THREAD_SAFE
7413 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7420 for (
int o = 0, k = 0; k <= p; k++)
7421 for (
int j = 0; j <= p; j++)
7422 for (
int i = 0; i <= p; i++)
7424 shape(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_z(k);
7431 const int p =
Order;
7433 #ifdef MFEM_THREAD_SAFE
7434 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7435 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
7442 for (
int o = 0, k = 0; k <= p; k++)
7443 for (
int j = 0; j <= p; j++)
7444 for (
int i = 0; i <= p; i++)
7446 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j)* shape_z(k);
7447 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j)* shape_z(k);
7448 dshape(dof_map[o],2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
7465 #ifndef MFEM_THREAD_SAFE
7475 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
7485 for (
int i = 1; i < p; i++)
7489 for (
int i = 1; i < p; i++)
7493 for (
int i = 1; i < p; i++)
7499 for (
int j = 1; j < p; j++)
7500 for (
int i = 1; i + j < p; i++)
7502 const double w = cp[i] + cp[j] + cp[p-i-j];
7507 for (
int k = 0; k <
Dof; k++)
7515 for (
int j = 0; j <= p; j++)
7516 for (
int i = 0; i + j <= p; i++)
7518 T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
7529 const int p =
Order;
7531 #ifdef MFEM_THREAD_SAFE
7532 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(
Dof);
7539 for (
int o = 0, j = 0; j <= p; j++)
7540 for (
int i = 0; i + j <= p; i++)
7542 u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
7551 const int p =
Order;
7553 #ifdef MFEM_THREAD_SAFE
7554 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
7555 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
7563 for (
int o = 0, j = 0; j <= p; j++)
7564 for (
int i = 0; i + j <= p; i++)
7567 du(o,0) = ((dshape_x(i)* shape_l(k)) -
7568 ( shape_x(i)*dshape_l(k)))*shape_y(j);
7569 du(o,1) = ((dshape_y(j)* shape_l(k)) -
7570 ( shape_y(j)*dshape_l(k)))*shape_x(i);
7574 Ti.
Mult(du, dshape);
7584 #ifndef MFEM_THREAD_SAFE
7596 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7607 for (
int i = 1; i < p; i++)
7611 for (
int i = 1; i < p; i++)
7615 for (
int i = 1; i < p; i++)
7619 for (
int i = 1; i < p; i++)
7623 for (
int i = 1; i < p; i++)
7627 for (
int i = 1; i < p; i++)
7633 for (
int j = 1; j < p; j++)
7634 for (
int i = 1; i + j < p; i++)
7636 double w = cp[i] + cp[j] + cp[p-i-j];
7639 for (
int j = 1; j < p; j++)
7640 for (
int i = 1; i + j < p; i++)
7642 double w = cp[i] + cp[j] + cp[p-i-j];
7645 for (
int j = 1; j < p; j++)
7646 for (
int i = 1; i + j < p; i++)
7648 double w = cp[i] + cp[j] + cp[p-i-j];
7651 for (
int j = 1; j < p; j++)
7652 for (
int i = 1; i + j < p; i++)
7654 double w = cp[i] + cp[j] + cp[p-i-j];
7659 for (
int k = 1; k < p; k++)
7660 for (
int j = 1; j + k < p; j++)
7661 for (
int i = 1; i + j + k < p; i++)
7663 double w = cp[i] + cp[j] + cp[k] + cp[p-i-j-k];
7668 for (
int m = 0; m <
Dof; m++)
7677 for (
int k = 0; k <= p; k++)
7678 for (
int j = 0; j + k <= p; j++)
7679 for (
int i = 0; i + j + k <= p; i++)
7681 T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
7692 const int p =
Order;
7694 #ifdef MFEM_THREAD_SAFE
7695 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7704 for (
int o = 0, k = 0; k <= p; k++)
7705 for (
int j = 0; j + k <= p; j++)
7706 for (
int i = 0; i + j + k <= p; i++)
7708 u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
7717 const int p =
Order;
7719 #ifdef MFEM_THREAD_SAFE
7720 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7721 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
7730 for (
int o = 0, k = 0; k <= p; k++)
7731 for (
int j = 0; j + k <= p; j++)
7732 for (
int i = 0; i + j + k <= p; i++)
7734 int l = p - i - j - k;
7735 du(o,0) = ((dshape_x(i)* shape_l(l)) -
7736 ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
7737 du(o,1) = ((dshape_y(j)* shape_l(l)) -
7738 ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
7739 du(o,2) = ((dshape_z(k)* shape_l(l)) -
7740 ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
7744 Ti.
Mult(du, dshape);
7752 #ifndef MFEM_THREAD_SAFE
7762 Index(
int p) { p2p3 = 2*p + 3; }
7763 int operator()(
int i,
int j) {
return ((p2p3-j)*j)/2+i; }
7777 for (
int i = 1; i < p; i++)
7782 for (
int i = 1; i < p; i++)
7787 for (
int i = 1; i < p; i++)
7794 for (
int j = 1; j < p; j++)
7795 for (
int i = 1; i + j < p; i++)
7804 const int p,
const double l1,
const double l2,
double *shape)
7806 const double l3 = 1. - l1 - l2;
7816 for (
int o = 0, j = 0; j <= p; j++)
7820 for (
int i = 0; i <= p - j; i++)
7830 const int p,
const double l1,
const double l2,
7831 double *dshape_1d,
double *dshape)
7833 const int dof = ((p + 1)*(p + 2))/2;
7834 const double l3 = 1. - l1 - l2;
7838 for (
int o = 0, j = 0; j <= p; j++)
7842 for (
int i = 0; i <= p - j; i++)
7844 dshape[o++] = s*dshape_1d[i];
7849 for (
int i = 0; i <= p; i++)
7853 for (
int o = i, j = 0; j <= p - i; j++)
7855 dshape[dof + o] = s*dshape_1d[j];
7865 #ifdef MFEM_THREAD_SAFE
7869 for (
int i = 0; i <
Dof; i++)
7878 #ifdef MFEM_THREAD_SAFE
7883 for (
int d = 0; d < 2; d++)
7885 for (
int i = 0; i <
Dof; i++)
7897 #ifndef MFEM_THREAD_SAFE
7907 int tri(
int k) {
return (k*(k + 1))/2; }
7908 int tet(
int k) {
return (k*(k + 1)*(k + 2))/6; }
7909 Index(
int p_) { p = p_; dof = tet(p + 1); }
7910 int operator()(
int i,
int j,
int k)
7911 {
return dof - tet(p - k) - tri(p + 1 - k - j) + i; }
7927 for (
int i = 1; i < p; i++)
7932 for (
int i = 1; i < p; i++)
7937 for (
int i = 1; i < p; i++)
7942 for (
int i = 1; i < p; i++)
7947 for (
int i = 1; i < p; i++)
7952 for (
int i = 1; i < p; i++)
7959 for (
int j = 1; j < p; j++)
7960 for (
int i = 1; i + j < p; i++)
7965 for (
int j = 1; j < p; j++)
7966 for (
int i = 1; i + j < p; i++)
7971 for (
int j = 1; j < p; j++)
7972 for (
int i = 1; i + j < p; i++)
7977 for (
int j = 1; j < p; j++)
7978 for (
int i = 1; i + j < p; i++)
7985 for (
int k = 1; k < p; k++)
7986 for (
int j = 1; j + k < p; j++)
7987 for (
int i = 1; i + j + k < p; i++)
7996 const int p,
const double l1,
const double l2,
const double l3,
7999 const double l4 = 1. - l1 - l2 - l3;
8008 for (
int o = 0, k = 0; k <= p; k++)
8011 const double ek = bp[k]*l3k;
8013 for (
int j = 0; j <= p - k; j++)
8016 double ekj = ek*bpk[j]*l2j;
8017 for (
int i = 0; i <= p - k - j; i++)
8029 const int p,
const double l1,
const double l2,
const double l3,
8030 double *dshape_1d,
double *dshape)
8032 const int dof = ((p + 1)*(p + 2)*(p + 3))/6;
8033 const double l4 = 1. - l1 - l2 - l3;
8041 for (
int o = 0, k = 0; k <= p; k++)
8044 const double ek = bp[k]*l3k;
8046 for (
int j = 0; j <= p - k; j++)
8049 double ekj = ek*bpk[j]*l2j;
8050 for (
int i = 0; i <= p - k - j; i++)
8052 dshape[o++] = dshape_1d[i]*ekj;
8063 for (
int ok = 0, k = 0; k <= p; k++)
8066 const double ek = bp[k]*l3k;
8068 for (
int i = 0; i <= p - k; i++)
8071 double eki = ek*bpk[i]*l1i;
8073 for (
int j = 0; j <= p - k - i; j++)
8075 dshape[dof + o] = dshape_1d[j]*eki;
8081 ok += ((p - k + 2)*(p - k + 1))/2;
8088 for (
int j = 0; j <= p; j++)
8091 const double ej = bp[j]*l2j;
8093 for (
int i = 0; i <= p - j; i++)
8096 double eji = ej*bpj[i]*l1i;
8097 int m = ((p + 2)*(p + 1))/2;
8098 int n = ((p - j + 2)*(p - j + 1))/2;
8099 for (
int o = i, k = 0; k <= p - j - i; k++)
8104 dshape[2*dof + o - n] = dshape_1d[k]*eji;
8117 #ifdef MFEM_THREAD_SAFE
8121 for (
int i = 0; i <
Dof; i++)
8130 #ifdef MFEM_THREAD_SAFE
8135 for (
int d = 0; d < 3; d++)
8137 for (
int i = 0; i <
Dof; i++)
8163 #ifndef MFEM_THREAD_SAFE
8168 for (
int i = 0; i <= p; i++)
8177 basis1d->
Eval(ip.
x, shape);
8183 #ifdef MFEM_THREAD_SAFE
8188 basis1d->
Eval(ip.
x, shape_x, dshape_x);
8194 const int p =
Order;
8207 for (
int i = 0; i <= p; i++)
8214 for (
int i = 0; i <= p; i++)
8226 #ifndef MFEM_THREAD_SAFE
8237 for (
int i = 0; i <= p; i++)
8253 #ifdef MFEM_THREAD_SAFE
8264 dofs[vertex*
Order] = 1.0;
8287 #ifndef MFEM_THREAD_SAFE
8294 for (
int o = 0, j = 0; j <= p; j++)
8295 for (
int i = 0; i <= p; i++)
8304 const int p =
Order;
8306 #ifdef MFEM_THREAD_SAFE
8307 Vector shape_x(p+1), shape_y(p+1);
8310 basis1d->
Eval(ip.
x, shape_x);
8311 basis1d->
Eval(ip.
y, shape_y);
8313 for (
int o = 0, j = 0; j <= p; j++)
8314 for (
int i = 0; i <= p; i++)
8316 shape(o++) = shape_x(i)*shape_y(j);
8323 const int p =
Order;
8325 #ifdef MFEM_THREAD_SAFE
8326 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
8329 basis1d->
Eval(ip.
x, shape_x, dshape_x);
8330 basis1d->
Eval(ip.
y, shape_y, dshape_y);
8332 for (
int o = 0, j = 0; j <= p; j++)
8333 for (
int i = 0; i <= p; i++)
8335 dshape(o,0) = dshape_x(i)* shape_y(j);
8336 dshape(o,1) = shape_x(i)*dshape_y(j); o++;
8342 const int p =
Order;
8352 #ifdef MFEM_THREAD_SAFE
8353 Vector shape_x(p+1), shape_y(p+1);
8356 for (
int i = 0; i <= p; i++)
8365 for (
int o = 0, j = 0; j <= p; j++)
8366 for (
int i = 0; i <= p; i++)
8368 dofs[o++] = shape_x(i)*shape_x(j);
8372 for (
int o = 0, j = 0; j <= p; j++)
8373 for (
int i = 0; i <= p; i++)
8375 dofs[o++] = shape_y(i)*shape_x(j);
8379 for (
int o = 0, j = 0; j <= p; j++)
8380 for (
int i = 0; i <= p; i++)
8382 dofs[o++] = shape_y(i)*shape_y(j);
8386 for (
int o = 0, j = 0; j <= p; j++)
8387 for (
int i = 0; i <= p; i++)
8389 dofs[o++] = shape_x(i)*shape_y(j);
8400 #ifndef MFEM_THREAD_SAFE
8413 for (
int o = 0, j = 0; j <= p; j++)
8414 for (
int i = 0; i <= p; i++)
8424 const int p =
Order;
8426 #ifdef MFEM_THREAD_SAFE
8427 Vector shape_x(p+1), shape_y(p+1);
8433 for (
int o = 0, j = 0; j <= p; j++)
8434 for (
int i = 0; i <= p; i++)
8436 shape(o++) = shape_x(i)*shape_y(j);
8443 const int p =
Order;
8445 #ifdef MFEM_THREAD_SAFE
8446 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
8452 for (
int o = 0, j = 0; j <= p; j++)
8453 for (
int i = 0; i <= p; i++)
8455 dshape(o,0) = dshape_x(i)* shape_y(j);
8456 dshape(o,1) = shape_x(i)*dshape_y(j); o++;
8462 const int p =
Order;
8467 case 0: dofs[0] = 1.0;
break;
8468 case 1: dofs[p] = 1.0;
break;
8469 case 2: dofs[p*(p + 2)] = 1.0;
break;
8470 case 3: dofs[p*(p + 1)] = 1.0;
break;
8494 #ifndef MFEM_THREAD_SAFE
8503 for (
int o = 0, k = 0; k <= p; k++)
8504 for (
int j = 0; j <= p; j++)
8505 for (
int i = 0; i <= p; i++)
8514 const int p =
Order;
8516 #ifdef MFEM_THREAD_SAFE
8517 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
8520 basis1d->
Eval(ip.
x, shape_x);
8521 basis1d->
Eval(ip.
y, shape_y);
8522 basis1d->
Eval(ip.
z, shape_z);
8524 for (
int o = 0, k = 0; k <= p; k++)
8525 for (
int j = 0; j <= p; j++)
8526 for (
int i = 0; i <= p; i++)
8528 shape(o++) = shape_x(i)*shape_y(j)*shape_z(k);
8535 const int p =
Order;
8537 #ifdef MFEM_THREAD_SAFE
8538 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
8539 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
8542 basis1d->
Eval(ip.
x, shape_x, dshape_x);
8543 basis1d->
Eval(ip.
y, shape_y, dshape_y);
8544 basis1d->
Eval(ip.
z, shape_z, dshape_z);
8546 for (
int o = 0, k = 0; k <= p; k++)
8547 for (
int j = 0; j <= p; j++)
8548 for (
int i = 0; i <= p; i++)
8550 dshape(o,0) = dshape_x(i)* shape_y(j)* shape_z(k);
8551 dshape(o,1) = shape_x(i)*dshape_y(j)* shape_z(k);
8552 dshape(o,2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
8558 const int p =
Order;
8568 #ifdef MFEM_THREAD_SAFE
8569 Vector shape_x(p+1), shape_y(p+1);
8572 for (
int i = 0; i <= p; i++)
8581 for (
int o = 0, k = 0; k <= p; k++)
8582 for (
int j = 0; j <= p; j++)
8583 for (
int i = 0; i <= p; i++)
8585 dofs[o++] = shape_x(i)*shape_x(j)*shape_x(k);
8589 for (
int o = 0, k = 0; k <= p; k++)
8590 for (
int j = 0; j <= p; j++)
8591 for (
int i = 0; i <= p; i++)
8593 dofs[o++] = shape_y(i)*shape_x(j)*shape_x(k);
8597 for (
int o = 0, k = 0; k <= p; k++)
8598 for (
int j = 0; j <= p; j++)
8599 for (
int i = 0; i <= p; i++)
8601 dofs[o++] = shape_y(i)*shape_y(j)*shape_x(k);
8605 for (
int o = 0, k = 0; k <= p; k++)
8606 for (
int j = 0; j <= p; j++)
8607 for (
int i = 0; i <= p; i++)
8609 dofs[o++] = shape_x(i)*shape_y(j)*shape_x(k);
8613 for (
int o = 0, k = 0; k <= p; k++)
8614 for (
int j = 0; j <= p; j++)
8615 for (
int i = 0; i <= p; i++)
8617 dofs[o++] = shape_x(i)*shape_x(j)*shape_y(k);
8621 for (
int o = 0, k = 0; k <= p; k++)
8622 for (
int j = 0; j <= p; j++)
8623 for (
int i = 0; i <= p; i++)
8625 dofs[o++] = shape_y(i)*shape_x(j)*shape_y(k);
8629 for (
int o = 0, k = 0; k <= p; k++)
8630 for (
int j = 0; j <= p; j++)
8631 for (
int i = 0; i <= p; i++)
8633 dofs[o++] = shape_y(i)*shape_y(j)*shape_y(k);
8637 for (
int o = 0, k = 0; k <= p; k++)
8638 for (
int j = 0; j <= p; j++)
8639 for (
int i = 0; i <= p; i++)
8641 dofs[o++] = shape_x(i)*shape_y(j)*shape_y(k);
8652 #ifndef MFEM_THREAD_SAFE
8667 for (
int o = 0, k = 0; k <= p; k++)
8668 for (
int j = 0; j <= p; j++)
8669 for (
int i = 0; i <= p; i++)
8679 const int p =
Order;
8681 #ifdef MFEM_THREAD_SAFE
8682 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
8689 for (
int o = 0, k = 0; k <= p; k++)
8690 for (
int j = 0; j <= p; j++)
8691 for (
int i = 0; i <= p; i++)
8693 shape(o++) = shape_x(i)*shape_y(j)*shape_z(k);
8700 const int p =
Order;
8702 #ifdef MFEM_THREAD_SAFE
8703 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
8704 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
8711 for (
int o = 0, k = 0; k <= p; k++)
8712 for (
int j = 0; j <= p; j++)
8713 for (
int i = 0; i <= p; i++)
8715 dshape(o,0) = dshape_x(i)* shape_y(j)* shape_z(k);
8716 dshape(o,1) = shape_x(i)*dshape_y(j)* shape_z(k);
8717 dshape(o,2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
8723 const int p =
Order;
8728 case 0: dofs[0] = 1.0;
break;
8729 case 1: dofs[p] = 1.0;
break;
8730 case 2: dofs[p*(p + 2)] = 1.0;
break;
8731 case 3: dofs[p*(p + 1)] = 1.0;
break;
8732 case 4: dofs[p*(p + 1)*(p + 1)] = 1.0;
break;
8733 case 5: dofs[p + p*(p + 1)*(p + 1)] = 1.0;
break;
8734 case 6: dofs[
Dof - 1] = 1.0;
break;
8735 case 7: dofs[
Dof - p - 1] = 1.0;
break;
8754 #ifndef MFEM_THREAD_SAFE
8764 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
8767 for (
int o = 0, j = 0; j <= p; j++)
8768 for (
int i = 0; i + j <= p; i++)
8770 double w = op[i] + op[j] + op[p-i-j];
8775 for (
int k = 0; k <
Dof; k++)
8782 for (
int o = 0, j = 0; j <= p; j++)
8783 for (
int i = 0; i + j <= p; i++)
8785 T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
8796 const int p =
Order;
8798 #ifdef MFEM_THREAD_SAFE
8799 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(
Dof);
8806 for (
int o = 0, j = 0; j <= p; j++)
8807 for (
int i = 0; i + j <= p; i++)
8809 u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
8818 const int p =
Order;
8820 #ifdef MFEM_THREAD_SAFE
8821 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
8822 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
8830 for (
int o = 0, j = 0; j <= p; j++)
8831 for (
int i = 0; i + j <= p; i++)
8834 du(o,0) = ((dshape_x(i)* shape_l(k)) -
8835 ( shape_x(i)*dshape_l(k)))*shape_y(j);
8836 du(o,1) = ((dshape_y(j)* shape_l(k)) -
8837 ( shape_y(j)*dshape_l(k)))*shape_x(i);
8841 Ti.
Mult(du, dshape);
8849 for (
int i = 0; i <
Dof; i++)
8852 dofs[i] = pow(1.0 - ip.
x - ip.
y,
Order);
8856 for (
int i = 0; i <
Dof; i++)
8859 dofs[i] = pow(ip.
x,
Order);
8863 for (
int i = 0; i <
Dof; i++)
8866 dofs[i] = pow(ip.
y,
Order);
8877 #ifndef MFEM_THREAD_SAFE
8887 for (
int o = 0, j = 0; j <= p; j++)
8888 for (
int i = 0; i + j <= p; i++)
8904 #ifdef MFEM_THREAD_SAFE
8917 case 0: dofs[0] = 1.0;
break;
8918 case 1: dofs[
Order] = 1.0;
break;
8919 case 2: dofs[
Dof-1] = 1.0;
break;
8938 #ifndef MFEM_THREAD_SAFE
8950 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8953 for (
int o = 0, k = 0; k <= p; k++)
8954 for (
int j = 0; j + k <= p; j++)
8955 for (
int i = 0; i + j + k <= p; i++)
8957 double w = op[i] + op[j] + op[k] + op[p-i-j-k];
8962 for (
int m = 0; m <
Dof; m++)
8970 for (
int o = 0, k = 0; k <= p; k++)
8971 for (
int j = 0; j + k <= p; j++)
8972 for (
int i = 0; i + j + k <= p; i++)
8974 T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
8985 const int p =
Order;
8987 #ifdef MFEM_THREAD_SAFE
8988 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8997 for (
int o = 0, k = 0; k <= p; k++)
8998 for (
int j = 0; j + k <= p; j++)
8999 for (
int i = 0; i + j + k <= p; i++)
9001 u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
9010 const int p =
Order;
9012 #ifdef MFEM_THREAD_SAFE
9013 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
9014 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
9023 for (
int o = 0, k = 0; k <= p; k++)
9024 for (
int j = 0; j + k <= p; j++)
9025 for (
int i = 0; i + j + k <= p; i++)
9027 int l = p - i - j - k;
9028 du(o,0) = ((dshape_x(i)* shape_l(l)) -
9029 ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
9030 du(o,1) = ((dshape_y(j)* shape_l(l)) -
9031 ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
9032 du(o,2) = ((dshape_z(k)* shape_l(l)) -
9033 ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
9037 Ti.
Mult(du, dshape);
9045 for (
int i = 0; i <
Dof; i++)
9048 dofs[i] = pow(1.0 - ip.
x - ip.
y - ip.
z,
Order);
9052 for (
int i = 0; i <
Dof; i++)
9055 dofs[i] = pow(ip.
x,
Order);
9059 for (
int i = 0; i <
Dof; i++)
9062 dofs[i] = pow(ip.
y,
Order);
9065 for (
int i = 0; i <
Dof; i++)
9068 dofs[i] = pow(ip.
z,
Order);
9079 #ifndef MFEM_THREAD_SAFE
9089 for (
int o = 0, k = 0; k <= p; k++)
9090 for (
int j = 0; j + k <= p; j++)
9091 for (
int i = 0; i + j + k <= p; i++)
9108 #ifdef MFEM_THREAD_SAFE
9121 case 0: dofs[0] = 1.0;
break;
9122 case 1: dofs[
Order] = 1.0;
break;
9123 case 2: dofs[(
Order*(
Order+3))/2] = 1.0;
break;
9124 case 3: dofs[
Dof-1] = 1.0;
break;
9129 const double RT_QuadrilateralElement::nk[8] =
9130 { 0., -1., 1., 0., 0., 1., -1., 0. };
9135 cbasis1d(
poly1d.ClosedBasis(p + 1)), obasis1d(
poly1d.OpenBasis(p)),
9136 dof_map(Dof), dof2nk(Dof)
9140 const int dof2 =
Dof/2;
9142 #ifndef MFEM_THREAD_SAFE
9153 for (
int i = 0; i <= p; i++)
9155 dof_map[1*dof2 + i + 0*(p + 1)] = o++;
9157 for (
int i = 0; i <= p; i++)
9159 dof_map[0*dof2 + (p + 1) + i*(p + 2)] = o++;
9161 for (
int i = 0; i <= p; i++)
9163 dof_map[1*dof2 + (p - i) + (p + 1)*(p + 1)] = o++;
9165 for (
int i = 0; i <= p; i++)
9167 dof_map[0*dof2 + 0 + (p - i)*(p + 2)] = o++;
9171 for (
int j = 0; j <= p; j++)
9172 for (
int i = 1; i <= p; i++)
9174 dof_map[0*dof2 + i + j*(p + 2)] = o++;
9176 for (
int j = 1; j <= p; j++)
9177 for (
int i = 0; i <= p; i++)
9179 dof_map[1*dof2 + i + j*(p + 1)] = o++;
9184 for (
int j = 0; j <= p; j++)
9185 for (
int i = 0; i <= p/2; i++)
9187 int idx = 0*dof2 + i + j*(p + 2);
9188 dof_map[idx] = -1 - dof_map[idx];
9191 for (
int j = p/2 + 1; j <= p; j++)
9193 int idx = 0*dof2 + (p/2 + 1) + j*(p + 2);
9194 dof_map[idx] = -1 - dof_map[idx];
9197 for (
int j = 0; j <= p/2; j++)
9198 for (
int i = 0; i <= p; i++)
9200 int idx = 1*dof2 + i + j*(p + 1);
9201 dof_map[idx] = -1 - dof_map[idx];
9204 for (
int i = 0; i <= p/2; i++)
9206 int idx = 1*dof2 + i + (p/2 + 1)*(p + 1);
9207 dof_map[idx] = -1 - dof_map[idx];
9211 for (
int j = 0; j <= p; j++)
9212 for (
int i = 0; i <= p + 1; i++)
9215 if ((idx = dof_map[o++]) < 0)
9226 for (
int j = 0; j <= p + 1; j++)
9227 for (
int i = 0; i <= p; i++)
9230 if ((idx = dof_map[o++]) < 0)
9246 const int pp1 =
Order;
9248 #ifdef MFEM_THREAD_SAFE
9249 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9252 cbasis1d.
Eval(ip.
x, shape_cx);
9253 obasis1d.
Eval(ip.
x, shape_ox);
9254 cbasis1d.
Eval(ip.
y, shape_cy);
9255 obasis1d.
Eval(ip.
y, shape_oy);
9258 for (
int j = 0; j < pp1; j++)
9259 for (
int i = 0; i <= pp1; i++)
9262 if ((idx = dof_map[o++]) < 0)
9264 idx = -1 - idx, s = -1;
9270 shape(idx,0) = s*shape_cx(i)*shape_oy(j);
9273 for (
int j = 0; j <= pp1; j++)
9274 for (
int i = 0; i < pp1; i++)
9277 if ((idx = dof_map[o++]) < 0)
9279 idx = -1 - idx, s = -1;
9286 shape(idx,1) = s*shape_ox(i)*shape_cy(j);
9293 const int pp1 =
Order;
9295 #ifdef MFEM_THREAD_SAFE
9296 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9297 Vector dshape_cx(pp1 + 1), dshape_cy(pp1 + 1);
9300 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
9301 obasis1d.
Eval(ip.
x, shape_ox);
9302 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
9303 obasis1d.
Eval(ip.
y, shape_oy);
9306 for (
int j = 0; j < pp1; j++)
9307 for (
int i = 0; i <= pp1; i++)
9310 if ((idx = dof_map[o++]) < 0)
9312 idx = -1 - idx, s = -1;
9318 divshape(idx) = s*dshape_cx(i)*shape_oy(j);
9320 for (
int j = 0; j <= pp1; j++)
9321 for (
int i = 0; i < pp1; i++)
9324 if ((idx = dof_map[o++]) < 0)
9326 idx = -1 - idx, s = -1;
9332 divshape(idx) = s*shape_ox(i)*dshape_cy(j);
9337 const double RT_HexahedronElement::nk[18] =
9338 { 0.,0.,-1., 0.,-1.,0., 1.,0.,0., 0.,1.,0., -1.,0.,0., 0.,0.,1. };
9343 cbasis1d(
poly1d.ClosedBasis(p + 1)), obasis1d(
poly1d.OpenBasis(p)),
9344 dof_map(Dof), dof2nk(Dof)
9348 const int dof3 =
Dof/3;
9350 #ifndef MFEM_THREAD_SAFE
9364 for (
int j = 0; j <= p; j++)
9365 for (
int i = 0; i <= p; i++)
9367 dof_map[2*dof3 + i + ((p - j) + 0*(p + 1))*(p + 1)] = o++;
9369 for (
int j = 0; j <= p; j++)
9370 for (
int i = 0; i <= p; i++)
9372 dof_map[1*dof3 + i + (0 + j*(p + 2))*(p + 1)] = o++;
9374 for (
int j = 0; j <= p; j++)
9375 for (
int i = 0; i <= p; i++)
9377 dof_map[0*dof3 + (p + 1) + (i + j*(p + 1))*(p + 2)] = o++;
9379 for (
int j = 0; j <= p; j++)
9380 for (
int i = 0; i <= p; i++)
9382 dof_map[1*dof3 + (p - i) + ((p + 1) + j*(p + 2))*(p + 1)] = o++;
9384 for (
int j = 0; j <= p; j++)
9385 for (
int i = 0; i <= p; i++)
9387 dof_map[0*dof3 + 0 + ((p - i) + j*(p + 1))*(p + 2)] = o++;
9389 for (
int j = 0; j <= p; j++)
9390 for (
int i = 0; i <= p; i++)
9392 dof_map[2*dof3 + i + (j + (p + 1)*(p + 1))*(p + 1)] = o++;
9397 for (
int k = 0; k <= p; k++)
9398 for (
int j = 0; j <= p; j++)
9399 for (
int i = 1; i <= p; i++)
9401 dof_map[0*dof3 + i + (j + k*(p + 1))*(p + 2)] = o++;
9404 for (
int k = 0; k <= p; k++)
9405 for (
int j = 1; j <= p; j++)
9406 for (
int i = 0; i <= p; i++)
9408 dof_map[1*dof3 + i + (j + k*(p + 2))*(p + 1)] = o++;
9411 for (
int k = 1; k <= p; k++)
9412 for (
int j = 0; j <= p; j++)
9413 for (
int i = 0; i <= p; i++)
9415 dof_map[2*dof3 + i + (j + k*(p + 1))*(p + 1)] = o++;
9423 for (
int k = 0; k <= p; k++)
9424 for (
int j = 0; j <= p; j++)
9425 for (
int i = 0; i <= p/2; i++)
9427 int idx = 0*dof3 + i + (j + k*(p + 1))*(p + 2);
9428 dof_map[idx] = -1 - dof_map[idx];
9431 for (
int k = 0; k <= p; k++)
9432 for (
int j = 0; j <= p/2; j++)
9433 for (
int i = 0; i <= p; i++)
9435 int idx = 1*dof3 + i + (j + k*(p + 2))*(p + 1);
9436 dof_map[idx] = -1 - dof_map[idx];
9439 for (
int k = 0; k <= p/2; k++)
9440 for (
int j = 0; j <= p; j++)
9441 for (
int i = 0; i <= p; i++)
9443 int idx = 2*dof3 + i + (j + k*(p + 1))*(p + 1);
9444 dof_map[idx] = -1 - dof_map[idx];
9449 for (
int k = 0; k <= p; k++)
9450 for (
int j = 0; j <= p; j++)
9451 for (
int i = 0; i <= p + 1; i++)
9454 if ((idx = dof_map[o++]) < 0)
9466 for (
int k = 0; k <= p; k++)
9467 for (
int j = 0; j <= p + 1; j++)
9468 for (
int i = 0; i <= p; i++)
9471 if ((idx = dof_map[o++]) < 0)
9483 for (
int k = 0; k <= p + 1; k++)
9484 for (
int j = 0; j <= p; j++)
9485 for (
int i = 0; i <= p; i++)
9488 if ((idx = dof_map[o++]) < 0)
9504 const int pp1 =
Order;
9506 #ifdef MFEM_THREAD_SAFE
9507 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9508 Vector shape_cz(pp1 + 1), shape_oz(pp1);
9511 cbasis1d.
Eval(ip.
x, shape_cx);
9512 obasis1d.
Eval(ip.
x, shape_ox);
9513 cbasis1d.
Eval(ip.
y, shape_cy);
9514 obasis1d.
Eval(ip.
y, shape_oy);
9515 cbasis1d.
Eval(ip.
z, shape_cz);
9516 obasis1d.
Eval(ip.
z, shape_oz);
9520 for (
int k = 0; k < pp1; k++)
9521 for (
int j = 0; j < pp1; j++)
9522 for (
int i = 0; i <= pp1; i++)
9525 if ((idx = dof_map[o++]) < 0)
9527 idx = -1 - idx, s = -1;
9533 shape(idx,0) = s*shape_cx(i)*shape_oy(j)*shape_oz(k);
9538 for (
int k = 0; k < pp1; k++)
9539 for (
int j = 0; j <= pp1; j++)
9540 for (
int i = 0; i < pp1; i++)
9543 if ((idx = dof_map[o++]) < 0)
9545 idx = -1 - idx, s = -1;
9552 shape(idx,1) = s*shape_ox(i)*shape_cy(j)*shape_oz(k);
9556 for (
int k = 0; k <= pp1; k++)
9557 for (
int j = 0; j < pp1; j++)
9558 for (
int i = 0; i < pp1; i++)
9561 if ((idx = dof_map[o++]) < 0)
9563 idx = -1 - idx, s = -1;
9571 shape(idx,2) = s*shape_ox(i)*shape_oy(j)*shape_cz(k);
9578 const int pp1 =
Order;
9580 #ifdef MFEM_THREAD_SAFE
9581 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9582 Vector shape_cz(pp1 + 1), shape_oz(pp1);
9583 Vector dshape_cx(pp1 + 1), dshape_cy(pp1 + 1), dshape_cz(pp1 + 1);
9586 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
9587 obasis1d.
Eval(ip.
x, shape_ox);
9588 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
9589 obasis1d.
Eval(ip.
y, shape_oy);
9590 cbasis1d.
Eval(ip.
z, shape_cz, dshape_cz);
9591 obasis1d.
Eval(ip.
z, shape_oz);
9595 for (
int k = 0; k < pp1; k++)
9596 for (
int j = 0; j < pp1; j++)
9597 for (
int i = 0; i <= pp1; i++)
9600 if ((idx = dof_map[o++]) < 0)
9602 idx = -1 - idx, s = -1;
9608 divshape(idx) = s*dshape_cx(i)*shape_oy(j)*shape_oz(k);
9611 for (
int k = 0; k < pp1; k++)
9612 for (
int j = 0; j <= pp1; j++)
9613 for (
int i = 0; i < pp1; i++)
9616 if ((idx = dof_map[o++]) < 0)
9618 idx = -1 - idx, s = -1;
9624 divshape(idx) = s*shape_ox(i)*dshape_cy(j)*shape_oz(k);
9627 for (
int k = 0; k <= pp1; k++)
9628 for (
int j = 0; j < pp1; j++)
9629 for (
int i = 0; i < pp1; i++)
9632 if ((idx = dof_map[o++]) < 0)
9634 idx = -1 - idx, s = -1;
9640 divshape(idx) = s*shape_ox(i)*shape_oy(j)*dshape_cz(k);
9645 const double RT_TriangleElement::nk[6] =
9646 { 0., -1., 1., 1., -1., 0. };
9648 const double RT_TriangleElement::c = 1./3.;
9657 #ifndef MFEM_THREAD_SAFE
9667 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
9672 for (
int i = 0; i <= p; i++)
9677 for (
int i = 0; i <= p; i++)
9682 for (
int i = 0; i <= p; i++)
9689 for (
int j = 0; j < p; j++)
9690 for (
int i = 0; i + j < p; i++)
9692 double w = iop[i] + iop[j] + iop[p-1-i-j];
9700 for (
int k = 0; k <
Dof; k++)
9706 const double *n_k = nk + 2*dof2nk[k];
9709 for (
int j = 0; j <= p; j++)
9710 for (
int i = 0; i + j <= p; i++)
9712 double s = shape_x(i)*shape_y(j)*shape_l(p-i-j);
9713 T(o++, k) = s*n_k[0];
9714 T(o++, k) = s*n_k[1];
9716 for (
int i = 0; i <= p; i++)
9718 double s = shape_x(i)*shape_y(p-i);
9719 T(o++, k) = s*((ip.
x - c)*n_k[0] + (ip.
y - c)*n_k[1]);
9730 const int p =
Order - 1;
9732 #ifdef MFEM_THREAD_SAFE
9733 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
9742 for (
int j = 0; j <= p; j++)
9743 for (
int i = 0; i + j <= p; i++)
9745 double s = shape_x(i)*shape_y(j)*shape_l(p-i-j);
9746 u(o,0) = s; u(o,1) = 0; o++;
9747 u(o,0) = 0; u(o,1) = s; o++;
9749 for (
int i = 0; i <= p; i++)
9751 double s = shape_x(i)*shape_y(p-i);
9752 u(o,0) = (ip.
x - c)*s;
9753 u(o,1) = (ip.
y - c)*s;
9763 const int p =
Order - 1;
9765 #ifdef MFEM_THREAD_SAFE
9766 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
9767 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
9776 for (
int j = 0; j <= p; j++)
9777 for (
int i = 0; i + j <= p; i++)
9780 divu(o++) = (dshape_x(i)*shape_l(k) -
9781 shape_x(i)*dshape_l(k))*shape_y(j);
9782 divu(o++) = (dshape_y(j)*shape_l(k) -
9783 shape_y(j)*dshape_l(k))*shape_x(i);
9785 for (
int i = 0; i <= p; i++)
9788 divu(o++) = ((shape_x(i) + (ip.
x - c)*dshape_x(i))*shape_y(j) +
9789 (shape_y(j) + (ip.
y - c)*dshape_y(j))*shape_x(i));
9792 Ti.
Mult(divu, divshape);
9796 const double RT_TetrahedronElement::nk[12] =
9797 { 1,1,1, -1,0,0, 0,-1,0, 0,0,-1 };
9800 const double RT_TetrahedronElement::c = 1./4.;
9809 #ifndef MFEM_THREAD_SAFE
9821 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
9827 for (
int j = 0; j <= p; j++)
9828 for (
int i = 0; i + j <= p; i++)
9830 double w = bop[i] + bop[j] + bop[p-i-j];
9834 for (
int j = 0; j <= p; j++)
9835 for (
int i = 0; i + j <= p; i++)
9837 double w = bop[i] + bop[j] + bop[p-i-j];
9841 for (
int j = 0; j <= p; j++)
9842 for (
int i = 0; i + j <= p; i++)
9844 double w = bop[i] + bop[j] + bop[p-i-j];
9848 for (
int j = 0; j <= p; j++)
9849 for (
int i = 0; i + j <= p; i++)
9851 double w = bop[i] + bop[j] + bop[p-i-j];
9857 for (
int k = 0; k < p; k++)
9858 for (
int j = 0; j + k < p; j++)
9859 for (
int i = 0; i + j + k < p; i++)
9861 double w = iop[i] + iop[j] + iop[k] + iop[p-1-i-j-k];
9871 for (
int m = 0; m <
Dof; m++)
9878 const double *nm = nk + 3*dof2nk[m];
9881 for (
int k = 0; k <= p; k++)
9882 for (
int j = 0; j + k <= p; j++)
9883 for (
int i = 0; i + j + k <= p; i++)
9885 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
9886 T(o++, m) = s * nm[0];
9887 T(o++, m) = s * nm[1];
9888 T(o++, m) = s * nm[2];
9890 for (
int j = 0; j <= p; j++)
9891 for (
int i = 0; i + j <= p; i++)
9893 double s = shape_x(i)*shape_y(j)*shape_z(p-i-j);
9894 T(o++, m) = s*((ip.
x - c)*nm[0] + (ip.
y - c)*nm[1] +
9906 const int p =
Order - 1;
9908 #ifdef MFEM_THREAD_SAFE
9909 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
9919 for (
int k = 0; k <= p; k++)
9920 for (
int j = 0; j + k <= p; j++)
9921 for (
int i = 0; i + j + k <= p; i++)
9923 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
9924 u(o,0) = s; u(o,1) = 0; u(o,2) = 0; o++;
9925 u(o,0) = 0; u(o,1) = s; u(o,2) = 0; o++;
9926 u(o,0) = 0; u(o,1) = 0; u(o,2) = s; o++;
9928 for (
int j = 0; j <= p; j++)
9929 for (
int i = 0; i + j <= p; i++)
9931 double s = shape_x(i)*shape_y(j)*shape_z(p-i-j);
9932 u(o,0) = (ip.
x - c)*s; u(o,1) = (ip.
y - c)*s; u(o,2) = (ip.
z - c)*s;
9942 const int p =
Order - 1;
9944 #ifdef MFEM_THREAD_SAFE
9945 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
9946 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
9956 for (
int k = 0; k <= p; k++)
9957 for (
int j = 0; j + k <= p; j++)
9958 for (
int i = 0; i + j + k <= p; i++)
9960 int l = p - i - j - k;
9961 divu(o++) = (dshape_x(i)*shape_l(l) -
9962 shape_x(i)*dshape_l(l))*shape_y(j)*shape_z(k);
9963 divu(o++) = (dshape_y(j)*shape_l(l) -
9964 shape_y(j)*dshape_l(l))*shape_x(i)*shape_z(k);
9965 divu(o++) = (dshape_z(k)*shape_l(l) -
9966 shape_z(k)*dshape_l(l))*shape_x(i)*shape_y(j);
9968 for (
int j = 0; j <= p; j++)
9969 for (
int i = 0; i + j <= p; i++)
9973 (shape_x(i) + (ip.
x - c)*dshape_x(i))*shape_y(j)*shape_z(k) +
9974 (shape_y(j) + (ip.
y - c)*dshape_y(j))*shape_x(i)*shape_z(k) +
9975 (shape_z(k) + (ip.
z - c)*dshape_z(k))*shape_x(i)*shape_y(j);
9978 Ti.
Mult(divu, divshape);
9982 const double ND_HexahedronElement::tk[18] =
9983 { 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,0.,0., 0.,-1.,0., 0.,0.,-1. };
9988 cbasis1d(
poly1d.ClosedBasis(p)), obasis1d(
poly1d.OpenBasis(p - 1)),
9989 dof_map(Dof), dof2tk(Dof)
9993 const int dof3 =
Dof/3;
9995 #ifndef MFEM_THREAD_SAFE
10009 for (
int i = 0; i < p; i++)
10011 dof_map[0*dof3 + i + (0 + 0*(p + 1))*p] = o++;
10013 for (
int i = 0; i < p; i++)
10015 dof_map[1*dof3 + p + (i + 0*p)*(p + 1)] = o++;
10017 for (
int i = 0; i < p; i++)
10019 dof_map[0*dof3 + i + (p + 0*(p + 1))*p] = o++;
10021 for (
int i = 0; i < p; i++)
10023 dof_map[1*dof3 + 0 + (i + 0*p)*(p + 1)] = o++;
10025 for (
int i = 0; i < p; i++)
10027 dof_map[0*dof3 + i + (0 + p*(p + 1))*p] = o++;
10029 for (
int i = 0; i < p; i++)
10031 dof_map[1*dof3 + p + (i + p*p)*(p + 1)] = o++;
10033 for (
int i = 0; i < p; i++)
10035 dof_map[0*dof3 + i + (p + p*(p + 1))*p] = o++;
10037 for (
int i = 0; i < p; i++)
10039 dof_map[1*dof3 + 0 + (i + p*p)*(p + 1)] = o++;
10041 for (
int i = 0; i < p; i++)
10043 dof_map[2*dof3 + 0 + (0 + i*(p + 1))*(p + 1)] = o++;
10045 for (
int i = 0; i < p; i++)
10047 dof_map[2*dof3 + p + (0 + i*(p + 1))*(p + 1)] = o++;
10049 for (
int i = 0; i < p; i++)
10051 dof_map[2*dof3 + p + (p + i*(p + 1))*(p + 1)] = o++;
10053 for (
int i = 0; i < p; i++)
10055 dof_map[2*dof3 + 0 + (p + i*(p + 1))*(p + 1)] = o++;
10060 for (
int j = 1; j < p; j++)
10061 for (
int i = 0; i < p; i++)
10063 dof_map[0*dof3 + i + ((p - j) + 0*(p + 1))*p] = o++;
10065 for (
int j = 0; j < p; j++)
10066 for (
int i = 1; i < p; i++)
10068 dof_map[1*dof3 + i + ((p - 1 - j) + 0*p)*(p + 1)] = -1 - (o++);
10071 for (
int k = 1; k < p; k++)
10072 for (
int i = 0; i < p; i++)
10074 dof_map[0*dof3 + i + (0 + k*(p + 1))*p] = o++;
10076 for (
int k = 0; k < p; k++)
10077 for (
int i = 1; i < p; i++ )
10079 dof_map[2*dof3 + i + (0 + k*(p + 1))*(p + 1)] = o++;
10082 for (
int k = 1; k < p; k++)
10083 for (
int j = 0; j < p; j++)
10085 dof_map[1*dof3 + p + (j + k*p)*(p + 1)] = o++;
10087 for (
int k = 0; k < p; k++)
10088 for (
int j = 1; j < p; j++)
10090 dof_map[2*dof3 + p + (j + k*(p + 1))*(p + 1)] = o++;
10093 for (
int k = 1; k < p; k++)
10094 for (
int i = 0; i < p; i++)
10096 dof_map[0*dof3 + (p - 1 - i) + (p + k*(p + 1))*p] = -1 - (o++);
10098 for (
int k = 0; k < p; k++)
10099 for (
int i = 1; i < p; i++)
10101 dof_map[2*dof3 + (p - i) + (p + k*(p + 1))*(p + 1)] = o++;
10104 for (
int k = 1; k < p; k++)
10105 for (
int j = 0; j < p; j++)
10107 dof_map[1*dof3 + 0 + ((p - 1 - j) + k*p)*(p + 1)] = -1 - (o++);
10109 for (
int k = 0; k < p; k++)
10110 for (
int j = 1; j < p; j++)
10112 dof_map[2*dof3 + 0 + ((p - j) + k*(p + 1))*(p + 1)] = o++;
10115 for (
int j = 1; j < p; j++)
10116 for (
int i = 0; i < p; i++)
10118 dof_map[0*dof3 + i + (j + p*(p + 1))*p] = o++;
10120 for (
int j = 0; j < p; j++)
10121 for (
int i = 1; i < p; i++)
10123 dof_map[1*dof3 + i + (j + p*p)*(p + 1)] = o++;
10128 for (
int k = 1; k < p; k++)
10129 for (
int j = 1; j < p; j++)
10130 for (
int i = 0; i < p; i++)
10132 dof_map[0*dof3 + i + (j + k*(p + 1))*p] = o++;
10135 for (
int k = 1; k < p; k++)
10136 for (
int j = 0; j < p; j++)
10137 for (
int i = 1; i < p; i++)
10139 dof_map[1*dof3 + i + (j + k*p)*(p + 1)] = o++;
10142 for (
int k = 0; k < p; k++)
10143 for (
int j = 1; j < p; j++)
10144 for (
int i = 1; i < p; i++)
10146 dof_map[2*dof3 + i + (j + k*(p + 1))*(p + 1)] = o++;
10152 for (
int k = 0; k <= p; k++)
10153 for (
int j = 0; j <= p; j++)
10154 for (
int i = 0; i < p; i++)
10157 if ((idx = dof_map[o++]) < 0)
10159 dof2tk[idx = -1 - idx] = 3;
10168 for (
int k = 0; k <= p; k++)
10169 for (
int j = 0; j < p; j++)
10170 for (
int i = 0; i <= p; i++)
10173 if ((idx = dof_map[o++]) < 0)
10175 dof2tk[idx = -1 - idx] = 4;
10184 for (
int k = 0; k < p; k++)
10185 for (
int j = 0; j <= p; j++)
10186 for (
int i = 0; i <= p; i++)
10189 if ((idx = dof_map[o++]) < 0)
10191 dof2tk[idx = -1 - idx] = 5;
10204 const int p =
Order;
10206 #ifdef MFEM_THREAD_SAFE
10207 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10208 Vector shape_cz(p + 1), shape_oz(p);
10211 cbasis1d.
Eval(ip.
x, shape_cx);
10212 obasis1d.
Eval(ip.
x, shape_ox);
10213 cbasis1d.
Eval(ip.
y, shape_cy);
10214 obasis1d.
Eval(ip.
y, shape_oy);
10215 cbasis1d.
Eval(ip.
z, shape_cz);
10216 obasis1d.
Eval(ip.
z, shape_oz);
10220 for (
int k = 0; k <= p; k++)
10221 for (
int j = 0; j <= p; j++)
10222 for (
int i = 0; i < p; i++)
10225 if ((idx = dof_map[o++]) < 0)
10227 idx = -1 - idx, s = -1;
10233 shape(idx,0) = s*shape_ox(i)*shape_cy(j)*shape_cz(k);
10238 for (
int k = 0; k <= p; k++)
10239 for (
int j = 0; j < p; j++)
10240 for (
int i = 0; i <= p; i++)
10243 if ((idx = dof_map[o++]) < 0)
10245 idx = -1 - idx, s = -1;
10252 shape(idx,1) = s*shape_cx(i)*shape_oy(j)*shape_cz(k);
10256 for (
int k = 0; k < p; k++)
10257 for (
int j = 0; j <= p; j++)
10258 for (
int i = 0; i <= p; i++)
10261 if ((idx = dof_map[o++]) < 0)
10263 idx = -1 - idx, s = -1;
10271 shape(idx,2) = s*shape_cx(i)*shape_cy(j)*shape_oz(k);
10278 const int p =
Order;
10280 #ifdef MFEM_THREAD_SAFE
10281 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10282 Vector shape_cz(p + 1), shape_oz(p);
10283 Vector dshape_cx(p + 1), dshape_cy(p + 1), dshape_cz(p + 1);
10286 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
10287 obasis1d.
Eval(ip.
x, shape_ox);
10288 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
10289 obasis1d.
Eval(ip.
y, shape_oy);
10290 cbasis1d.
Eval(ip.
z, shape_cz, dshape_cz);
10291 obasis1d.
Eval(ip.
z, shape_oz);
10295 for (
int k = 0; k <= p; k++)
10296 for (
int j = 0; j <= p; j++)
10297 for (
int i = 0; i < p; i++)
10300 if ((idx = dof_map[o++]) < 0)
10302 idx = -1 - idx, s = -1;
10308 curl_shape(idx,0) = 0.;
10309 curl_shape(idx,1) = s*shape_ox(i)* shape_cy(j)*dshape_cz(k);
10310 curl_shape(idx,2) = -s*shape_ox(i)*dshape_cy(j)* shape_cz(k);
10313 for (
int k = 0; k <= p; k++)
10314 for (
int j = 0; j < p; j++)
10315 for (
int i = 0; i <= p; i++)
10318 if ((idx = dof_map[o++]) < 0)
10320 idx = -1 - idx, s = -1;
10326 curl_shape(idx,0) = -s* shape_cx(i)*shape_oy(j)*dshape_cz(k);
10327 curl_shape(idx,1) = 0.;
10328 curl_shape(idx,2) = s*dshape_cx(i)*shape_oy(j)* shape_cz(k);
10331 for (
int k = 0; k < p; k++)
10332 for (
int j = 0; j <= p; j++)
10333 for (
int i = 0; i <= p; i++)
10336 if ((idx = dof_map[o++]) < 0)
10338 idx = -1 - idx, s = -1;
10344 curl_shape(idx,0) = s* shape_cx(i)*dshape_cy(j)*shape_oz(k);
10345 curl_shape(idx,1) = -s*dshape_cx(i)* shape_cy(j)*shape_oz(k);
10346 curl_shape(idx,2) = 0.;
10351 const double ND_QuadrilateralElement::tk[8] =
10352 { 1.,0., 0.,1., -1.,0., 0.,-1. };
10357 cbasis1d(
poly1d.ClosedBasis(p)), obasis1d(
poly1d.OpenBasis(p - 1)),
10358 dof_map(Dof), dof2tk(Dof)
10362 const int dof2 =
Dof/2;
10364 #ifndef MFEM_THREAD_SAFE
10375 for (
int i = 0; i < p; i++)
10377 dof_map[0*dof2 + i + 0*p] = o++;
10379 for (
int j = 0; j < p; j++)
10381 dof_map[1*dof2 + p + j*(p + 1)] = o++;
10383 for (
int i = 0; i < p; i++)
10385 dof_map[0*dof2 + (p - 1 - i) + p*p] = -1 - (o++);
10387 for (
int j = 0; j < p; j++)
10389 dof_map[1*dof2 + 0 + (p - 1 - j)*(p + 1)] = -1 - (o++);
10394 for (
int j = 1; j < p; j++)
10395 for (
int i = 0; i < p; i++)
10397 dof_map[0*dof2 + i + j*p] = o++;
10400 for (
int j = 0; j < p; j++)
10401 for (
int i = 1; i < p; i++)
10403 dof_map[1*dof2 + i + j*(p + 1)] = o++;
10409 for (
int j = 0; j <= p; j++)
10410 for (
int i = 0; i < p; i++)
10413 if ((idx = dof_map[o++]) < 0)
10415 dof2tk[idx = -1 - idx] = 2;
10424 for (
int j = 0; j < p; j++)
10425 for (
int i = 0; i <= p; i++)
10428 if ((idx = dof_map[o++]) < 0)
10430 dof2tk[idx = -1 - idx] = 3;
10443 const int p =
Order;
10445 #ifdef MFEM_THREAD_SAFE
10446 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10449 cbasis1d.
Eval(ip.
x, shape_cx);
10450 obasis1d.
Eval(ip.
x, shape_ox);
10451 cbasis1d.
Eval(ip.
y, shape_cy);
10452 obasis1d.
Eval(ip.
y, shape_oy);
10456 for (
int j = 0; j <= p; j++)
10457 for (
int i = 0; i < p; i++)
10460 if ((idx = dof_map[o++]) < 0)
10462 idx = -1 - idx, s = -1;
10468 shape(idx,0) = s*shape_ox(i)*shape_cy(j);
10472 for (
int j = 0; j < p; j++)
10473 for (
int i = 0; i <= p; i++)
10476 if ((idx = dof_map[o++]) < 0)
10478 idx = -1 - idx, s = -1;
10485 shape(idx,1) = s*shape_cx(i)*shape_oy(j);
10492 const int p =
Order;
10494 #ifdef MFEM_THREAD_SAFE
10495 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10496 Vector dshape_cx(p + 1), dshape_cy(p + 1);
10499 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
10500 obasis1d.
Eval(ip.
x, shape_ox);
10501 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
10502 obasis1d.
Eval(ip.
y, shape_oy);
10506 for (
int j = 0; j <= p; j++)
10507 for (
int i = 0; i < p; i++)
10510 if ((idx = dof_map[o++]) < 0)
10512 idx = -1 - idx, s = -1;
10518 curl_shape(idx,0) = -s*shape_ox(i)*dshape_cy(j);
10521 for (
int j = 0; j < p; j++)
10522 for (
int i = 0; i <= p; i++)
10525 if ((idx = dof_map[o++]) < 0)
10527 idx = -1 - idx, s = -1;
10533 curl_shape(idx,0) = s*dshape_cx(i)*shape_oy(j);
10538 const double ND_TetrahedronElement::tk[18] =
10539 { 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,1.,0., -1.,0.,1., 0.,-1.,1. };
10541 const double ND_TetrahedronElement::c = 1./4.;
10551 const int pm1 = p - 1, pm2 = p - 2, pm3 = p - 3;
10553 #ifndef MFEM_THREAD_SAFE
10564 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
10569 for (
int i = 0; i < p; i++)
10574 for (
int i = 0; i < p; i++)
10579 for (
int i = 0; i < p; i++)
10584 for (
int i = 0; i < p; i++)
10589 for (
int i = 0; i < p; i++)
10594 for (
int i = 0; i < p; i++)
10601 for (
int j = 0; j <= pm2; j++)
10602 for (
int i = 0; i + j <= pm2; i++)
10604 double w = fop[i] + fop[j] + fop[pm2-i-j];
10610 for (
int j = 0; j <= pm2; j++)
10611 for (
int i = 0; i + j <= pm2; i++)
10613 double w = fop[i] + fop[j] + fop[pm2-i-j];
10619 for (
int j = 0; j <= pm2; j++)
10620 for (
int i = 0; i + j <= pm2; i++)
10622 double w = fop[i] + fop[j] + fop[pm2-i-j];
10628 for (
int j = 0; j <= pm2; j++)
10629 for (
int i = 0; i + j <= pm2; i++)
10631 double w = fop[i] + fop[j] + fop[pm2-i-j];
10639 for (
int k = 0; k <= pm3; k++)
10640 for (
int j = 0; j + k <= pm3; j++)
10641 for (
int i = 0; i + j + k <= pm3; i++)
10643 double w = iop[i] + iop[j] + iop[k] + iop[pm3-i-j-k];
10653 for (
int m = 0; m <
Dof; m++)
10656 const double *tm = tk + 3*dof2tk[m];
10664 for (
int k = 0; k <= pm1; k++)
10665 for (
int j = 0; j + k <= pm1; j++)
10666 for (
int i = 0; i + j + k <= pm1; i++)
10668 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
10669 T(o++, m) = s * tm[0];
10670 T(o++, m) = s * tm[1];
10671 T(o++, m) = s * tm[2];
10673 for (
int k = 0; k <= pm1; k++)
10674 for (
int j = 0; j + k <= pm1; j++)
10676 double s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
10677 T(o++, m) = s*((ip.
y - c)*tm[0] - (ip.
x - c)*tm[1]);
10678 T(o++, m) = s*((ip.
z - c)*tm[0] - (ip.
x - c)*tm[2]);
10680 for (
int k = 0; k <= pm1; k++)
10683 shape_y(pm1-k)*shape_z(k)*((ip.
z - c)*tm[1] - (ip.
y - c)*tm[2]);
10694 const int pm1 =
Order - 1;
10696 #ifdef MFEM_THREAD_SAFE
10697 const int p =
Order;
10698 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
10708 for (
int k = 0; k <= pm1; k++)
10709 for (
int j = 0; j + k <= pm1; j++)
10710 for (
int i = 0; i + j + k <= pm1; i++)
10712 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
10713 u(n,0) = s; u(n,1) = 0.; u(n,2) = 0.; n++;
10714 u(n,0) = 0.; u(n,1) = s; u(n,2) = 0.; n++;
10715 u(n,0) = 0.; u(n,1) = 0.; u(n,2) = s; n++;
10717 for (
int k = 0; k <= pm1; k++)
10718 for (
int j = 0; j + k <= pm1; j++)
10720 double s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
10721 u(n,0) = s*(ip.
y - c); u(n,1) = -s*(ip.
x - c); u(n,2) = 0.; n++;
10722 u(n,0) = s*(ip.
z - c); u(n,1) = 0.; u(n,2) = -s*(ip.
x - c); n++;
10724 for (
int k = 0; k <= pm1; k++)
10726 double s = shape_y(pm1-k)*shape_z(k);
10727 u(n,0) = 0.; u(n,1) = s*(ip.
z - c); u(n,2) = -s*(ip.
y - c); n++;
10736 const int pm1 =
Order - 1;
10738 #ifdef MFEM_THREAD_SAFE
10739 const int p =
Order;
10740 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
10741 Vector dshape_x(p), dshape_y(p), dshape_z(p), dshape_l(p);
10751 for (
int k = 0; k <= pm1; k++)
10752 for (
int j = 0; j + k <= pm1; j++)
10753 for (
int i = 0; i + j + k <= pm1; i++)
10756 const double dx = (dshape_x(i)*shape_l(l) -
10757 shape_x(i)*dshape_l(l))*shape_y(j)*shape_z(k);
10758 const double dy = (dshape_y(j)*shape_l(l) -
10759 shape_y(j)*dshape_l(l))*shape_x(i)*shape_z(k);
10760 const double dz = (dshape_z(k)*shape_l(l) -
10761 shape_z(k)*dshape_l(l))*shape_x(i)*shape_y(j);
10763 u(n,0) = 0.; u(n,1) = dz; u(n,2) = -dy; n++;
10764 u(n,0) = -dz; u(n,1) = 0.; u(n,2) = dx; n++;
10765 u(n,0) = dy; u(n,1) = -dx; u(n,2) = 0.; n++;
10767 for (
int k = 0; k <= pm1; k++)
10768 for (
int j = 0; j + k <= pm1; j++)
10770 int i = pm1 - j - k;
10773 u(n,0) = shape_x(i)*(ip.
x - c)*shape_y(j)*dshape_z(k);
10774 u(n,1) = shape_x(i)*shape_y(j)*(ip.
y - c)*dshape_z(k);
10776 -((dshape_x(i)*(ip.
x - c) + shape_x(i))*shape_y(j)*shape_z(k) +
10777 (dshape_y(j)*(ip.
y - c) + shape_y(j))*shape_x(i)*shape_z(k));
10780 u(n,0) = -shape_x(i)*(ip.
x - c)*dshape_y(j)*shape_z(k);
10781 u(n,1) = (shape_x(i)*shape_y(j)*(dshape_z(k)*(ip.
z - c) + shape_z(k)) +
10782 (dshape_x(i)*(ip.
x - c) + shape_x(i))*shape_y(j)*shape_z(k));
10783 u(n,2) = -shape_x(i)*dshape_y(j)*shape_z(k)*(ip.
z - c);
10786 for (
int k = 0; k <= pm1; k++)
10790 u(n,0) = -((dshape_y(j)*(ip.
y - c) + shape_y(j))*shape_z(k) +
10791 shape_y(j)*(dshape_z(k)*(ip.
z - c) + shape_z(k)));
10796 Ti.
Mult(u, curl_shape);
10800 const double ND_TriangleElement::tk[8] =
10801 { 1.,0., -1.,1., 0.,-1., 0.,1. };
10803 const double ND_TriangleElement::c = 1./3.;
10812 const int pm1 = p - 1, pm2 = p - 2;
10814 #ifndef MFEM_THREAD_SAFE
10824 Vector shape_x(p), shape_y(p), shape_l(p);
10829 for (
int i = 0; i < p; i++)
10834 for (
int i = 0; i < p; i++)
10839 for (
int i = 0; i < p; i++)
10846 for (
int j = 0; j <= pm2; j++)
10847 for (
int i = 0; i + j <= pm2; i++)
10849 double w = iop[i] + iop[j] + iop[pm2-i-j];
10857 for (
int m = 0; m <
Dof; m++)
10860 const double *tm = tk + 2*dof2tk[m];
10867 for (
int j = 0; j <= pm1; j++)
10868 for (
int i = 0; i + j <= pm1; i++)
10870 double s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
10871 T(n++, m) = s * tm[0];
10872 T(n++, m) = s * tm[1];
10874 for (
int j = 0; j <= pm1; j++)
10877 shape_x(pm1-j)*shape_y(j)*((ip.
y - c)*tm[0] - (ip.
x - c)*tm[1]);
10888 const int pm1 =
Order - 1;
10890 #ifdef MFEM_THREAD_SAFE
10891 const int p =
Order;
10892 Vector shape_x(p), shape_y(p), shape_l(p);
10901 for (
int j = 0; j <= pm1; j++)
10902 for (
int i = 0; i + j <= pm1; i++)
10904 double s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
10905 u(n,0) = s; u(n,1) = 0; n++;
10906 u(n,0) = 0; u(n,1) = s; n++;
10908 for (
int j = 0; j <= pm1; j++)
10910 double s = shape_x(pm1-j)*shape_y(j);
10911 u(n,0) = s*(ip.
y - c);
10912 u(n,1) = -s*(ip.
x - c);
10922 const int pm1 =
Order - 1;
10924 #ifdef MFEM_THREAD_SAFE
10925 const int p =
Order;
10926 Vector shape_x(p), shape_y(p), shape_l(p);
10927 Vector dshape_x(p), dshape_y(p), dshape_l(p);
10936 for (
int j = 0; j <= pm1; j++)
10937 for (
int i = 0; i + j <= pm1; i++)
10940 const double dx = (dshape_x(i)*shape_l(l) -
10941 shape_x(i)*dshape_l(l)) * shape_y(j);
10942 const double dy = (dshape_y(j)*shape_l(l) -
10943 shape_y(j)*dshape_l(l)) * shape_x(i);
10949 for (
int j = 0; j <= pm1; j++)
10953 curlu(n++) = -((dshape_x(i)*(ip.
x - c) + shape_x(i)) * shape_y(j) +
10954 (dshape_y(j)*(ip.
y - c) + shape_y(j)) * shape_x(i));
10958 Ti.
Mult(curlu, curl2d);
10962 const double ND_SegmentElement::tk[1] = { 1. };
10967 obasis1d(
poly1d.OpenBasis(p - 1)), dof2tk(Dof)
10972 for (
int i = 0; i < p; i++)
10991 kv[0]->CalcShape(shape,
ijk[0], ip.
x);
10994 for (
int i = 0; i <=
Order; i++)
10996 sum += (shape(i) *=
weights(i));
11008 kv[0]->CalcDShape(grad,
ijk[0], ip.
x);
11010 double sum = 0.0, dsum = 0.0;
11011 for (
int i = 0; i <=
Order; i++)
11014 dsum += ( grad(i) *=
weights(i));
11018 add(sum, grad, -dsum*sum*sum,
shape_x, grad);
11028 for (
int o = 0, j = 0; j <=
Order; j++)
11030 const double sy =
shape_y(j);
11031 for (
int i = 0; i <=
Order; i++, o++)
11043 double sum, dsum[2];
11051 sum = dsum[0] = dsum[1] = 0.0;
11052 for (
int o = 0, j = 0; j <=
Order; j++)
11055 for (
int i = 0; i <=
Order; i++, o++)
11065 dsum[0] *= sum*sum;
11066 dsum[1] *= sum*sum;
11068 for (
int o = 0; o <
Dof; o++)
11070 dshape(o,0) = dshape(o,0)*sum -
u(o)*dsum[0];
11071 dshape(o,1) = dshape(o,1)*sum -
u(o)*dsum[1];
11083 for (
int o = 0, k = 0; k <=
Order; k++)
11085 const double sz =
shape_z(k);
11086 for (
int j = 0; j <=
Order; j++)
11088 const double sy_sz =
shape_y(j)*sz;
11089 for (
int i = 0; i <=
Order; i++, o++)
11102 double sum, dsum[3];
11112 sum = dsum[0] = dsum[1] = dsum[2] = 0.0;
11113 for (
int o = 0, k = 0; k <=
Order; k++)
11116 for (
int j = 0; j <=
Order; j++)
11118 const double sy_sz =
shape_y(j)* sz;
11119 const double dsy_sz =
dshape_y(j)* sz;
11120 const double sy_dsz =
shape_y(j)*dsz;
11121 for (
int i = 0; i <=
Order; i++, o++)
11133 dsum[0] *= sum*sum;
11134 dsum[1] *= sum*sum;
11135 dsum[2] *= sum*sum;
11137 for (
int o = 0; o <
Dof; o++)
11139 dshape(o,0) = dshape(o,0)*sum -
u(o)*dsum[0];
11140 dshape(o,1) = dshape(o,1)*sum -
u(o)*dsum[1];
11141 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)
Resize 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
H1Pos_TetrahedronElement(const int p)
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
Data type dense matrix using column-major storage.
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
void ProjectCurl_ND(const double *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) 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
void Factor()
Factor the current DenseMatrix, *a.
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
static void CalcShape(const int p, const double x, const double y, const double z, double *shape)
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 CalcDShape(const int p, const double x, const double y, const double z, double *dshape_1d, double *dshape)
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
H1Pos_TriangleElement(const int p)
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 Mult(const Vector &x, Vector &y) const
Matrix vector multiplication with the inverse of dense matrix.
virtual void ProjectCurl(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
CrouzeixRaviartFiniteElement()
static void CalcShape(const int p, const double x, const double y, double *shape)
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 MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
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
ND_SegmentElement(const int p)
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)
static void CalcDShape(const int p, const double x, const double y, double *dshape_1d, double *dshape)
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
void ProjectCurl_2D(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) 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 CalcVShape(const IntegrationPoint &ip, DenseMatrix &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
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
void Mult(const double *x, double *y) const
Matrix vector multiplication.
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)
Change the size of the DenseMatrix to s x s.
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