16 #include "../mesh/nurbs.hpp"
35 #ifndef MFEM_THREAD_SAFE
43 mfem_error (
"FiniteElement::CalcVShape (ip, ...)\n"
44 " is not implemented for this class!");
50 mfem_error (
"FiniteElement::CalcVShape (trans, ...)\n"
51 " is not implemented for this class!");
57 mfem_error (
"FiniteElement::CalcDivShape (ip, ...)\n"
58 " is not implemented for this class!");
65 div_shape *= (1.0 / Trans.
Weight());
71 mfem_error (
"FiniteElement::CalcCurlShape (ip, ...)\n"
72 " is not implemented for this class!");
82 #ifdef MFEM_THREAD_SAFE
87 curl_shape *= (1.0 / Trans.
Weight());
93 curl_shape *= (1.0 / Trans.
Weight());
96 MFEM_ABORT(
"Invalid dimension, Dim = " <<
Dim);
102 mfem_error (
"FiniteElement::GetFaceDofs (...)");
108 mfem_error (
"FiniteElement::CalcHessian (...) is not overloaded !");
114 mfem_error (
"GetLocalInterpolation (...) is not overloaded !");
121 MFEM_ABORT(
"method is not overloaded !");
127 mfem_error (
"FiniteElement::Project (...) is not overloaded !");
133 mfem_error (
"FiniteElement::Project (...) (vector) is not overloaded !");
139 mfem_error(
"FiniteElement::ProjectMatrixCoefficient() is not overloaded !");
144 mfem_error(
"FiniteElement::ProjectDelta(...) is not implemented for "
151 mfem_error(
"FiniteElement::Project(...) (fe version) is not implemented "
152 "for this element!");
159 mfem_error(
"FiniteElement::ProjectGrad(...) is not implemented for "
167 mfem_error(
"FiniteElement::ProjectCurl(...) is not implemented for "
175 mfem_error(
"FiniteElement::ProjectDiv(...) is not implemented for "
193 #ifdef MFEM_THREAD_SAFE
209 #ifdef MFEM_THREAD_SAFE
216 for (
int i = 0; i < fine_fe.
Dof; i++)
221 for (
int j = 0; j <
Dof; j++)
222 if (fabs(I(i,j) =
c_shape(j)) < 1.0e-12)
247 Vector fine_shape(fs), coarse_shape(cs);
248 DenseMatrix fine_mass(fs), fine_coarse_mass(fs, cs);
265 fine_mass_inv.
Mult(fine_coarse_mass, I);
285 for (
int i = 0; i <
Dof; i++)
288 for (
int j = 0; j < fe.
GetDof(); j++)
290 curl(i,j) = curl_shape(j,0);
298 for (
int i = 0; i <
Dof; i++)
304 dofs(i) = coeff.
Eval (Trans, ip);
307 dofs(i) *= Trans.
Weight();
318 for (
int i = 0; i <
Dof; i++)
322 vc.
Eval (x, Trans, ip);
327 for (
int j = 0; j < x.Size(); j++)
329 dofs(Dof*j+i) = x(j);
341 for (
int k = 0; k <
Dof; k++)
346 for (
int r = 0; r < MQ.Height(); r++)
348 for (
int d = 0; d < MQ.Width(); d++)
350 dofs(k+Dof*(d+MQ.Width()*r)) = MQ(r,d);
366 for (
int k = 0; k <
Dof; k++)
369 for (
int j = 0; j < shape.Size(); j++)
371 I(k,j) = (fabs(shape(j)) < 1e-12) ? 0.0 : shape(j);
380 for (
int k = 0; k <
Dof; k++)
391 I(k+d*Dof,j) =
vshape(j,d);
407 for (
int k = 0; k <
Dof; k++)
413 Mult(dshape, Jinv, grad_k);
418 for (
int j = 0; j < grad_k.Height(); j++)
419 for (
int d = 0; d <
Dim; d++)
421 grad(k+d*Dof,j) = grad_k(j,d);
434 for (
int k = 0; k <
Dof; k++)
442 for (
int j = 0; j < div_shape.Size(); j++)
444 div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j)/detJ;
449 for (
int j = 0; j < div_shape.Size(); j++)
451 div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j);
461 for (
int i = 0; i <
Dof; i++)
465 dofs(i) = coeff.
Eval(Trans, ip);
491 pos_mass_inv.
Mult(mixed_mass, I);
496 void VectorFiniteElement::CalcShape (
499 mfem_error (
"Error: Cannot use scalar CalcShape(...) function with\n"
500 " VectorFiniteElements!");
503 void VectorFiniteElement::CalcDShape (
504 const IntegrationPoint &ip, DenseMatrix &dshape )
const
506 mfem_error (
"Error: Cannot use scalar CalcDShape(...) function with\n"
507 " VectorFiniteElements!");
539 MFEM_ABORT(
"Invalid dimension, Dim = " <<
Dim);
543 MFEM_ABORT(
"Invalid MapType = " <<
MapType);
551 #ifdef MFEM_THREAD_SAFE
556 shape *= (1.0 / Trans.
Weight());
563 #ifdef MFEM_THREAD_SAFE
576 MFEM_ASSERT(vc.
GetVDim() == sdim,
"");
578 const bool square_J = (
Dim == sdim);
580 for (
int k = 0; k <
Dof; k++)
586 if (!square_J) { dofs(k) /= Trans.
Weight(); }
597 MFEM_ASSERT(mc.
GetWidth() == sdim,
"");
598 const bool square_J = (
Dim == sdim);
600 Vector nk_phys(sdim), dofs_k(MQ.Height());
601 MFEM_ASSERT(dofs.
Size() ==
Dof*MQ.Height(),
"");
603 for (
int k = 0; k <
Dof; k++)
609 if (!square_J) { nk_phys /= T.
Weight(); }
610 MQ.Mult(nk_phys, dofs_k);
611 for (
int r = 0; r < MQ.Height(); r++)
613 dofs(k+Dof*r) = dofs_k(r);
629 for (
int k = 0; k <
Dof; k++)
638 double w = 1.0/Trans.
Weight();
639 for (
int d = 0; d <
Dim; d++)
645 for (
int j = 0; j < shape.Size(); j++)
652 for (
int d = 0; d < sdim; d++)
654 I(k,j+d*shape.Size()) = s*vk[d];
661 mfem_error(
"VectorFiniteElement::Project_RT (fe version)");
671 mfem_error(
"VectorFiniteElement::ProjectGrad_RT works only in 2D!");
679 for (
int k = 0; k <
Dof; k++)
682 tk[0] = nk[d2n[k]*
Dim+1];
683 tk[1] = -nk[d2n[k]*
Dim];
684 dshape.Mult(tk, grad_k);
685 for (
int j = 0; j < grad_k.Size(); j++)
687 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
696 #ifdef MFEM_THREAD_SAFE
709 for (
int k = 0; k <
Dof; k++)
716 J *= 1.0 / Trans.
Weight();
723 for (
int j = 0; j < curl_k.Size(); j++)
725 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
738 for (
int k = 0; k <
Dof; k++)
741 curl_shape.Mult(nk + d2n[k]*
Dim, curl_k);
742 for (
int j = 0; j < curl_k.Size(); j++)
744 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
756 for (
int k = 0; k <
Dof; k++)
773 MFEM_ASSERT(mc.
GetWidth() == sdim,
"");
775 Vector tk_phys(sdim), dofs_k(MQ.Height());
776 MFEM_ASSERT(dofs.
Size() ==
Dof*MQ.Height(),
"");
778 for (
int k = 0; k <
Dof; k++)
784 MQ.Mult(tk_phys, dofs_k);
785 for (
int r = 0; r < MQ.Height(); r++)
787 dofs(k+Dof*r) = dofs_k(r);
803 for (
int k = 0; k <
Dof; k++)
812 double w = 1.0/Trans.
Weight();
813 for (
int d = 0; d < sdim; d++)
819 for (
int j = 0; j < shape.Size(); j++)
826 for (
int d = 0; d < sdim; d++)
828 I(k, j + d*shape.Size()) = s*vk[d];
835 mfem_error(
"VectorFiniteElement::Project_ND (fe version)");
849 for (
int k = 0; k <
Dof; k++)
852 dshape.Mult(tk + d2t[k]*
Dim, grad_k);
853 for (
int j = 0; j < grad_k.Size(); j++)
855 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
869 #ifdef MFEM_THREAD_SAFE
879 for (
int k = 0; k <
Dof; k++)
890 for (
int i = 0; i <
Dim; i++)
892 Ikj +=
vshape(j, i) * vk[i];
894 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
906 #ifdef MFEM_THREAD_SAFE
916 for (
int k = 0; k <
Dof; k++)
927 for (
int i = 0; i <
Dim; i++)
929 Ikj +=
vshape(j, i) * vk[i];
931 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
965 shape(0) = 1. - ip.
x;
990 shape(0) = 1. - ip.
x - ip.
y;
998 dshape(0,0) = -1.; dshape(0,1) = -1.;
999 dshape(1,0) = 1.; dshape(1,1) = 0.;
1000 dshape(2,0) = 0.; dshape(2,1) = 1.;
1019 shape(0) = (1. - ip.
x) * (1. - ip.
y) ;
1020 shape(1) = ip.
x * (1. - ip.
y) ;
1021 shape(2) = ip.
x * ip.
y ;
1022 shape(3) = (1. - ip.
x) * ip.
y ;
1028 dshape(0,0) = -1. + ip.
y; dshape(0,1) = -1. + ip.
x ;
1029 dshape(1,0) = 1. - ip.
y; dshape(1,1) = -ip.
x ;
1030 dshape(2,0) = ip.
y ; dshape(2,1) = ip.
x ;
1031 dshape(3,0) = -ip.
y ; dshape(3,1) = 1. - ip.
x ;
1037 h(0,0) = 0.; h(0,1) = 1.; h(0,2) = 0.;
1038 h(1,0) = 0.; h(1,1) = -1.; h(1,2) = 0.;
1039 h(2,0) = 0.; h(2,1) = 1.; h(2,2) = 0.;
1040 h(3,0) = 0.; h(3,1) = -1.; h(3,2) = 0.;
1058 const double x = ip.
x, y = ip.
y;
1060 shape(0) = 5./3. - 2. * (x + y);
1061 shape(1) = 2. * (x - 1./6.);
1062 shape(2) = 2. * (y - 1./6.);
1068 dshape(0,0) = -2.; dshape(0,1) = -2.;
1069 dshape(1,0) = 2.; dshape(1,1) = 0.;
1070 dshape(2,0) = 0.; dshape(2,1) = 2.;
1075 dofs(vertex) = 2./3.;
1076 dofs((vertex+1)%3) = 1./6.;
1077 dofs((vertex+2)%3) = 1./6.;
1082 const double GaussBiLinear2DFiniteElement::p[] =
1083 { 0.2113248654051871177454256, 0.7886751345948128822545744 };
1101 const double x = ip.
x, y = ip.
y;
1103 shape(0) = 3. * (p[1] - x) * (p[1] - y);
1104 shape(1) = 3. * (x - p[0]) * (p[1] - y);
1105 shape(2) = 3. * (x - p[0]) * (y - p[0]);
1106 shape(3) = 3. * (p[1] - x) * (y - p[0]);
1112 const double x = ip.
x, y = ip.
y;
1114 dshape(0,0) = 3. * (y - p[1]); dshape(0,1) = 3. * (x - p[1]);
1115 dshape(1,0) = 3. * (p[1] - y); dshape(1,1) = 3. * (p[0] - x);
1116 dshape(2,0) = 3. * (y - p[0]); dshape(2,1) = 3. * (x - p[0]);
1117 dshape(3,0) = 3. * (p[0] - y); dshape(3,1) = 3. * (p[1] - x);
1123 dofs(vertex) = p[1]*p[1];
1124 dofs((vertex+1)%4) = p[0]*p[1];
1125 dofs((vertex+2)%4) = p[0]*p[0];
1126 dofs((vertex+3)%4) = p[0]*p[1];
1147 shape(0) = 1. - ip.
x - ip.
y;
1155 dshape(0,0) = -1.; dshape(0,1) = -1.;
1156 dshape(1,0) = 1.; dshape(1,1) = 0.;
1157 dshape(2,0) = 0.; dshape(2,1) = 1.;
1173 double l1 = 1.0 - x, l2 = x, l3 = 2. * x - 1.;
1175 shape(0) = l1 * (-l3);
1177 shape(2) = 4. * l1 * l2;
1185 dshape(0,0) = 4. * x - 3.;
1186 dshape(1,0) = 4. * x - 1.;
1187 dshape(2,0) = 4. - 8. * x;
1202 const double x = ip.
x, x1 = 1. - x;
1206 shape(2) = 2. * x * x1;
1212 const double x = ip.
x;
1214 dshape(0,0) = 2. * x - 2.;
1215 dshape(1,0) = 2. * x;
1216 dshape(2,0) = 2. - 4. * x;
1239 double x = ip.
x, y = ip.
y;
1240 double l1 = 1.-x-y, l2 = x, l3 = y;
1242 shape(0) = l1 * (2. * l1 - 1.);
1243 shape(1) = l2 * (2. * l2 - 1.);
1244 shape(2) = l3 * (2. * l3 - 1.);
1245 shape(3) = 4. * l1 * l2;
1246 shape(4) = 4. * l2 * l3;
1247 shape(5) = 4. * l3 * l1;
1253 double x = ip.
x, y = ip.
y;
1256 dshape(0,1) = 4. * (x + y) - 3.;
1258 dshape(1,0) = 4. * x - 1.;
1262 dshape(2,1) = 4. * y - 1.;
1264 dshape(3,0) = -4. * (2. * x + y - 1.);
1265 dshape(3,1) = -4. * x;
1267 dshape(4,0) = 4. * y;
1268 dshape(4,1) = 4. * x;
1270 dshape(5,0) = -4. * y;
1271 dshape(5,1) = -4. * (x + 2. * y - 1.);
1311 case 0: dofs(3) = 0.25; dofs(5) = 0.25;
break;
1312 case 1: dofs(3) = 0.25; dofs(4) = 0.25;
break;
1313 case 2: dofs(4) = 0.25; dofs(5) = 0.25;
break;
1319 const double GaussQuad2DFiniteElement::p[] =
1320 { 0.0915762135097707434595714634022015, 0.445948490915964886318329253883051 };
1338 for (
int i = 0; i < 6; i++)
1355 const double x = ip.
x, y = ip.
y;
1369 const double x = ip.
x, y = ip.
y;
1370 D(0,0) = 0.; D(0,1) = 0.;
1371 D(1,0) = 1.; D(1,1) = 0.;
1372 D(2,0) = 0.; D(2,1) = 1.;
1373 D(3,0) = 2. * x; D(3,1) = 0.;
1374 D(4,0) = y; D(4,1) = x;
1375 D(5,0) = 0.; D(5,1) = 2. * y;
1407 double x = ip.
x, y = ip.
y;
1408 double l1x, l2x, l3x, l1y, l2y, l3y;
1410 l1x = (x - 1.) * (2. * x - 1);
1411 l2x = 4. * x * (1. - x);
1412 l3x = x * (2. * x - 1.);
1413 l1y = (y - 1.) * (2. * y - 1);
1414 l2y = 4. * y * (1. - y);
1415 l3y = y * (2. * y - 1.);
1417 shape(0) = l1x * l1y;
1418 shape(4) = l2x * l1y;
1419 shape(1) = l3x * l1y;
1420 shape(7) = l1x * l2y;
1421 shape(8) = l2x * l2y;
1422 shape(5) = l3x * l2y;
1423 shape(3) = l1x * l3y;
1424 shape(6) = l2x * l3y;
1425 shape(2) = l3x * l3y;
1431 double x = ip.
x, y = ip.
y;
1432 double l1x, l2x, l3x, l1y, l2y, l3y;
1433 double d1x, d2x, d3x, d1y, d2y, d3y;
1435 l1x = (x - 1.) * (2. * x - 1);
1436 l2x = 4. * x * (1. - x);
1437 l3x = x * (2. * x - 1.);
1438 l1y = (y - 1.) * (2. * y - 1);
1439 l2y = 4. * y * (1. - y);
1440 l3y = y * (2. * y - 1.);
1449 dshape(0,0) = d1x * l1y;
1450 dshape(0,1) = l1x * d1y;
1452 dshape(4,0) = d2x * l1y;
1453 dshape(4,1) = l2x * d1y;
1455 dshape(1,0) = d3x * l1y;
1456 dshape(1,1) = l3x * d1y;
1458 dshape(7,0) = d1x * l2y;
1459 dshape(7,1) = l1x * d2y;
1461 dshape(8,0) = d2x * l2y;
1462 dshape(8,1) = l2x * d2y;
1464 dshape(5,0) = d3x * l2y;
1465 dshape(5,1) = l3x * d2y;
1467 dshape(3,0) = d1x * l3y;
1468 dshape(3,1) = l1x * d3y;
1470 dshape(6,0) = d2x * l3y;
1471 dshape(6,1) = l2x * d3y;
1473 dshape(2,0) = d3x * l3y;
1474 dshape(2,1) = l3x * d3y;
1486 case 0: dofs(4) = 0.25; dofs(7) = 0.25;
break;
1487 case 1: dofs(4) = 0.25; dofs(5) = 0.25;
break;
1488 case 2: dofs(5) = 0.25; dofs(6) = 0.25;
break;
1489 case 3: dofs(6) = 0.25; dofs(7) = 0.25;
break;
1521 double x = ip.
x, y = ip.
y;
1522 double l1x, l2x, l3x, l1y, l2y, l3y;
1524 l1x = (1. - x) * (1. - x);
1525 l2x = 2. * x * (1. - x);
1527 l1y = (1. - y) * (1. - y);
1528 l2y = 2. * y * (1. - y);
1531 shape(0) = l1x * l1y;
1532 shape(4) = l2x * l1y;
1533 shape(1) = l3x * l1y;
1534 shape(7) = l1x * l2y;
1535 shape(8) = l2x * l2y;
1536 shape(5) = l3x * l2y;
1537 shape(3) = l1x * l3y;
1538 shape(6) = l2x * l3y;
1539 shape(2) = l3x * l3y;
1545 double x = ip.
x, y = ip.
y;
1546 double l1x, l2x, l3x, l1y, l2y, l3y;
1547 double d1x, d2x, d3x, d1y, d2y, d3y;
1549 l1x = (1. - x) * (1. - x);
1550 l2x = 2. * x * (1. - x);
1552 l1y = (1. - y) * (1. - y);
1553 l2y = 2. * y * (1. - y);
1563 dshape(0,0) = d1x * l1y;
1564 dshape(0,1) = l1x * d1y;
1566 dshape(4,0) = d2x * l1y;
1567 dshape(4,1) = l2x * d1y;
1569 dshape(1,0) = d3x * l1y;
1570 dshape(1,1) = l3x * d1y;
1572 dshape(7,0) = d1x * l2y;
1573 dshape(7,1) = l1x * d2y;
1575 dshape(8,0) = d2x * l2y;
1576 dshape(8,1) = l2x * d2y;
1578 dshape(5,0) = d3x * l2y;
1579 dshape(5,1) = l3x * d2y;
1581 dshape(3,0) = d1x * l3y;
1582 dshape(3,1) = l1x * d3y;
1584 dshape(6,0) = d2x * l3y;
1585 dshape(6,1) = l2x * d3y;
1587 dshape(2,0) = d3x * l3y;
1588 dshape(2,1) = l3x * d3y;
1596 Vector xx(&tr_ip.
x, 2), shape(s, 9);
1598 for (
int i = 0; i < 9; i++)
1602 for (
int j = 0; j < 9; j++)
1603 if (fabs(I(i,j) = s[j]) < 1.0e-12)
1608 for (
int i = 0; i < 9; i++)
1610 double *d = &I(0,i);
1611 d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
1612 d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
1613 d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
1614 d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
1615 d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
1616 0.25 * (d[0] + d[1] + d[2] + d[3]);
1625 for (
int i = 0; i < 9; i++)
1629 d[i] = coeff.
Eval(Trans, ip);
1631 d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
1632 d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
1633 d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
1634 d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
1635 d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
1636 0.25 * (d[0] + d[1] + d[2] + d[3]);
1646 for (
int i = 0; i < 9; i++)
1650 vc.
Eval (x, Trans, ip);
1651 for (
int j = 0; j < x.Size(); j++)
1656 for (
int j = 0; j < x.Size(); j++)
1658 double *d = &dofs(9*j);
1660 d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
1661 d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
1662 d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
1663 d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
1664 d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
1665 0.25 * (d[0] + d[1] + d[2] + d[3]);
1673 const double p1 = 0.5*(1.-sqrt(3./5.));
1698 const double a = sqrt(5./3.);
1699 const double p1 = 0.5*(1.-sqrt(3./5.));
1701 double x = a*(ip.
x-p1), y = a*(ip.
y-p1);
1702 double l1x, l2x, l3x, l1y, l2y, l3y;
1704 l1x = (x - 1.) * (2. * x - 1);
1705 l2x = 4. * x * (1. - x);
1706 l3x = x * (2. * x - 1.);
1707 l1y = (y - 1.) * (2. * y - 1);
1708 l2y = 4. * y * (1. - y);
1709 l3y = y * (2. * y - 1.);
1711 shape(0) = l1x * l1y;
1712 shape(4) = l2x * l1y;
1713 shape(1) = l3x * l1y;
1714 shape(7) = l1x * l2y;
1715 shape(8) = l2x * l2y;
1716 shape(5) = l3x * l2y;
1717 shape(3) = l1x * l3y;
1718 shape(6) = l2x * l3y;
1719 shape(2) = l3x * l3y;
1725 const double a = sqrt(5./3.);
1726 const double p1 = 0.5*(1.-sqrt(3./5.));
1728 double x = a*(ip.
x-p1), y = a*(ip.
y-p1);
1729 double l1x, l2x, l3x, l1y, l2y, l3y;
1730 double d1x, d2x, d3x, d1y, d2y, d3y;
1732 l1x = (x - 1.) * (2. * x - 1);
1733 l2x = 4. * x * (1. - x);
1734 l3x = x * (2. * x - 1.);
1735 l1y = (y - 1.) * (2. * y - 1);
1736 l2y = 4. * y * (1. - y);
1737 l3y = y * (2. * y - 1.);
1739 d1x = a * (4. * x - 3.);
1740 d2x = a * (4. - 8. * x);
1741 d3x = a * (4. * x - 1.);
1742 d1y = a * (4. * y - 3.);
1743 d2y = a * (4. - 8. * y);
1744 d3y = a * (4. * y - 1.);
1746 dshape(0,0) = d1x * l1y;
1747 dshape(0,1) = l1x * d1y;
1749 dshape(4,0) = d2x * l1y;
1750 dshape(4,1) = l2x * d1y;
1752 dshape(1,0) = d3x * l1y;
1753 dshape(1,1) = l3x * d1y;
1755 dshape(7,0) = d1x * l2y;
1756 dshape(7,1) = l1x * d2y;
1758 dshape(8,0) = d2x * l2y;
1759 dshape(8,1) = l2x * d2y;
1761 dshape(5,0) = d3x * l2y;
1762 dshape(5,1) = l3x * d2y;
1764 dshape(3,0) = d1x * l3y;
1765 dshape(3,1) = l1x * d3y;
1767 dshape(6,0) = d2x * l3y;
1768 dshape(6,1) = l2x * d3y;
1770 dshape(2,0) = d3x * l3y;
1771 dshape(2,1) = l3x * d3y;
1814 double x = ip.
x, y = ip.
y;
1816 double w1x, w2x, w3x, w1y, w2y, w3y;
1817 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
1819 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
1820 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
1822 l0x = (- 4.5) * w1x * w2x * w3x;
1823 l1x = ( 13.5) * x * w2x * w3x;
1824 l2x = (-13.5) * x * w1x * w3x;
1825 l3x = ( 4.5) * x * w1x * w2x;
1827 l0y = (- 4.5) * w1y * w2y * w3y;
1828 l1y = ( 13.5) * y * w2y * w3y;
1829 l2y = (-13.5) * y * w1y * w3y;
1830 l3y = ( 4.5) * y * w1y * w2y;
1832 shape(0) = l0x * l0y;
1833 shape(1) = l3x * l0y;
1834 shape(2) = l3x * l3y;
1835 shape(3) = l0x * l3y;
1836 shape(4) = l1x * l0y;
1837 shape(5) = l2x * l0y;
1838 shape(6) = l3x * l1y;
1839 shape(7) = l3x * l2y;
1840 shape(8) = l2x * l3y;
1841 shape(9) = l1x * l3y;
1842 shape(10) = l0x * l2y;
1843 shape(11) = l0x * l1y;
1844 shape(12) = l1x * l1y;
1845 shape(13) = l2x * l1y;
1846 shape(14) = l1x * l2y;
1847 shape(15) = l2x * l2y;
1853 double x = ip.
x, y = ip.
y;
1855 double w1x, w2x, w3x, w1y, w2y, w3y;
1856 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
1857 double d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
1859 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
1860 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
1862 l0x = (- 4.5) * w1x * w2x * w3x;
1863 l1x = ( 13.5) * x * w2x * w3x;
1864 l2x = (-13.5) * x * w1x * w3x;
1865 l3x = ( 4.5) * x * w1x * w2x;
1867 l0y = (- 4.5) * w1y * w2y * w3y;
1868 l1y = ( 13.5) * y * w2y * w3y;
1869 l2y = (-13.5) * y * w1y * w3y;
1870 l3y = ( 4.5) * y * w1y * w2y;
1872 d0x = -5.5 + ( 18. - 13.5 * x) * x;
1873 d1x = 9. + (-45. + 40.5 * x) * x;
1874 d2x = -4.5 + ( 36. - 40.5 * x) * x;
1875 d3x = 1. + (- 9. + 13.5 * x) * x;
1877 d0y = -5.5 + ( 18. - 13.5 * y) * y;
1878 d1y = 9. + (-45. + 40.5 * y) * y;
1879 d2y = -4.5 + ( 36. - 40.5 * y) * y;
1880 d3y = 1. + (- 9. + 13.5 * y) * y;
1882 dshape( 0,0) = d0x * l0y; dshape( 0,1) = l0x * d0y;
1883 dshape( 1,0) = d3x * l0y; dshape( 1,1) = l3x * d0y;
1884 dshape( 2,0) = d3x * l3y; dshape( 2,1) = l3x * d3y;
1885 dshape( 3,0) = d0x * l3y; dshape( 3,1) = l0x * d3y;
1886 dshape( 4,0) = d1x * l0y; dshape( 4,1) = l1x * d0y;
1887 dshape( 5,0) = d2x * l0y; dshape( 5,1) = l2x * d0y;
1888 dshape( 6,0) = d3x * l1y; dshape( 6,1) = l3x * d1y;
1889 dshape( 7,0) = d3x * l2y; dshape( 7,1) = l3x * d2y;
1890 dshape( 8,0) = d2x * l3y; dshape( 8,1) = l2x * d3y;
1891 dshape( 9,0) = d1x * l3y; dshape( 9,1) = l1x * d3y;
1892 dshape(10,0) = d0x * l2y; dshape(10,1) = l0x * d2y;
1893 dshape(11,0) = d0x * l1y; dshape(11,1) = l0x * d1y;
1894 dshape(12,0) = d1x * l1y; dshape(12,1) = l1x * d1y;
1895 dshape(13,0) = d2x * l1y; dshape(13,1) = l2x * d1y;
1896 dshape(14,0) = d1x * l2y; dshape(14,1) = l1x * d2y;
1897 dshape(15,0) = d2x * l2y; dshape(15,1) = l2x * d2y;
1903 double x = ip.
x, y = ip.
y;
1905 double w1x, w2x, w3x, w1y, w2y, w3y;
1906 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
1907 double d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
1908 double h0x, h1x, h2x, h3x, h0y, h1y, h2y, h3y;
1910 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
1911 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
1913 l0x = (- 4.5) * w1x * w2x * w3x;
1914 l1x = ( 13.5) * x * w2x * w3x;
1915 l2x = (-13.5) * x * w1x * w3x;
1916 l3x = ( 4.5) * x * w1x * w2x;
1918 l0y = (- 4.5) * w1y * w2y * w3y;
1919 l1y = ( 13.5) * y * w2y * w3y;
1920 l2y = (-13.5) * y * w1y * w3y;
1921 l3y = ( 4.5) * y * w1y * w2y;
1923 d0x = -5.5 + ( 18. - 13.5 * x) * x;
1924 d1x = 9. + (-45. + 40.5 * x) * x;
1925 d2x = -4.5 + ( 36. - 40.5 * x) * x;
1926 d3x = 1. + (- 9. + 13.5 * x) * x;
1928 d0y = -5.5 + ( 18. - 13.5 * y) * y;
1929 d1y = 9. + (-45. + 40.5 * y) * y;
1930 d2y = -4.5 + ( 36. - 40.5 * y) * y;
1931 d3y = 1. + (- 9. + 13.5 * y) * y;
1933 h0x = -27. * x + 18.;
1934 h1x = 81. * x - 45.;
1935 h2x = -81. * x + 36.;
1938 h0y = -27. * y + 18.;
1939 h1y = 81. * y - 45.;
1940 h2y = -81. * y + 36.;
1943 h( 0,0) = h0x * l0y; h( 0,1) = d0x * d0y; h( 0,2) = l0x * h0y;
1944 h( 1,0) = h3x * l0y; h( 1,1) = d3x * d0y; h( 1,2) = l3x * h0y;
1945 h( 2,0) = h3x * l3y; h( 2,1) = d3x * d3y; h( 2,2) = l3x * h3y;
1946 h( 3,0) = h0x * l3y; h( 3,1) = d0x * d3y; h( 3,2) = l0x * h3y;
1947 h( 4,0) = h1x * l0y; h( 4,1) = d1x * d0y; h( 4,2) = l1x * h0y;
1948 h( 5,0) = h2x * l0y; h( 5,1) = d2x * d0y; h( 5,2) = l2x * h0y;
1949 h( 6,0) = h3x * l1y; h( 6,1) = d3x * d1y; h( 6,2) = l3x * h1y;
1950 h( 7,0) = h3x * l2y; h( 7,1) = d3x * d2y; h( 7,2) = l3x * h2y;
1951 h( 8,0) = h2x * l3y; h( 8,1) = d2x * d3y; h( 8,2) = l2x * h3y;
1952 h( 9,0) = h1x * l3y; h( 9,1) = d1x * d3y; h( 9,2) = l1x * h3y;
1953 h(10,0) = h0x * l2y; h(10,1) = d0x * d2y; h(10,2) = l0x * h2y;
1954 h(11,0) = h0x * l1y; h(11,1) = d0x * d1y; h(11,2) = l0x * h1y;
1955 h(12,0) = h1x * l1y; h(12,1) = d1x * d1y; h(12,2) = l1x * h1y;
1956 h(13,0) = h2x * l1y; h(13,1) = d2x * d1y; h(13,2) = l2x * h1y;
1957 h(14,0) = h1x * l2y; h(14,1) = d1x * d2y; h(14,2) = l1x * h2y;
1958 h(15,0) = h2x * l2y; h(15,1) = d2x * d2y; h(15,2) = l2x * h2y;
1977 l3 = (0.33333333333333333333-x),
1978 l4 = (0.66666666666666666667-x);
1980 shape(0) = 4.5 * l2 * l3 * l4;
1981 shape(1) = 4.5 * l1 * l3 * l4;
1982 shape(2) = 13.5 * l1 * l2 * l4;
1983 shape(3) = -13.5 * l1 * l2 * l3;
1991 dshape(0,0) = -5.5 + x * (18. - 13.5 * x);
1992 dshape(1,0) = 1. - x * (9. - 13.5 * x);
1993 dshape(2,0) = 9. - x * (45. - 40.5 * x);
1994 dshape(3,0) = -4.5 + x * (36. - 40.5 * x);
2026 double x = ip.
x, y = ip.
y;
2027 double l1 = (-1. + x + y),
2031 shape(0) = -0.5*l1*(3.*l1 + 1.)*(3.*l1 + 2.);
2032 shape(1) = 0.5*x*(lx - 1.)*lx;
2033 shape(2) = 0.5*y*(-1. + ly)*ly;
2034 shape(3) = 4.5*x*l1*(3.*l1 + 1.);
2035 shape(4) = -4.5*x*lx*l1;
2036 shape(5) = 4.5*x*lx*y;
2037 shape(6) = 4.5*x*y*ly;
2038 shape(7) = -4.5*y*l1*ly;
2039 shape(8) = 4.5*y*l1*(1. + 3.*l1);
2040 shape(9) = -27.*x*y*l1;
2046 double x = ip.
x, y = ip.
y;
2048 dshape(0,0) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
2049 dshape(1,0) = 1. + 4.5*x*(-2. + 3.*x);
2051 dshape(3,0) = 4.5*(2. + 9.*x*x - 5.*y + 3.*y*y + 2.*x*(-5. + 6.*y));
2052 dshape(4,0) = -4.5*(1. - 1.*y + x*(-8. + 9.*x + 6.*y));
2053 dshape(5,0) = 4.5*(-1. + 6.*x)*y;
2054 dshape(6,0) = 4.5*y*(-1. + 3.*y);
2055 dshape(7,0) = 4.5*(1. - 3.*y)*y;
2056 dshape(8,0) = 4.5*y*(-5. + 6.*x + 6.*y);
2057 dshape(9,0) = -27.*y*(-1. + 2.*x + y);
2059 dshape(0,1) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
2061 dshape(2,1) = 1. + 4.5*y*(-2. + 3.*y);
2062 dshape(3,1) = 4.5*x*(-5. + 6.*x + 6.*y);
2063 dshape(4,1) = 4.5*(1. - 3.*x)*x;
2064 dshape(5,1) = 4.5*x*(-1. + 3.*x);
2065 dshape(6,1) = 4.5*x*(-1. + 6.*y);
2066 dshape(7,1) = -4.5*(1. + x*(-1. + 6.*y) + y*(-8. + 9.*y));
2067 dshape(8,1) = 4.5*(2. + 3.*x*x + y*(-10. + 9.*y) + x*(-5. + 12.*y));
2068 dshape(9,1) = -27.*x*(-1. + x + 2.*y);
2074 double x = ip.
x, y = ip.
y;
2076 h(0,0) = 18.-27.*(x+y);
2077 h(0,1) = 18.-27.*(x+y);
2078 h(0,2) = 18.-27.*(x+y);
2088 h(3,0) = -45.+81.*x+54.*y;
2089 h(3,1) = -22.5+54.*x+27.*y;
2092 h(4,0) = 36.-81.*x-27.*y;
2097 h(5,1) = -4.5+27.*x;
2101 h(6,1) = -4.5+27.*y;
2106 h(7,2) = 36.-27.*x-81.*y;
2109 h(8,1) = -22.5+27.*x+54.*y;
2110 h(8,2) = -45.+54.*x+81.*y;
2113 h(9,1) = 27.-54.*(x+y);
2186 double x = ip.
x, y = ip.
y, z = ip.
z;
2188 shape(0) = -((-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z)*
2189 (-1 + 3*x + 3*y + 3*z))/2.;
2190 shape(4) = (9*x*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
2191 shape(5) = (-9*x*(-1 + 3*x)*(-1 + x + y + z))/2.;
2192 shape(1) = (x*(2 + 9*(-1 + x)*x))/2.;
2193 shape(6) = (9*y*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
2194 shape(19) = -27*x*y*(-1 + x + y + z);
2195 shape(10) = (9*x*(-1 + 3*x)*y)/2.;
2196 shape(7) = (-9*y*(-1 + 3*y)*(-1 + x + y + z))/2.;
2197 shape(11) = (9*x*y*(-1 + 3*y))/2.;
2198 shape(2) = (y*(2 + 9*(-1 + y)*y))/2.;
2199 shape(8) = (9*z*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
2200 shape(18) = -27*x*z*(-1 + x + y + z);
2201 shape(12) = (9*x*(-1 + 3*x)*z)/2.;
2202 shape(17) = -27*y*z*(-1 + x + y + z);
2203 shape(16) = 27*x*y*z;
2204 shape(14) = (9*y*(-1 + 3*y)*z)/2.;
2205 shape(9) = (-9*z*(-1 + x + y + z)*(-1 + 3*z))/2.;
2206 shape(13) = (9*x*z*(-1 + 3*z))/2.;
2207 shape(15) = (9*y*z*(-1 + 3*z))/2.;
2208 shape(3) = (z*(2 + 9*(-1 + z)*z))/2.;
2214 double x = ip.
x, y = ip.
y, z = ip.
z;
2216 dshape(0,0) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
2217 x*(-4 + 6*y + 6*z)))/2.;
2218 dshape(0,1) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
2219 x*(-4 + 6*y + 6*z)))/2.;
2220 dshape(0,2) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
2221 x*(-4 + 6*y + 6*z)))/2.;
2222 dshape(4,0) = (9*(9*pow(x,2) + (-1 + y + z)*(-2 + 3*y + 3*z) +
2223 2*x*(-5 + 6*y + 6*z)))/2.;
2224 dshape(4,1) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
2225 dshape(4,2) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
2226 dshape(5,0) = (-9*(1 - y - z + x*(-8 + 9*x + 6*y + 6*z)))/2.;
2227 dshape(5,1) = (9*(1 - 3*x)*x)/2.;
2228 dshape(5,2) = (9*(1 - 3*x)*x)/2.;
2229 dshape(1,0) = 1 + (9*x*(-2 + 3*x))/2.;
2232 dshape(6,0) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
2233 dshape(6,1) = (9*(2 + 3*pow(x,2) - 10*y - 5*z + 3*(y + z)*(3*y + z) +
2234 x*(-5 + 12*y + 6*z)))/2.;
2235 dshape(6,2) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
2236 dshape(19,0) = -27*y*(-1 + 2*x + y + z);
2237 dshape(19,1) = -27*x*(-1 + x + 2*y + z);
2238 dshape(19,2) = -27*x*y;
2239 dshape(10,0) = (9*(-1 + 6*x)*y)/2.;
2240 dshape(10,1) = (9*x*(-1 + 3*x))/2.;
2242 dshape(7,0) = (9*(1 - 3*y)*y)/2.;
2243 dshape(7,1) = (-9*(1 + x*(-1 + 6*y) - z + y*(-8 + 9*y + 6*z)))/2.;
2244 dshape(7,2) = (9*(1 - 3*y)*y)/2.;
2245 dshape(11,0) = (9*y*(-1 + 3*y))/2.;
2246 dshape(11,1) = (9*x*(-1 + 6*y))/2.;
2249 dshape(2,1) = 1 + (9*y*(-2 + 3*y))/2.;
2251 dshape(8,0) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
2252 dshape(8,1) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
2253 dshape(8,2) = (9*(2 + 3*pow(x,2) - 5*y - 10*z + 3*(y + z)*(y + 3*z) +
2254 x*(-5 + 6*y + 12*z)))/2.;
2255 dshape(18,0) = -27*z*(-1 + 2*x + y + z);
2256 dshape(18,1) = -27*x*z;
2257 dshape(18,2) = -27*x*(-1 + x + y + 2*z);
2258 dshape(12,0) = (9*(-1 + 6*x)*z)/2.;
2260 dshape(12,2) = (9*x*(-1 + 3*x))/2.;
2261 dshape(17,0) = -27*y*z;
2262 dshape(17,1) = -27*z*(-1 + x + 2*y + z);
2263 dshape(17,2) = -27*y*(-1 + x + y + 2*z);
2264 dshape(16,0) = 27*y*z;
2265 dshape(16,1) = 27*x*z;
2266 dshape(16,2) = 27*x*y;
2268 dshape(14,1) = (9*(-1 + 6*y)*z)/2.;
2269 dshape(14,2) = (9*y*(-1 + 3*y))/2.;
2270 dshape(9,0) = (9*(1 - 3*z)*z)/2.;
2271 dshape(9,1) = (9*(1 - 3*z)*z)/2.;
2272 dshape(9,2) = (9*(-1 + x + y + 8*z - 6*(x + y)*z - 9*pow(z,2)))/2.;
2273 dshape(13,0) = (9*z*(-1 + 3*z))/2.;
2275 dshape(13,2) = (9*x*(-1 + 6*z))/2.;
2277 dshape(15,1) = (9*z*(-1 + 3*z))/2.;
2278 dshape(15,2) = (9*y*(-1 + 6*z))/2.;
2281 dshape(3,2) = 1 + (9*z*(-2 + 3*z))/2.;
2347 shape(0) = 1. - ip.
x - ip.
y - ip.
z;
2356 if (dshape.
Height() == 4)
2358 double *A = &dshape(0,0);
2359 A[0] = -1.; A[4] = -1.; A[8] = -1.;
2360 A[1] = 1.; A[5] = 0.; A[9] = 0.;
2361 A[2] = 0.; A[6] = 1.; A[10] = 0.;
2362 A[3] = 0.; A[7] = 0.; A[11] = 1.;
2366 dshape(0,0) = -1.; dshape(0,1) = -1.; dshape(0,2) = -1.;
2367 dshape(1,0) = 1.; dshape(1,1) = 0.; dshape(1,2) = 0.;
2368 dshape(2,0) = 0.; dshape(2,1) = 1.; dshape(2,2) = 0.;
2369 dshape(3,0) = 0.; dshape(3,1) = 0.; dshape(3,2) = 1.;
2376 static int face_dofs[4][3] = {{1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {0, 1, 2}};
2379 *dofs = face_dofs[face];
2421 double L0, L1, L2, L3;
2423 L0 = 1. - ip.
x - ip.
y - ip.
z;
2428 shape(0) = L0 * ( 2.0 * L0 - 1.0 );
2429 shape(1) = L1 * ( 2.0 * L1 - 1.0 );
2430 shape(2) = L2 * ( 2.0 * L2 - 1.0 );
2431 shape(3) = L3 * ( 2.0 * L3 - 1.0 );
2432 shape(4) = 4.0 * L0 * L1;
2433 shape(5) = 4.0 * L0 * L2;
2434 shape(6) = 4.0 * L0 * L3;
2435 shape(7) = 4.0 * L1 * L2;
2436 shape(8) = 4.0 * L1 * L3;
2437 shape(9) = 4.0 * L2 * L3;
2448 L0 = 1.0 - x - y - z;
2450 dshape(0,0) = dshape(0,1) = dshape(0,2) = 1.0 - 4.0 * L0;
2451 dshape(1,0) = -1.0 + 4.0 * x; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
2452 dshape(2,0) = 0.0; dshape(2,1) = -1.0 + 4.0 * y; dshape(2,2) = 0.0;
2453 dshape(3,0) = dshape(3,1) = 0.0; dshape(3,2) = -1.0 + 4.0 * z;
2454 dshape(4,0) = 4.0 * (L0 - x); dshape(4,1) = dshape(4,2) = -4.0 * x;
2455 dshape(5,0) = dshape(5,2) = -4.0 * y; dshape(5,1) = 4.0 * (L0 - y);
2456 dshape(6,0) = dshape(6,1) = -4.0 * z; dshape(6,2) = 4.0 * (L0 - z);
2457 dshape(7,0) = 4.0 * y; dshape(7,1) = 4.0 * x; dshape(7,2) = 0.0;
2458 dshape(8,0) = 4.0 * z; dshape(8,1) = 0.0; dshape(8,2) = 4.0 * x;
2459 dshape(9,0) = 0.0; dshape(9,1) = 4.0 * z; dshape(9,2) = 4.0 * y;
2501 double x = ip.
x, y = ip.
y, z = ip.
z;
2502 double ox = 1.-x, oy = 1.-y, oz = 1.-z;
2504 shape(0) = ox * oy * oz;
2505 shape(1) = x * oy * oz;
2506 shape(2) = x * y * oz;
2507 shape(3) = ox * y * oz;
2508 shape(4) = ox * oy * z;
2509 shape(5) = x * oy * z;
2510 shape(6) = x * y * z;
2511 shape(7) = ox * y * z;
2517 double x = ip.
x, y = ip.
y, z = ip.
z;
2518 double ox = 1.-x, oy = 1.-y, oz = 1.-z;
2520 dshape(0,0) = - oy * oz;
2521 dshape(0,1) = - ox * oz;
2522 dshape(0,2) = - ox * oy;
2524 dshape(1,0) = oy * oz;
2525 dshape(1,1) = - x * oz;
2526 dshape(1,2) = - x * oy;
2528 dshape(2,0) = y * oz;
2529 dshape(2,1) = x * oz;
2530 dshape(2,2) = - x * y;
2532 dshape(3,0) = - y * oz;
2533 dshape(3,1) = ox * oz;
2534 dshape(3,2) = - ox * y;
2536 dshape(4,0) = - oy * z;
2537 dshape(4,1) = - ox * z;
2538 dshape(4,2) = ox * oy;
2540 dshape(5,0) = oy * z;
2541 dshape(5,1) = - x * z;
2542 dshape(5,2) = x * oy;
2544 dshape(6,0) = y * z;
2545 dshape(6,1) = x * z;
2546 dshape(6,2) = x * y;
2548 dshape(7,0) = - y * z;
2549 dshape(7,1) = ox * z;
2550 dshape(7,2) = ox * y;
2585 shape(0) = 1.0 - 2.0 * ip.
y;
2586 shape(1) = -1.0 + 2.0 * ( ip.
x + ip.
y );
2587 shape(2) = 1.0 - 2.0 * ip.
x;
2593 dshape(0,0) = 0.0; dshape(0,1) = -2.0;
2594 dshape(1,0) = 2.0; dshape(1,1) = 2.0;
2595 dshape(2,0) = -2.0; dshape(2,1) = 0.0;
2616 const double l1 = ip.
x+ip.
y-0.5, l2 = 1.-l1, l3 = ip.
x-ip.
y+0.5, l4 = 1.-l3;
2627 const double x2 = 2.*ip.
x, y2 = 2.*ip.
y;
2629 dshape(0,0) = 1. - x2; dshape(0,1) = -2. + y2;
2630 dshape(1,0) = x2; dshape(1,1) = 1. - y2;
2631 dshape(2,0) = 1. - x2; dshape(2,1) = y2;
2632 dshape(3,0) = -2. + x2; dshape(3,1) = 1. - y2;
2650 double x = ip.
x, y = ip.
y;
2653 shape(0,1) = y - 1.;
2656 shape(2,0) = x - 1.;
2668 const double RT0TriangleFiniteElement::nk[3][2] =
2669 { {0, -1}, {1, 1}, {-1, 0} };
2675 #ifdef MFEM_THREAD_SAFE
2681 for (k = 0; k < 3; k++)
2684 for (j = 0; j < 3; j++)
2687 if (j == k) { d -= 1.0; }
2688 if (fabs(d) > 1.0e-12)
2690 mfem::err <<
"RT0TriangleFiniteElement::GetLocalInterpolation (...)\n"
2691 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2707 for (k = 0; k < 3; k++)
2710 ip.
x = vk[0]; ip.
y = vk[1];
2713 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2714 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2715 for (j = 0; j < 3; j++)
2716 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2729 #ifdef MFEM_THREAD_SAFE
2733 for (
int k = 0; k < 3; k++)
2741 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2742 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
2762 double x = ip.
x, y = ip.
y;
2765 shape(0,1) = y - 1.;
2770 shape(3,0) = x - 1.;
2783 const double RT0QuadFiniteElement::nk[4][2] =
2784 { {0, -1}, {1, 0}, {0, 1}, {-1, 0} };
2790 #ifdef MFEM_THREAD_SAFE
2796 for (k = 0; k < 4; k++)
2799 for (j = 0; j < 4; j++)
2802 if (j == k) { d -= 1.0; }
2803 if (fabs(d) > 1.0e-12)
2805 mfem::err <<
"RT0QuadFiniteElement::GetLocalInterpolation (...)\n"
2806 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2822 for (k = 0; k < 4; k++)
2825 ip.
x = vk[0]; ip.
y = vk[1];
2828 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2829 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2830 for (j = 0; j < 4; j++)
2831 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2844 #ifdef MFEM_THREAD_SAFE
2848 for (
int k = 0; k < 4; k++)
2856 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2857 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
2885 double x = ip.
x, y = ip.
y;
2887 shape(0,0) = -2 * x * (-1 + x + 2 * y);
2888 shape(0,1) = -2 * (-1 + y) * (-1 + x + 2 * y);
2889 shape(1,0) = 2 * x * (x - y);
2890 shape(1,1) = 2 * (x - y) * (-1 + y);
2891 shape(2,0) = 2 * x * (-1 + 2 * x + y);
2892 shape(2,1) = 2 * y * (-1 + 2 * x + y);
2893 shape(3,0) = 2 * x * (-1 + x + 2 * y);
2894 shape(3,1) = 2 * y * (-1 + x + 2 * y);
2895 shape(4,0) = -2 * (-1 + x) * (x - y);
2896 shape(4,1) = 2 * y * (-x + y);
2897 shape(5,0) = -2 * (-1 + x) * (-1 + 2 * x + y);
2898 shape(5,1) = -2 * y * (-1 + 2 * x + y);
2899 shape(6,0) = -3 * x * (-2 + 2 * x + y);
2900 shape(6,1) = -3 * y * (-1 + 2 * x + y);
2901 shape(7,0) = -3 * x * (-1 + x + 2 * y);
2902 shape(7,1) = -3 * y * (-2 + x + 2 * y);
2908 double x = ip.
x, y = ip.
y;
2910 divshape(0) = -2 * (-4 + 3 * x + 6 * y);
2911 divshape(1) = 2 + 6 * x - 6 * y;
2912 divshape(2) = -4 + 12 * x + 6 * y;
2913 divshape(3) = -4 + 6 * x + 12 * y;
2914 divshape(4) = 2 - 6 * x + 6 * y;
2915 divshape(5) = -2 * (-4 + 6 * x + 3 * y);
2916 divshape(6) = -9 * (-1 + 2 * x + y);
2917 divshape(7) = -9 * (-1 + x + 2 * y);
2920 const double RT1TriangleFiniteElement::nk[8][2] =
2932 #ifdef MFEM_THREAD_SAFE
2938 for (k = 0; k < 8; k++)
2941 for (j = 0; j < 8; j++)
2944 if (j == k) { d -= 1.0; }
2945 if (fabs(d) > 1.0e-12)
2947 mfem::err <<
"RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
2948 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2964 for (k = 0; k < 8; k++)
2967 ip.
x = vk[0]; ip.
y = vk[1];
2970 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2971 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2972 for (j = 0; j < 8; j++)
2973 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2985 #ifdef MFEM_THREAD_SAFE
2989 for (
int k = 0; k < 8; k++)
2997 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2998 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
3041 double x = ip.
x, y = ip.
y;
3045 shape(0,1) = -( 1. - 3.*y + 2.*y*y)*( 2. - 3.*x);
3047 shape(1,1) = -( 1. - 3.*y + 2.*y*y)*(-1. + 3.*x);
3049 shape(2,0) = (-x + 2.*x*x)*( 2. - 3.*y);
3051 shape(3,0) = (-x + 2.*x*x)*(-1. + 3.*y);
3055 shape(4,1) = (-y + 2.*y*y)*(-1. + 3.*x);
3057 shape(5,1) = (-y + 2.*y*y)*( 2. - 3.*x);
3059 shape(6,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y);
3061 shape(7,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y);
3064 shape(8,0) = (4.*x - 4.*x*x)*( 2. - 3.*y);
3066 shape(9,0) = (4.*x - 4.*x*x)*(-1. + 3.*y);
3070 shape(10,1) = (4.*y - 4.*y*y)*( 2. - 3.*x);
3072 shape(11,1) = (4.*y - 4.*y*y)*(-1. + 3.*x);
3078 double x = ip.
x, y = ip.
y;
3080 divshape(0) = -(-3. + 4.*y)*( 2. - 3.*x);
3081 divshape(1) = -(-3. + 4.*y)*(-1. + 3.*x);
3082 divshape(2) = (-1. + 4.*x)*( 2. - 3.*y);
3083 divshape(3) = (-1. + 4.*x)*(-1. + 3.*y);
3084 divshape(4) = (-1. + 4.*y)*(-1. + 3.*x);
3085 divshape(5) = (-1. + 4.*y)*( 2. - 3.*x);
3086 divshape(6) = -(-3. + 4.*x)*(-1. + 3.*y);
3087 divshape(7) = -(-3. + 4.*x)*( 2. - 3.*y);
3088 divshape(8) = ( 4. - 8.*x)*( 2. - 3.*y);
3089 divshape(9) = ( 4. - 8.*x)*(-1. + 3.*y);
3090 divshape(10) = ( 4. - 8.*y)*( 2. - 3.*x);
3091 divshape(11) = ( 4. - 8.*y)*(-1. + 3.*x);
3094 const double RT1QuadFiniteElement::nk[12][2] =
3114 #ifdef MFEM_THREAD_SAFE
3120 for (k = 0; k < 12; k++)
3123 for (j = 0; j < 12; j++)
3126 if (j == k) { d -= 1.0; }
3127 if (fabs(d) > 1.0e-12)
3129 mfem::err <<
"RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
3130 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
3146 for (k = 0; k < 12; k++)
3149 ip.
x = vk[0]; ip.
y = vk[1];
3152 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
3153 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
3154 for (j = 0; j < 12; j++)
3155 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
3167 #ifdef MFEM_THREAD_SAFE
3171 for (
int k = 0; k < 12; k++)
3179 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
3180 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
3184 const double RT2TriangleFiniteElement::M[15][15] =
3187 0, -5.3237900077244501311, 5.3237900077244501311, 16.647580015448900262,
3188 0, 24.442740046346700787, -16.647580015448900262, -12.,
3189 -19.118950038622250656, -47.237900077244501311, 0, -34.414110069520051180,
3190 12., 30.590320061795601049, 15.295160030897800524
3193 0, 1.5, -1.5, -15., 0, 2.625, 15., 15., -4.125, 30., 0, -14.625, -15.,
3197 0, -0.67620999227554986889, 0.67620999227554986889, 7.3524199845510997378,
3198 0, -3.4427400463467007866, -7.3524199845510997378, -12.,
3199 4.1189500386222506555, -0.76209992275549868892, 0, 7.4141100695200511800,
3200 12., -6.5903200617956010489, -3.2951600308978005244
3203 0, 0, 1.5, 0, 0, 1.5, -11.471370023173350393, 0, 2.4713700231733503933,
3204 -11.471370023173350393, 0, 2.4713700231733503933, 15.295160030897800524,
3205 0, -3.2951600308978005244
3208 0, 0, 4.875, 0, 0, 4.875, -16.875, 0, -16.875, -16.875, 0, -16.875, 10.5,
3212 0, 0, 1.5, 0, 0, 1.5, 2.4713700231733503933, 0, -11.471370023173350393,
3213 2.4713700231733503933, 0, -11.471370023173350393, -3.2951600308978005244,
3214 0, 15.295160030897800524
3217 -0.67620999227554986889, 0, -3.4427400463467007866, 0,
3218 7.3524199845510997378, 0.67620999227554986889, 7.4141100695200511800, 0,
3219 -0.76209992275549868892, 4.1189500386222506555, -12.,
3220 -7.3524199845510997378, -3.2951600308978005244, -6.5903200617956010489,
3224 1.5, 0, 2.625, 0, -15., -1.5, -14.625, 0, 30., -4.125, 15., 15., 10.5,
3228 -5.3237900077244501311, 0, 24.442740046346700787, 0, 16.647580015448900262,
3229 5.3237900077244501311, -34.414110069520051180, 0, -47.237900077244501311,
3230 -19.118950038622250656, -12., -16.647580015448900262, 15.295160030897800524,
3231 30.590320061795601049, 12.
3233 { 0, 0, 18., 0, 0, 6., -42., 0, -30., -26., 0, -14., 24., 32., 8.},
3234 { 0, 0, 6., 0, 0, 18., -14., 0, -26., -30., 0, -42., 8., 32., 24.},
3235 { 0, 0, -6., 0, 0, -4., 30., 0, 4., 22., 0, 4., -24., -16., 0},
3236 { 0, 0, -4., 0, 0, -8., 20., 0, 8., 36., 0, 8., -16., -32., 0},
3237 { 0, 0, -8., 0, 0, -4., 8., 0, 36., 8., 0, 20., 0, -32., -16.},
3238 { 0, 0, -4., 0, 0, -6., 4., 0, 22., 4., 0, 30., 0, -16., -24.}
3244 const double p = 0.11270166537925831148;
3281 double x = ip.
x, y = ip.
y;
3283 double Bx[15] = {1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y, 0., x*x*x,
3286 double By[15] = {0., 1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y,
3290 for (
int i = 0; i < 15; i++)
3292 double cx = 0.0, cy = 0.0;
3293 for (
int j = 0; j < 15; j++)
3295 cx += M[i][j] * Bx[j];
3296 cy += M[i][j] * By[j];
3306 double x = ip.
x, y = ip.
y;
3308 double DivB[15] = {0., 0., 1., 0., 0., 1., 2.*x, 0., y, x, 0., 2.*y,
3309 4.*x*x, 4.*x*y, 4.*y*y
3312 for (
int i = 0; i < 15; i++)
3315 for (
int j = 0; j < 15; j++)
3317 div += M[i][j] * DivB[j];
3323 const double RT2QuadFiniteElement::pt[4] = {0.,1./3.,2./3.,1.};
3325 const double RT2QuadFiniteElement::dpt[3] = {0.25,0.5,0.75};
3367 double x = ip.
x, y = ip.
y;
3369 double ax0 = pt[0] - x;
3370 double ax1 = pt[1] - x;
3371 double ax2 = pt[2] - x;
3372 double ax3 = pt[3] - x;
3374 double by0 = dpt[0] - y;
3375 double by1 = dpt[1] - y;
3376 double by2 = dpt[2] - y;
3378 double ay0 = pt[0] - y;
3379 double ay1 = pt[1] - y;
3380 double ay2 = pt[2] - y;
3381 double ay3 = pt[3] - y;
3383 double bx0 = dpt[0] - x;
3384 double bx1 = dpt[1] - x;
3385 double bx2 = dpt[2] - x;
3387 double A01 = pt[0] - pt[1];
3388 double A02 = pt[0] - pt[2];
3389 double A12 = pt[1] - pt[2];
3390 double A03 = pt[0] - pt[3];
3391 double A13 = pt[1] - pt[3];
3392 double A23 = pt[2] - pt[3];
3394 double B01 = dpt[0] - dpt[1];
3395 double B02 = dpt[0] - dpt[2];
3396 double B12 = dpt[1] - dpt[2];
3398 double tx0 = (bx1*bx2)/(B01*B02);
3399 double tx1 = -(bx0*bx2)/(B01*B12);
3400 double tx2 = (bx0*bx1)/(B02*B12);
3402 double ty0 = (by1*by2)/(B01*B02);
3403 double ty1 = -(by0*by2)/(B01*B12);
3404 double ty2 = (by0*by1)/(B02*B12);
3408 shape(0, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx0;
3410 shape(1, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx1;
3412 shape(2, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx2;
3414 shape(3, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty0;
3416 shape(4, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty1;
3418 shape(5, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty2;
3422 shape(6, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx2;
3424 shape(7, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx1;
3426 shape(8, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx0;
3428 shape(9, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty2;
3430 shape(10, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty1;
3432 shape(11, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty0;
3435 shape(12, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty0;
3437 shape(13, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty1;
3439 shape(14, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty2;
3442 shape(15, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty0;
3444 shape(16, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty1;
3446 shape(17, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty2;
3450 shape(18, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx0;
3452 shape(19, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx1;
3454 shape(20, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx2;
3457 shape(21, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx0;
3459 shape(22, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx1;
3461 shape(23, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx2;
3467 double x = ip.
x, y = ip.
y;
3469 double a01 = pt[0]*pt[1];
3470 double a02 = pt[0]*pt[2];
3471 double a12 = pt[1]*pt[2];
3472 double a03 = pt[0]*pt[3];
3473 double a13 = pt[1]*pt[3];
3474 double a23 = pt[2]*pt[3];
3476 double bx0 = dpt[0] - x;
3477 double bx1 = dpt[1] - x;
3478 double bx2 = dpt[2] - x;
3480 double by0 = dpt[0] - y;
3481 double by1 = dpt[1] - y;
3482 double by2 = dpt[2] - y;
3484 double A01 = pt[0] - pt[1];
3485 double A02 = pt[0] - pt[2];
3486 double A12 = pt[1] - pt[2];
3487 double A03 = pt[0] - pt[3];
3488 double A13 = pt[1] - pt[3];
3489 double A23 = pt[2] - pt[3];
3491 double A012 = pt[0] + pt[1] + pt[2];
3492 double A013 = pt[0] + pt[1] + pt[3];
3493 double A023 = pt[0] + pt[2] + pt[3];
3494 double A123 = pt[1] + pt[2] + pt[3];
3496 double B01 = dpt[0] - dpt[1];
3497 double B02 = dpt[0] - dpt[2];
3498 double B12 = dpt[1] - dpt[2];
3500 double tx0 = (bx1*bx2)/(B01*B02);
3501 double tx1 = -(bx0*bx2)/(B01*B12);
3502 double tx2 = (bx0*bx1)/(B02*B12);
3504 double ty0 = (by1*by2)/(B01*B02);
3505 double ty1 = -(by0*by2)/(B01*B12);
3506 double ty2 = (by0*by1)/(B02*B12);
3509 divshape(0) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx0;
3510 divshape(1) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx1;
3511 divshape(2) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx2;
3513 divshape(3) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty0;
3514 divshape(4) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty1;
3515 divshape(5) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty2;
3517 divshape(6) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx2;
3518 divshape(7) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx1;
3519 divshape(8) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx0;
3521 divshape(9) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty2;
3522 divshape(10) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty1;
3523 divshape(11) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty0;
3525 divshape(12) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty0;
3526 divshape(13) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty1;
3527 divshape(14) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty2;
3529 divshape(15) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty0;
3530 divshape(16) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty1;
3531 divshape(17) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty2;
3533 divshape(18) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx0;
3534 divshape(19) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx1;
3535 divshape(20) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx2;
3537 divshape(21) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx0;
3538 divshape(22) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx1;
3539 divshape(23) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx2;
3542 const double RT2QuadFiniteElement::nk[24][2] =
3545 {0,-1}, {0,-1}, {0,-1},
3547 {1, 0}, {1, 0}, {1, 0},
3549 {0, 1}, {0, 1}, {0, 1},
3551 {-1,0}, {-1,0}, {-1,0},
3553 {1, 0}, {1, 0}, {1, 0},
3555 {1, 0}, {1, 0}, {1, 0},
3557 {0, 1}, {0, 1}, {0, 1},
3559 {0, 1}, {0, 1}, {0, 1}
3566 #ifdef MFEM_THREAD_SAFE
3572 for (k = 0; k < 24; k++)
3575 for (j = 0; j < 24; j++)
3578 if (j == k) { d -= 1.0; }
3579 if (fabs(d) > 1.0e-12)
3581 mfem::err <<
"RT2QuadFiniteElement::GetLocalInterpolation (...)\n"
3582 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
3598 for (k = 0; k < 24; k++)
3601 ip.
x = vk[0]; ip.
y = vk[1];
3604 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
3605 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
3606 for (j = 0; j < 24; j++)
3607 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
3619 #ifdef MFEM_THREAD_SAFE
3623 for (
int k = 0; k < 24; k++)
3631 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
3632 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
3648 shape(0) = 2. - 3. * x;
3649 shape(1) = 3. * x - 1.;
3663 const double p = 0.11270166537925831148;
3673 const double p = 0.11270166537925831148;
3674 const double w = 1./((1-2*p)*(1-2*p));
3677 shape(0) = (2*x-1)*(x-1+p)*w;
3678 shape(1) = 4*(x-1+p)*(p-x)*w;
3679 shape(2) = (2*x-1)*(x-p)*w;
3685 const double p = 0.11270166537925831148;
3686 const double w = 1./((1-2*p)*(1-2*p));
3689 dshape(0,0) = (-3+4*x+2*p)*w;
3690 dshape(1,0) = (4-8*x)*w;
3691 dshape(2,0) = (-1+4*x-2*p)*w;
3702 for (i = 1; i < m; i++)
3708 #ifndef MFEM_THREAD_SAFE
3713 for (i = 1; i <= m; i++)
3715 rwk(i) = rwk(i-1) * ( (double)(m) / (double)(i) );
3717 for (i = 0; i < m/2+1; i++)
3719 rwk(m-i) = ( rwk(i) *= rwk(m-i) );
3721 for (i = m-1; i >= 0; i -= 2)
3730 double w, wk, x = ip.
x;
3733 #ifdef MFEM_THREAD_SAFE
3737 k = (int) floor ( m * x + 0.5 );
3738 k = k > m ? m : k < 0 ? 0 : k;
3741 for (i = 0; i <= m; i++)
3744 wk *= ( rxxk(i) = x - (double)(i) / m );
3746 w = wk * ( rxxk(k) = x - (double)(k) / m );
3750 shape(0) = w * rwk(0) / rxxk(0);
3754 shape(0) = wk * rwk(0);
3758 shape(1) = w * rwk(m) / rxxk(m);
3762 shape(1) = wk * rwk(k);
3764 for (i = 1; i < m; i++)
3767 shape(i+1) = w * rwk(i) / rxxk(i);
3771 shape(k+1) = wk * rwk(k);
3778 double s, srx, w, wk, x = ip.
x;
3781 #ifdef MFEM_THREAD_SAFE
3785 k = (int) floor ( m * x + 0.5 );
3786 k = k > m ? m : k < 0 ? 0 : k;
3789 for (i = 0; i <= m; i++)
3792 wk *= ( rxxk(i) = x - (double)(i) / m );
3794 w = wk * ( rxxk(k) = x - (double)(k) / m );
3796 for (i = 0; i <= m; i++)
3798 rxxk(i) = 1.0 / rxxk(i);
3801 for (i = 0; i <= m; i++)
3810 dshape(0,0) = (s - w * rxxk(0)) * rwk(0) * rxxk(0);
3814 dshape(0,0) = wk * srx * rwk(0);
3818 dshape(1,0) = (s - w * rxxk(m)) * rwk(m) * rxxk(m);
3822 dshape(1,0) = wk * srx * rwk(k);
3824 for (i = 1; i < m; i++)
3827 dshape(i+1,0) = (s - w * rxxk(i)) * rwk(i) * rxxk(i);
3831 dshape(k+1,0) = wk * srx * rwk(k);
3860 double L0, L1, L2, L3;
3862 L1 = ip.
x; L2 = ip.
y; L3 = ip.
z; L0 = 1.0 - L1 - L2 - L3;
3863 shape(0) = 1.0 - 3.0 * L0;
3864 shape(1) = 1.0 - 3.0 * L1;
3865 shape(2) = 1.0 - 3.0 * L2;
3866 shape(3) = 1.0 - 3.0 * L3;
3872 dshape(0,0) = 3.0; dshape(0,1) = 3.0; dshape(0,2) = 3.0;
3873 dshape(1,0) = -3.0; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
3874 dshape(2,0) = 0.0; dshape(2,1) = -3.0; dshape(2,2) = 0.0;
3875 dshape(3,0) = 0.0; dshape(3,1) = 0.0; dshape(3,2) = -3.0;
3896 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3917 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3931 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3932 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3933 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3934 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3935 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3936 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3937 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3938 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3940 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3941 I[ 9] = 1; J[ 9] = 2; K[ 9] = 0;
3942 I[10] = 2; J[10] = 1; K[10] = 0;
3943 I[11] = 0; J[11] = 2; K[11] = 0;
3944 I[12] = 2; J[12] = 0; K[12] = 1;
3945 I[13] = 1; J[13] = 2; K[13] = 1;
3946 I[14] = 2; J[14] = 1; K[14] = 1;
3947 I[15] = 0; J[15] = 2; K[15] = 1;
3948 I[16] = 0; J[16] = 0; K[16] = 2;
3949 I[17] = 1; J[17] = 0; K[17] = 2;
3950 I[18] = 1; J[18] = 1; K[18] = 2;
3951 I[19] = 0; J[19] = 1; K[19] = 2;
3953 I[20] = 2; J[20] = 2; K[20] = 0;
3954 I[21] = 2; J[21] = 0; K[21] = 2;
3955 I[22] = 1; J[22] = 2; K[22] = 2;
3956 I[23] = 2; J[23] = 1; K[23] = 2;
3957 I[24] = 0; J[24] = 2; K[24] = 2;
3958 I[25] = 2; J[25] = 2; K[25] = 1;
3960 I[26] = 2; J[26] = 2; K[26] = 2;
3962 else if (degree == 3)
3968 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3969 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3970 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3971 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3972 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3973 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3974 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3975 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3977 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3978 I[ 9] = 3; J[ 9] = 0; K[ 9] = 0;
3979 I[10] = 1; J[10] = 2; K[10] = 0;
3980 I[11] = 1; J[11] = 3; K[11] = 0;
3981 I[12] = 2; J[12] = 1; K[12] = 0;
3982 I[13] = 3; J[13] = 1; K[13] = 0;
3983 I[14] = 0; J[14] = 2; K[14] = 0;
3984 I[15] = 0; J[15] = 3; K[15] = 0;
3985 I[16] = 2; J[16] = 0; K[16] = 1;
3986 I[17] = 3; J[17] = 0; K[17] = 1;
3987 I[18] = 1; J[18] = 2; K[18] = 1;
3988 I[19] = 1; J[19] = 3; K[19] = 1;
3989 I[20] = 2; J[20] = 1; K[20] = 1;
3990 I[21] = 3; J[21] = 1; K[21] = 1;
3991 I[22] = 0; J[22] = 2; K[22] = 1;
3992 I[23] = 0; J[23] = 3; K[23] = 1;
3993 I[24] = 0; J[24] = 0; K[24] = 2;
3994 I[25] = 0; J[25] = 0; K[25] = 3;
3995 I[26] = 1; J[26] = 0; K[26] = 2;
3996 I[27] = 1; J[27] = 0; K[27] = 3;
3997 I[28] = 1; J[28] = 1; K[28] = 2;
3998 I[29] = 1; J[29] = 1; K[29] = 3;
3999 I[30] = 0; J[30] = 1; K[30] = 2;
4000 I[31] = 0; J[31] = 1; K[31] = 3;
4002 I[32] = 2; J[32] = 3; K[32] = 0;
4003 I[33] = 3; J[33] = 3; K[33] = 0;
4004 I[34] = 2; J[34] = 2; K[34] = 0;
4005 I[35] = 3; J[35] = 2; K[35] = 0;
4006 I[36] = 2; J[36] = 0; K[36] = 2;
4007 I[37] = 3; J[37] = 0; K[37] = 2;
4008 I[38] = 2; J[38] = 0; K[38] = 3;
4009 I[39] = 3; J[39] = 0; K[39] = 3;
4010 I[40] = 1; J[40] = 2; K[40] = 2;
4011 I[41] = 1; J[41] = 3; K[41] = 2;
4012 I[42] = 1; J[42] = 2; K[42] = 3;
4013 I[43] = 1; J[43] = 3; K[43] = 3;
4014 I[44] = 3; J[44] = 1; K[44] = 2;
4015 I[45] = 2; J[45] = 1; K[45] = 2;
4016 I[46] = 3; J[46] = 1; K[46] = 3;
4017 I[47] = 2; J[47] = 1; K[47] = 3;
4018 I[48] = 0; J[48] = 3; K[48] = 2;
4019 I[49] = 0; J[49] = 2; K[49] = 2;
4020 I[50] = 0; J[50] = 3; K[50] = 3;
4021 I[51] = 0; J[51] = 2; K[51] = 3;
4022 I[52] = 2; J[52] = 2; K[52] = 1;
4023 I[53] = 3; J[53] = 2; K[53] = 1;
4024 I[54] = 2; J[54] = 3; K[54] = 1;
4025 I[55] = 3; J[55] = 3; K[55] = 1;
4027 I[56] = 2; J[56] = 2; K[56] = 2;
4028 I[57] = 3; J[57] = 2; K[57] = 2;
4029 I[58] = 3; J[58] = 3; K[58] = 2;
4030 I[59] = 2; J[59] = 3; K[59] = 2;
4031 I[60] = 2; J[60] = 2; K[60] = 3;
4032 I[61] = 3; J[61] = 2; K[61] = 3;
4033 I[62] = 3; J[62] = 3; K[62] = 3;
4034 I[63] = 2; J[63] = 3; K[63] = 3;
4038 mfem_error (
"LagrangeHexFiniteElement::LagrangeHexFiniteElement");
4042 dof1d = fe1d ->
GetDof();
4044 #ifndef MFEM_THREAD_SAFE
4054 for (
int n = 0; n <
Dof; n++)
4069 #ifdef MFEM_THREAD_SAFE
4070 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
4077 for (
int n = 0; n <
Dof; n++)
4079 shape(n) = shape1dx(I[n]) * shape1dy(J[n]) * shape1dz(K[n]);
4090 #ifdef MFEM_THREAD_SAFE
4091 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
4092 DenseMatrix dshape1dx(dof1d,1), dshape1dy(dof1d,1), dshape1dz(dof1d,1);
4103 for (
int n = 0; n <
Dof; n++)
4105 dshape(n,0) = dshape1dx(I[n],0) * shape1dy(J[n]) * shape1dz(K[n]);
4106 dshape(n,1) = shape1dx(I[n]) * dshape1dy(J[n],0) * shape1dz(K[n]);
4107 dshape(n,2) = shape1dx(I[n]) * shape1dy(J[n]) * dshape1dz(K[n],0);
4136 shape(0) = 1.0 - 2.0 * x;
4143 shape(1) = 2.0 * x - 1.0;
4144 shape(2) = 2.0 - 2.0 * x;
4155 dshape(0,0) = - 2.0;
4163 dshape(2,0) = - 2.0;
4190 L0 = 2.0 * ( 1. - ip.
x - ip.
y );
4191 L1 = 2.0 * ( ip.
x );
4192 L2 = 2.0 * ( ip.
y );
4201 for (i = 0; i < 6; i++)
4208 shape(0) = L0 - 1.0;
4215 shape(1) = L1 - 1.0;
4222 shape(2) = L2 - 1.0;
4226 shape(3) = 1.0 - L2;
4227 shape(4) = 1.0 - L0;
4228 shape(5) = 1.0 - L1;
4238 L0 = 2.0 * ( 1. - ip.
x - ip.
y );
4239 L1 = 2.0 * ( ip.
x );
4240 L2 = 2.0 * ( ip.
y );
4242 double DL0[2], DL1[2], DL2[2];
4243 DL0[0] = -2.0; DL0[1] = -2.0;
4244 DL1[0] = 2.0; DL1[1] = 0.0;
4245 DL2[0] = 0.0; DL2[1] = 2.0;
4247 for (i = 0; i < 6; i++)
4248 for (j = 0; j < 2; j++)
4255 for (j = 0; j < 2; j++)
4257 dshape(0,j) = DL0[j];
4258 dshape(3,j) = DL1[j];
4259 dshape(5,j) = DL2[j];
4264 for (j = 0; j < 2; j++)
4266 dshape(3,j) = DL0[j];
4267 dshape(1,j) = DL1[j];
4268 dshape(4,j) = DL2[j];
4273 for (j = 0; j < 2; j++)
4275 dshape(5,j) = DL0[j];
4276 dshape(4,j) = DL1[j];
4277 dshape(2,j) = DL2[j];
4282 for (j = 0; j < 2; j++)
4284 dshape(3,j) = - DL2[j];
4285 dshape(4,j) = - DL0[j];
4286 dshape(5,j) = - DL1[j];
4331 double L0, L1, L2, L3, L4, L5;
4332 L0 = 2.0 * ( 1. - ip.
x - ip.
y - ip.
z );
4333 L1 = 2.0 * ( ip.
x );
4334 L2 = 2.0 * ( ip.
y );
4335 L3 = 2.0 * ( ip.
z );
4336 L4 = 2.0 * ( ip.
x + ip.
y );
4337 L5 = 2.0 * ( ip.
y + ip.
z );
4350 for (i = 0; i < 10; i++)
4357 shape(0) = L0 - 1.0;
4365 shape(1) = L1 - 1.0;
4373 shape(2) = L2 - 1.0;
4381 shape(3) = L3 - 1.0;
4383 else if ((L4 <= 1.0) && (L5 <= 1.0))
4385 shape(4) = 1.0 - L5;
4387 shape(6) = 1.0 - L4;
4388 shape(8) = 1.0 - L0;
4390 else if ((L4 >= 1.0) && (L5 <= 1.0))
4392 shape(4) = 1.0 - L5;
4393 shape(5) = 1.0 - L1;
4394 shape(7) = L4 - 1.0;
4397 else if ((L4 <= 1.0) && (L5 >= 1.0))
4399 shape(5) = 1.0 - L3;
4400 shape(6) = 1.0 - L4;
4402 shape(9) = L5 - 1.0;
4404 else if ((L4 >= 1.0) && (L5 >= 1.0))
4407 shape(7) = L4 - 1.0;
4408 shape(8) = 1.0 - L2;
4409 shape(9) = L5 - 1.0;
4418 double L0, L1, L2, L3, L4, L5;
4419 L0 = 2.0 * ( 1. - ip.
x - ip.
y - ip.
z );
4420 L1 = 2.0 * ( ip.
x );
4421 L2 = 2.0 * ( ip.
y );
4422 L3 = 2.0 * ( ip.
z );
4423 L4 = 2.0 * ( ip.
x + ip.
y );
4424 L5 = 2.0 * ( ip.
y + ip.
z );
4426 double DL0[3], DL1[3], DL2[3], DL3[3], DL4[3], DL5[3];
4427 DL0[0] = -2.0; DL0[1] = -2.0; DL0[2] = -2.0;
4428 DL1[0] = 2.0; DL1[1] = 0.0; DL1[2] = 0.0;
4429 DL2[0] = 0.0; DL2[1] = 2.0; DL2[2] = 0.0;
4430 DL3[0] = 0.0; DL3[1] = 0.0; DL3[2] = 2.0;
4431 DL4[0] = 2.0; DL4[1] = 2.0; DL4[2] = 0.0;
4432 DL5[0] = 0.0; DL5[1] = 2.0; DL5[2] = 2.0;
4434 for (i = 0; i < 10; i++)
4435 for (j = 0; j < 3; j++)
4442 for (j = 0; j < 3; j++)
4444 dshape(0,j) = DL0[j];
4445 dshape(4,j) = DL1[j];
4446 dshape(5,j) = DL2[j];
4447 dshape(6,j) = DL3[j];
4452 for (j = 0; j < 3; j++)
4454 dshape(4,j) = DL0[j];
4455 dshape(1,j) = DL1[j];
4456 dshape(7,j) = DL2[j];
4457 dshape(8,j) = DL3[j];
4462 for (j = 0; j < 3; j++)
4464 dshape(5,j) = DL0[j];
4465 dshape(7,j) = DL1[j];
4466 dshape(2,j) = DL2[j];
4467 dshape(9,j) = DL3[j];
4472 for (j = 0; j < 3; j++)
4474 dshape(6,j) = DL0[j];
4475 dshape(8,j) = DL1[j];
4476 dshape(9,j) = DL2[j];
4477 dshape(3,j) = DL3[j];
4480 else if ((L4 <= 1.0) && (L5 <= 1.0))
4482 for (j = 0; j < 3; j++)
4484 dshape(4,j) = - DL5[j];
4485 dshape(5,j) = DL2[j];
4486 dshape(6,j) = - DL4[j];
4487 dshape(8,j) = - DL0[j];
4490 else if ((L4 >= 1.0) && (L5 <= 1.0))
4492 for (j = 0; j < 3; j++)
4494 dshape(4,j) = - DL5[j];
4495 dshape(5,j) = - DL1[j];
4496 dshape(7,j) = DL4[j];
4497 dshape(8,j) = DL3[j];
4500 else if ((L4 <= 1.0) && (L5 >= 1.0))
4502 for (j = 0; j < 3; j++)
4504 dshape(5,j) = - DL3[j];
4505 dshape(6,j) = - DL4[j];
4506 dshape(8,j) = DL1[j];
4507 dshape(9,j) = DL5[j];
4510 else if ((L4 >= 1.0) && (L5 >= 1.0))
4512 for (j = 0; j < 3; j++)
4514 dshape(5,j) = DL0[j];
4515 dshape(7,j) = DL4[j];
4516 dshape(8,j) = - DL2[j];
4517 dshape(9,j) = DL5[j];
4550 double x = ip.
x, y = ip.
y;
4552 Lx = 2.0 * ( 1. - x );
4553 Ly = 2.0 * ( 1. - y );
4562 for (i = 0; i < 9; i++)
4567 if ((x <= 0.5) && (y <= 0.5))
4569 shape(0) = (Lx - 1.0) * (Ly - 1.0);
4570 shape(4) = (2.0 - Lx) * (Ly - 1.0);
4571 shape(8) = (2.0 - Lx) * (2.0 - Ly);
4572 shape(7) = (Lx - 1.0) * (2.0 - Ly);
4574 else if ((x >= 0.5) && (y <= 0.5))
4576 shape(4) = Lx * (Ly - 1.0);
4577 shape(1) = (1.0 - Lx) * (Ly - 1.0);
4578 shape(5) = (1.0 - Lx) * (2.0 - Ly);
4579 shape(8) = Lx * (2.0 - Ly);
4581 else if ((x >= 0.5) && (y >= 0.5))
4583 shape(8) = Lx * Ly ;
4584 shape(5) = (1.0 - Lx) * Ly ;
4585 shape(2) = (1.0 - Lx) * (1.0 - Ly);
4586 shape(6) = Lx * (1.0 - Ly);
4588 else if ((x <= 0.5) && (y >= 0.5))
4590 shape(7) = (Lx - 1.0) * Ly ;
4591 shape(8) = (2.0 - Lx) * Ly ;
4592 shape(6) = (2.0 - Lx) * (1.0 - Ly);
4593 shape(3) = (Lx - 1.0) * (1.0 - Ly);
4601 double x = ip.
x, y = ip.
y;
4603 Lx = 2.0 * ( 1. - x );
4604 Ly = 2.0 * ( 1. - y );
4606 for (i = 0; i < 9; i++)
4607 for (j = 0; j < 2; j++)
4612 if ((x <= 0.5) && (y <= 0.5))
4614 dshape(0,0) = 2.0 * (1.0 - Ly);
4615 dshape(0,1) = 2.0 * (1.0 - Lx);
4617 dshape(4,0) = 2.0 * (Ly - 1.0);
4618 dshape(4,1) = -2.0 * (2.0 - Lx);
4620 dshape(8,0) = 2.0 * (2.0 - Ly);
4621 dshape(8,1) = 2.0 * (2.0 - Lx);
4623 dshape(7,0) = -2.0 * (2.0 - Ly);
4624 dshape(7,0) = 2.0 * (Lx - 1.0);
4626 else if ((x >= 0.5) && (y <= 0.5))
4628 dshape(4,0) = -2.0 * (Ly - 1.0);
4629 dshape(4,1) = -2.0 * Lx;
4631 dshape(1,0) = 2.0 * (Ly - 1.0);
4632 dshape(1,1) = -2.0 * (1.0 - Lx);
4634 dshape(5,0) = 2.0 * (2.0 - Ly);
4635 dshape(5,1) = 2.0 * (1.0 - Lx);
4637 dshape(8,0) = -2.0 * (2.0 - Ly);
4638 dshape(8,1) = 2.0 * Lx;
4640 else if ((x >= 0.5) && (y >= 0.5))
4642 dshape(8,0) = -2.0 * Ly;
4643 dshape(8,1) = -2.0 * Lx;
4645 dshape(5,0) = 2.0 * Ly;
4646 dshape(5,1) = -2.0 * (1.0 - Lx);
4648 dshape(2,0) = 2.0 * (1.0 - Ly);
4649 dshape(2,1) = 2.0 * (1.0 - Lx);
4651 dshape(6,0) = -2.0 * (1.0 - Ly);
4652 dshape(6,1) = 2.0 * Lx;
4654 else if ((x <= 0.5) && (y >= 0.5))
4656 dshape(7,0) = -2.0 * Ly;
4657 dshape(7,1) = -2.0 * (Lx - 1.0);
4659 dshape(8,0) = 2.0 * Ly ;
4660 dshape(8,1) = -2.0 * (2.0 - Lx);
4662 dshape(6,0) = 2.0 * (1.0 - Ly);
4663 dshape(6,1) = 2.0 * (2.0 - Lx);
4665 dshape(3,0) = -2.0 * (1.0 - Ly);
4666 dshape(3,1) = 2.0 * (Lx - 1.0);
4677 I[ 0] = 0.0; J[ 0] = 0.0; K[ 0] = 0.0;
4678 I[ 1] = 1.0; J[ 1] = 0.0; K[ 1] = 0.0;
4679 I[ 2] = 1.0; J[ 2] = 1.0; K[ 2] = 0.0;
4680 I[ 3] = 0.0; J[ 3] = 1.0; K[ 3] = 0.0;
4681 I[ 4] = 0.0; J[ 4] = 0.0; K[ 4] = 1.0;
4682 I[ 5] = 1.0; J[ 5] = 0.0; K[ 5] = 1.0;
4683 I[ 6] = 1.0; J[ 6] = 1.0; K[ 6] = 1.0;
4684 I[ 7] = 0.0; J[ 7] = 1.0; K[ 7] = 1.0;
4686 I[ 8] = 0.5; J[ 8] = 0.0; K[ 8] = 0.0;
4687 I[ 9] = 1.0; J[ 9] = 0.5; K[ 9] = 0.0;
4688 I[10] = 0.5; J[10] = 1.0; K[10] = 0.0;
4689 I[11] = 0.0; J[11] = 0.5; K[11] = 0.0;
4690 I[12] = 0.5; J[12] = 0.0; K[12] = 1.0;
4691 I[13] = 1.0; J[13] = 0.5; K[13] = 1.0;
4692 I[14] = 0.5; J[14] = 1.0; K[14] = 1.0;
4693 I[15] = 0.0; J[15] = 0.5; K[15] = 1.0;
4694 I[16] = 0.0; J[16] = 0.0; K[16] = 0.5;
4695 I[17] = 1.0; J[17] = 0.0; K[17] = 0.5;
4696 I[18] = 1.0; J[18] = 1.0; K[18] = 0.5;
4697 I[19] = 0.0; J[19] = 1.0; K[19] = 0.5;
4699 I[20] = 0.5; J[20] = 0.5; K[20] = 0.0;
4700 I[21] = 0.5; J[21] = 0.0; K[21] = 0.5;
4701 I[22] = 1.0; J[22] = 0.5; K[22] = 0.5;
4702 I[23] = 0.5; J[23] = 1.0; K[23] = 0.5;
4703 I[24] = 0.0; J[24] = 0.5; K[24] = 0.5;
4704 I[25] = 0.5; J[25] = 0.5; K[25] = 1.0;
4706 I[26] = 0.5; J[26] = 0.5; K[26] = 0.5;
4708 for (
int n = 0; n < 27; n++)
4721 double x = ip.
x, y = ip.
y, z = ip.
z;
4723 for (i = 0; i < 27; i++)
4728 if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5))
4743 else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5))
4758 else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5))
4773 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5))
4788 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5))
4803 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5))
4818 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5))
4849 shape(N[0]) = Lx * Ly * Lz;
4850 shape(N[1]) = (1 - Lx) * Ly * Lz;
4851 shape(N[2]) = (1 - Lx) * (1 - Ly) * Lz;
4852 shape(N[3]) = Lx * (1 - Ly) * Lz;
4853 shape(N[4]) = Lx * Ly * (1 - Lz);
4854 shape(N[5]) = (1 - Lx) * Ly * (1 - Lz);
4855 shape(N[6]) = (1 - Lx) * (1 - Ly) * (1 - Lz);
4856 shape(N[7]) = Lx * (1 - Ly) * (1 - Lz);
4864 double x = ip.
x, y = ip.
y, z = ip.
z;
4866 for (i = 0; i < 27; i++)
4867 for (j = 0; j < 3; j++)
4872 if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5))
4887 else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5))
4902 else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5))
4917 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5))
4932 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5))
4947 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5))
4962 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5))
4993 dshape(N[0],0) = -2.0 * Ly * Lz ;
4994 dshape(N[0],1) = -2.0 * Lx * Lz ;
4995 dshape(N[0],2) = -2.0 * Lx * Ly ;
4997 dshape(N[1],0) = 2.0 * Ly * Lz ;
4998 dshape(N[1],1) = -2.0 * (1 - Lx) * Lz ;
4999 dshape(N[1],2) = -2.0 * (1 - Lx) * Ly ;
5001 dshape(N[2],0) = 2.0 * (1 - Ly) * Lz ;
5002 dshape(N[2],1) = 2.0 * (1 - Lx) * Lz ;
5003 dshape(N[2],2) = -2.0 * (1 - Lx) * (1 - Ly);
5005 dshape(N[3],0) = -2.0 * (1 - Ly) * Lz ;
5006 dshape(N[3],1) = 2.0 * Lx * Lz ;
5007 dshape(N[3],2) = -2.0 * Lx * (1 - Ly);
5009 dshape(N[4],0) = -2.0 * Ly * (1 - Lz);
5010 dshape(N[4],1) = -2.0 * Lx * (1 - Lz);
5011 dshape(N[4],2) = 2.0 * Lx * Ly ;
5013 dshape(N[5],0) = 2.0 * Ly * (1 - Lz);
5014 dshape(N[5],1) = -2.0 * (1 - Lx) * (1 - Lz);
5015 dshape(N[5],2) = 2.0 * (1 - Lx) * Ly ;
5017 dshape(N[6],0) = 2.0 * (1 - Ly) * (1 - Lz);
5018 dshape(N[6],1) = 2.0 * (1 - Lx) * (1 - Lz);
5019 dshape(N[6],2) = 2.0 * (1 - Lx) * (1 - Ly);
5021 dshape(N[7],0) = -2.0 * (1 - Ly) * (1 - Lz);
5022 dshape(N[7],1) = 2.0 * Lx * (1 - Lz);
5023 dshape(N[7],2) = 2.0 * Lx * (1 - Ly);
5083 double x = ip.
x, y = ip.
y, z = ip.
z;
5085 shape(0,0) = (1. - y) * (1. - z);
5089 shape(2,0) = y * (1. - z);
5093 shape(4,0) = z * (1. - y);
5102 shape(1,1) = x * (1. - z);
5106 shape(3,1) = (1. - x) * (1. - z);
5114 shape(7,1) = (1. - x) * z;
5119 shape(8,2) = (1. - x) * (1. - y);
5123 shape(9,2) = x * (1. - y);
5127 shape(10,2) = x * y;
5131 shape(11,2) = y * (1. - x);
5139 double x = ip.
x, y = ip.
y, z = ip.
z;
5141 curl_shape(0,0) = 0.;
5142 curl_shape(0,1) = y - 1.;
5143 curl_shape(0,2) = 1. - z;
5145 curl_shape(2,0) = 0.;
5146 curl_shape(2,1) = -y;
5147 curl_shape(2,2) = z - 1.;
5149 curl_shape(4,0) = 0;
5150 curl_shape(4,1) = 1. - y;
5151 curl_shape(4,2) = z;
5153 curl_shape(6,0) = 0.;
5154 curl_shape(6,1) = y;
5155 curl_shape(6,2) = -z;
5157 curl_shape(1,0) = x;
5158 curl_shape(1,1) = 0.;
5159 curl_shape(1,2) = 1. - z;
5161 curl_shape(3,0) = 1. - x;
5162 curl_shape(3,1) = 0.;
5163 curl_shape(3,2) = z - 1.;
5165 curl_shape(5,0) = -x;
5166 curl_shape(5,1) = 0.;
5167 curl_shape(5,2) = z;
5169 curl_shape(7,0) = x - 1.;
5170 curl_shape(7,1) = 0.;
5171 curl_shape(7,2) = -z;
5173 curl_shape(8,0) = x - 1.;
5174 curl_shape(8,1) = 1. - y;
5175 curl_shape(8,2) = 0.;
5177 curl_shape(9,0) = -x;
5178 curl_shape(9,1) = y - 1.;
5179 curl_shape(9,2) = 0;
5181 curl_shape(10,0) = x;
5182 curl_shape(10,1) = -y;
5183 curl_shape(10,2) = 0.;
5185 curl_shape(11,0) = 1. - x;
5186 curl_shape(11,1) = y;
5187 curl_shape(11,2) = 0.;
5190 const double Nedelec1HexFiniteElement::tk[12][3] =
5192 {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
5193 {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
5194 {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}
5201 #ifdef MFEM_THREAD_SAFE
5206 for (k = 0; k < 12; k++)
5209 for (j = 0; j < 12; j++)
5211 double d = (
vshape(j,0)*tk[k][0] +
vshape(j,1)*tk[k][1] +
5213 if (j == k) { d -= 1.0; }
5214 if (fabs(d) > 1.0e-12)
5216 mfem::err <<
"Nedelec1HexFiniteElement::GetLocalInterpolation (...)\n"
5217 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5225 ip.
x = ip.
y = ip.
z = 0.0;
5232 for (k = 0; k < 12; k++)
5235 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5238 vk[0] =
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2];
5239 vk[1] =
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2];
5240 vk[2] =
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2];
5241 for (j = 0; j < 12; j++)
5242 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5243 vshape(j,2)*vk[2])) < 1.0e-12)
5257 for (
int k = 0; k < 12; k++)
5265 vk[0] * (
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2] ) +
5266 vk[1] * (
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2] ) +
5267 vk[2] * (
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2] );
5304 double x = ip.
x, y = ip.
y, z = ip.
z;
5306 shape(0,0) = 1. - y - z;
5311 shape(1,1) = 1. - x - z;
5316 shape(2,2) = 1. - x - y;
5335 curl_shape(0,0) = 0.;
5336 curl_shape(0,1) = -2.;
5337 curl_shape(0,2) = 2.;
5339 curl_shape(1,0) = 2.;
5340 curl_shape(1,1) = 0.;
5341 curl_shape(1,2) = -2.;
5343 curl_shape(2,0) = -2.;
5344 curl_shape(2,1) = 2.;
5345 curl_shape(2,2) = 0.;
5347 curl_shape(3,0) = 0.;
5348 curl_shape(3,1) = 0.;
5349 curl_shape(3,2) = 2.;
5351 curl_shape(4,0) = 0.;
5352 curl_shape(4,1) = -2.;
5353 curl_shape(4,2) = 0.;
5355 curl_shape(5,0) = 2.;
5356 curl_shape(5,1) = 0.;
5357 curl_shape(5,2) = 0.;
5360 const double Nedelec1TetFiniteElement::tk[6][3] =
5361 {{1,0,0}, {0,1,0}, {0,0,1}, {-1,1,0}, {-1,0,1}, {0,-1,1}};
5367 #ifdef MFEM_THREAD_SAFE
5372 for (k = 0; k < 6; k++)
5375 for (j = 0; j < 6; j++)
5377 double d = (
vshape(j,0)*tk[k][0] +
vshape(j,1)*tk[k][1] +
5379 if (j == k) { d -= 1.0; }
5380 if (fabs(d) > 1.0e-12)
5382 mfem::err <<
"Nedelec1TetFiniteElement::GetLocalInterpolation (...)\n"
5383 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5391 ip.
x = ip.
y = ip.
z = 0.0;
5398 for (k = 0; k < 6; k++)
5401 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5404 vk[0] =
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2];
5405 vk[1] =
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2];
5406 vk[2] =
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2];
5407 for (j = 0; j < 6; j++)
5408 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5409 vshape(j,2)*vk[2])) < 1.0e-12)
5423 for (
int k = 0; k < 6; k++)
5431 vk[0] * (
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2] ) +
5432 vk[1] * (
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2] ) +
5433 vk[2] * (
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2] );
5470 double x = ip.
x, y = ip.
y, z = ip.
z;
5474 shape(0,2) = z - 1.;
5477 shape(1,1) = y - 1.;
5488 shape(4,0) = x - 1.;
5508 const double RT0HexFiniteElement::nk[6][3] =
5509 {{0,0,-1}, {0,-1,0}, {1,0,0}, {0,1,0}, {-1,0,0}, {0,0,1}};
5515 #ifdef MFEM_THREAD_SAFE
5521 for (k = 0; k < 6; k++)
5524 for (j = 0; j < 6; j++)
5526 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5528 if (j == k) { d -= 1.0; }
5529 if (fabs(d) > 1.0e-12)
5531 mfem::err <<
"RT0HexFiniteElement::GetLocalInterpolation (...)\n"
5532 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5540 ip.
x = ip.
y = ip.
z = 0.0;
5548 for (k = 0; k < 6; k++)
5551 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5554 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2];
5555 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2];
5556 vk[2] =
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2];
5557 for (j = 0; j < 6; j++)
5558 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5559 vshape(j,2)*vk[2])) < 1.0e-12)
5572 #ifdef MFEM_THREAD_SAFE
5576 for (
int k = 0; k < 6; k++)
5585 vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2] ) +
5586 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2] ) +
5587 vk[2] * (
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2] );
5716 double x = ip.
x, y = ip.
y, z = ip.
z;
5720 shape(2,2) = -(1. - 3.*z + 2.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5723 shape(3,2) = -(1. - 3.*z + 2.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5726 shape(0,2) = -(1. - 3.*z + 2.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5729 shape(1,2) = -(1. - 3.*z + 2.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5732 shape(4,1) = -(1. - 3.*y + 2.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5735 shape(5,1) = -(1. - 3.*y + 2.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5738 shape(6,1) = -(1. - 3.*y + 2.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5741 shape(7,1) = -(1. - 3.*y + 2.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5744 shape(8,0) = (-x + 2.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5747 shape(9,0) = (-x + 2.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5750 shape(10,0) = (-x + 2.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5753 shape(11,0) = (-x + 2.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5758 shape(13,1) = (-y + 2.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5761 shape(12,1) = (-y + 2.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5764 shape(15,1) = (-y + 2.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5767 shape(14,1) = (-y + 2.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5770 shape(17,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5773 shape(16,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5776 shape(19,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5779 shape(18,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5785 shape(20,2) = (-z + 2.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5788 shape(21,2) = (-z + 2.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5791 shape(22,2) = (-z + 2.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5794 shape(23,2) = (-z + 2.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5796 shape(24,0) = (4.*x - 4.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5799 shape(25,0) = (4.*x - 4.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5802 shape(26,0) = (4.*x - 4.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5805 shape(27,0) = (4.*x - 4.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5810 shape(28,1) = (4.*y - 4.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5813 shape(29,1) = (4.*y - 4.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5816 shape(30,1) = (4.*y - 4.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5819 shape(31,1) = (4.*y - 4.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5824 shape(32,2) = (4.*z - 4.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5827 shape(33,2) = (4.*z - 4.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5830 shape(34,2) = (4.*z - 4.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5833 shape(35,2) = (4.*z - 4.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5839 double x = ip.
x, y = ip.
y, z = ip.
z;
5841 divshape(2) = -(-3. + 4.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5842 divshape(3) = -(-3. + 4.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5843 divshape(0) = -(-3. + 4.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5844 divshape(1) = -(-3. + 4.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5846 divshape(4) = -(-3. + 4.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5847 divshape(5) = -(-3. + 4.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5848 divshape(6) = -(-3. + 4.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5849 divshape(7) = -(-3. + 4.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5851 divshape(8) = (-1. + 4.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5852 divshape(9) = (-1. + 4.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5853 divshape(10) = (-1. + 4.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5854 divshape(11) = (-1. + 4.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5856 divshape(13) = (-1. + 4.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5857 divshape(12) = (-1. + 4.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5858 divshape(15) = (-1. + 4.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5859 divshape(14) = (-1. + 4.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5861 divshape(17) = -(-3. + 4.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5862 divshape(16) = -(-3. + 4.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5863 divshape(19) = -(-3. + 4.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5864 divshape(18) = -(-3. + 4.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5866 divshape(20) = (-1. + 4.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5867 divshape(21) = (-1. + 4.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5868 divshape(22) = (-1. + 4.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5869 divshape(23) = (-1. + 4.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5871 divshape(24) = ( 4. - 8.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5872 divshape(25) = ( 4. - 8.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5873 divshape(26) = ( 4. - 8.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5874 divshape(27) = ( 4. - 8.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5876 divshape(28) = ( 4. - 8.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5877 divshape(29) = ( 4. - 8.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5878 divshape(30) = ( 4. - 8.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5879 divshape(31) = ( 4. - 8.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5881 divshape(32) = ( 4. - 8.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5882 divshape(33) = ( 4. - 8.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5883 divshape(34) = ( 4. - 8.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5884 divshape(35) = ( 4. - 8.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5887 const double RT1HexFiniteElement::nk[36][3] =
5889 {0, 0,-1}, {0, 0,-1}, {0, 0,-1}, {0, 0,-1},
5890 {0,-1, 0}, {0,-1, 0}, {0,-1, 0}, {0,-1, 0},
5891 {1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0},
5892 {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0},
5893 {-1,0, 0}, {-1,0, 0}, {-1,0, 0}, {-1,0, 0},
5894 {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1},
5895 {1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0},
5896 {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0},
5897 {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1}
5904 #ifdef MFEM_THREAD_SAFE
5910 for (k = 0; k < 36; k++)
5913 for (j = 0; j < 36; j++)
5915 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5917 if (j == k) { d -= 1.0; }
5918 if (fabs(d) > 1.0e-12)
5920 mfem::err <<
"RT0HexFiniteElement::GetLocalInterpolation (...)\n"
5921 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5929 ip.
x = ip.
y = ip.
z = 0.0;
5937 for (k = 0; k < 36; k++)
5940 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5943 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2];
5944 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2];
5945 vk[2] =
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2];
5946 for (j = 0; j < 36; j++)
5947 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5948 vshape(j,2)*vk[2])) < 1.0e-12)
5961 #ifdef MFEM_THREAD_SAFE
5965 for (
int k = 0; k < 36; k++)
5974 vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2] ) +
5975 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2] ) +
5976 vk[2] * (
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2] );
6004 double x2 = 2.0*ip.
x, y2 = 2.0*ip.
y, z2 = 2.0*ip.
z;
6010 shape(1,0) = x2 - 2.0;
6015 shape(2,1) = y2 - 2.0;
6020 shape(3,2) = z2 - 2.0;
6032 const double RT0TetFiniteElement::nk[4][3] =
6033 {{.5,.5,.5}, {-.5,0,0}, {0,-.5,0}, {0,0,-.5}};
6039 #ifdef MFEM_THREAD_SAFE
6045 for (k = 0; k < 4; k++)
6048 for (j = 0; j < 4; j++)
6050 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
6052 if (j == k) { d -= 1.0; }
6053 if (fabs(d) > 1.0e-12)
6055 mfem::err <<
"RT0TetFiniteElement::GetLocalInterpolation (...)\n"
6056 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
6064 ip.
x = ip.
y = ip.
z = 0.0;
6072 for (k = 0; k < 4; k++)
6075 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
6078 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2];
6079 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2];
6080 vk[2] =
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2];
6081 for (j = 0; j < 4; j++)
6082 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
6083 vshape(j,2)*vk[2])) < 1.0e-12)
6096 #ifdef MFEM_THREAD_SAFE
6100 for (
int k = 0; k < 4; k++)
6109 vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2] ) +
6110 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2] ) +
6111 vk[2] * (
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2] );
6146 double x = 2. * ip.
x - 1.;
6147 double y = 2. * ip.
y - 1.;
6148 double z = 2. * ip.
z - 1.;
6149 double f5 = x * x - y * y;
6150 double f6 = y * y - z * z;
6152 shape(0) = (1./6.) * (1. - 3. * z - f5 - 2. * f6);
6153 shape(1) = (1./6.) * (1. - 3. * y - f5 + f6);
6154 shape(2) = (1./6.) * (1. + 3. * x + 2. * f5 + f6);
6155 shape(3) = (1./6.) * (1. + 3. * y - f5 + f6);
6156 shape(4) = (1./6.) * (1. - 3. * x + 2. * f5 + f6);
6157 shape(5) = (1./6.) * (1. + 3. * z - f5 - 2. * f6);
6163 const double a = 2./3.;
6165 double xt = a * (1. - 2. * ip.
x);
6166 double yt = a * (1. - 2. * ip.
y);
6167 double zt = a * (1. - 2. * ip.
z);
6171 dshape(0,2) = -1. - 2. * zt;
6174 dshape(1,1) = -1. - 2. * yt;
6177 dshape(2,0) = 1. - 2. * xt;
6182 dshape(3,1) = 1. - 2. * yt;
6185 dshape(4,0) = -1. - 2. * xt;
6191 dshape(5,2) = 1. - 2. * zt;
6205 for (
int i = 0; i <= p; i++)
6219 for (
int i = 0; i <= p; i++)
6221 for (
int j = 0; j < i; j++)
6223 double xij = x(i) - x(j);
6228 for (
int i = 0; i <= p; i++)
6235 for (
int i = 0; i < p; i++)
6239 mfem_error(
"Poly_1D::Basis::Basis : nodes are not increasing!");
6265 int i, k, p = x.Size() - 1;
6275 for (k = 0; k < p; k++)
6277 if (y >= (x(k) + x(k+1))/2)
6283 for (i = k+1; i <= p; i++)
6290 l = lk * (y - x(k));
6292 for (i = 0; i < k; i++)
6294 u(i) = l * w(i) / (y - x(i));
6297 for (i++; i <= p; i++)
6299 u(i) = l * w(i) / (y - x(i));
6324 int i, k, p = x.Size() - 1;
6325 double l, lp, lk, sk, si;
6335 for (k = 0; k < p; k++)
6337 if (y >= (x(k) + x(k+1))/2)
6343 for (i = k+1; i <= p; i++)
6350 l = lk * (y - x(k));
6353 for (i = 0; i < k; i++)
6355 si = 1.0/(y - x(i));
6357 u(i) = l * si * w(i);
6360 for (i++; i <= p; i++)
6362 si = 1.0/(y - x(i));
6364 u(i) = l * si * w(i);
6368 for (i = 0; i < k; i++)
6370 d(i) = (lp * w(i) - u(i))/(y - x(i));
6373 for (i++; i <= p; i++)
6375 d(i) = (lp * w(i) - u(i))/(y - x(i));
6392 for (
int i = 0; i <= p; i++)
6394 binom(i,0) = binom(i,i) = 1;
6395 for (
int j = 1; j < i; j++)
6397 binom(i,j) = binom(i-1,j) + binom(i-1,j-1);
6406 for (
int i = 0; i <= p; i++)
6409 double s = sin(M_PI_2*(i + 0.5)/(p + 1));
6414 void Poly_1D::CalcMono(
const int p,
const double x,
double *u)
6418 for (
int n = 1; n <= p; n++)
6424 void Poly_1D::CalcMono(
const int p,
const double x,
double *u,
double *d)
6429 for (
int n = 1; n <= p; n++)
6446 const int *b =
Binom(p);
6449 for (i = 1; i < p; i++)
6456 for (i--; i > 0; i--)
6466 double *u,
double *d)
6476 const int *b =
Binom(p);
6477 const double xpy = x + y, ptx = p*x;
6480 for (i = 1; i < p; i++)
6482 d[i] = b[i]*z*(i*xpy - ptx);
6489 for (i--; i > 0; i--)
6510 const int *b =
Binom(p);
6511 const double xpy = x + y, ptx = p*x;
6514 for (i = 1; i < p; i++)
6516 d[i] = b[i]*z*(i*xpy - ptx);
6521 for (i--; i > 0; i--)
6530 void Poly_1D::CalcLegendre(
const int p,
const double x,
double *u)
6536 if (p == 0) {
return; }
6537 u[1] = z = 2.*x - 1.;
6538 for (
int n = 1; n < p; n++)
6540 u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
6544 void Poly_1D::CalcLegendre(
const int p,
const double x,
double *u,
double *d)
6553 if (p == 0) {
return; }
6554 u[1] = z = 2.*x - 1.;
6556 for (
int n = 1; n < p; n++)
6558 u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
6559 d[n+1] = (4*n + 2)*u[n] + d[n-1];
6563 void Poly_1D::CalcChebyshev(
const int p,
const double x,
double *u)
6570 if (p == 0) {
return; }
6571 u[1] = z = 2.*x - 1.;
6572 for (
int n = 1; n < p; n++)
6574 u[n+1] = 2*z*u[n] - u[n-1];
6578 void Poly_1D::CalcChebyshev(
const int p,
const double x,
double *u,
double *d)
6591 if (p == 0) {
return; }
6592 u[1] = z = 2.*x - 1.;
6594 for (
int n = 1; n < p; n++)
6596 u[n+1] = 2*z*u[n] - u[n-1];
6597 d[n+1] = (n + 1)*(z*d[n]/n + 2*u[n]);
6608 if (points_container.find(btype) == points_container.end())
6613 if (pts.
Size() <= p)
6619 pts[p] =
new double[p + 1];
6629 if ( bases_container.find(btype) == bases_container.end() )
6635 if (bases.
Size() <= p)
6639 if (bases[p] == NULL)
6649 for (PointsMap::iterator it = points_container.begin();
6650 it != points_container.end() ; ++it)
6653 for (
int i = 0 ; i < pts.
Size() ; ++i )
6660 for (BasisMap::iterator it = bases_container.begin();
6661 it != bases_container.end() ; ++it)
6664 for (
int i = 0 ; i < bases.
Size() ; ++i )
6690 for (
int i = 1; i < p; i++)
6698 const int p1 = p + 1;
6709 for (
int i = 1; i < p; i++)
6713 for (
int i = 1; i < p; i++)
6717 for (
int i = 1; i < p; i++)
6721 for (
int i = 1; i < p; i++)
6727 for (
int j = 1; j < p; j++)
6729 for (
int i = 1; i < p; i++)
6738 const int p1 = p + 1;
6742 dof_map[0 + (0 + 0*p1)*p1] = 0;
6743 dof_map[p + (0 + 0*p1)*p1] = 1;
6744 dof_map[p + (p + 0*p1)*p1] = 2;
6745 dof_map[0 + (p + 0*p1)*p1] = 3;
6746 dof_map[0 + (0 + p*p1)*p1] = 4;
6747 dof_map[p + (0 + p*p1)*p1] = 5;
6748 dof_map[p + (p + p*p1)*p1] = 6;
6749 dof_map[0 + (p + p*p1)*p1] = 7;
6753 for (
int i = 1; i < p; i++)
6755 dof_map[i + (0 + 0*p1)*p1] = o++;
6757 for (
int i = 1; i < p; i++)
6759 dof_map[p + (i + 0*p1)*p1] = o++;
6761 for (
int i = 1; i < p; i++)
6763 dof_map[i + (p + 0*p1)*p1] = o++;
6765 for (
int i = 1; i < p; i++)
6767 dof_map[0 + (i + 0*p1)*p1] = o++;
6769 for (
int i = 1; i < p; i++)
6771 dof_map[i + (0 + p*p1)*p1] = o++;
6773 for (
int i = 1; i < p; i++)
6775 dof_map[p + (i + p*p1)*p1] = o++;
6777 for (
int i = 1; i < p; i++)
6779 dof_map[i + (p + p*p1)*p1] = o++;
6781 for (
int i = 1; i < p; i++)
6783 dof_map[0 + (i + p*p1)*p1] = o++;
6785 for (
int i = 1; i < p; i++)
6787 dof_map[0 + (0 + i*p1)*p1] = o++;
6789 for (
int i = 1; i < p; i++)
6791 dof_map[p + (0 + i*p1)*p1] = o++;
6793 for (
int i = 1; i < p; i++)
6795 dof_map[p + (p + i*p1)*p1] = o++;
6797 for (
int i = 1; i < p; i++)
6799 dof_map[0 + (p + i*p1)*p1] = o++;
6803 for (
int j = 1; j < p; j++)
6805 for (
int i = 1; i < p; i++)
6807 dof_map[i + ((p-j) + 0*p1)*p1] = o++;
6810 for (
int j = 1; j < p; j++)
6812 for (
int i = 1; i < p; i++)
6814 dof_map[i + (0 + j*p1)*p1] = o++;
6817 for (
int j = 1; j < p; j++)
6819 for (
int i = 1; i < p; i++)
6821 dof_map[p + (i + j*p1)*p1] = o++;
6824 for (
int j = 1; j < p; j++)
6826 for (
int i = 1; i < p; i++)
6828 dof_map[(p-i) + (p + j*p1)*p1] = o++;
6831 for (
int j = 1; j < p; j++)
6833 for (
int i = 1; i < p; i++)
6835 dof_map[0 + ((p-i) + j*p1)*p1] = o++;
6838 for (
int j = 1; j < p; j++)
6840 for (
int i = 1; i < p; i++)
6842 dof_map[i + (j + p*p1)*p1] = o++;
6847 for (
int k = 1; k < p; k++)
6849 for (
int j = 1; j < p; j++)
6851 for (
int i = 1; i < p; i++)
6853 dof_map[i + (j + k*p1)*p1] = o++;
6860 MFEM_ABORT(
"invalid dimension: " << dims);
6871 MFEM_ABORT(
"invalid DofMapType: " << dmtype);
6886 const int dims,
const int p,
const DofMapType dmtype)
6888 Pow(p + 1, dims), p,
6898 #ifndef MFEM_THREAD_SAFE
6905 for (
int i = 1; i < p; i++)
6914 const int p =
Order;
6916 #ifdef MFEM_THREAD_SAFE
6922 shape(0) = shape_x(0);
6923 shape(1) = shape_x(p);
6924 for (
int i = 1; i < p; i++)
6926 shape(i+1) = shape_x(i);
6933 const int p =
Order;
6935 #ifdef MFEM_THREAD_SAFE
6936 Vector shape_x(p+1), dshape_x(p+1);
6941 dshape(0,0) = dshape_x(0);
6942 dshape(1,0) = dshape_x(p);
6943 for (
int i = 1; i < p; i++)
6945 dshape(i+1,0) = dshape_x(i);
6951 const int p =
Order;
6959 for (
int i = 1; i < p; i++)
6968 for (
int i = 1; i < p; i++)
6982 #ifndef MFEM_THREAD_SAFE
6983 const int p1 = p + 1;
6992 for (
int j = 0; j <= p; j++)
6994 for (
int i = 0; i <= p; i++)
7004 const int p =
Order;
7006 #ifdef MFEM_THREAD_SAFE
7007 Vector shape_x(p+1), shape_y(p+1);
7013 for (
int o = 0, j = 0; j <= p; j++)
7014 for (
int i = 0; i <= p; i++)
7016 shape(
dof_map[o++]) = shape_x(i)*shape_y(j);
7023 const int p =
Order;
7025 #ifdef MFEM_THREAD_SAFE
7026 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
7032 for (
int o = 0, j = 0; j <= p; j++)
7034 for (
int i = 0; i <= p; i++)
7036 dshape(
dof_map[o],0) = dshape_x(i)* shape_y(j);
7037 dshape(
dof_map[o],1) = shape_x(i)*dshape_y(j); o++;
7044 const int p =
Order;
7047 #ifdef MFEM_THREAD_SAFE
7048 Vector shape_x(p+1), shape_y(p+1);
7051 for (
int i = 0; i <= p; i++)
7060 for (
int o = 0, j = 0; j <= p; j++)
7061 for (
int i = 0; i <= p; i++)
7063 dofs(
dof_map[o++]) = shape_x(i)*shape_x(j);
7067 for (
int o = 0, j = 0; j <= p; j++)
7068 for (
int i = 0; i <= p; i++)
7070 dofs(
dof_map[o++]) = shape_y(i)*shape_x(j);
7074 for (
int o = 0, j = 0; j <= p; j++)
7075 for (
int i = 0; i <= p; i++)
7077 dofs(
dof_map[o++]) = shape_y(i)*shape_y(j);
7081 for (
int o = 0, j = 0; j <= p; j++)
7082 for (
int i = 0; i <= p; i++)
7084 dofs(
dof_map[o++]) = shape_x(i)*shape_y(j);
7096 #ifndef MFEM_THREAD_SAFE
7097 const int p1 = p + 1;
7108 for (
int k = 0; k <= p; k++)
7109 for (
int j = 0; j <= p; j++)
7110 for (
int i = 0; i <= p; i++)
7119 const int p =
Order;
7121 #ifdef MFEM_THREAD_SAFE
7122 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7129 for (
int o = 0, k = 0; k <= p; k++)
7130 for (
int j = 0; j <= p; j++)
7131 for (
int i = 0; i <= p; i++)
7133 shape(
dof_map[o++]) = shape_x(i)*shape_y(j)*shape_z(k);
7140 const int p =
Order;
7142 #ifdef MFEM_THREAD_SAFE
7143 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7144 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
7151 for (
int o = 0, k = 0; k <= p; k++)
7152 for (
int j = 0; j <= p; j++)
7153 for (
int i = 0; i <= p; i++)
7155 dshape(
dof_map[o],0) = dshape_x(i)* shape_y(j)* shape_z(k);
7156 dshape(
dof_map[o],1) = shape_x(i)*dshape_y(j)* shape_z(k);
7157 dshape(
dof_map[o],2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
7163 const int p =
Order;
7166 #ifdef MFEM_THREAD_SAFE
7167 Vector shape_x(p+1), shape_y(p+1);
7170 for (
int i = 0; i <= p; i++)
7179 for (
int o = 0, k = 0; k <= p; k++)
7180 for (
int j = 0; j <= p; j++)
7181 for (
int i = 0; i <= p; i++)
7183 dofs(
dof_map[o++]) = shape_x(i)*shape_x(j)*shape_x(k);
7187 for (
int o = 0, k = 0; k <= p; k++)
7188 for (
int j = 0; j <= p; j++)
7189 for (
int i = 0; i <= p; i++)
7191 dofs(
dof_map[o++]) = shape_y(i)*shape_x(j)*shape_x(k);
7195 for (
int o = 0, k = 0; k <= p; k++)
7196 for (
int j = 0; j <= p; j++)
7197 for (
int i = 0; i <= p; i++)
7199 dofs(
dof_map[o++]) = shape_y(i)*shape_y(j)*shape_x(k);
7203 for (
int o = 0, k = 0; k <= p; k++)
7204 for (
int j = 0; j <= p; j++)
7205 for (
int i = 0; i <= p; i++)
7207 dofs(
dof_map[o++]) = shape_x(i)*shape_y(j)*shape_x(k);
7211 for (
int o = 0, k = 0; k <= p; k++)
7212 for (
int j = 0; j <= p; j++)
7213 for (
int i = 0; i <= p; i++)
7215 dofs(
dof_map[o++]) = shape_x(i)*shape_x(j)*shape_y(k);
7219 for (
int o = 0, k = 0; k <= p; k++)
7220 for (
int j = 0; j <= p; j++)
7221 for (
int i = 0; i <= p; i++)
7223 dofs(
dof_map[o++]) = shape_y(i)*shape_x(j)*shape_y(k);
7227 for (
int o = 0, k = 0; k <= p; k++)
7228 for (
int j = 0; j <= p; j++)
7229 for (
int i = 0; i <= p; i++)
7231 dofs(
dof_map[o++]) = shape_y(i)*shape_y(j)*shape_y(k);
7235 for (
int o = 0, k = 0; k <= p; k++)
7236 for (
int j = 0; j <= p; j++)
7237 for (
int i = 0; i <= p; i++)
7239 dofs(
dof_map[o++]) = shape_x(i)*shape_y(j)*shape_y(k);
7249 #ifndef MFEM_THREAD_SAFE
7258 for (
int i = 1; i < p; i++)
7267 const int p =
Order;
7269 #ifdef MFEM_THREAD_SAFE
7276 shape(0) = shape_x(0);
7277 shape(1) = shape_x(p);
7278 for (
int i = 1; i < p; i++)
7280 shape(i+1) = shape_x(i);
7287 const int p =
Order;
7289 #ifdef MFEM_THREAD_SAFE
7290 Vector shape_x(p+1), dshape_x(p+1);
7296 dshape(0,0) = dshape_x(0);
7297 dshape(1,0) = dshape_x(p);
7298 for (
int i = 1; i < p; i++)
7300 dshape(i+1,0) = dshape_x(i);
7314 #ifndef MFEM_THREAD_SAFE
7315 const int p1 = p + 1;
7324 for (
int j = 0; j <= p; j++)
7325 for (
int i = 0; i <= p; i++)
7334 const int p =
Order;
7336 #ifdef MFEM_THREAD_SAFE
7337 Vector shape_x(p+1), shape_y(p+1);
7344 for (
int o = 0, j = 0; j <= p; j++)
7345 for (
int i = 0; i <= p; i++)
7347 shape(
dof_map[o++]) = shape_x(i)*shape_y(j);
7354 const int p =
Order;
7356 #ifdef MFEM_THREAD_SAFE
7357 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
7364 for (
int o = 0, j = 0; j <= p; j++)
7365 for (
int i = 0; i <= p; i++)
7367 dshape(
dof_map[o],0) = dshape_x(i)* shape_y(j);
7368 dshape(
dof_map[o],1) = shape_x(i)*dshape_y(j); o++;
7382 #ifndef MFEM_THREAD_SAFE
7383 const int p1 = p + 1;
7394 for (
int k = 0; k <= p; k++)
7395 for (
int j = 0; j <= p; j++)
7396 for (
int i = 0; i <= p; i++)
7404 const int p =
Order;
7406 #ifdef MFEM_THREAD_SAFE
7407 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7414 for (
int o = 0, k = 0; k <= p; k++)
7415 for (
int j = 0; j <= p; j++)
7416 for (
int i = 0; i <= p; i++)
7418 shape(
dof_map[o++]) = shape_x(i)*shape_y(j)*shape_z(k);
7425 const int p =
Order;
7427 #ifdef MFEM_THREAD_SAFE
7428 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7429 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
7436 for (
int o = 0, k = 0; k <= p; k++)
7437 for (
int j = 0; j <= p; j++)
7438 for (
int i = 0; i <= p; i++)
7440 dshape(
dof_map[o],0) = dshape_x(i)* shape_y(j)* shape_z(k);
7441 dshape(
dof_map[o],1) = shape_x(i)*dshape_y(j)* shape_z(k);
7442 dshape(
dof_map[o],2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
7459 #ifndef MFEM_THREAD_SAFE
7469 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
7479 for (
int i = 1; i < p; i++)
7483 for (
int i = 1; i < p; i++)
7487 for (
int i = 1; i < p; i++)
7493 for (
int j = 1; j < p; j++)
7494 for (
int i = 1; i + j < p; i++)
7496 const double w = cp[i] + cp[j] + cp[p-i-j];
7501 for (
int k = 0; k <
Dof; k++)
7509 for (
int j = 0; j <= p; j++)
7510 for (
int i = 0; i + j <= p; i++)
7512 T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
7523 const int p =
Order;
7525 #ifdef MFEM_THREAD_SAFE
7526 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(
Dof);
7533 for (
int o = 0, j = 0; j <= p; j++)
7534 for (
int i = 0; i + j <= p; i++)
7536 u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
7545 const int p =
Order;
7547 #ifdef MFEM_THREAD_SAFE
7548 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
7549 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
7557 for (
int o = 0, j = 0; j <= p; j++)
7558 for (
int i = 0; i + j <= p; i++)
7561 du(o,0) = ((dshape_x(i)* shape_l(k)) -
7562 ( shape_x(i)*dshape_l(k)))*shape_y(j);
7563 du(o,1) = ((dshape_y(j)* shape_l(k)) -
7564 ( shape_y(j)*dshape_l(k)))*shape_x(i);
7568 Ti.
Mult(du, dshape);
7578 #ifndef MFEM_THREAD_SAFE
7590 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7601 for (
int i = 1; i < p; i++)
7605 for (
int i = 1; i < p; i++)
7609 for (
int i = 1; i < p; i++)
7613 for (
int i = 1; i < p; i++)
7617 for (
int i = 1; i < p; i++)
7621 for (
int i = 1; i < p; i++)
7627 for (
int j = 1; j < p; j++)
7628 for (
int i = 1; i + j < p; i++)
7630 double w = cp[i] + cp[j] + cp[p-i-j];
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];
7653 for (
int k = 1; k < p; k++)
7654 for (
int j = 1; j + k < p; j++)
7655 for (
int i = 1; i + j + k < p; i++)
7657 double w = cp[i] + cp[j] + cp[k] + cp[p-i-j-k];
7662 for (
int m = 0; m <
Dof; m++)
7671 for (
int k = 0; k <= p; k++)
7672 for (
int j = 0; j + k <= p; j++)
7673 for (
int i = 0; i + j + k <= p; i++)
7675 T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
7686 const int p =
Order;
7688 #ifdef MFEM_THREAD_SAFE
7689 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7698 for (
int o = 0, k = 0; k <= p; k++)
7699 for (
int j = 0; j + k <= p; j++)
7700 for (
int i = 0; i + j + k <= p; i++)
7702 u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
7711 const int p =
Order;
7713 #ifdef MFEM_THREAD_SAFE
7714 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7715 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
7724 for (
int o = 0, k = 0; k <= p; k++)
7725 for (
int j = 0; j + k <= p; j++)
7726 for (
int i = 0; i + j + k <= p; i++)
7728 int l = p - i - j - k;
7729 du(o,0) = ((dshape_x(i)* shape_l(l)) -
7730 ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
7731 du(o,1) = ((dshape_y(j)* shape_l(l)) -
7732 ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
7733 du(o,2) = ((dshape_z(k)* shape_l(l)) -
7734 ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
7738 Ti.
Mult(du, dshape);
7746 #ifndef MFEM_THREAD_SAFE
7756 Index(
int p) { p2p3 = 2*p + 3; }
7757 int operator()(
int i,
int j) {
return ((p2p3-j)*j)/2+i; }
7771 for (
int i = 1; i < p; i++)
7776 for (
int i = 1; i < p; i++)
7781 for (
int i = 1; i < p; i++)
7788 for (
int j = 1; j < p; j++)
7789 for (
int i = 1; i + j < p; i++)
7798 const int p,
const double l1,
const double l2,
double *shape)
7800 const double l3 = 1. - l1 - l2;
7810 for (
int o = 0, j = 0; j <= p; j++)
7814 for (
int i = 0; i <= p - j; i++)
7824 const int p,
const double l1,
const double l2,
7825 double *dshape_1d,
double *dshape)
7827 const int dof = ((p + 1)*(p + 2))/2;
7828 const double l3 = 1. - l1 - l2;
7832 for (
int o = 0, j = 0; j <= p; j++)
7836 for (
int i = 0; i <= p - j; i++)
7838 dshape[o++] = s*dshape_1d[i];
7843 for (
int i = 0; i <= p; i++)
7847 for (
int o = i, j = 0; j <= p - i; j++)
7849 dshape[dof + o] = s*dshape_1d[j];
7859 #ifdef MFEM_THREAD_SAFE
7863 for (
int i = 0; i <
Dof; i++)
7872 #ifdef MFEM_THREAD_SAFE
7877 for (
int d = 0; d < 2; d++)
7879 for (
int i = 0; i <
Dof; i++)
7891 #ifndef MFEM_THREAD_SAFE
7901 int tri(
int k) {
return (k*(k + 1))/2; }
7902 int tet(
int k) {
return (k*(k + 1)*(k + 2))/6; }
7903 Index(
int p_) { p = p_; dof = tet(p + 1); }
7904 int operator()(
int i,
int j,
int k)
7905 {
return dof - tet(p - k) - tri(p + 1 - k - j) + i; }
7921 for (
int i = 1; i < p; i++)
7926 for (
int i = 1; i < p; i++)
7931 for (
int i = 1; i < p; i++)
7936 for (
int i = 1; i < p; i++)
7941 for (
int i = 1; i < p; i++)
7946 for (
int i = 1; i < p; i++)
7953 for (
int j = 1; j < p; j++)
7954 for (
int i = 1; i + j < 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++)
7979 for (
int k = 1; k < p; k++)
7980 for (
int j = 1; j + k < p; j++)
7981 for (
int i = 1; i + j + k < p; i++)
7990 const int p,
const double l1,
const double l2,
const double l3,
7993 const double l4 = 1. - l1 - l2 - l3;
8002 for (
int o = 0, k = 0; k <= p; k++)
8005 const double ek = bp[k]*l3k;
8007 for (
int j = 0; j <= p - k; j++)
8010 double ekj = ek*bpk[j]*l2j;
8011 for (
int i = 0; i <= p - k - j; i++)
8023 const int p,
const double l1,
const double l2,
const double l3,
8024 double *dshape_1d,
double *dshape)
8026 const int dof = ((p + 1)*(p + 2)*(p + 3))/6;
8027 const double l4 = 1. - l1 - l2 - l3;
8035 for (
int o = 0, k = 0; k <= p; k++)
8038 const double ek = bp[k]*l3k;
8040 for (
int j = 0; j <= p - k; j++)
8043 double ekj = ek*bpk[j]*l2j;
8044 for (
int i = 0; i <= p - k - j; i++)
8046 dshape[o++] = dshape_1d[i]*ekj;
8057 for (
int ok = 0, k = 0; k <= p; k++)
8060 const double ek = bp[k]*l3k;
8062 for (
int i = 0; i <= p - k; i++)
8065 double eki = ek*bpk[i]*l1i;
8067 for (
int j = 0; j <= p - k - i; j++)
8069 dshape[dof + o] = dshape_1d[j]*eki;
8075 ok += ((p - k + 2)*(p - k + 1))/2;
8082 for (
int j = 0; j <= p; j++)
8085 const double ej = bp[j]*l2j;
8087 for (
int i = 0; i <= p - j; i++)
8090 double eji = ej*bpj[i]*l1i;
8091 int m = ((p + 2)*(p + 1))/2;
8092 int n = ((p - j + 2)*(p - j + 1))/2;
8093 for (
int o = i, k = 0; k <= p - j - i; k++)
8098 dshape[2*dof + o - n] = dshape_1d[k]*eji;
8111 #ifdef MFEM_THREAD_SAFE
8115 for (
int i = 0; i <
Dof; i++)
8124 #ifdef MFEM_THREAD_SAFE
8129 for (
int d = 0; d < 3; d++)
8131 for (
int i = 0; i <
Dof; i++)
8144 #ifndef MFEM_THREAD_SAFE
8149 for (
int i = 0; i <= p; i++)
8164 #ifdef MFEM_THREAD_SAFE
8174 const int p =
Order;
8180 for (
int i = 0; i <= p; i++)
8187 for (
int i = 0; i <= p; i++)
8199 #ifndef MFEM_THREAD_SAFE
8210 for (
int i = 0; i <= p; i++)
8226 #ifdef MFEM_THREAD_SAFE
8237 dofs[vertex*
Order] = 1.0;
8246 #ifndef MFEM_THREAD_SAFE
8253 for (
int o = 0, j = 0; j <= p; j++)
8254 for (
int i = 0; i <= p; i++)
8263 const int p =
Order;
8265 #ifdef MFEM_THREAD_SAFE
8266 Vector shape_x(p+1), shape_y(p+1);
8272 for (
int o = 0, j = 0; j <= p; j++)
8273 for (
int i = 0; i <= p; i++)
8275 shape(o++) = shape_x(i)*shape_y(j);
8282 const int p =
Order;
8284 #ifdef MFEM_THREAD_SAFE
8285 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
8291 for (
int o = 0, j = 0; j <= p; j++)
8292 for (
int i = 0; i <= p; i++)
8294 dshape(o,0) = dshape_x(i)* shape_y(j);
8295 dshape(o,1) = shape_x(i)*dshape_y(j); o++;
8301 const int p =
Order;
8304 #ifdef MFEM_THREAD_SAFE
8305 Vector shape_x(p+1), shape_y(p+1);
8308 for (
int i = 0; i <= p; i++)
8317 for (
int o = 0, j = 0; j <= p; j++)
8318 for (
int i = 0; i <= p; i++)
8320 dofs[o++] = shape_x(i)*shape_x(j);
8324 for (
int o = 0, j = 0; j <= p; j++)
8325 for (
int i = 0; i <= p; i++)
8327 dofs[o++] = shape_y(i)*shape_x(j);
8331 for (
int o = 0, j = 0; j <= p; j++)
8332 for (
int i = 0; i <= p; i++)
8334 dofs[o++] = shape_y(i)*shape_y(j);
8338 for (
int o = 0, j = 0; j <= p; j++)
8339 for (
int i = 0; i <= p; i++)
8341 dofs[o++] = shape_x(i)*shape_y(j);
8351 #ifndef MFEM_THREAD_SAFE
8364 for (
int o = 0, j = 0; j <= p; j++)
8365 for (
int i = 0; i <= p; i++)
8375 const int p =
Order;
8377 #ifdef MFEM_THREAD_SAFE
8378 Vector shape_x(p+1), shape_y(p+1);
8384 for (
int o = 0, j = 0; j <= p; j++)
8385 for (
int i = 0; i <= p; i++)
8387 shape(o++) = shape_x(i)*shape_y(j);
8394 const int p =
Order;
8396 #ifdef MFEM_THREAD_SAFE
8397 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
8403 for (
int o = 0, j = 0; j <= p; j++)
8404 for (
int i = 0; i <= p; i++)
8406 dshape(o,0) = dshape_x(i)* shape_y(j);
8407 dshape(o,1) = shape_x(i)*dshape_y(j); o++;
8413 const int p =
Order;
8418 case 0: dofs[0] = 1.0;
break;
8419 case 1: dofs[p] = 1.0;
break;
8420 case 2: dofs[p*(p + 2)] = 1.0;
break;
8421 case 3: dofs[p*(p + 1)] = 1.0;
break;
8431 #ifndef MFEM_THREAD_SAFE
8440 for (
int o = 0, k = 0; k <= p; k++)
8441 for (
int j = 0; j <= p; j++)
8442 for (
int i = 0; i <= p; i++)
8451 const int p =
Order;
8453 #ifdef MFEM_THREAD_SAFE
8454 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
8461 for (
int o = 0, k = 0; k <= p; k++)
8462 for (
int j = 0; j <= p; j++)
8463 for (
int i = 0; i <= p; i++)
8465 shape(o++) = shape_x(i)*shape_y(j)*shape_z(k);
8472 const int p =
Order;
8474 #ifdef MFEM_THREAD_SAFE
8475 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
8476 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
8483 for (
int o = 0, k = 0; k <= p; k++)
8484 for (
int j = 0; j <= p; j++)
8485 for (
int i = 0; i <= p; i++)
8487 dshape(o,0) = dshape_x(i)* shape_y(j)* shape_z(k);
8488 dshape(o,1) = shape_x(i)*dshape_y(j)* shape_z(k);
8489 dshape(o,2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
8495 const int p =
Order;
8498 #ifdef MFEM_THREAD_SAFE
8499 Vector shape_x(p+1), shape_y(p+1);
8502 for (
int i = 0; i <= p; i++)
8511 for (
int o = 0, k = 0; k <= p; k++)
8512 for (
int j = 0; j <= p; j++)
8513 for (
int i = 0; i <= p; i++)
8515 dofs[o++] = shape_x(i)*shape_x(j)*shape_x(k);
8519 for (
int o = 0, k = 0; k <= p; k++)
8520 for (
int j = 0; j <= p; j++)
8521 for (
int i = 0; i <= p; i++)
8523 dofs[o++] = shape_y(i)*shape_x(j)*shape_x(k);
8527 for (
int o = 0, k = 0; k <= p; k++)
8528 for (
int j = 0; j <= p; j++)
8529 for (
int i = 0; i <= p; i++)
8531 dofs[o++] = shape_y(i)*shape_y(j)*shape_x(k);
8535 for (
int o = 0, k = 0; k <= p; k++)
8536 for (
int j = 0; j <= p; j++)
8537 for (
int i = 0; i <= p; i++)
8539 dofs[o++] = shape_x(i)*shape_y(j)*shape_x(k);
8543 for (
int o = 0, k = 0; k <= p; k++)
8544 for (
int j = 0; j <= p; j++)
8545 for (
int i = 0; i <= p; i++)
8547 dofs[o++] = shape_x(i)*shape_x(j)*shape_y(k);
8551 for (
int o = 0, k = 0; k <= p; k++)
8552 for (
int j = 0; j <= p; j++)
8553 for (
int i = 0; i <= p; i++)
8555 dofs[o++] = shape_y(i)*shape_x(j)*shape_y(k);
8559 for (
int o = 0, k = 0; k <= p; k++)
8560 for (
int j = 0; j <= p; j++)
8561 for (
int i = 0; i <= p; i++)
8563 dofs[o++] = shape_y(i)*shape_y(j)*shape_y(k);
8567 for (
int o = 0, k = 0; k <= p; k++)
8568 for (
int j = 0; j <= p; j++)
8569 for (
int i = 0; i <= p; i++)
8571 dofs[o++] = shape_x(i)*shape_y(j)*shape_y(k);
8581 #ifndef MFEM_THREAD_SAFE
8596 for (
int o = 0, k = 0; k <= p; k++)
8597 for (
int j = 0; j <= p; j++)
8598 for (
int i = 0; i <= p; i++)
8608 const int p =
Order;
8610 #ifdef MFEM_THREAD_SAFE
8611 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
8618 for (
int o = 0, k = 0; k <= p; k++)
8619 for (
int j = 0; j <= p; j++)
8620 for (
int i = 0; i <= p; i++)
8622 shape(o++) = shape_x(i)*shape_y(j)*shape_z(k);
8629 const int p =
Order;
8631 #ifdef MFEM_THREAD_SAFE
8632 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
8633 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
8640 for (
int o = 0, k = 0; k <= p; k++)
8641 for (
int j = 0; j <= p; j++)
8642 for (
int i = 0; i <= p; i++)
8644 dshape(o,0) = dshape_x(i)* shape_y(j)* shape_z(k);
8645 dshape(o,1) = shape_x(i)*dshape_y(j)* shape_z(k);
8646 dshape(o,2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
8652 const int p =
Order;
8657 case 0: dofs[0] = 1.0;
break;
8658 case 1: dofs[p] = 1.0;
break;
8659 case 2: dofs[p*(p + 2)] = 1.0;
break;
8660 case 3: dofs[p*(p + 1)] = 1.0;
break;
8661 case 4: dofs[p*(p + 1)*(p + 1)] = 1.0;
break;
8662 case 5: dofs[p + p*(p + 1)*(p + 1)] = 1.0;
break;
8663 case 6: dofs[
Dof - 1] = 1.0;
break;
8664 case 7: dofs[
Dof - p - 1] = 1.0;
break;
8675 #ifndef MFEM_THREAD_SAFE
8685 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
8688 for (
int o = 0, j = 0; j <= p; j++)
8689 for (
int i = 0; i + j <= p; i++)
8691 double w = op[i] + op[j] + op[p-i-j];
8696 for (
int k = 0; k <
Dof; k++)
8703 for (
int o = 0, j = 0; j <= p; j++)
8704 for (
int i = 0; i + j <= p; i++)
8706 T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
8717 const int p =
Order;
8719 #ifdef MFEM_THREAD_SAFE
8720 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(
Dof);
8727 for (
int o = 0, j = 0; j <= p; j++)
8728 for (
int i = 0; i + j <= p; i++)
8730 u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
8739 const int p =
Order;
8741 #ifdef MFEM_THREAD_SAFE
8742 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
8743 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
8751 for (
int o = 0, j = 0; j <= p; j++)
8752 for (
int i = 0; i + j <= p; i++)
8755 du(o,0) = ((dshape_x(i)* shape_l(k)) -
8756 ( shape_x(i)*dshape_l(k)))*shape_y(j);
8757 du(o,1) = ((dshape_y(j)* shape_l(k)) -
8758 ( shape_y(j)*dshape_l(k)))*shape_x(i);
8762 Ti.
Mult(du, dshape);
8770 for (
int i = 0; i <
Dof; i++)
8773 dofs[i] = pow(1.0 - ip.
x - ip.
y,
Order);
8777 for (
int i = 0; i <
Dof; i++)
8780 dofs[i] = pow(ip.
x,
Order);
8784 for (
int i = 0; i <
Dof; i++)
8787 dofs[i] = pow(ip.
y,
Order);
8798 #ifndef MFEM_THREAD_SAFE
8808 for (
int o = 0, j = 0; j <= p; j++)
8809 for (
int i = 0; i + j <= p; i++)
8825 #ifdef MFEM_THREAD_SAFE
8838 case 0: dofs[0] = 1.0;
break;
8839 case 1: dofs[
Order] = 1.0;
break;
8840 case 2: dofs[
Dof-1] = 1.0;
break;
8851 #ifndef MFEM_THREAD_SAFE
8863 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8866 for (
int o = 0, k = 0; k <= p; k++)
8867 for (
int j = 0; j + k <= p; j++)
8868 for (
int i = 0; i + j + k <= p; i++)
8870 double w = op[i] + op[j] + op[k] + op[p-i-j-k];
8875 for (
int m = 0; m <
Dof; m++)
8883 for (
int o = 0, k = 0; k <= p; k++)
8884 for (
int j = 0; j + k <= p; j++)
8885 for (
int i = 0; i + j + k <= p; i++)
8887 T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
8898 const int p =
Order;
8900 #ifdef MFEM_THREAD_SAFE
8901 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8910 for (
int o = 0, k = 0; k <= p; k++)
8911 for (
int j = 0; j + k <= p; j++)
8912 for (
int i = 0; i + j + k <= p; i++)
8914 u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
8923 const int p =
Order;
8925 #ifdef MFEM_THREAD_SAFE
8926 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8927 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
8936 for (
int o = 0, k = 0; k <= p; k++)
8937 for (
int j = 0; j + k <= p; j++)
8938 for (
int i = 0; i + j + k <= p; i++)
8940 int l = p - i - j - k;
8941 du(o,0) = ((dshape_x(i)* shape_l(l)) -
8942 ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
8943 du(o,1) = ((dshape_y(j)* shape_l(l)) -
8944 ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
8945 du(o,2) = ((dshape_z(k)* shape_l(l)) -
8946 ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
8950 Ti.
Mult(du, dshape);
8958 for (
int i = 0; i <
Dof; i++)
8961 dofs[i] = pow(1.0 - ip.
x - ip.
y - ip.
z,
Order);
8965 for (
int i = 0; i <
Dof; i++)
8968 dofs[i] = pow(ip.
x,
Order);
8972 for (
int i = 0; i <
Dof; i++)
8975 dofs[i] = pow(ip.
y,
Order);
8978 for (
int i = 0; i <
Dof; i++)
8981 dofs[i] = pow(ip.
z,
Order);
8992 #ifndef MFEM_THREAD_SAFE
9002 for (
int o = 0, k = 0; k <= p; k++)
9003 for (
int j = 0; j + k <= p; j++)
9004 for (
int i = 0; i + j + k <= p; i++)
9021 #ifdef MFEM_THREAD_SAFE
9034 case 0: dofs[0] = 1.0;
break;
9035 case 1: dofs[
Order] = 1.0;
break;
9036 case 2: dofs[(
Order*(
Order+3))/2] = 1.0;
break;
9037 case 3: dofs[
Dof-1] = 1.0;
break;
9042 const double RT_QuadrilateralElement::nk[8] =
9043 { 0., -1., 1., 0., 0., 1., -1., 0. };
9050 cbasis1d(
poly1d.GetBasis(p + 1, VerifyClosed(cb_type))),
9051 obasis1d(
poly1d.GetBasis(p, VerifyOpen(ob_type))),
9052 dof_map(Dof), dof2nk(Dof)
9056 const int dof2 =
Dof/2;
9058 #ifndef MFEM_THREAD_SAFE
9069 for (
int i = 0; i <= p; i++)
9071 dof_map[1*dof2 + i + 0*(p + 1)] = o++;
9073 for (
int i = 0; i <= p; i++)
9075 dof_map[0*dof2 + (p + 1) + i*(p + 2)] = o++;
9077 for (
int i = 0; i <= p; i++)
9079 dof_map[1*dof2 + (p - i) + (p + 1)*(p + 1)] = o++;
9081 for (
int i = 0; i <= p; i++)
9083 dof_map[0*dof2 + 0 + (p - i)*(p + 2)] = o++;
9087 for (
int j = 0; j <= p; j++)
9088 for (
int i = 1; i <= p; i++)
9090 dof_map[0*dof2 + i + j*(p + 2)] = o++;
9092 for (
int j = 1; j <= p; j++)
9093 for (
int i = 0; i <= p; i++)
9095 dof_map[1*dof2 + i + j*(p + 1)] = o++;
9100 for (
int j = 0; j <= p; j++)
9101 for (
int i = 0; i <= p/2; i++)
9103 int idx = 0*dof2 + i + j*(p + 2);
9104 dof_map[idx] = -1 - dof_map[idx];
9107 for (
int j = p/2 + 1; j <= p; j++)
9109 int idx = 0*dof2 + (p/2 + 1) + j*(p + 2);
9110 dof_map[idx] = -1 - dof_map[idx];
9113 for (
int j = 0; j <= p/2; j++)
9114 for (
int i = 0; i <= p; i++)
9116 int idx = 1*dof2 + i + j*(p + 1);
9117 dof_map[idx] = -1 - dof_map[idx];
9120 for (
int i = 0; i <= p/2; i++)
9122 int idx = 1*dof2 + i + (p/2 + 1)*(p + 1);
9123 dof_map[idx] = -1 - dof_map[idx];
9127 for (
int j = 0; j <= p; j++)
9128 for (
int i = 0; i <= p + 1; i++)
9131 if ((idx = dof_map[o++]) < 0)
9142 for (
int j = 0; j <= p + 1; j++)
9143 for (
int i = 0; i <= p; i++)
9146 if ((idx = dof_map[o++]) < 0)
9162 const int pp1 =
Order;
9164 #ifdef MFEM_THREAD_SAFE
9165 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9168 cbasis1d.
Eval(ip.
x, shape_cx);
9169 obasis1d.
Eval(ip.
x, shape_ox);
9170 cbasis1d.
Eval(ip.
y, shape_cy);
9171 obasis1d.
Eval(ip.
y, shape_oy);
9174 for (
int j = 0; j < pp1; j++)
9175 for (
int i = 0; i <= pp1; i++)
9178 if ((idx = dof_map[o++]) < 0)
9180 idx = -1 - idx, s = -1;
9186 shape(idx,0) = s*shape_cx(i)*shape_oy(j);
9189 for (
int j = 0; j <= pp1; j++)
9190 for (
int i = 0; i < pp1; i++)
9193 if ((idx = dof_map[o++]) < 0)
9195 idx = -1 - idx, s = -1;
9202 shape(idx,1) = s*shape_ox(i)*shape_cy(j);
9209 const int pp1 =
Order;
9211 #ifdef MFEM_THREAD_SAFE
9212 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9213 Vector dshape_cx(pp1 + 1), dshape_cy(pp1 + 1);
9216 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
9217 obasis1d.
Eval(ip.
x, shape_ox);
9218 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
9219 obasis1d.
Eval(ip.
y, shape_oy);
9222 for (
int j = 0; j < pp1; j++)
9223 for (
int i = 0; i <= pp1; i++)
9226 if ((idx = dof_map[o++]) < 0)
9228 idx = -1 - idx, s = -1;
9234 divshape(idx) = s*dshape_cx(i)*shape_oy(j);
9236 for (
int j = 0; j <= pp1; j++)
9237 for (
int i = 0; i < pp1; i++)
9240 if ((idx = dof_map[o++]) < 0)
9242 idx = -1 - idx, s = -1;
9248 divshape(idx) = s*shape_ox(i)*dshape_cy(j);
9253 const double RT_HexahedronElement::nk[18] =
9254 { 0.,0.,-1., 0.,-1.,0., 1.,0.,0., 0.,1.,0., -1.,0.,0., 0.,0.,1. };
9261 cbasis1d(
poly1d.GetBasis(p + 1, VerifyClosed(cb_type))),
9262 obasis1d(
poly1d.GetBasis(p, VerifyOpen(ob_type))),
9263 dof_map(Dof), dof2nk(Dof)
9267 const int dof3 =
Dof/3;
9269 #ifndef MFEM_THREAD_SAFE
9283 for (
int j = 0; j <= p; j++)
9284 for (
int i = 0; i <= p; i++)
9286 dof_map[2*dof3 + i + ((p - j) + 0*(p + 1))*(p + 1)] = o++;
9288 for (
int j = 0; j <= p; j++)
9289 for (
int i = 0; i <= p; i++)
9291 dof_map[1*dof3 + i + (0 + j*(p + 2))*(p + 1)] = o++;
9293 for (
int j = 0; j <= p; j++)
9294 for (
int i = 0; i <= p; i++)
9296 dof_map[0*dof3 + (p + 1) + (i + j*(p + 1))*(p + 2)] = o++;
9298 for (
int j = 0; j <= p; j++)
9299 for (
int i = 0; i <= p; i++)
9301 dof_map[1*dof3 + (p - i) + ((p + 1) + j*(p + 2))*(p + 1)] = o++;
9303 for (
int j = 0; j <= p; j++)
9304 for (
int i = 0; i <= p; i++)
9306 dof_map[0*dof3 + 0 + ((p - i) + j*(p + 1))*(p + 2)] = o++;
9308 for (
int j = 0; j <= p; j++)
9309 for (
int i = 0; i <= p; i++)
9311 dof_map[2*dof3 + i + (j + (p + 1)*(p + 1))*(p + 1)] = o++;
9316 for (
int k = 0; k <= p; k++)
9317 for (
int j = 0; j <= p; j++)
9318 for (
int i = 1; i <= p; i++)
9320 dof_map[0*dof3 + i + (j + k*(p + 1))*(p + 2)] = o++;
9323 for (
int k = 0; k <= p; k++)
9324 for (
int j = 1; j <= p; j++)
9325 for (
int i = 0; i <= p; i++)
9327 dof_map[1*dof3 + i + (j + k*(p + 2))*(p + 1)] = o++;
9330 for (
int k = 1; k <= p; k++)
9331 for (
int j = 0; j <= p; j++)
9332 for (
int i = 0; i <= p; i++)
9334 dof_map[2*dof3 + i + (j + k*(p + 1))*(p + 1)] = o++;
9342 for (
int k = 0; k <= p; k++)
9343 for (
int j = 0; j <= p; j++)
9344 for (
int i = 0; i <= p/2; i++)
9346 int idx = 0*dof3 + i + (j + k*(p + 1))*(p + 2);
9347 dof_map[idx] = -1 - dof_map[idx];
9350 for (
int k = 0; k <= p; k++)
9351 for (
int j = 0; j <= p/2; j++)
9352 for (
int i = 0; i <= p; i++)
9354 int idx = 1*dof3 + i + (j + k*(p + 2))*(p + 1);
9355 dof_map[idx] = -1 - dof_map[idx];
9358 for (
int k = 0; k <= p/2; k++)
9359 for (
int j = 0; j <= p; j++)
9360 for (
int i = 0; i <= p; i++)
9362 int idx = 2*dof3 + i + (j + k*(p + 1))*(p + 1);
9363 dof_map[idx] = -1 - dof_map[idx];
9368 for (
int k = 0; k <= p; k++)
9369 for (
int j = 0; j <= p; j++)
9370 for (
int i = 0; i <= p + 1; i++)
9373 if ((idx = dof_map[o++]) < 0)
9385 for (
int k = 0; k <= p; k++)
9386 for (
int j = 0; j <= p + 1; j++)
9387 for (
int i = 0; i <= p; i++)
9390 if ((idx = dof_map[o++]) < 0)
9402 for (
int k = 0; k <= p + 1; k++)
9403 for (
int j = 0; j <= p; j++)
9404 for (
int i = 0; i <= p; i++)
9407 if ((idx = dof_map[o++]) < 0)
9423 const int pp1 =
Order;
9425 #ifdef MFEM_THREAD_SAFE
9426 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9427 Vector shape_cz(pp1 + 1), shape_oz(pp1);
9430 cbasis1d.
Eval(ip.
x, shape_cx);
9431 obasis1d.
Eval(ip.
x, shape_ox);
9432 cbasis1d.
Eval(ip.
y, shape_cy);
9433 obasis1d.
Eval(ip.
y, shape_oy);
9434 cbasis1d.
Eval(ip.
z, shape_cz);
9435 obasis1d.
Eval(ip.
z, shape_oz);
9439 for (
int k = 0; k < pp1; k++)
9440 for (
int j = 0; j < pp1; j++)
9441 for (
int i = 0; i <= pp1; i++)
9444 if ((idx = dof_map[o++]) < 0)
9446 idx = -1 - idx, s = -1;
9452 shape(idx,0) = s*shape_cx(i)*shape_oy(j)*shape_oz(k);
9457 for (
int k = 0; k < pp1; k++)
9458 for (
int j = 0; j <= pp1; j++)
9459 for (
int i = 0; i < pp1; i++)
9462 if ((idx = dof_map[o++]) < 0)
9464 idx = -1 - idx, s = -1;
9471 shape(idx,1) = s*shape_ox(i)*shape_cy(j)*shape_oz(k);
9475 for (
int k = 0; k <= pp1; k++)
9476 for (
int j = 0; j < pp1; j++)
9477 for (
int i = 0; i < pp1; i++)
9480 if ((idx = dof_map[o++]) < 0)
9482 idx = -1 - idx, s = -1;
9490 shape(idx,2) = s*shape_ox(i)*shape_oy(j)*shape_cz(k);
9497 const int pp1 =
Order;
9499 #ifdef MFEM_THREAD_SAFE
9500 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9501 Vector shape_cz(pp1 + 1), shape_oz(pp1);
9502 Vector dshape_cx(pp1 + 1), dshape_cy(pp1 + 1), dshape_cz(pp1 + 1);
9505 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
9506 obasis1d.
Eval(ip.
x, shape_ox);
9507 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
9508 obasis1d.
Eval(ip.
y, shape_oy);
9509 cbasis1d.
Eval(ip.
z, shape_cz, dshape_cz);
9510 obasis1d.
Eval(ip.
z, shape_oz);
9514 for (
int k = 0; k < pp1; k++)
9515 for (
int j = 0; j < pp1; j++)
9516 for (
int i = 0; i <= pp1; i++)
9519 if ((idx = dof_map[o++]) < 0)
9521 idx = -1 - idx, s = -1;
9527 divshape(idx) = s*dshape_cx(i)*shape_oy(j)*shape_oz(k);
9530 for (
int k = 0; k < pp1; k++)
9531 for (
int j = 0; j <= pp1; j++)
9532 for (
int i = 0; i < pp1; i++)
9535 if ((idx = dof_map[o++]) < 0)
9537 idx = -1 - idx, s = -1;
9543 divshape(idx) = s*shape_ox(i)*dshape_cy(j)*shape_oz(k);
9546 for (
int k = 0; k <= pp1; k++)
9547 for (
int j = 0; j < pp1; j++)
9548 for (
int i = 0; i < pp1; i++)
9551 if ((idx = dof_map[o++]) < 0)
9553 idx = -1 - idx, s = -1;
9559 divshape(idx) = s*shape_ox(i)*shape_oy(j)*dshape_cz(k);
9564 const double RT_TriangleElement::nk[6] =
9565 { 0., -1., 1., 1., -1., 0. };
9567 const double RT_TriangleElement::c = 1./3.;
9577 #ifndef MFEM_THREAD_SAFE
9587 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
9592 for (
int i = 0; i <= p; i++)
9597 for (
int i = 0; i <= p; i++)
9602 for (
int i = 0; i <= p; i++)
9609 for (
int j = 0; j < p; j++)
9610 for (
int i = 0; i + j < p; i++)
9612 double w = iop[i] + iop[j] + iop[p-1-i-j];
9620 for (
int k = 0; k <
Dof; k++)
9626 const double *n_k = nk + 2*dof2nk[k];
9629 for (
int j = 0; j <= p; j++)
9630 for (
int i = 0; i + j <= p; i++)
9632 double s = shape_x(i)*shape_y(j)*shape_l(p-i-j);
9633 T(o++, k) = s*n_k[0];
9634 T(o++, k) = s*n_k[1];
9636 for (
int i = 0; i <= p; i++)
9638 double s = shape_x(i)*shape_y(p-i);
9639 T(o++, k) = s*((ip.
x - c)*n_k[0] + (ip.
y - c)*n_k[1]);
9650 const int p =
Order - 1;
9652 #ifdef MFEM_THREAD_SAFE
9653 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
9662 for (
int j = 0; j <= p; j++)
9663 for (
int i = 0; i + j <= p; i++)
9665 double s = shape_x(i)*shape_y(j)*shape_l(p-i-j);
9666 u(o,0) = s; u(o,1) = 0; o++;
9667 u(o,0) = 0; u(o,1) = s; o++;
9669 for (
int i = 0; i <= p; i++)
9671 double s = shape_x(i)*shape_y(p-i);
9672 u(o,0) = (ip.
x - c)*s;
9673 u(o,1) = (ip.
y - c)*s;
9683 const int p =
Order - 1;
9685 #ifdef MFEM_THREAD_SAFE
9686 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
9687 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
9696 for (
int j = 0; j <= p; j++)
9697 for (
int i = 0; i + j <= p; i++)
9700 divu(o++) = (dshape_x(i)*shape_l(k) -
9701 shape_x(i)*dshape_l(k))*shape_y(j);
9702 divu(o++) = (dshape_y(j)*shape_l(k) -
9703 shape_y(j)*dshape_l(k))*shape_x(i);
9705 for (
int i = 0; i <= p; i++)
9708 divu(o++) = ((shape_x(i) + (ip.
x - c)*dshape_x(i))*shape_y(j) +
9709 (shape_y(j) + (ip.
y - c)*dshape_y(j))*shape_x(i));
9712 Ti.
Mult(divu, divshape);
9716 const double RT_TetrahedronElement::nk[12] =
9717 { 1,1,1, -1,0,0, 0,-1,0, 0,0,-1 };
9720 const double RT_TetrahedronElement::c = 1./4.;
9730 #ifndef MFEM_THREAD_SAFE
9742 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
9748 for (
int j = 0; j <= p; j++)
9749 for (
int i = 0; i + j <= p; i++)
9751 double w = bop[i] + bop[j] + bop[p-i-j];
9755 for (
int j = 0; j <= p; j++)
9756 for (
int i = 0; i + j <= p; i++)
9758 double w = bop[i] + bop[j] + bop[p-i-j];
9762 for (
int j = 0; j <= p; j++)
9763 for (
int i = 0; i + j <= p; i++)
9765 double w = bop[i] + bop[j] + bop[p-i-j];
9769 for (
int j = 0; j <= p; j++)
9770 for (
int i = 0; i + j <= p; i++)
9772 double w = bop[i] + bop[j] + bop[p-i-j];
9778 for (
int k = 0; k < p; k++)
9779 for (
int j = 0; j + k < p; j++)
9780 for (
int i = 0; i + j + k < p; i++)
9782 double w = iop[i] + iop[j] + iop[k] + iop[p-1-i-j-k];
9792 for (
int m = 0; m <
Dof; m++)
9799 const double *nm = nk + 3*dof2nk[m];
9802 for (
int k = 0; k <= p; k++)
9803 for (
int j = 0; j + k <= p; j++)
9804 for (
int i = 0; i + j + k <= p; i++)
9806 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
9807 T(o++, m) = s * nm[0];
9808 T(o++, m) = s * nm[1];
9809 T(o++, m) = s * nm[2];
9811 for (
int j = 0; j <= p; j++)
9812 for (
int i = 0; i + j <= p; i++)
9814 double s = shape_x(i)*shape_y(j)*shape_z(p-i-j);
9815 T(o++, m) = s*((ip.
x - c)*nm[0] + (ip.
y - c)*nm[1] +
9827 const int p =
Order - 1;
9829 #ifdef MFEM_THREAD_SAFE
9830 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
9840 for (
int k = 0; k <= p; k++)
9841 for (
int j = 0; j + k <= p; j++)
9842 for (
int i = 0; i + j + k <= p; i++)
9844 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
9845 u(o,0) = s; u(o,1) = 0; u(o,2) = 0; o++;
9846 u(o,0) = 0; u(o,1) = s; u(o,2) = 0; o++;
9847 u(o,0) = 0; u(o,1) = 0; u(o,2) = s; o++;
9849 for (
int j = 0; j <= p; j++)
9850 for (
int i = 0; i + j <= p; i++)
9852 double s = shape_x(i)*shape_y(j)*shape_z(p-i-j);
9853 u(o,0) = (ip.
x - c)*s; u(o,1) = (ip.
y - c)*s; u(o,2) = (ip.
z - c)*s;
9863 const int p =
Order - 1;
9865 #ifdef MFEM_THREAD_SAFE
9866 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
9867 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
9877 for (
int k = 0; k <= p; k++)
9878 for (
int j = 0; j + k <= p; j++)
9879 for (
int i = 0; i + j + k <= p; i++)
9881 int l = p - i - j - k;
9882 divu(o++) = (dshape_x(i)*shape_l(l) -
9883 shape_x(i)*dshape_l(l))*shape_y(j)*shape_z(k);
9884 divu(o++) = (dshape_y(j)*shape_l(l) -
9885 shape_y(j)*dshape_l(l))*shape_x(i)*shape_z(k);
9886 divu(o++) = (dshape_z(k)*shape_l(l) -
9887 shape_z(k)*dshape_l(l))*shape_x(i)*shape_y(j);
9889 for (
int j = 0; j <= p; j++)
9890 for (
int i = 0; i + j <= p; i++)
9894 (shape_x(i) + (ip.
x - c)*dshape_x(i))*shape_y(j)*shape_z(k) +
9895 (shape_y(j) + (ip.
y - c)*dshape_y(j))*shape_x(i)*shape_z(k) +
9896 (shape_z(k) + (ip.
z - c)*dshape_z(k))*shape_x(i)*shape_y(j);
9899 Ti.
Mult(divu, divshape);
9903 const double ND_HexahedronElement::tk[18] =
9904 { 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,0.,0., 0.,-1.,0., 0.,0.,-1. };
9907 const int cb_type,
const int ob_type)
9910 cbasis1d(
poly1d.GetBasis(p, VerifyClosed(cb_type))),
9911 obasis1d(
poly1d.GetBasis(p - 1, VerifyOpen(ob_type))),
9912 dof_map(Dof), dof2tk(Dof)
9916 const int dof3 =
Dof/3;
9918 #ifndef MFEM_THREAD_SAFE
9932 for (
int i = 0; i < p; i++)
9934 dof_map[0*dof3 + i + (0 + 0*(p + 1))*p] = o++;
9936 for (
int i = 0; i < p; i++)
9938 dof_map[1*dof3 + p + (i + 0*p)*(p + 1)] = o++;
9940 for (
int i = 0; i < p; i++)
9942 dof_map[0*dof3 + i + (p + 0*(p + 1))*p] = o++;
9944 for (
int i = 0; i < p; i++)
9946 dof_map[1*dof3 + 0 + (i + 0*p)*(p + 1)] = o++;
9948 for (
int i = 0; i < p; i++)
9950 dof_map[0*dof3 + i + (0 + p*(p + 1))*p] = o++;
9952 for (
int i = 0; i < p; i++)
9954 dof_map[1*dof3 + p + (i + p*p)*(p + 1)] = o++;
9956 for (
int i = 0; i < p; i++)
9958 dof_map[0*dof3 + i + (p + p*(p + 1))*p] = o++;
9960 for (
int i = 0; i < p; i++)
9962 dof_map[1*dof3 + 0 + (i + p*p)*(p + 1)] = o++;
9964 for (
int i = 0; i < p; i++)
9966 dof_map[2*dof3 + 0 + (0 + i*(p + 1))*(p + 1)] = o++;
9968 for (
int i = 0; i < p; i++)
9970 dof_map[2*dof3 + p + (0 + i*(p + 1))*(p + 1)] = o++;
9972 for (
int i = 0; i < p; i++)
9974 dof_map[2*dof3 + p + (p + i*(p + 1))*(p + 1)] = o++;
9976 for (
int i = 0; i < p; i++)
9978 dof_map[2*dof3 + 0 + (p + i*(p + 1))*(p + 1)] = o++;
9983 for (
int j = 1; j < p; j++)
9984 for (
int i = 0; i < p; i++)
9986 dof_map[0*dof3 + i + ((p - j) + 0*(p + 1))*p] = o++;
9988 for (
int j = 0; j < p; j++)
9989 for (
int i = 1; i < p; i++)
9991 dof_map[1*dof3 + i + ((p - 1 - j) + 0*p)*(p + 1)] = -1 - (o++);
9994 for (
int k = 1; k < p; k++)
9995 for (
int i = 0; i < p; i++)
9997 dof_map[0*dof3 + i + (0 + k*(p + 1))*p] = o++;
9999 for (
int k = 0; k < p; k++)
10000 for (
int i = 1; i < p; i++ )
10002 dof_map[2*dof3 + i + (0 + k*(p + 1))*(p + 1)] = o++;
10005 for (
int k = 1; k < p; k++)
10006 for (
int j = 0; j < p; j++)
10008 dof_map[1*dof3 + p + (j + k*p)*(p + 1)] = o++;
10010 for (
int k = 0; k < p; k++)
10011 for (
int j = 1; j < p; j++)
10013 dof_map[2*dof3 + p + (j + k*(p + 1))*(p + 1)] = o++;
10016 for (
int k = 1; k < p; k++)
10017 for (
int i = 0; i < p; i++)
10019 dof_map[0*dof3 + (p - 1 - i) + (p + k*(p + 1))*p] = -1 - (o++);
10021 for (
int k = 0; k < p; k++)
10022 for (
int i = 1; i < p; i++)
10024 dof_map[2*dof3 + (p - i) + (p + k*(p + 1))*(p + 1)] = o++;
10027 for (
int k = 1; k < p; k++)
10028 for (
int j = 0; j < p; j++)
10030 dof_map[1*dof3 + 0 + ((p - 1 - j) + k*p)*(p + 1)] = -1 - (o++);
10032 for (
int k = 0; k < p; k++)
10033 for (
int j = 1; j < p; j++)
10035 dof_map[2*dof3 + 0 + ((p - j) + k*(p + 1))*(p + 1)] = o++;
10038 for (
int j = 1; j < p; j++)
10039 for (
int i = 0; i < p; i++)
10041 dof_map[0*dof3 + i + (j + p*(p + 1))*p] = o++;
10043 for (
int j = 0; j < p; j++)
10044 for (
int i = 1; i < p; i++)
10046 dof_map[1*dof3 + i + (j + p*p)*(p + 1)] = o++;
10051 for (
int k = 1; k < p; k++)
10052 for (
int j = 1; j < p; j++)
10053 for (
int i = 0; i < p; i++)
10055 dof_map[0*dof3 + i + (j + k*(p + 1))*p] = o++;
10058 for (
int k = 1; k < p; k++)
10059 for (
int j = 0; j < p; j++)
10060 for (
int i = 1; i < p; i++)
10062 dof_map[1*dof3 + i + (j + k*p)*(p + 1)] = o++;
10065 for (
int k = 0; k < p; k++)
10066 for (
int j = 1; j < p; j++)
10067 for (
int i = 1; i < p; i++)
10069 dof_map[2*dof3 + i + (j + k*(p + 1))*(p + 1)] = o++;
10075 for (
int k = 0; k <= p; k++)
10076 for (
int j = 0; j <= p; j++)
10077 for (
int i = 0; i < p; i++)
10080 if ((idx = dof_map[o++]) < 0)
10082 dof2tk[idx = -1 - idx] = 3;
10091 for (
int k = 0; k <= p; k++)
10092 for (
int j = 0; j < p; j++)
10093 for (
int i = 0; i <= p; i++)
10096 if ((idx = dof_map[o++]) < 0)
10098 dof2tk[idx = -1 - idx] = 4;
10107 for (
int k = 0; k < p; k++)
10108 for (
int j = 0; j <= p; j++)
10109 for (
int i = 0; i <= p; i++)
10112 if ((idx = dof_map[o++]) < 0)
10114 dof2tk[idx = -1 - idx] = 5;
10127 const int p =
Order;
10129 #ifdef MFEM_THREAD_SAFE
10130 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10131 Vector shape_cz(p + 1), shape_oz(p);
10134 cbasis1d.
Eval(ip.
x, shape_cx);
10135 obasis1d.
Eval(ip.
x, shape_ox);
10136 cbasis1d.
Eval(ip.
y, shape_cy);
10137 obasis1d.
Eval(ip.
y, shape_oy);
10138 cbasis1d.
Eval(ip.
z, shape_cz);
10139 obasis1d.
Eval(ip.
z, shape_oz);
10143 for (
int k = 0; k <= p; k++)
10144 for (
int j = 0; j <= p; j++)
10145 for (
int i = 0; i < p; i++)
10148 if ((idx = dof_map[o++]) < 0)
10150 idx = -1 - idx, s = -1;
10156 shape(idx,0) = s*shape_ox(i)*shape_cy(j)*shape_cz(k);
10161 for (
int k = 0; k <= p; k++)
10162 for (
int j = 0; j < p; j++)
10163 for (
int i = 0; i <= p; i++)
10166 if ((idx = dof_map[o++]) < 0)
10168 idx = -1 - idx, s = -1;
10175 shape(idx,1) = s*shape_cx(i)*shape_oy(j)*shape_cz(k);
10179 for (
int k = 0; k < p; k++)
10180 for (
int j = 0; j <= p; j++)
10181 for (
int i = 0; i <= p; i++)
10184 if ((idx = dof_map[o++]) < 0)
10186 idx = -1 - idx, s = -1;
10194 shape(idx,2) = s*shape_cx(i)*shape_cy(j)*shape_oz(k);
10201 const int p =
Order;
10203 #ifdef MFEM_THREAD_SAFE
10204 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10205 Vector shape_cz(p + 1), shape_oz(p);
10206 Vector dshape_cx(p + 1), dshape_cy(p + 1), dshape_cz(p + 1);
10209 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
10210 obasis1d.
Eval(ip.
x, shape_ox);
10211 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
10212 obasis1d.
Eval(ip.
y, shape_oy);
10213 cbasis1d.
Eval(ip.
z, shape_cz, dshape_cz);
10214 obasis1d.
Eval(ip.
z, shape_oz);
10218 for (
int k = 0; k <= p; k++)
10219 for (
int j = 0; j <= p; j++)
10220 for (
int i = 0; i < p; i++)
10223 if ((idx = dof_map[o++]) < 0)
10225 idx = -1 - idx, s = -1;
10231 curl_shape(idx,0) = 0.;
10232 curl_shape(idx,1) = s*shape_ox(i)* shape_cy(j)*dshape_cz(k);
10233 curl_shape(idx,2) = -s*shape_ox(i)*dshape_cy(j)* shape_cz(k);
10236 for (
int k = 0; k <= p; k++)
10237 for (
int j = 0; j < p; j++)
10238 for (
int i = 0; i <= p; i++)
10241 if ((idx = dof_map[o++]) < 0)
10243 idx = -1 - idx, s = -1;
10249 curl_shape(idx,0) = -s* shape_cx(i)*shape_oy(j)*dshape_cz(k);
10250 curl_shape(idx,1) = 0.;
10251 curl_shape(idx,2) = s*dshape_cx(i)*shape_oy(j)* shape_cz(k);
10254 for (
int k = 0; k < p; k++)
10255 for (
int j = 0; j <= p; j++)
10256 for (
int i = 0; i <= p; i++)
10259 if ((idx = dof_map[o++]) < 0)
10261 idx = -1 - idx, s = -1;
10267 curl_shape(idx,0) = s* shape_cx(i)*dshape_cy(j)*shape_oz(k);
10268 curl_shape(idx,1) = -s*dshape_cx(i)* shape_cy(j)*shape_oz(k);
10269 curl_shape(idx,2) = 0.;
10274 const double ND_QuadrilateralElement::tk[8] =
10275 { 1.,0., 0.,1., -1.,0., 0.,-1. };
10282 cbasis1d(
poly1d.GetBasis(p, VerifyClosed(cb_type))),
10283 obasis1d(
poly1d.GetBasis(p - 1, VerifyOpen(ob_type))),
10284 dof_map(Dof), dof2tk(Dof)
10288 const int dof2 =
Dof/2;
10290 #ifndef MFEM_THREAD_SAFE
10301 for (
int i = 0; i < p; i++)
10303 dof_map[0*dof2 + i + 0*p] = o++;
10305 for (
int j = 0; j < p; j++)
10307 dof_map[1*dof2 + p + j*(p + 1)] = o++;
10309 for (
int i = 0; i < p; i++)
10311 dof_map[0*dof2 + (p - 1 - i) + p*p] = -1 - (o++);
10313 for (
int j = 0; j < p; j++)
10315 dof_map[1*dof2 + 0 + (p - 1 - j)*(p + 1)] = -1 - (o++);
10320 for (
int j = 1; j < p; j++)
10321 for (
int i = 0; i < p; i++)
10323 dof_map[0*dof2 + i + j*p] = o++;
10326 for (
int j = 0; j < p; j++)
10327 for (
int i = 1; i < p; i++)
10329 dof_map[1*dof2 + i + j*(p + 1)] = o++;
10335 for (
int j = 0; j <= p; j++)
10336 for (
int i = 0; i < p; i++)
10339 if ((idx = dof_map[o++]) < 0)
10341 dof2tk[idx = -1 - idx] = 2;
10350 for (
int j = 0; j < p; j++)
10351 for (
int i = 0; i <= p; i++)
10354 if ((idx = dof_map[o++]) < 0)
10356 dof2tk[idx = -1 - idx] = 3;
10369 const int p =
Order;
10371 #ifdef MFEM_THREAD_SAFE
10372 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10375 cbasis1d.
Eval(ip.
x, shape_cx);
10376 obasis1d.
Eval(ip.
x, shape_ox);
10377 cbasis1d.
Eval(ip.
y, shape_cy);
10378 obasis1d.
Eval(ip.
y, shape_oy);
10382 for (
int j = 0; j <= p; j++)
10383 for (
int i = 0; i < p; i++)
10386 if ((idx = dof_map[o++]) < 0)
10388 idx = -1 - idx, s = -1;
10394 shape(idx,0) = s*shape_ox(i)*shape_cy(j);
10398 for (
int j = 0; j < p; j++)
10399 for (
int i = 0; i <= p; i++)
10402 if ((idx = dof_map[o++]) < 0)
10404 idx = -1 - idx, s = -1;
10411 shape(idx,1) = s*shape_cx(i)*shape_oy(j);
10418 const int p =
Order;
10420 #ifdef MFEM_THREAD_SAFE
10421 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10422 Vector dshape_cx(p + 1), dshape_cy(p + 1);
10425 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
10426 obasis1d.
Eval(ip.
x, shape_ox);
10427 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
10428 obasis1d.
Eval(ip.
y, shape_oy);
10432 for (
int j = 0; j <= p; j++)
10433 for (
int i = 0; i < p; i++)
10436 if ((idx = dof_map[o++]) < 0)
10438 idx = -1 - idx, s = -1;
10444 curl_shape(idx,0) = -s*shape_ox(i)*dshape_cy(j);
10447 for (
int j = 0; j < p; j++)
10448 for (
int i = 0; i <= p; i++)
10451 if ((idx = dof_map[o++]) < 0)
10453 idx = -1 - idx, s = -1;
10459 curl_shape(idx,0) = s*dshape_cx(i)*shape_oy(j);
10464 const double ND_TetrahedronElement::tk[18] =
10465 { 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,1.,0., -1.,0.,1., 0.,-1.,1. };
10467 const double ND_TetrahedronElement::c = 1./4.;
10477 const int pm1 = p - 1, pm2 = p - 2, pm3 = p - 3;
10479 #ifndef MFEM_THREAD_SAFE
10490 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
10495 for (
int i = 0; i < p; i++)
10500 for (
int i = 0; i < p; i++)
10505 for (
int i = 0; i < p; i++)
10510 for (
int i = 0; i < p; i++)
10515 for (
int i = 0; i < p; i++)
10520 for (
int i = 0; i < p; i++)
10527 for (
int j = 0; j <= pm2; j++)
10528 for (
int i = 0; i + j <= pm2; i++)
10530 double w = fop[i] + fop[j] + fop[pm2-i-j];
10536 for (
int j = 0; j <= pm2; j++)
10537 for (
int i = 0; i + j <= pm2; i++)
10539 double w = fop[i] + fop[j] + fop[pm2-i-j];
10545 for (
int j = 0; j <= pm2; j++)
10546 for (
int i = 0; i + j <= pm2; i++)
10548 double w = fop[i] + fop[j] + fop[pm2-i-j];
10554 for (
int j = 0; j <= pm2; j++)
10555 for (
int i = 0; i + j <= pm2; i++)
10557 double w = fop[i] + fop[j] + fop[pm2-i-j];
10565 for (
int k = 0; k <= pm3; k++)
10566 for (
int j = 0; j + k <= pm3; j++)
10567 for (
int i = 0; i + j + k <= pm3; i++)
10569 double w = iop[i] + iop[j] + iop[k] + iop[pm3-i-j-k];
10579 for (
int m = 0; m <
Dof; m++)
10582 const double *tm = tk + 3*dof2tk[m];
10590 for (
int k = 0; k <= pm1; k++)
10591 for (
int j = 0; j + k <= pm1; j++)
10592 for (
int i = 0; i + j + k <= pm1; i++)
10594 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
10595 T(o++, m) = s * tm[0];
10596 T(o++, m) = s * tm[1];
10597 T(o++, m) = s * tm[2];
10599 for (
int k = 0; k <= pm1; k++)
10600 for (
int j = 0; j + k <= pm1; j++)
10602 double s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
10603 T(o++, m) = s*((ip.
y - c)*tm[0] - (ip.
x - c)*tm[1]);
10604 T(o++, m) = s*((ip.
z - c)*tm[0] - (ip.
x - c)*tm[2]);
10606 for (
int k = 0; k <= pm1; k++)
10609 shape_y(pm1-k)*shape_z(k)*((ip.
z - c)*tm[1] - (ip.
y - c)*tm[2]);
10620 const int pm1 =
Order - 1;
10622 #ifdef MFEM_THREAD_SAFE
10623 const int p =
Order;
10624 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
10634 for (
int k = 0; k <= pm1; k++)
10635 for (
int j = 0; j + k <= pm1; j++)
10636 for (
int i = 0; i + j + k <= pm1; i++)
10638 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
10639 u(n,0) = s; u(n,1) = 0.; u(n,2) = 0.; n++;
10640 u(n,0) = 0.; u(n,1) = s; u(n,2) = 0.; n++;
10641 u(n,0) = 0.; u(n,1) = 0.; u(n,2) = s; n++;
10643 for (
int k = 0; k <= pm1; k++)
10644 for (
int j = 0; j + k <= pm1; j++)
10646 double s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
10647 u(n,0) = s*(ip.
y - c); u(n,1) = -s*(ip.
x - c); u(n,2) = 0.; n++;
10648 u(n,0) = s*(ip.
z - c); u(n,1) = 0.; u(n,2) = -s*(ip.
x - c); n++;
10650 for (
int k = 0; k <= pm1; k++)
10652 double s = shape_y(pm1-k)*shape_z(k);
10653 u(n,0) = 0.; u(n,1) = s*(ip.
z - c); u(n,2) = -s*(ip.
y - c); n++;
10662 const int pm1 =
Order - 1;
10664 #ifdef MFEM_THREAD_SAFE
10665 const int p =
Order;
10666 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
10667 Vector dshape_x(p), dshape_y(p), dshape_z(p), dshape_l(p);
10677 for (
int k = 0; k <= pm1; k++)
10678 for (
int j = 0; j + k <= pm1; j++)
10679 for (
int i = 0; i + j + k <= pm1; i++)
10682 const double dx = (dshape_x(i)*shape_l(l) -
10683 shape_x(i)*dshape_l(l))*shape_y(j)*shape_z(k);
10684 const double dy = (dshape_y(j)*shape_l(l) -
10685 shape_y(j)*dshape_l(l))*shape_x(i)*shape_z(k);
10686 const double dz = (dshape_z(k)*shape_l(l) -
10687 shape_z(k)*dshape_l(l))*shape_x(i)*shape_y(j);
10689 u(n,0) = 0.; u(n,1) = dz; u(n,2) = -dy; n++;
10690 u(n,0) = -dz; u(n,1) = 0.; u(n,2) = dx; n++;
10691 u(n,0) = dy; u(n,1) = -dx; u(n,2) = 0.; n++;
10693 for (
int k = 0; k <= pm1; k++)
10694 for (
int j = 0; j + k <= pm1; j++)
10696 int i = pm1 - j - k;
10699 u(n,0) = shape_x(i)*(ip.
x - c)*shape_y(j)*dshape_z(k);
10700 u(n,1) = shape_x(i)*shape_y(j)*(ip.
y - c)*dshape_z(k);
10702 -((dshape_x(i)*(ip.
x - c) + shape_x(i))*shape_y(j)*shape_z(k) +
10703 (dshape_y(j)*(ip.
y - c) + shape_y(j))*shape_x(i)*shape_z(k));
10706 u(n,0) = -shape_x(i)*(ip.
x - c)*dshape_y(j)*shape_z(k);
10707 u(n,1) = (shape_x(i)*shape_y(j)*(dshape_z(k)*(ip.
z - c) + shape_z(k)) +
10708 (dshape_x(i)*(ip.
x - c) + shape_x(i))*shape_y(j)*shape_z(k));
10709 u(n,2) = -shape_x(i)*dshape_y(j)*shape_z(k)*(ip.
z - c);
10712 for (
int k = 0; k <= pm1; k++)
10716 u(n,0) = -((dshape_y(j)*(ip.
y - c) + shape_y(j))*shape_z(k) +
10717 shape_y(j)*(dshape_z(k)*(ip.
z - c) + shape_z(k)));
10722 Ti.
Mult(u, curl_shape);
10726 const double ND_TriangleElement::tk[8] =
10727 { 1.,0., -1.,1., 0.,-1., 0.,1. };
10729 const double ND_TriangleElement::c = 1./3.;
10739 const int pm1 = p - 1, pm2 = p - 2;
10741 #ifndef MFEM_THREAD_SAFE
10751 Vector shape_x(p), shape_y(p), shape_l(p);
10756 for (
int i = 0; i < p; i++)
10761 for (
int i = 0; i < p; i++)
10766 for (
int i = 0; i < p; i++)
10773 for (
int j = 0; j <= pm2; j++)
10774 for (
int i = 0; i + j <= pm2; i++)
10776 double w = iop[i] + iop[j] + iop[pm2-i-j];
10784 for (
int m = 0; m <
Dof; m++)
10787 const double *tm = tk + 2*dof2tk[m];
10794 for (
int j = 0; j <= pm1; j++)
10795 for (
int i = 0; i + j <= pm1; i++)
10797 double s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
10798 T(n++, m) = s * tm[0];
10799 T(n++, m) = s * tm[1];
10801 for (
int j = 0; j <= pm1; j++)
10804 shape_x(pm1-j)*shape_y(j)*((ip.
y - c)*tm[0] - (ip.
x - c)*tm[1]);
10815 const int pm1 =
Order - 1;
10817 #ifdef MFEM_THREAD_SAFE
10818 const int p =
Order;
10819 Vector shape_x(p), shape_y(p), shape_l(p);
10828 for (
int j = 0; j <= pm1; j++)
10829 for (
int i = 0; i + j <= pm1; i++)
10831 double s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
10832 u(n,0) = s; u(n,1) = 0; n++;
10833 u(n,0) = 0; u(n,1) = s; n++;
10835 for (
int j = 0; j <= pm1; j++)
10837 double s = shape_x(pm1-j)*shape_y(j);
10838 u(n,0) = s*(ip.
y - c);
10839 u(n,1) = -s*(ip.
x - c);
10849 const int pm1 =
Order - 1;
10851 #ifdef MFEM_THREAD_SAFE
10852 const int p =
Order;
10853 Vector shape_x(p), shape_y(p), shape_l(p);
10854 Vector dshape_x(p), dshape_y(p), dshape_l(p);
10863 for (
int j = 0; j <= pm1; j++)
10864 for (
int i = 0; i + j <= pm1; i++)
10867 const double dx = (dshape_x(i)*shape_l(l) -
10868 shape_x(i)*dshape_l(l)) * shape_y(j);
10869 const double dy = (dshape_y(j)*shape_l(l) -
10870 shape_y(j)*dshape_l(l)) * shape_x(i);
10876 for (
int j = 0; j <= pm1; j++)
10880 curlu(n++) = -((dshape_x(i)*(ip.
x - c) + shape_x(i)) * shape_y(j) +
10881 (dshape_y(j)*(ip.
y - c) + shape_y(j)) * shape_x(i));
10885 Ti.
Mult(curlu, curl2d);
10889 const double ND_SegmentElement::tk[1] = { 1. };
10894 obasis1d(
poly1d.GetBasis(p - 1, VerifyOpen(ob_type))),
10900 for (
int i = 0; i < p; i++)
10927 kv[0]->CalcShape(shape,
ijk[0], ip.
x);
10930 for (
int i = 0; i <=
Order; i++)
10932 sum += (shape(i) *=
weights(i));
10944 kv[0]->CalcDShape(grad,
ijk[0], ip.
x);
10946 double sum = 0.0, dsum = 0.0;
10947 for (
int i = 0; i <=
Order; i++)
10950 dsum += ( grad(i) *=
weights(i));
10954 add(sum, grad, -dsum*sum*sum,
shape_x, grad);
10979 for (
int o = 0, j = 0; j <=
Orders[1]; j++)
10981 const double sy =
shape_y(j);
10982 for (
int i = 0; i <=
Orders[0]; i++, o++)
10994 double sum, dsum[2];
11002 sum = dsum[0] = dsum[1] = 0.0;
11003 for (
int o = 0, j = 0; j <=
Orders[1]; j++)
11006 for (
int i = 0; i <=
Orders[0]; i++, o++)
11016 dsum[0] *= sum*sum;
11017 dsum[1] *= sum*sum;
11019 for (
int o = 0; o <
Dof; o++)
11021 dshape(o,0) = dshape(o,0)*sum -
u(o)*dsum[0];
11022 dshape(o,1) = dshape(o,1)*sum -
u(o)*dsum[1];
11054 for (
int o = 0, k = 0; k <=
Orders[2]; k++)
11056 const double sz =
shape_z(k);
11057 for (
int j = 0; j <=
Orders[1]; j++)
11059 const double sy_sz =
shape_y(j)*sz;
11060 for (
int i = 0; i <=
Orders[0]; i++, o++)
11073 double sum, dsum[3];
11083 sum = dsum[0] = dsum[1] = dsum[2] = 0.0;
11084 for (
int o = 0, k = 0; k <=
Orders[2]; k++)
11087 for (
int j = 0; j <=
Orders[1]; j++)
11089 const double sy_sz =
shape_y(j)* sz;
11090 const double dsy_sz =
dshape_y(j)* sz;
11091 const double sy_dsz =
shape_y(j)*dsz;
11092 for (
int i = 0; i <=
Orders[0]; i++, o++)
11104 dsum[0] *= sum*sum;
11105 dsum[1] *= sum*sum;
11106 dsum[2] *= sum*sum;
11108 for (
int o = 0; o <
Dof; o++)
11110 dshape(o,0) = dshape(o,0)*sum -
u(o)*dsum[0];
11111 dshape(o,1) = dshape(o,1)*sum -
u(o)*dsum[1];
11112 dshape(o,2) = dshape(o,2)*sum -
u(o)*dsum[2];
int GetNPoints() const
Returns the number of the points in the integration rule.
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
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
Multiply a matrix A with the transpose of a matrix B: A*Bt.
RT0TriangleFiniteElement()
ND_SegmentElement(const int p, const int ob_type=BasisType::GaussLegendre)
int Size() const
Logical size of the array.
For scalar fields; preserves point values.
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
RT2TriangleFiniteElement()
int GetDim() const
Returns the reference space dimension for the finite element.
GaussQuad2DFiniteElement()
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Class for an integration rule - an Array of IntegrationPoint.
CrouzeixRaviartQuadFiniteElement()
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Linear3DFiniteElement()
Construct a linear FE on tetrahedron.
static int Check(int b_type)
If the input does not represents a valid BasisType, abort with an error; otherwise return the input...
static double CalcDelta(const int p, const double x)
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
RT1TriangleFiniteElement()
No derivatives implemented.
virtual void ProjectDelta(int vertex, Vector &dofs) const
void ProjectMatrixCoefficient_RT(const double *nk, const Array< int > &d2n, MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
void ProjectMatrixCoefficient_ND(const double *tk, const Array< int > &d2t, MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
RefinedBiLinear2DFiniteElement()
Construct a biquadratic FE on quadrilateral.
NodalTensorFiniteElement(const int dims, const int p, const int btype, const DofMapType dmtype)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
void LocalInterpolation_RT(const VectorFiniteElement &cfe, const double *nk, const Array< int > &d2n, ElementTransformation &Trans, DenseMatrix &I) const
L2Pos_TriangleElement(const int p)
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
H1Pos_SegmentElement(const int p)
Basis(const int p, const double *nodes, EvalType etype=Barycentric)
Create a nodal or positive (Bernstein) basis.
virtual void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Array< const KnotVector * > kv
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
void CalcPhysDivShape(ElementTransformation &Trans, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in physical space at the po...
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
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
PositiveTensorFiniteElement(const int dims, const int p, const DofMapType dmtype)
void GivePolyPoints(const int np, double *pts, const int type)
L2Pos_TetrahedronElement(const int p)
virtual void ProjectMatrixCoefficient(MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
void CalcVShape_RT(ElementTransformation &Trans, DenseMatrix &shape) const
void SetSize(int s)
Resize the vector to size s.
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &h) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
evaluate derivatives of shape function - constant 0
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)=0
H1Pos_TetrahedronElement(const int p)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
double InnerProduct(const double *x, const double *y) const
Compute y^t A x.
LagrangeHexFiniteElement(int degree)
TensorBasisElement(const int dims, const int p, const int btype, const DofMapType dmtype)
L2_SegmentElement(const int p, const int btype=BasisType::GaussLegendre)
BiQuadPos2DFiniteElement()
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Data type dense matrix using column-major storage.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
GaussBiQuad2DFiniteElement()
int Size() const
Returns the size of the vector.
virtual void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
static int GetQuadrature1D(int b_type)
Get the corresponding Quadrature1D constant, when that makes sense; otherwise return Quadrature1D::In...
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
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
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void SetOrder() const
Update the NURBSFiniteElement according to the currently set knot vectors.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
RefinedTriLinear3DFiniteElement()
Construct a biquadratic FE on quadrilateral.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Quadratic3DFiniteElement()
Construct a quadratic FE on tetrahedron.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
void Factor()
Factor the current DenseMatrix, *a.
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void ProjectDelta(int vertex, Vector &dofs) const
void NodalLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const ScalarFiniteElement &fine_fe) const
Nodal interpolation.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
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
double * GetData() const
Return a pointer to the beginning of the Vector data.
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
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Implements CalcDivShape methods.
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &h) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Linear2DFiniteElement()
Construct a linear FE on triangle.
const double * ClosedPoints(const int p, const int btype=BasisType::GaussLobatto)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Possible basis types. Note that not all elements can use all BasisType(s).
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
void add(const Vector &v1, const Vector &v2, Vector &v)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
H1_HexahedronElement(const int p, const int btype=BasisType::GaussLobatto)
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
void CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in physical space at the point de...
L2Pos_HexahedronElement(const int p)
ND_TriangleElement(const int p)
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
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
const double * GetPoints(const int p, const int btype)
Get the coordinates of the points of the given BasisType, btype.
H1Pos_TriangleElement(const int p)
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Cubic3DFiniteElement()
Construct 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
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
virtual void 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
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
CrouzeixRaviartFiniteElement()
static void CalcShape(const int p, const double x, const double y, double *shape)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
void SetSize(int m, int n)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Implements CalcCurlShape methods.
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)
int Dof
Number of degrees of freedom.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
RefinedLinear1DFiniteElement()
Construct a quadratic FE on interval.
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
ND_HexahedronElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
For scalar fields; preserves volume integrals.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
static void ChebyshevPoints(const int p, double *x)
void AddMult_a_VWt(const double a, const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += a * v w^t.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
ND_TetrahedronElement(const int p)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Nedelec1TetFiniteElement()
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
H1_QuadrilateralElement(const int p, const int btype=BasisType::GaussLobatto)
GaussLinear2DFiniteElement()
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
int GeomType
Geometry::Type of the reference element.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
void Invert()
Replaces the current matrix with its inverse.
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
static int VerifyOpen(int b_type)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
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
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
VVt += a * v v^t.
H1_SegmentElement(const int p, const int btype=BasisType::GaussLobatto)
Nedelec1HexFiniteElement()
void Set2(const double x1, const double x2)
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
void CalcPhysDShape(ElementTransformation &Trans, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in physical space at the poi...
int GetVDim()
Returns dimension of the vector.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
L2_QuadrilateralElement(const int p, const int btype=BasisType::GaussLegendre)
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
Basis & GetBasis(const int p, const int btype)
Get a Poly_1D::Basis object of the given degree and BasisType, btype.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
void GetColumn(int c, Vector &col) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
const double * OpenPoints(const int p, const int btype=BasisType::GaussLegendre)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &h) const
void ScalarLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const ScalarFiniteElement &fine_fe) const
"Interpolation" defined through local L2-projection.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
void ProjectGrad_ND(const double *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
L2_TetrahedronElement(const int p, const int btype=BasisType::GaussLegendre)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
int GetGeomType() const
Returns the Geometry::Type of the reference element.
int Dim
Dimension of reference space.
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
void Set3(const double x1, const double x2, const double x3)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
int Orders[Geometry::MaxDim]
Anisotropic orders.
double * Data() const
Returns the matrix data array.
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
void LocalInterpolation_ND(const VectorFiniteElement &cfe, const double *tk, const Array< int > &d2t, ElementTransformation &Trans, DenseMatrix &I) const
ND_QuadrilateralElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
~LagrangeHexFiniteElement()
void CalcVShape_ND(ElementTransformation &Trans, DenseMatrix &shape) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
int GetDof() const
Returns the number of degrees of freedom in the finite element.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
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
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
void CalcPhysShape(ElementTransformation &Trans, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in physical space at the point ...
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
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
static const int * Binom(const int p)
Get a pointer to an array containing the binomial coefficients "p choose k" for k=0,...,p for the given p.
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
H1_TetrahedronElement(const int p, const int btype=BasisType::GaussLobatto)
virtual void SetOrder() const
Update the NURBSFiniteElement according to the currently set knot vectors.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
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
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
H1Pos_QuadrilateralElement(const int p)
H1_TriangleElement(const int p, const int btype=BasisType::GaussLobatto)
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 ProjectMatrixCoefficient(MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
void 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
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
static int VerifyClosed(int b_type)
virtual void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
virtual void 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.
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
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
RT_TetrahedronElement(const int p)
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
BiQuad2DFiniteElement()
Construct a biquadratic FE on quadrilateral.
virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
P0TriangleFiniteElement()
Construct P0 triangle finite element.
static int VerifyNodal(int b_type)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
RT_TriangleElement(const int p)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
void Eval(const double x, Vector &u) const
L2Pos_SegmentElement(const int p)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
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
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
void Mult(const double *x, double *y) const
Matrix vector multiplication.
static void CalcBernstein(const int p, const double x, double *u)
virtual void SetOrder() const
Update the NURBSFiniteElement according to the currently set knot vectors.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Quad2DFiniteElement()
Construct a quadratic FE on triangle.
Describes the space on each element.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void GetTransferMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const
Return interpolation matrix, I, which maps dofs from a coarse element, fe, to the fine dofs on this f...
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
GaussBiLinear2DFiniteElement()
RT_HexahedronElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void 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
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
RT_QuadrilateralElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
RotTriLinearHexFiniteElement()
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
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
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
int Order
Order/degree of the shape functions.
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
L2_HexahedronElement(const int p, const int btype=BasisType::GaussLegendre)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
static void CalcDBinomTerms(const int p, const double x, const double y, double *d)
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
TriLinear3DFiniteElement()
Construct a tri-linear FE on cube.
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
L2_TriangleElement(const int p, const int btype=BasisType::GaussLegendre)
BiLinear2DFiniteElement()
Construct a bilinear FE on quadrilateral.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Lagrange1DFiniteElement(int degree)
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
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
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
void 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