15 #include "../mesh/nurbs.hpp"
33 #ifndef MFEM_THREAD_SAFE
41 mfem_error (
"FiniteElement::CalcVShape (ip, ...)\n"
42 " is not implemented for this class!");
48 mfem_error (
"FiniteElement::CalcVShape (trans, ...)\n"
49 " is not implemented for this class!");
55 mfem_error (
"FiniteElement::CalcDivShape (ip, ...)\n"
56 " is not implemented for this class!");
63 div_shape *= (1.0 / Trans.
Weight());
69 mfem_error (
"FiniteElement::CalcCurlShape (ip, ...)\n"
70 " is not implemented for this class!");
80 #ifdef MFEM_THREAD_SAFE
85 curl_shape *= (1.0 / Trans.
Weight());
91 curl_shape *= (1.0 / Trans.
Weight());
94 MFEM_ABORT(
"Invalid dimension, Dim = " <<
Dim);
100 mfem_error (
"FiniteElement::GetFaceDofs (...)");
106 mfem_error (
"FiniteElement::CalcHessian (...) is not overloaded !");
112 mfem_error (
"GetLocalInterpolation (...) is not overloaded !");
118 mfem_error (
"FiniteElement::Project (...) is not overloaded !");
124 mfem_error (
"FiniteElement::Project (...) (vector) is not overloaded !");
130 mfem_error(
"FiniteElement::ProjectMatrixCoefficient() is not overloaded !");
135 mfem_error(
"FiniteElement::ProjectDelta(...) is not implemented for "
142 mfem_error(
"FiniteElement::Project(...) (fe version) is not implemented "
143 "for this element!");
150 mfem_error(
"FiniteElement::ProjectGrad(...) is not implemented for "
158 mfem_error(
"FiniteElement::ProjectCurl(...) is not implemented for "
166 mfem_error(
"FiniteElement::ProjectDiv(...) is not implemented for "
184 #ifdef MFEM_THREAD_SAFE
199 #ifdef MFEM_THREAD_SAFE
205 for (
int i = 0; i < fine_fe.
Dof; i++)
210 for (
int j = 0; j <
Dof; j++)
211 if (fabs (I (i,j) = c_shape (j)) < 1.0e-12)
233 for (
int i = 0; i <
Dof; i++)
236 for (
int j = 0; j < fe.
GetDof(); j++)
238 curl(i,j) = curl_shape(j,0);
246 for (
int i = 0; i <
Dof; i++)
252 dofs(i) = coeff.
Eval (Trans, ip);
255 dofs(i) *= Trans.
Weight();
266 for (
int i = 0; i <
Dof; i++)
270 vc.
Eval (x, Trans, ip);
275 for (
int j = 0; j < x.Size(); j++)
277 dofs(Dof*j+i) = x(j);
289 for (
int k = 0; k <
Dof; k++)
294 for (
int r = 0; r < MQ.Height(); r++)
296 for (
int d = 0; d < MQ.Width(); d++)
298 dofs(k+Dof*(d+MQ.Width()*r)) = MQ(r,d);
314 for (
int k = 0; k <
Dof; k++)
317 for (
int j = 0; j < shape.Size(); j++)
319 I(k,j) = (fabs(shape(j)) < 1e-12) ? 0.0 : shape(j);
328 for (
int k = 0; k <
Dof; k++)
339 I(k+d*Dof,j) =
vshape(j,d);
355 for (
int k = 0; k <
Dof; k++)
361 Mult(dshape, Jinv, grad_k);
366 for (
int j = 0; j < grad_k.Height(); j++)
367 for (
int d = 0; d <
Dim; d++)
369 grad(k+d*Dof,j) = grad_k(j,d);
382 for (
int k = 0; k <
Dof; k++)
390 for (
int j = 0; j < div_shape.Size(); j++)
392 div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j)/detJ;
397 for (
int j = 0; j < div_shape.Size(); j++)
399 div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j);
417 Vector fine_shape(fs), coarse_shape(cs);
418 DenseMatrix fine_mass(fs), fine_coarse_mass(fs, cs);
435 fine_mass_inv.
Mult(fine_coarse_mass, I);
448 for (
int i = 0; i <
Dof; i++)
452 dofs(i) = coeff.
Eval(Trans, ip);
478 pos_mass_inv.
Mult(mixed_mass, I);
483 void VectorFiniteElement::CalcShape (
486 mfem_error (
"Error: Cannot use scalar CalcShape(...) function with\n"
487 " VectorFiniteElements!");
490 void VectorFiniteElement::CalcDShape (
491 const IntegrationPoint &ip, DenseMatrix &dshape )
const
493 mfem_error (
"Error: Cannot use scalar CalcDShape(...) function with\n"
494 " VectorFiniteElements!");
526 MFEM_ABORT(
"Invalid dimension, Dim = " <<
Dim);
530 MFEM_ABORT(
"Invalid MapType = " <<
MapType);
538 #ifdef MFEM_THREAD_SAFE
543 shape *= (1.0 / Trans.
Weight());
550 #ifdef MFEM_THREAD_SAFE
563 MFEM_ASSERT(vc.
GetVDim() == sdim,
"");
565 const bool square_J = (
Dim == sdim);
567 for (
int k = 0; k <
Dof; k++)
573 if (!square_J) { dofs(k) /= Trans.
Weight(); }
584 MFEM_ASSERT(mc.
GetWidth() == sdim,
"");
585 const bool square_J = (
Dim == sdim);
587 Vector nk_phys(sdim), dofs_k(MQ.Height());
588 MFEM_ASSERT(dofs.
Size() ==
Dof*MQ.Height(),
"");
590 for (
int k = 0; k <
Dof; k++)
596 if (!square_J) { nk_phys /= T.
Weight(); }
597 MQ.Mult(nk_phys, dofs_k);
598 for (
int r = 0; r < MQ.Height(); r++)
600 dofs(k+Dof*r) = dofs_k(r);
616 for (
int k = 0; k <
Dof; k++)
625 double w = 1.0/Trans.
Weight();
626 for (
int d = 0; d <
Dim; d++)
632 for (
int j = 0; j < shape.Size(); j++)
639 for (
int d = 0; d < sdim; d++)
641 I(k,j+d*shape.Size()) = s*vk[d];
648 mfem_error(
"VectorFiniteElement::Project_RT (fe version)");
658 mfem_error(
"VectorFiniteElement::ProjectGrad_RT works only in 2D!");
666 for (
int k = 0; k <
Dof; k++)
669 tk[0] = nk[d2n[k]*
Dim+1];
670 tk[1] = -nk[d2n[k]*
Dim];
671 dshape.Mult(tk, grad_k);
672 for (
int j = 0; j < grad_k.Size(); j++)
674 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
683 #ifdef MFEM_THREAD_SAFE
696 for (
int k = 0; k <
Dof; k++)
703 J *= 1.0 / Trans.
Weight();
710 for (
int j = 0; j < curl_k.Size(); j++)
712 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
725 for (
int k = 0; k <
Dof; k++)
728 curl_shape.Mult(nk + d2n[k]*
Dim, curl_k);
729 for (
int j = 0; j < curl_k.Size(); j++)
731 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
743 for (
int k = 0; k <
Dof; k++)
760 MFEM_ASSERT(mc.
GetWidth() == sdim,
"");
762 Vector tk_phys(sdim), dofs_k(MQ.Height());
763 MFEM_ASSERT(dofs.
Size() ==
Dof*MQ.Height(),
"");
765 for (
int k = 0; k <
Dof; k++)
771 MQ.Mult(tk_phys, dofs_k);
772 for (
int r = 0; r < MQ.Height(); r++)
774 dofs(k+Dof*r) = dofs_k(r);
790 for (
int k = 0; k <
Dof; k++)
799 double w = 1.0/Trans.
Weight();
800 for (
int d = 0; d < sdim; d++)
806 for (
int j = 0; j < shape.Size(); j++)
813 for (
int d = 0; d < sdim; d++)
815 I(k, j + d*shape.Size()) = s*vk[d];
822 mfem_error(
"VectorFiniteElement::Project_ND (fe version)");
836 for (
int k = 0; k <
Dof; k++)
839 dshape.Mult(tk + d2t[k]*
Dim, grad_k);
840 for (
int j = 0; j < grad_k.Size(); j++)
842 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
854 #ifdef MFEM_THREAD_SAFE
861 for (
int k = 0; k <
Dof; k++)
869 for (
int j = 0; j <
Dof; j++)
872 for (
int i = 0; i <
Dim; i++)
874 Ikj +=
vshape(j, i) * vk[i];
876 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
888 #ifdef MFEM_THREAD_SAFE
895 for (
int k = 0; k <
Dof; k++)
903 for (
int j = 0; j <
Dof; j++)
906 for (
int i = 0; i <
Dim; i++)
908 Ikj +=
vshape(j, i) * vk[i];
910 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
944 shape(0) = 1. - ip.
x;
969 shape(0) = 1. - ip.
x - ip.
y;
977 dshape(0,0) = -1.; dshape(0,1) = -1.;
978 dshape(1,0) = 1.; dshape(1,1) = 0.;
979 dshape(2,0) = 0.; dshape(2,1) = 1.;
998 shape(0) = (1. - ip.
x) * (1. - ip.
y) ;
999 shape(1) = ip.
x * (1. - ip.
y) ;
1000 shape(2) = ip.
x * ip.
y ;
1001 shape(3) = (1. - ip.
x) * ip.
y ;
1007 dshape(0,0) = -1. + ip.
y; dshape(0,1) = -1. + ip.
x ;
1008 dshape(1,0) = 1. - ip.
y; dshape(1,1) = -ip.
x ;
1009 dshape(2,0) = ip.
y ; dshape(2,1) = ip.
x ;
1010 dshape(3,0) = -ip.
y ; dshape(3,1) = 1. - ip.
x ;
1016 h(0,0) = 0.; h(0,1) = 1.; h(0,2) = 0.;
1017 h(1,0) = 0.; h(1,1) = -1.; h(1,2) = 0.;
1018 h(2,0) = 0.; h(2,1) = 1.; h(2,2) = 0.;
1019 h(3,0) = 0.; h(3,1) = -1.; h(3,2) = 0.;
1037 const double x = ip.
x, y = ip.
y;
1039 shape(0) = 5./3. - 2. * (x + y);
1040 shape(1) = 2. * (x - 1./6.);
1041 shape(2) = 2. * (y - 1./6.);
1047 dshape(0,0) = -2.; dshape(0,1) = -2.;
1048 dshape(1,0) = 2.; dshape(1,1) = 0.;
1049 dshape(2,0) = 0.; dshape(2,1) = 2.;
1054 dofs(vertex) = 2./3.;
1055 dofs((vertex+1)%3) = 1./6.;
1056 dofs((vertex+2)%3) = 1./6.;
1061 const double GaussBiLinear2DFiniteElement::p[] =
1062 { 0.2113248654051871177454256, 0.7886751345948128822545744 };
1080 const double x = ip.
x, y = ip.
y;
1082 shape(0) = 3. * (p[1] - x) * (p[1] - y);
1083 shape(1) = 3. * (x - p[0]) * (p[1] - y);
1084 shape(2) = 3. * (x - p[0]) * (y - p[0]);
1085 shape(3) = 3. * (p[1] - x) * (y - p[0]);
1091 const double x = ip.
x, y = ip.
y;
1093 dshape(0,0) = 3. * (y - p[1]); dshape(0,1) = 3. * (x - p[1]);
1094 dshape(1,0) = 3. * (p[1] - y); dshape(1,1) = 3. * (p[0] - x);
1095 dshape(2,0) = 3. * (y - p[0]); dshape(2,1) = 3. * (x - p[0]);
1096 dshape(3,0) = 3. * (p[0] - y); dshape(3,1) = 3. * (p[1] - x);
1102 dofs(vertex) = p[1]*p[1];
1103 dofs((vertex+1)%4) = p[0]*p[1];
1104 dofs((vertex+2)%4) = p[0]*p[0];
1105 dofs((vertex+3)%4) = p[0]*p[1];
1126 shape(0) = 1. - ip.
x - ip.
y;
1134 dshape(0,0) = -1.; dshape(0,1) = -1.;
1135 dshape(1,0) = 1.; dshape(1,1) = 0.;
1136 dshape(2,0) = 0.; dshape(2,1) = 1.;
1152 double l1 = 1.0 - x, l2 = x, l3 = 2. * x - 1.;
1154 shape(0) = l1 * (-l3);
1156 shape(2) = 4. * l1 * l2;
1164 dshape(0,0) = 4. * x - 3.;
1165 dshape(1,0) = 4. * x - 1.;
1166 dshape(2,0) = 4. - 8. * x;
1181 const double x = ip.
x, x1 = 1. - x;
1185 shape(2) = 2. * x * x1;
1191 const double x = ip.
x;
1193 dshape(0,0) = 2. * x - 2.;
1194 dshape(1,0) = 2. * x;
1195 dshape(2,0) = 2. - 4. * x;
1218 double x = ip.
x, y = ip.
y;
1219 double l1 = 1.-x-y, l2 = x, l3 = y;
1221 shape(0) = l1 * (2. * l1 - 1.);
1222 shape(1) = l2 * (2. * l2 - 1.);
1223 shape(2) = l3 * (2. * l3 - 1.);
1224 shape(3) = 4. * l1 * l2;
1225 shape(4) = 4. * l2 * l3;
1226 shape(5) = 4. * l3 * l1;
1232 double x = ip.
x, y = ip.
y;
1235 dshape(0,1) = 4. * (x + y) - 3.;
1237 dshape(1,0) = 4. * x - 1.;
1241 dshape(2,1) = 4. * y - 1.;
1243 dshape(3,0) = -4. * (2. * x + y - 1.);
1244 dshape(3,1) = -4. * x;
1246 dshape(4,0) = 4. * y;
1247 dshape(4,1) = 4. * x;
1249 dshape(5,0) = -4. * y;
1250 dshape(5,1) = -4. * (x + 2. * y - 1.);
1290 case 0: dofs(3) = 0.25; dofs(5) = 0.25;
break;
1291 case 1: dofs(3) = 0.25; dofs(4) = 0.25;
break;
1292 case 2: dofs(4) = 0.25; dofs(5) = 0.25;
break;
1298 const double GaussQuad2DFiniteElement::p[] =
1299 { 0.0915762135097707434595714634022015, 0.445948490915964886318329253883051 };
1317 for (
int i = 0; i < 6; i++)
1334 const double x = ip.
x, y = ip.
y;
1348 const double x = ip.
x, y = ip.
y;
1349 D(0,0) = 0.; D(0,1) = 0.;
1350 D(1,0) = 1.; D(1,1) = 0.;
1351 D(2,0) = 0.; D(2,1) = 1.;
1352 D(3,0) = 2. * x; D(3,1) = 0.;
1353 D(4,0) = y; D(4,1) = x;
1354 D(5,0) = 0.; D(5,1) = 2. * y;
1386 double x = ip.
x, y = ip.
y;
1387 double l1x, l2x, l3x, l1y, l2y, l3y;
1389 l1x = (x - 1.) * (2. * x - 1);
1390 l2x = 4. * x * (1. - x);
1391 l3x = x * (2. * x - 1.);
1392 l1y = (y - 1.) * (2. * y - 1);
1393 l2y = 4. * y * (1. - y);
1394 l3y = y * (2. * y - 1.);
1396 shape(0) = l1x * l1y;
1397 shape(4) = l2x * l1y;
1398 shape(1) = l3x * l1y;
1399 shape(7) = l1x * l2y;
1400 shape(8) = l2x * l2y;
1401 shape(5) = l3x * l2y;
1402 shape(3) = l1x * l3y;
1403 shape(6) = l2x * l3y;
1404 shape(2) = l3x * l3y;
1410 double x = ip.
x, y = ip.
y;
1411 double l1x, l2x, l3x, l1y, l2y, l3y;
1412 double d1x, d2x, d3x, d1y, d2y, d3y;
1414 l1x = (x - 1.) * (2. * x - 1);
1415 l2x = 4. * x * (1. - x);
1416 l3x = x * (2. * x - 1.);
1417 l1y = (y - 1.) * (2. * y - 1);
1418 l2y = 4. * y * (1. - y);
1419 l3y = y * (2. * y - 1.);
1428 dshape(0,0) = d1x * l1y;
1429 dshape(0,1) = l1x * d1y;
1431 dshape(4,0) = d2x * l1y;
1432 dshape(4,1) = l2x * d1y;
1434 dshape(1,0) = d3x * l1y;
1435 dshape(1,1) = l3x * d1y;
1437 dshape(7,0) = d1x * l2y;
1438 dshape(7,1) = l1x * d2y;
1440 dshape(8,0) = d2x * l2y;
1441 dshape(8,1) = l2x * d2y;
1443 dshape(5,0) = d3x * l2y;
1444 dshape(5,1) = l3x * d2y;
1446 dshape(3,0) = d1x * l3y;
1447 dshape(3,1) = l1x * d3y;
1449 dshape(6,0) = d2x * l3y;
1450 dshape(6,1) = l2x * d3y;
1452 dshape(2,0) = d3x * l3y;
1453 dshape(2,1) = l3x * d3y;
1465 case 0: dofs(4) = 0.25; dofs(7) = 0.25;
break;
1466 case 1: dofs(4) = 0.25; dofs(5) = 0.25;
break;
1467 case 2: dofs(5) = 0.25; dofs(6) = 0.25;
break;
1468 case 3: dofs(6) = 0.25; dofs(7) = 0.25;
break;
1500 double x = ip.
x, y = ip.
y;
1501 double l1x, l2x, l3x, l1y, l2y, l3y;
1503 l1x = (1. - x) * (1. - x);
1504 l2x = 2. * x * (1. - x);
1506 l1y = (1. - y) * (1. - y);
1507 l2y = 2. * y * (1. - y);
1510 shape(0) = l1x * l1y;
1511 shape(4) = l2x * l1y;
1512 shape(1) = l3x * l1y;
1513 shape(7) = l1x * l2y;
1514 shape(8) = l2x * l2y;
1515 shape(5) = l3x * l2y;
1516 shape(3) = l1x * l3y;
1517 shape(6) = l2x * l3y;
1518 shape(2) = l3x * l3y;
1524 double x = ip.
x, y = ip.
y;
1525 double l1x, l2x, l3x, l1y, l2y, l3y;
1526 double d1x, d2x, d3x, d1y, d2y, d3y;
1528 l1x = (1. - x) * (1. - x);
1529 l2x = 2. * x * (1. - x);
1531 l1y = (1. - y) * (1. - y);
1532 l2y = 2. * y * (1. - y);
1542 dshape(0,0) = d1x * l1y;
1543 dshape(0,1) = l1x * d1y;
1545 dshape(4,0) = d2x * l1y;
1546 dshape(4,1) = l2x * d1y;
1548 dshape(1,0) = d3x * l1y;
1549 dshape(1,1) = l3x * d1y;
1551 dshape(7,0) = d1x * l2y;
1552 dshape(7,1) = l1x * d2y;
1554 dshape(8,0) = d2x * l2y;
1555 dshape(8,1) = l2x * d2y;
1557 dshape(5,0) = d3x * l2y;
1558 dshape(5,1) = l3x * d2y;
1560 dshape(3,0) = d1x * l3y;
1561 dshape(3,1) = l1x * d3y;
1563 dshape(6,0) = d2x * l3y;
1564 dshape(6,1) = l2x * d3y;
1566 dshape(2,0) = d3x * l3y;
1567 dshape(2,1) = l3x * d3y;
1575 Vector xx(&tr_ip.
x, 2), shape(s, 9);
1577 for (
int i = 0; i < 9; i++)
1581 for (
int j = 0; j < 9; j++)
1582 if (fabs(I(i,j) = s[j]) < 1.0e-12)
1587 for (
int i = 0; i < 9; i++)
1589 double *d = &I(0,i);
1590 d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
1591 d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
1592 d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
1593 d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
1594 d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
1595 0.25 * (d[0] + d[1] + d[2] + d[3]);
1604 for (
int i = 0; i < 9; i++)
1608 d[i] = coeff.
Eval(Trans, ip);
1610 d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
1611 d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
1612 d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
1613 d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
1614 d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
1615 0.25 * (d[0] + d[1] + d[2] + d[3]);
1625 for (
int i = 0; i < 9; i++)
1629 vc.
Eval (x, Trans, ip);
1630 for (
int j = 0; j < x.Size(); j++)
1635 for (
int j = 0; j < x.Size(); j++)
1637 double *d = &dofs(9*j);
1639 d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
1640 d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
1641 d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
1642 d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
1643 d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
1644 0.25 * (d[0] + d[1] + d[2] + d[3]);
1652 const double p1 = 0.5*(1.-sqrt(3./5.));
1677 const double a = sqrt(5./3.);
1678 const double p1 = 0.5*(1.-sqrt(3./5.));
1680 double x = a*(ip.
x-p1), y = a*(ip.
y-p1);
1681 double l1x, l2x, l3x, l1y, l2y, l3y;
1683 l1x = (x - 1.) * (2. * x - 1);
1684 l2x = 4. * x * (1. - x);
1685 l3x = x * (2. * x - 1.);
1686 l1y = (y - 1.) * (2. * y - 1);
1687 l2y = 4. * y * (1. - y);
1688 l3y = y * (2. * y - 1.);
1690 shape(0) = l1x * l1y;
1691 shape(4) = l2x * l1y;
1692 shape(1) = l3x * l1y;
1693 shape(7) = l1x * l2y;
1694 shape(8) = l2x * l2y;
1695 shape(5) = l3x * l2y;
1696 shape(3) = l1x * l3y;
1697 shape(6) = l2x * l3y;
1698 shape(2) = l3x * l3y;
1704 const double a = sqrt(5./3.);
1705 const double p1 = 0.5*(1.-sqrt(3./5.));
1707 double x = a*(ip.
x-p1), y = a*(ip.
y-p1);
1708 double l1x, l2x, l3x, l1y, l2y, l3y;
1709 double d1x, d2x, d3x, d1y, d2y, d3y;
1711 l1x = (x - 1.) * (2. * x - 1);
1712 l2x = 4. * x * (1. - x);
1713 l3x = x * (2. * x - 1.);
1714 l1y = (y - 1.) * (2. * y - 1);
1715 l2y = 4. * y * (1. - y);
1716 l3y = y * (2. * y - 1.);
1718 d1x = a * (4. * x - 3.);
1719 d2x = a * (4. - 8. * x);
1720 d3x = a * (4. * x - 1.);
1721 d1y = a * (4. * y - 3.);
1722 d2y = a * (4. - 8. * y);
1723 d3y = a * (4. * y - 1.);
1725 dshape(0,0) = d1x * l1y;
1726 dshape(0,1) = l1x * d1y;
1728 dshape(4,0) = d2x * l1y;
1729 dshape(4,1) = l2x * d1y;
1731 dshape(1,0) = d3x * l1y;
1732 dshape(1,1) = l3x * d1y;
1734 dshape(7,0) = d1x * l2y;
1735 dshape(7,1) = l1x * d2y;
1737 dshape(8,0) = d2x * l2y;
1738 dshape(8,1) = l2x * d2y;
1740 dshape(5,0) = d3x * l2y;
1741 dshape(5,1) = l3x * d2y;
1743 dshape(3,0) = d1x * l3y;
1744 dshape(3,1) = l1x * d3y;
1746 dshape(6,0) = d2x * l3y;
1747 dshape(6,1) = l2x * d3y;
1749 dshape(2,0) = d3x * l3y;
1750 dshape(2,1) = l3x * d3y;
1793 double x = ip.
x, y = ip.
y;
1795 double w1x, w2x, w3x, w1y, w2y, w3y;
1796 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
1798 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
1799 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
1801 l0x = (- 4.5) * w1x * w2x * w3x;
1802 l1x = ( 13.5) * x * w2x * w3x;
1803 l2x = (-13.5) * x * w1x * w3x;
1804 l3x = ( 4.5) * x * w1x * w2x;
1806 l0y = (- 4.5) * w1y * w2y * w3y;
1807 l1y = ( 13.5) * y * w2y * w3y;
1808 l2y = (-13.5) * y * w1y * w3y;
1809 l3y = ( 4.5) * y * w1y * w2y;
1811 shape(0) = l0x * l0y;
1812 shape(1) = l3x * l0y;
1813 shape(2) = l3x * l3y;
1814 shape(3) = l0x * l3y;
1815 shape(4) = l1x * l0y;
1816 shape(5) = l2x * l0y;
1817 shape(6) = l3x * l1y;
1818 shape(7) = l3x * l2y;
1819 shape(8) = l2x * l3y;
1820 shape(9) = l1x * l3y;
1821 shape(10) = l0x * l2y;
1822 shape(11) = l0x * l1y;
1823 shape(12) = l1x * l1y;
1824 shape(13) = l2x * l1y;
1825 shape(14) = l1x * l2y;
1826 shape(15) = l2x * l2y;
1832 double x = ip.
x, y = ip.
y;
1834 double w1x, w2x, w3x, w1y, w2y, w3y;
1835 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
1836 double d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
1838 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
1839 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
1841 l0x = (- 4.5) * w1x * w2x * w3x;
1842 l1x = ( 13.5) * x * w2x * w3x;
1843 l2x = (-13.5) * x * w1x * w3x;
1844 l3x = ( 4.5) * x * w1x * w2x;
1846 l0y = (- 4.5) * w1y * w2y * w3y;
1847 l1y = ( 13.5) * y * w2y * w3y;
1848 l2y = (-13.5) * y * w1y * w3y;
1849 l3y = ( 4.5) * y * w1y * w2y;
1851 d0x = -5.5 + ( 18. - 13.5 * x) * x;
1852 d1x = 9. + (-45. + 40.5 * x) * x;
1853 d2x = -4.5 + ( 36. - 40.5 * x) * x;
1854 d3x = 1. + (- 9. + 13.5 * x) * x;
1856 d0y = -5.5 + ( 18. - 13.5 * y) * y;
1857 d1y = 9. + (-45. + 40.5 * y) * y;
1858 d2y = -4.5 + ( 36. - 40.5 * y) * y;
1859 d3y = 1. + (- 9. + 13.5 * y) * y;
1861 dshape( 0,0) = d0x * l0y; dshape( 0,1) = l0x * d0y;
1862 dshape( 1,0) = d3x * l0y; dshape( 1,1) = l3x * d0y;
1863 dshape( 2,0) = d3x * l3y; dshape( 2,1) = l3x * d3y;
1864 dshape( 3,0) = d0x * l3y; dshape( 3,1) = l0x * d3y;
1865 dshape( 4,0) = d1x * l0y; dshape( 4,1) = l1x * d0y;
1866 dshape( 5,0) = d2x * l0y; dshape( 5,1) = l2x * d0y;
1867 dshape( 6,0) = d3x * l1y; dshape( 6,1) = l3x * d1y;
1868 dshape( 7,0) = d3x * l2y; dshape( 7,1) = l3x * d2y;
1869 dshape( 8,0) = d2x * l3y; dshape( 8,1) = l2x * d3y;
1870 dshape( 9,0) = d1x * l3y; dshape( 9,1) = l1x * d3y;
1871 dshape(10,0) = d0x * l2y; dshape(10,1) = l0x * d2y;
1872 dshape(11,0) = d0x * l1y; dshape(11,1) = l0x * d1y;
1873 dshape(12,0) = d1x * l1y; dshape(12,1) = l1x * d1y;
1874 dshape(13,0) = d2x * l1y; dshape(13,1) = l2x * d1y;
1875 dshape(14,0) = d1x * l2y; dshape(14,1) = l1x * d2y;
1876 dshape(15,0) = d2x * l2y; dshape(15,1) = l2x * d2y;
1882 double x = ip.
x, y = ip.
y;
1884 double w1x, w2x, w3x, w1y, w2y, w3y;
1885 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
1886 double d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
1887 double h0x, h1x, h2x, h3x, h0y, h1y, h2y, h3y;
1889 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
1890 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
1892 l0x = (- 4.5) * w1x * w2x * w3x;
1893 l1x = ( 13.5) * x * w2x * w3x;
1894 l2x = (-13.5) * x * w1x * w3x;
1895 l3x = ( 4.5) * x * w1x * w2x;
1897 l0y = (- 4.5) * w1y * w2y * w3y;
1898 l1y = ( 13.5) * y * w2y * w3y;
1899 l2y = (-13.5) * y * w1y * w3y;
1900 l3y = ( 4.5) * y * w1y * w2y;
1902 d0x = -5.5 + ( 18. - 13.5 * x) * x;
1903 d1x = 9. + (-45. + 40.5 * x) * x;
1904 d2x = -4.5 + ( 36. - 40.5 * x) * x;
1905 d3x = 1. + (- 9. + 13.5 * x) * x;
1907 d0y = -5.5 + ( 18. - 13.5 * y) * y;
1908 d1y = 9. + (-45. + 40.5 * y) * y;
1909 d2y = -4.5 + ( 36. - 40.5 * y) * y;
1910 d3y = 1. + (- 9. + 13.5 * y) * y;
1912 h0x = -27. * x + 18.;
1913 h1x = 81. * x - 45.;
1914 h2x = -81. * x + 36.;
1917 h0y = -27. * y + 18.;
1918 h1y = 81. * y - 45.;
1919 h2y = -81. * y + 36.;
1922 h( 0,0) = h0x * l0y; h( 0,1) = d0x * d0y; h( 0,2) = l0x * h0y;
1923 h( 1,0) = h3x * l0y; h( 1,1) = d3x * d0y; h( 1,2) = l3x * h0y;
1924 h( 2,0) = h3x * l3y; h( 2,1) = d3x * d3y; h( 2,2) = l3x * h3y;
1925 h( 3,0) = h0x * l3y; h( 3,1) = d0x * d3y; h( 3,2) = l0x * h3y;
1926 h( 4,0) = h1x * l0y; h( 4,1) = d1x * d0y; h( 4,2) = l1x * h0y;
1927 h( 5,0) = h2x * l0y; h( 5,1) = d2x * d0y; h( 5,2) = l2x * h0y;
1928 h( 6,0) = h3x * l1y; h( 6,1) = d3x * d1y; h( 6,2) = l3x * h1y;
1929 h( 7,0) = h3x * l2y; h( 7,1) = d3x * d2y; h( 7,2) = l3x * h2y;
1930 h( 8,0) = h2x * l3y; h( 8,1) = d2x * d3y; h( 8,2) = l2x * h3y;
1931 h( 9,0) = h1x * l3y; h( 9,1) = d1x * d3y; h( 9,2) = l1x * h3y;
1932 h(10,0) = h0x * l2y; h(10,1) = d0x * d2y; h(10,2) = l0x * h2y;
1933 h(11,0) = h0x * l1y; h(11,1) = d0x * d1y; h(11,2) = l0x * h1y;
1934 h(12,0) = h1x * l1y; h(12,1) = d1x * d1y; h(12,2) = l1x * h1y;
1935 h(13,0) = h2x * l1y; h(13,1) = d2x * d1y; h(13,2) = l2x * h1y;
1936 h(14,0) = h1x * l2y; h(14,1) = d1x * d2y; h(14,2) = l1x * h2y;
1937 h(15,0) = h2x * l2y; h(15,1) = d2x * d2y; h(15,2) = l2x * h2y;
1956 l3 = (0.33333333333333333333-x),
1957 l4 = (0.66666666666666666667-x);
1959 shape(0) = 4.5 * l2 * l3 * l4;
1960 shape(1) = 4.5 * l1 * l3 * l4;
1961 shape(2) = 13.5 * l1 * l2 * l4;
1962 shape(3) = -13.5 * l1 * l2 * l3;
1970 dshape(0,0) = -5.5 + x * (18. - 13.5 * x);
1971 dshape(1,0) = 1. - x * (9. - 13.5 * x);
1972 dshape(2,0) = 9. - x * (45. - 40.5 * x);
1973 dshape(3,0) = -4.5 + x * (36. - 40.5 * x);
2005 double x = ip.
x, y = ip.
y;
2006 double l1 = (-1. + x + y),
2010 shape(0) = -0.5*l1*(3.*l1 + 1.)*(3.*l1 + 2.);
2011 shape(1) = 0.5*x*(lx - 1.)*lx;
2012 shape(2) = 0.5*y*(-1. + ly)*ly;
2013 shape(3) = 4.5*x*l1*(3.*l1 + 1.);
2014 shape(4) = -4.5*x*lx*l1;
2015 shape(5) = 4.5*x*lx*y;
2016 shape(6) = 4.5*x*y*ly;
2017 shape(7) = -4.5*y*l1*ly;
2018 shape(8) = 4.5*y*l1*(1. + 3.*l1);
2019 shape(9) = -27.*x*y*l1;
2025 double x = ip.
x, y = ip.
y;
2027 dshape(0,0) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
2028 dshape(1,0) = 1. + 4.5*x*(-2. + 3.*x);
2030 dshape(3,0) = 4.5*(2. + 9.*x*x - 5.*y + 3.*y*y + 2.*x*(-5. + 6.*y));
2031 dshape(4,0) = -4.5*(1. - 1.*y + x*(-8. + 9.*x + 6.*y));
2032 dshape(5,0) = 4.5*(-1. + 6.*x)*y;
2033 dshape(6,0) = 4.5*y*(-1. + 3.*y);
2034 dshape(7,0) = 4.5*(1. - 3.*y)*y;
2035 dshape(8,0) = 4.5*y*(-5. + 6.*x + 6.*y);
2036 dshape(9,0) = -27.*y*(-1. + 2.*x + y);
2038 dshape(0,1) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
2040 dshape(2,1) = 1. + 4.5*y*(-2. + 3.*y);
2041 dshape(3,1) = 4.5*x*(-5. + 6.*x + 6.*y);
2042 dshape(4,1) = 4.5*(1. - 3.*x)*x;
2043 dshape(5,1) = 4.5*x*(-1. + 3.*x);
2044 dshape(6,1) = 4.5*x*(-1. + 6.*y);
2045 dshape(7,1) = -4.5*(1. + x*(-1. + 6.*y) + y*(-8. + 9.*y));
2046 dshape(8,1) = 4.5*(2. + 3.*x*x + y*(-10. + 9.*y) + x*(-5. + 12.*y));
2047 dshape(9,1) = -27.*x*(-1. + x + 2.*y);
2053 double x = ip.
x, y = ip.
y;
2055 h(0,0) = 18.-27.*(x+y);
2056 h(0,1) = 18.-27.*(x+y);
2057 h(0,2) = 18.-27.*(x+y);
2067 h(3,0) = -45.+81.*x+54.*y;
2068 h(3,1) = -22.5+54.*x+27.*y;
2071 h(4,0) = 36.-81.*x-27.*y;
2076 h(5,1) = -4.5+27.*x;
2080 h(6,1) = -4.5+27.*y;
2085 h(7,2) = 36.-27.*x-81.*y;
2088 h(8,1) = -22.5+27.*x+54.*y;
2089 h(8,2) = -45.+54.*x+81.*y;
2092 h(9,1) = 27.-54.*(x+y);
2165 double x = ip.
x, y = ip.
y, z = ip.
z;
2167 shape(0) = -((-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z)*
2168 (-1 + 3*x + 3*y + 3*z))/2.;
2169 shape(4) = (9*x*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
2170 shape(5) = (-9*x*(-1 + 3*x)*(-1 + x + y + z))/2.;
2171 shape(1) = (x*(2 + 9*(-1 + x)*x))/2.;
2172 shape(6) = (9*y*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
2173 shape(19) = -27*x*y*(-1 + x + y + z);
2174 shape(10) = (9*x*(-1 + 3*x)*y)/2.;
2175 shape(7) = (-9*y*(-1 + 3*y)*(-1 + x + y + z))/2.;
2176 shape(11) = (9*x*y*(-1 + 3*y))/2.;
2177 shape(2) = (y*(2 + 9*(-1 + y)*y))/2.;
2178 shape(8) = (9*z*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
2179 shape(18) = -27*x*z*(-1 + x + y + z);
2180 shape(12) = (9*x*(-1 + 3*x)*z)/2.;
2181 shape(17) = -27*y*z*(-1 + x + y + z);
2182 shape(16) = 27*x*y*z;
2183 shape(14) = (9*y*(-1 + 3*y)*z)/2.;
2184 shape(9) = (-9*z*(-1 + x + y + z)*(-1 + 3*z))/2.;
2185 shape(13) = (9*x*z*(-1 + 3*z))/2.;
2186 shape(15) = (9*y*z*(-1 + 3*z))/2.;
2187 shape(3) = (z*(2 + 9*(-1 + z)*z))/2.;
2193 double x = ip.
x, y = ip.
y, z = ip.
z;
2195 dshape(0,0) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
2196 x*(-4 + 6*y + 6*z)))/2.;
2197 dshape(0,1) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
2198 x*(-4 + 6*y + 6*z)))/2.;
2199 dshape(0,2) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
2200 x*(-4 + 6*y + 6*z)))/2.;
2201 dshape(4,0) = (9*(9*pow(x,2) + (-1 + y + z)*(-2 + 3*y + 3*z) +
2202 2*x*(-5 + 6*y + 6*z)))/2.;
2203 dshape(4,1) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
2204 dshape(4,2) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
2205 dshape(5,0) = (-9*(1 - y - z + x*(-8 + 9*x + 6*y + 6*z)))/2.;
2206 dshape(5,1) = (9*(1 - 3*x)*x)/2.;
2207 dshape(5,2) = (9*(1 - 3*x)*x)/2.;
2208 dshape(1,0) = 1 + (9*x*(-2 + 3*x))/2.;
2211 dshape(6,0) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
2212 dshape(6,1) = (9*(2 + 3*pow(x,2) - 10*y - 5*z + 3*(y + z)*(3*y + z) +
2213 x*(-5 + 12*y + 6*z)))/2.;
2214 dshape(6,2) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
2215 dshape(19,0) = -27*y*(-1 + 2*x + y + z);
2216 dshape(19,1) = -27*x*(-1 + x + 2*y + z);
2217 dshape(19,2) = -27*x*y;
2218 dshape(10,0) = (9*(-1 + 6*x)*y)/2.;
2219 dshape(10,1) = (9*x*(-1 + 3*x))/2.;
2221 dshape(7,0) = (9*(1 - 3*y)*y)/2.;
2222 dshape(7,1) = (-9*(1 + x*(-1 + 6*y) - z + y*(-8 + 9*y + 6*z)))/2.;
2223 dshape(7,2) = (9*(1 - 3*y)*y)/2.;
2224 dshape(11,0) = (9*y*(-1 + 3*y))/2.;
2225 dshape(11,1) = (9*x*(-1 + 6*y))/2.;
2228 dshape(2,1) = 1 + (9*y*(-2 + 3*y))/2.;
2230 dshape(8,0) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
2231 dshape(8,1) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
2232 dshape(8,2) = (9*(2 + 3*pow(x,2) - 5*y - 10*z + 3*(y + z)*(y + 3*z) +
2233 x*(-5 + 6*y + 12*z)))/2.;
2234 dshape(18,0) = -27*z*(-1 + 2*x + y + z);
2235 dshape(18,1) = -27*x*z;
2236 dshape(18,2) = -27*x*(-1 + x + y + 2*z);
2237 dshape(12,0) = (9*(-1 + 6*x)*z)/2.;
2239 dshape(12,2) = (9*x*(-1 + 3*x))/2.;
2240 dshape(17,0) = -27*y*z;
2241 dshape(17,1) = -27*z*(-1 + x + 2*y + z);
2242 dshape(17,2) = -27*y*(-1 + x + y + 2*z);
2243 dshape(16,0) = 27*y*z;
2244 dshape(16,1) = 27*x*z;
2245 dshape(16,2) = 27*x*y;
2247 dshape(14,1) = (9*(-1 + 6*y)*z)/2.;
2248 dshape(14,2) = (9*y*(-1 + 3*y))/2.;
2249 dshape(9,0) = (9*(1 - 3*z)*z)/2.;
2250 dshape(9,1) = (9*(1 - 3*z)*z)/2.;
2251 dshape(9,2) = (9*(-1 + x + y + 8*z - 6*(x + y)*z - 9*pow(z,2)))/2.;
2252 dshape(13,0) = (9*z*(-1 + 3*z))/2.;
2254 dshape(13,2) = (9*x*(-1 + 6*z))/2.;
2256 dshape(15,1) = (9*z*(-1 + 3*z))/2.;
2257 dshape(15,2) = (9*y*(-1 + 6*z))/2.;
2260 dshape(3,2) = 1 + (9*z*(-2 + 3*z))/2.;
2326 shape(0) = 1. - ip.
x - ip.
y - ip.
z;
2335 if (dshape.
Height() == 4)
2337 double *A = &dshape(0,0);
2338 A[0] = -1.; A[4] = -1.; A[8] = -1.;
2339 A[1] = 1.; A[5] = 0.; A[9] = 0.;
2340 A[2] = 0.; A[6] = 1.; A[10] = 0.;
2341 A[3] = 0.; A[7] = 0.; A[11] = 1.;
2345 dshape(0,0) = -1.; dshape(0,1) = -1.; dshape(0,2) = -1.;
2346 dshape(1,0) = 1.; dshape(1,1) = 0.; dshape(1,2) = 0.;
2347 dshape(2,0) = 0.; dshape(2,1) = 1.; dshape(2,2) = 0.;
2348 dshape(3,0) = 0.; dshape(3,1) = 0.; dshape(3,2) = 1.;
2355 static int face_dofs[4][3] = {{1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {0, 1, 2}};
2358 *dofs = face_dofs[face];
2400 double L0, L1, L2, L3;
2402 L0 = 1. - ip.
x - ip.
y - ip.
z;
2407 shape(0) = L0 * ( 2.0 * L0 - 1.0 );
2408 shape(1) = L1 * ( 2.0 * L1 - 1.0 );
2409 shape(2) = L2 * ( 2.0 * L2 - 1.0 );
2410 shape(3) = L3 * ( 2.0 * L3 - 1.0 );
2411 shape(4) = 4.0 * L0 * L1;
2412 shape(5) = 4.0 * L0 * L2;
2413 shape(6) = 4.0 * L0 * L3;
2414 shape(7) = 4.0 * L1 * L2;
2415 shape(8) = 4.0 * L1 * L3;
2416 shape(9) = 4.0 * L2 * L3;
2427 L0 = 1.0 - x - y - z;
2429 dshape(0,0) = dshape(0,1) = dshape(0,2) = 1.0 - 4.0 * L0;
2430 dshape(1,0) = -1.0 + 4.0 * x; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
2431 dshape(2,0) = 0.0; dshape(2,1) = -1.0 + 4.0 * y; dshape(2,2) = 0.0;
2432 dshape(3,0) = dshape(3,1) = 0.0; dshape(3,2) = -1.0 + 4.0 * z;
2433 dshape(4,0) = 4.0 * (L0 - x); dshape(4,1) = dshape(4,2) = -4.0 * x;
2434 dshape(5,0) = dshape(5,2) = -4.0 * y; dshape(5,1) = 4.0 * (L0 - y);
2435 dshape(6,0) = dshape(6,1) = -4.0 * z; dshape(6,2) = 4.0 * (L0 - z);
2436 dshape(7,0) = 4.0 * y; dshape(7,1) = 4.0 * x; dshape(7,2) = 0.0;
2437 dshape(8,0) = 4.0 * z; dshape(8,1) = 0.0; dshape(8,2) = 4.0 * x;
2438 dshape(9,0) = 0.0; dshape(9,1) = 4.0 * z; dshape(9,2) = 4.0 * y;
2480 double x = ip.
x, y = ip.
y, z = ip.
z;
2481 double ox = 1.-x, oy = 1.-y, oz = 1.-z;
2483 shape(0) = ox * oy * oz;
2484 shape(1) = x * oy * oz;
2485 shape(2) = x * y * oz;
2486 shape(3) = ox * y * oz;
2487 shape(4) = ox * oy * z;
2488 shape(5) = x * oy * z;
2489 shape(6) = x * y * z;
2490 shape(7) = ox * y * z;
2496 double x = ip.
x, y = ip.
y, z = ip.
z;
2497 double ox = 1.-x, oy = 1.-y, oz = 1.-z;
2499 dshape(0,0) = - oy * oz;
2500 dshape(0,1) = - ox * oz;
2501 dshape(0,2) = - ox * oy;
2503 dshape(1,0) = oy * oz;
2504 dshape(1,1) = - x * oz;
2505 dshape(1,2) = - x * oy;
2507 dshape(2,0) = y * oz;
2508 dshape(2,1) = x * oz;
2509 dshape(2,2) = - x * y;
2511 dshape(3,0) = - y * oz;
2512 dshape(3,1) = ox * oz;
2513 dshape(3,2) = - ox * y;
2515 dshape(4,0) = - oy * z;
2516 dshape(4,1) = - ox * z;
2517 dshape(4,2) = ox * oy;
2519 dshape(5,0) = oy * z;
2520 dshape(5,1) = - x * z;
2521 dshape(5,2) = x * oy;
2523 dshape(6,0) = y * z;
2524 dshape(6,1) = x * z;
2525 dshape(6,2) = x * y;
2527 dshape(7,0) = - y * z;
2528 dshape(7,1) = ox * z;
2529 dshape(7,2) = ox * y;
2564 shape(0) = 1.0 - 2.0 * ip.
y;
2565 shape(1) = -1.0 + 2.0 * ( ip.
x + ip.
y );
2566 shape(2) = 1.0 - 2.0 * ip.
x;
2572 dshape(0,0) = 0.0; dshape(0,1) = -2.0;
2573 dshape(1,0) = 2.0; dshape(1,1) = 2.0;
2574 dshape(2,0) = -2.0; dshape(2,1) = 0.0;
2595 const double l1 = ip.
x+ip.
y-0.5, l2 = 1.-l1, l3 = ip.
x-ip.
y+0.5, l4 = 1.-l3;
2606 const double x2 = 2.*ip.
x, y2 = 2.*ip.
y;
2608 dshape(0,0) = 1. - x2; dshape(0,1) = -2. + y2;
2609 dshape(1,0) = x2; dshape(1,1) = 1. - y2;
2610 dshape(2,0) = 1. - x2; dshape(2,1) = y2;
2611 dshape(3,0) = -2. + x2; dshape(3,1) = 1. - y2;
2629 double x = ip.
x, y = ip.
y;
2632 shape(0,1) = y - 1.;
2635 shape(2,0) = x - 1.;
2647 const double RT0TriangleFiniteElement::nk[3][2] =
2648 { {0, -1}, {1, 1}, {-1, 0} };
2654 #ifdef MFEM_THREAD_SAFE
2660 for (k = 0; k < 3; k++)
2663 for (j = 0; j < 3; j++)
2666 if (j == k) { d -= 1.0; }
2667 if (fabs(d) > 1.0e-12)
2669 mfem::err <<
"RT0TriangleFiniteElement::GetLocalInterpolation (...)\n"
2670 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2686 for (k = 0; k < 3; k++)
2689 ip.
x = vk[0]; ip.
y = vk[1];
2692 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2693 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2694 for (j = 0; j < 3; j++)
2695 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2708 #ifdef MFEM_THREAD_SAFE
2712 for (
int k = 0; k < 3; k++)
2720 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2721 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
2741 double x = ip.
x, y = ip.
y;
2744 shape(0,1) = y - 1.;
2749 shape(3,0) = x - 1.;
2762 const double RT0QuadFiniteElement::nk[4][2] =
2763 { {0, -1}, {1, 0}, {0, 1}, {-1, 0} };
2769 #ifdef MFEM_THREAD_SAFE
2775 for (k = 0; k < 4; k++)
2778 for (j = 0; j < 4; j++)
2781 if (j == k) { d -= 1.0; }
2782 if (fabs(d) > 1.0e-12)
2784 mfem::err <<
"RT0QuadFiniteElement::GetLocalInterpolation (...)\n"
2785 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2801 for (k = 0; k < 4; k++)
2804 ip.
x = vk[0]; ip.
y = vk[1];
2807 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2808 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2809 for (j = 0; j < 4; j++)
2810 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2823 #ifdef MFEM_THREAD_SAFE
2827 for (
int k = 0; k < 4; k++)
2835 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2836 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
2864 double x = ip.
x, y = ip.
y;
2866 shape(0,0) = -2 * x * (-1 + x + 2 * y);
2867 shape(0,1) = -2 * (-1 + y) * (-1 + x + 2 * y);
2868 shape(1,0) = 2 * x * (x - y);
2869 shape(1,1) = 2 * (x - y) * (-1 + y);
2870 shape(2,0) = 2 * x * (-1 + 2 * x + y);
2871 shape(2,1) = 2 * y * (-1 + 2 * x + y);
2872 shape(3,0) = 2 * x * (-1 + x + 2 * y);
2873 shape(3,1) = 2 * y * (-1 + x + 2 * y);
2874 shape(4,0) = -2 * (-1 + x) * (x - y);
2875 shape(4,1) = 2 * y * (-x + y);
2876 shape(5,0) = -2 * (-1 + x) * (-1 + 2 * x + y);
2877 shape(5,1) = -2 * y * (-1 + 2 * x + y);
2878 shape(6,0) = -3 * x * (-2 + 2 * x + y);
2879 shape(6,1) = -3 * y * (-1 + 2 * x + y);
2880 shape(7,0) = -3 * x * (-1 + x + 2 * y);
2881 shape(7,1) = -3 * y * (-2 + x + 2 * y);
2887 double x = ip.
x, y = ip.
y;
2889 divshape(0) = -2 * (-4 + 3 * x + 6 * y);
2890 divshape(1) = 2 + 6 * x - 6 * y;
2891 divshape(2) = -4 + 12 * x + 6 * y;
2892 divshape(3) = -4 + 6 * x + 12 * y;
2893 divshape(4) = 2 - 6 * x + 6 * y;
2894 divshape(5) = -2 * (-4 + 6 * x + 3 * y);
2895 divshape(6) = -9 * (-1 + 2 * x + y);
2896 divshape(7) = -9 * (-1 + x + 2 * y);
2899 const double RT1TriangleFiniteElement::nk[8][2] =
2911 #ifdef MFEM_THREAD_SAFE
2917 for (k = 0; k < 8; k++)
2920 for (j = 0; j < 8; j++)
2923 if (j == k) { d -= 1.0; }
2924 if (fabs(d) > 1.0e-12)
2926 mfem::err <<
"RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
2927 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2943 for (k = 0; k < 8; k++)
2946 ip.
x = vk[0]; ip.
y = vk[1];
2949 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2950 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2951 for (j = 0; j < 8; j++)
2952 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2964 #ifdef MFEM_THREAD_SAFE
2968 for (
int k = 0; k < 8; k++)
2976 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2977 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
3020 double x = ip.
x, y = ip.
y;
3024 shape(0,1) = -( 1. - 3.*y + 2.*y*y)*( 2. - 3.*x);
3026 shape(1,1) = -( 1. - 3.*y + 2.*y*y)*(-1. + 3.*x);
3028 shape(2,0) = (-x + 2.*x*x)*( 2. - 3.*y);
3030 shape(3,0) = (-x + 2.*x*x)*(-1. + 3.*y);
3034 shape(4,1) = (-y + 2.*y*y)*(-1. + 3.*x);
3036 shape(5,1) = (-y + 2.*y*y)*( 2. - 3.*x);
3038 shape(6,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y);
3040 shape(7,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y);
3043 shape(8,0) = (4.*x - 4.*x*x)*( 2. - 3.*y);
3045 shape(9,0) = (4.*x - 4.*x*x)*(-1. + 3.*y);
3049 shape(10,1) = (4.*y - 4.*y*y)*( 2. - 3.*x);
3051 shape(11,1) = (4.*y - 4.*y*y)*(-1. + 3.*x);
3057 double x = ip.
x, y = ip.
y;
3059 divshape(0) = -(-3. + 4.*y)*( 2. - 3.*x);
3060 divshape(1) = -(-3. + 4.*y)*(-1. + 3.*x);
3061 divshape(2) = (-1. + 4.*x)*( 2. - 3.*y);
3062 divshape(3) = (-1. + 4.*x)*(-1. + 3.*y);
3063 divshape(4) = (-1. + 4.*y)*(-1. + 3.*x);
3064 divshape(5) = (-1. + 4.*y)*( 2. - 3.*x);
3065 divshape(6) = -(-3. + 4.*x)*(-1. + 3.*y);
3066 divshape(7) = -(-3. + 4.*x)*( 2. - 3.*y);
3067 divshape(8) = ( 4. - 8.*x)*( 2. - 3.*y);
3068 divshape(9) = ( 4. - 8.*x)*(-1. + 3.*y);
3069 divshape(10) = ( 4. - 8.*y)*( 2. - 3.*x);
3070 divshape(11) = ( 4. - 8.*y)*(-1. + 3.*x);
3073 const double RT1QuadFiniteElement::nk[12][2] =
3093 #ifdef MFEM_THREAD_SAFE
3099 for (k = 0; k < 12; k++)
3102 for (j = 0; j < 12; j++)
3105 if (j == k) { d -= 1.0; }
3106 if (fabs(d) > 1.0e-12)
3108 mfem::err <<
"RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
3109 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
3125 for (k = 0; k < 12; k++)
3128 ip.
x = vk[0]; ip.
y = vk[1];
3131 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
3132 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
3133 for (j = 0; j < 12; j++)
3134 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
3146 #ifdef MFEM_THREAD_SAFE
3150 for (
int k = 0; k < 12; k++)
3158 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
3159 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
3163 const double RT2TriangleFiniteElement::M[15][15] =
3166 0, -5.3237900077244501311, 5.3237900077244501311, 16.647580015448900262,
3167 0, 24.442740046346700787, -16.647580015448900262, -12.,
3168 -19.118950038622250656, -47.237900077244501311, 0, -34.414110069520051180,
3169 12., 30.590320061795601049, 15.295160030897800524
3172 0, 1.5, -1.5, -15., 0, 2.625, 15., 15., -4.125, 30., 0, -14.625, -15.,
3176 0, -0.67620999227554986889, 0.67620999227554986889, 7.3524199845510997378,
3177 0, -3.4427400463467007866, -7.3524199845510997378, -12.,
3178 4.1189500386222506555, -0.76209992275549868892, 0, 7.4141100695200511800,
3179 12., -6.5903200617956010489, -3.2951600308978005244
3182 0, 0, 1.5, 0, 0, 1.5, -11.471370023173350393, 0, 2.4713700231733503933,
3183 -11.471370023173350393, 0, 2.4713700231733503933, 15.295160030897800524,
3184 0, -3.2951600308978005244
3187 0, 0, 4.875, 0, 0, 4.875, -16.875, 0, -16.875, -16.875, 0, -16.875, 10.5,
3191 0, 0, 1.5, 0, 0, 1.5, 2.4713700231733503933, 0, -11.471370023173350393,
3192 2.4713700231733503933, 0, -11.471370023173350393, -3.2951600308978005244,
3193 0, 15.295160030897800524
3196 -0.67620999227554986889, 0, -3.4427400463467007866, 0,
3197 7.3524199845510997378, 0.67620999227554986889, 7.4141100695200511800, 0,
3198 -0.76209992275549868892, 4.1189500386222506555, -12.,
3199 -7.3524199845510997378, -3.2951600308978005244, -6.5903200617956010489,
3203 1.5, 0, 2.625, 0, -15., -1.5, -14.625, 0, 30., -4.125, 15., 15., 10.5,
3207 -5.3237900077244501311, 0, 24.442740046346700787, 0, 16.647580015448900262,
3208 5.3237900077244501311, -34.414110069520051180, 0, -47.237900077244501311,
3209 -19.118950038622250656, -12., -16.647580015448900262, 15.295160030897800524,
3210 30.590320061795601049, 12.
3212 { 0, 0, 18., 0, 0, 6., -42., 0, -30., -26., 0, -14., 24., 32., 8.},
3213 { 0, 0, 6., 0, 0, 18., -14., 0, -26., -30., 0, -42., 8., 32., 24.},
3214 { 0, 0, -6., 0, 0, -4., 30., 0, 4., 22., 0, 4., -24., -16., 0},
3215 { 0, 0, -4., 0, 0, -8., 20., 0, 8., 36., 0, 8., -16., -32., 0},
3216 { 0, 0, -8., 0, 0, -4., 8., 0, 36., 8., 0, 20., 0, -32., -16.},
3217 { 0, 0, -4., 0, 0, -6., 4., 0, 22., 4., 0, 30., 0, -16., -24.}
3223 const double p = 0.11270166537925831148;
3260 double x = ip.
x, y = ip.
y;
3262 double Bx[15] = {1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y, 0., x*x*x,
3265 double By[15] = {0., 1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y,
3269 for (
int i = 0; i < 15; i++)
3271 double cx = 0.0, cy = 0.0;
3272 for (
int j = 0; j < 15; j++)
3274 cx += M[i][j] * Bx[j];
3275 cy += M[i][j] * By[j];
3285 double x = ip.
x, y = ip.
y;
3287 double DivB[15] = {0., 0., 1., 0., 0., 1., 2.*x, 0., y, x, 0., 2.*y,
3288 4.*x*x, 4.*x*y, 4.*y*y
3291 for (
int i = 0; i < 15; i++)
3294 for (
int j = 0; j < 15; j++)
3296 div += M[i][j] * DivB[j];
3302 const double RT2QuadFiniteElement::pt[4] = {0.,1./3.,2./3.,1.};
3304 const double RT2QuadFiniteElement::dpt[3] = {0.25,0.5,0.75};
3346 double x = ip.
x, y = ip.
y;
3348 double ax0 = pt[0] - x;
3349 double ax1 = pt[1] - x;
3350 double ax2 = pt[2] - x;
3351 double ax3 = pt[3] - x;
3353 double by0 = dpt[0] - y;
3354 double by1 = dpt[1] - y;
3355 double by2 = dpt[2] - y;
3357 double ay0 = pt[0] - y;
3358 double ay1 = pt[1] - y;
3359 double ay2 = pt[2] - y;
3360 double ay3 = pt[3] - y;
3362 double bx0 = dpt[0] - x;
3363 double bx1 = dpt[1] - x;
3364 double bx2 = dpt[2] - x;
3366 double A01 = pt[0] - pt[1];
3367 double A02 = pt[0] - pt[2];
3368 double A12 = pt[1] - pt[2];
3369 double A03 = pt[0] - pt[3];
3370 double A13 = pt[1] - pt[3];
3371 double A23 = pt[2] - pt[3];
3373 double B01 = dpt[0] - dpt[1];
3374 double B02 = dpt[0] - dpt[2];
3375 double B12 = dpt[1] - dpt[2];
3377 double tx0 = (bx1*bx2)/(B01*B02);
3378 double tx1 = -(bx0*bx2)/(B01*B12);
3379 double tx2 = (bx0*bx1)/(B02*B12);
3381 double ty0 = (by1*by2)/(B01*B02);
3382 double ty1 = -(by0*by2)/(B01*B12);
3383 double ty2 = (by0*by1)/(B02*B12);
3387 shape(0, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx0;
3389 shape(1, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx1;
3391 shape(2, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx2;
3393 shape(3, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty0;
3395 shape(4, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty1;
3397 shape(5, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty2;
3401 shape(6, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx2;
3403 shape(7, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx1;
3405 shape(8, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx0;
3407 shape(9, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty2;
3409 shape(10, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty1;
3411 shape(11, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty0;
3414 shape(12, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty0;
3416 shape(13, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty1;
3418 shape(14, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty2;
3421 shape(15, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty0;
3423 shape(16, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty1;
3425 shape(17, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty2;
3429 shape(18, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx0;
3431 shape(19, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx1;
3433 shape(20, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx2;
3436 shape(21, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx0;
3438 shape(22, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx1;
3440 shape(23, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx2;
3446 double x = ip.
x, y = ip.
y;
3448 double a01 = pt[0]*pt[1];
3449 double a02 = pt[0]*pt[2];
3450 double a12 = pt[1]*pt[2];
3451 double a03 = pt[0]*pt[3];
3452 double a13 = pt[1]*pt[3];
3453 double a23 = pt[2]*pt[3];
3455 double bx0 = dpt[0] - x;
3456 double bx1 = dpt[1] - x;
3457 double bx2 = dpt[2] - x;
3459 double by0 = dpt[0] - y;
3460 double by1 = dpt[1] - y;
3461 double by2 = dpt[2] - y;
3463 double A01 = pt[0] - pt[1];
3464 double A02 = pt[0] - pt[2];
3465 double A12 = pt[1] - pt[2];
3466 double A03 = pt[0] - pt[3];
3467 double A13 = pt[1] - pt[3];
3468 double A23 = pt[2] - pt[3];
3470 double A012 = pt[0] + pt[1] + pt[2];
3471 double A013 = pt[0] + pt[1] + pt[3];
3472 double A023 = pt[0] + pt[2] + pt[3];
3473 double A123 = pt[1] + pt[2] + pt[3];
3475 double B01 = dpt[0] - dpt[1];
3476 double B02 = dpt[0] - dpt[2];
3477 double B12 = dpt[1] - dpt[2];
3479 double tx0 = (bx1*bx2)/(B01*B02);
3480 double tx1 = -(bx0*bx2)/(B01*B12);
3481 double tx2 = (bx0*bx1)/(B02*B12);
3483 double ty0 = (by1*by2)/(B01*B02);
3484 double ty1 = -(by0*by2)/(B01*B12);
3485 double ty2 = (by0*by1)/(B02*B12);
3488 divshape(0) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx0;
3489 divshape(1) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx1;
3490 divshape(2) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx2;
3492 divshape(3) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty0;
3493 divshape(4) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty1;
3494 divshape(5) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty2;
3496 divshape(6) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx2;
3497 divshape(7) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx1;
3498 divshape(8) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx0;
3500 divshape(9) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty2;
3501 divshape(10) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty1;
3502 divshape(11) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty0;
3504 divshape(12) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty0;
3505 divshape(13) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty1;
3506 divshape(14) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty2;
3508 divshape(15) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty0;
3509 divshape(16) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty1;
3510 divshape(17) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty2;
3512 divshape(18) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx0;
3513 divshape(19) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx1;
3514 divshape(20) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx2;
3516 divshape(21) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx0;
3517 divshape(22) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx1;
3518 divshape(23) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx2;
3521 const double RT2QuadFiniteElement::nk[24][2] =
3524 {0,-1}, {0,-1}, {0,-1},
3526 {1, 0}, {1, 0}, {1, 0},
3528 {0, 1}, {0, 1}, {0, 1},
3530 {-1,0}, {-1,0}, {-1,0},
3532 {1, 0}, {1, 0}, {1, 0},
3534 {1, 0}, {1, 0}, {1, 0},
3536 {0, 1}, {0, 1}, {0, 1},
3538 {0, 1}, {0, 1}, {0, 1}
3545 #ifdef MFEM_THREAD_SAFE
3551 for (k = 0; k < 24; k++)
3554 for (j = 0; j < 24; j++)
3557 if (j == k) { d -= 1.0; }
3558 if (fabs(d) > 1.0e-12)
3560 mfem::err <<
"RT2QuadFiniteElement::GetLocalInterpolation (...)\n"
3561 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
3577 for (k = 0; k < 24; k++)
3580 ip.
x = vk[0]; ip.
y = vk[1];
3583 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
3584 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
3585 for (j = 0; j < 24; j++)
3586 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
3598 #ifdef MFEM_THREAD_SAFE
3602 for (
int k = 0; k < 24; k++)
3610 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
3611 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
3627 shape(0) = 2. - 3. * x;
3628 shape(1) = 3. * x - 1.;
3642 const double p = 0.11270166537925831148;
3652 const double p = 0.11270166537925831148;
3653 const double w = 1./((1-2*p)*(1-2*p));
3656 shape(0) = (2*x-1)*(x-1+p)*w;
3657 shape(1) = 4*(x-1+p)*(p-x)*w;
3658 shape(2) = (2*x-1)*(x-p)*w;
3664 const double p = 0.11270166537925831148;
3665 const double w = 1./((1-2*p)*(1-2*p));
3668 dshape(0,0) = (-3+4*x+2*p)*w;
3669 dshape(1,0) = (4-8*x)*w;
3670 dshape(2,0) = (-1+4*x-2*p)*w;
3681 for (i = 1; i < m; i++)
3687 #ifndef MFEM_THREAD_SAFE
3692 for (i = 1; i <= m; i++)
3694 rwk(i) = rwk(i-1) * ( (double)(m) / (double)(i) );
3696 for (i = 0; i < m/2+1; i++)
3698 rwk(m-i) = ( rwk(i) *= rwk(m-i) );
3700 for (i = m-1; i >= 0; i -= 2)
3709 double w, wk, x = ip.
x;
3712 #ifdef MFEM_THREAD_SAFE
3716 k = (int) floor ( m * x + 0.5 );
3717 k = k > m ? m : k < 0 ? 0 : k;
3720 for (i = 0; i <= m; i++)
3723 wk *= ( rxxk(i) = x - (double)(i) / m );
3725 w = wk * ( rxxk(k) = x - (double)(k) / m );
3729 shape(0) = w * rwk(0) / rxxk(0);
3733 shape(0) = wk * rwk(0);
3737 shape(1) = w * rwk(m) / rxxk(m);
3741 shape(1) = wk * rwk(k);
3743 for (i = 1; i < m; i++)
3746 shape(i+1) = w * rwk(i) / rxxk(i);
3750 shape(k+1) = wk * rwk(k);
3757 double s, srx, w, wk, x = ip.
x;
3760 #ifdef MFEM_THREAD_SAFE
3764 k = (int) floor ( m * x + 0.5 );
3765 k = k > m ? m : k < 0 ? 0 : k;
3768 for (i = 0; i <= m; i++)
3771 wk *= ( rxxk(i) = x - (double)(i) / m );
3773 w = wk * ( rxxk(k) = x - (double)(k) / m );
3775 for (i = 0; i <= m; i++)
3777 rxxk(i) = 1.0 / rxxk(i);
3780 for (i = 0; i <= m; i++)
3789 dshape(0,0) = (s - w * rxxk(0)) * rwk(0) * rxxk(0);
3793 dshape(0,0) = wk * srx * rwk(0);
3797 dshape(1,0) = (s - w * rxxk(m)) * rwk(m) * rxxk(m);
3801 dshape(1,0) = wk * srx * rwk(k);
3803 for (i = 1; i < m; i++)
3806 dshape(i+1,0) = (s - w * rxxk(i)) * rwk(i) * rxxk(i);
3810 dshape(k+1,0) = wk * srx * rwk(k);
3839 double L0, L1, L2, L3;
3841 L1 = ip.
x; L2 = ip.
y; L3 = ip.
z; L0 = 1.0 - L1 - L2 - L3;
3842 shape(0) = 1.0 - 3.0 * L0;
3843 shape(1) = 1.0 - 3.0 * L1;
3844 shape(2) = 1.0 - 3.0 * L2;
3845 shape(3) = 1.0 - 3.0 * L3;
3851 dshape(0,0) = 3.0; dshape(0,1) = 3.0; dshape(0,2) = 3.0;
3852 dshape(1,0) = -3.0; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
3853 dshape(2,0) = 0.0; dshape(2,1) = -3.0; dshape(2,2) = 0.0;
3854 dshape(3,0) = 0.0; dshape(3,1) = 0.0; dshape(3,2) = -3.0;
3875 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3896 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3910 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3911 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3912 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3913 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3914 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3915 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3916 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3917 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3919 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3920 I[ 9] = 1; J[ 9] = 2; K[ 9] = 0;
3921 I[10] = 2; J[10] = 1; K[10] = 0;
3922 I[11] = 0; J[11] = 2; K[11] = 0;
3923 I[12] = 2; J[12] = 0; K[12] = 1;
3924 I[13] = 1; J[13] = 2; K[13] = 1;
3925 I[14] = 2; J[14] = 1; K[14] = 1;
3926 I[15] = 0; J[15] = 2; K[15] = 1;
3927 I[16] = 0; J[16] = 0; K[16] = 2;
3928 I[17] = 1; J[17] = 0; K[17] = 2;
3929 I[18] = 1; J[18] = 1; K[18] = 2;
3930 I[19] = 0; J[19] = 1; K[19] = 2;
3932 I[20] = 2; J[20] = 2; K[20] = 0;
3933 I[21] = 2; J[21] = 0; K[21] = 2;
3934 I[22] = 1; J[22] = 2; K[22] = 2;
3935 I[23] = 2; J[23] = 1; K[23] = 2;
3936 I[24] = 0; J[24] = 2; K[24] = 2;
3937 I[25] = 2; J[25] = 2; K[25] = 1;
3939 I[26] = 2; J[26] = 2; K[26] = 2;
3941 else if (degree == 3)
3947 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3948 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3949 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3950 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3951 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3952 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3953 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3954 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3956 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3957 I[ 9] = 3; J[ 9] = 0; K[ 9] = 0;
3958 I[10] = 1; J[10] = 2; K[10] = 0;
3959 I[11] = 1; J[11] = 3; K[11] = 0;
3960 I[12] = 2; J[12] = 1; K[12] = 0;
3961 I[13] = 3; J[13] = 1; K[13] = 0;
3962 I[14] = 0; J[14] = 2; K[14] = 0;
3963 I[15] = 0; J[15] = 3; K[15] = 0;
3964 I[16] = 2; J[16] = 0; K[16] = 1;
3965 I[17] = 3; J[17] = 0; K[17] = 1;
3966 I[18] = 1; J[18] = 2; K[18] = 1;
3967 I[19] = 1; J[19] = 3; K[19] = 1;
3968 I[20] = 2; J[20] = 1; K[20] = 1;
3969 I[21] = 3; J[21] = 1; K[21] = 1;
3970 I[22] = 0; J[22] = 2; K[22] = 1;
3971 I[23] = 0; J[23] = 3; K[23] = 1;
3972 I[24] = 0; J[24] = 0; K[24] = 2;
3973 I[25] = 0; J[25] = 0; K[25] = 3;
3974 I[26] = 1; J[26] = 0; K[26] = 2;
3975 I[27] = 1; J[27] = 0; K[27] = 3;
3976 I[28] = 1; J[28] = 1; K[28] = 2;
3977 I[29] = 1; J[29] = 1; K[29] = 3;
3978 I[30] = 0; J[30] = 1; K[30] = 2;
3979 I[31] = 0; J[31] = 1; K[31] = 3;
3981 I[32] = 2; J[32] = 3; K[32] = 0;
3982 I[33] = 3; J[33] = 3; K[33] = 0;
3983 I[34] = 2; J[34] = 2; K[34] = 0;
3984 I[35] = 3; J[35] = 2; K[35] = 0;
3985 I[36] = 2; J[36] = 0; K[36] = 2;
3986 I[37] = 3; J[37] = 0; K[37] = 2;
3987 I[38] = 2; J[38] = 0; K[38] = 3;
3988 I[39] = 3; J[39] = 0; K[39] = 3;
3989 I[40] = 1; J[40] = 2; K[40] = 2;
3990 I[41] = 1; J[41] = 3; K[41] = 2;
3991 I[42] = 1; J[42] = 2; K[42] = 3;
3992 I[43] = 1; J[43] = 3; K[43] = 3;
3993 I[44] = 3; J[44] = 1; K[44] = 2;
3994 I[45] = 2; J[45] = 1; K[45] = 2;
3995 I[46] = 3; J[46] = 1; K[46] = 3;
3996 I[47] = 2; J[47] = 1; K[47] = 3;
3997 I[48] = 0; J[48] = 3; K[48] = 2;
3998 I[49] = 0; J[49] = 2; K[49] = 2;
3999 I[50] = 0; J[50] = 3; K[50] = 3;
4000 I[51] = 0; J[51] = 2; K[51] = 3;
4001 I[52] = 2; J[52] = 2; K[52] = 1;
4002 I[53] = 3; J[53] = 2; K[53] = 1;
4003 I[54] = 2; J[54] = 3; K[54] = 1;
4004 I[55] = 3; J[55] = 3; K[55] = 1;
4006 I[56] = 2; J[56] = 2; K[56] = 2;
4007 I[57] = 3; J[57] = 2; K[57] = 2;
4008 I[58] = 3; J[58] = 3; K[58] = 2;
4009 I[59] = 2; J[59] = 3; K[59] = 2;
4010 I[60] = 2; J[60] = 2; K[60] = 3;
4011 I[61] = 3; J[61] = 2; K[61] = 3;
4012 I[62] = 3; J[62] = 3; K[62] = 3;
4013 I[63] = 2; J[63] = 3; K[63] = 3;
4017 mfem_error (
"LagrangeHexFiniteElement::LagrangeHexFiniteElement");
4021 dof1d = fe1d ->
GetDof();
4023 #ifndef MFEM_THREAD_SAFE
4033 for (
int n = 0; n <
Dof; n++)
4048 #ifdef MFEM_THREAD_SAFE
4049 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
4056 for (
int n = 0; n <
Dof; n++)
4058 shape(n) = shape1dx(I[n]) * shape1dy(J[n]) * shape1dz(K[n]);
4069 #ifdef MFEM_THREAD_SAFE
4070 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
4071 DenseMatrix dshape1dx(dof1d,1), dshape1dy(dof1d,1), dshape1dz(dof1d,1);
4082 for (
int n = 0; n <
Dof; n++)
4084 dshape(n,0) = dshape1dx(I[n],0) * shape1dy(J[n]) * shape1dz(K[n]);
4085 dshape(n,1) = shape1dx(I[n]) * dshape1dy(J[n],0) * shape1dz(K[n]);
4086 dshape(n,2) = shape1dx(I[n]) * shape1dy(J[n]) * dshape1dz(K[n],0);
4115 shape(0) = 1.0 - 2.0 * x;
4122 shape(1) = 2.0 * x - 1.0;
4123 shape(2) = 2.0 - 2.0 * x;
4134 dshape(0,0) = - 2.0;
4142 dshape(2,0) = - 2.0;
4169 L0 = 2.0 * ( 1. - ip.
x - ip.
y );
4170 L1 = 2.0 * ( ip.
x );
4171 L2 = 2.0 * ( ip.
y );
4180 for (i = 0; i < 6; i++)
4187 shape(0) = L0 - 1.0;
4194 shape(1) = L1 - 1.0;
4201 shape(2) = L2 - 1.0;
4205 shape(3) = 1.0 - L2;
4206 shape(4) = 1.0 - L0;
4207 shape(5) = 1.0 - L1;
4217 L0 = 2.0 * ( 1. - ip.
x - ip.
y );
4218 L1 = 2.0 * ( ip.
x );
4219 L2 = 2.0 * ( ip.
y );
4221 double DL0[2], DL1[2], DL2[2];
4222 DL0[0] = -2.0; DL0[1] = -2.0;
4223 DL1[0] = 2.0; DL1[1] = 0.0;
4224 DL2[0] = 0.0; DL2[1] = 2.0;
4226 for (i = 0; i < 6; i++)
4227 for (j = 0; j < 2; j++)
4234 for (j = 0; j < 2; j++)
4236 dshape(0,j) = DL0[j];
4237 dshape(3,j) = DL1[j];
4238 dshape(5,j) = DL2[j];
4243 for (j = 0; j < 2; j++)
4245 dshape(3,j) = DL0[j];
4246 dshape(1,j) = DL1[j];
4247 dshape(4,j) = DL2[j];
4252 for (j = 0; j < 2; j++)
4254 dshape(5,j) = DL0[j];
4255 dshape(4,j) = DL1[j];
4256 dshape(2,j) = DL2[j];
4261 for (j = 0; j < 2; j++)
4263 dshape(3,j) = - DL2[j];
4264 dshape(4,j) = - DL0[j];
4265 dshape(5,j) = - DL1[j];
4310 double L0, L1, L2, L3, L4, L5;
4311 L0 = 2.0 * ( 1. - ip.
x - ip.
y - ip.
z );
4312 L1 = 2.0 * ( ip.
x );
4313 L2 = 2.0 * ( ip.
y );
4314 L3 = 2.0 * ( ip.
z );
4315 L4 = 2.0 * ( ip.
x + ip.
y );
4316 L5 = 2.0 * ( ip.
y + ip.
z );
4329 for (i = 0; i < 10; i++)
4336 shape(0) = L0 - 1.0;
4344 shape(1) = L1 - 1.0;
4352 shape(2) = L2 - 1.0;
4360 shape(3) = L3 - 1.0;
4362 else if ((L4 <= 1.0) && (L5 <= 1.0))
4364 shape(4) = 1.0 - L5;
4366 shape(6) = 1.0 - L4;
4367 shape(8) = 1.0 - L0;
4369 else if ((L4 >= 1.0) && (L5 <= 1.0))
4371 shape(4) = 1.0 - L5;
4372 shape(5) = 1.0 - L1;
4373 shape(7) = L4 - 1.0;
4376 else if ((L4 <= 1.0) && (L5 >= 1.0))
4378 shape(5) = 1.0 - L3;
4379 shape(6) = 1.0 - L4;
4381 shape(9) = L5 - 1.0;
4383 else if ((L4 >= 1.0) && (L5 >= 1.0))
4386 shape(7) = L4 - 1.0;
4387 shape(8) = 1.0 - L2;
4388 shape(9) = L5 - 1.0;
4397 double L0, L1, L2, L3, L4, L5;
4398 L0 = 2.0 * ( 1. - ip.
x - ip.
y - ip.
z );
4399 L1 = 2.0 * ( ip.
x );
4400 L2 = 2.0 * ( ip.
y );
4401 L3 = 2.0 * ( ip.
z );
4402 L4 = 2.0 * ( ip.
x + ip.
y );
4403 L5 = 2.0 * ( ip.
y + ip.
z );
4405 double DL0[3], DL1[3], DL2[3], DL3[3], DL4[3], DL5[3];
4406 DL0[0] = -2.0; DL0[1] = -2.0; DL0[2] = -2.0;
4407 DL1[0] = 2.0; DL1[1] = 0.0; DL1[2] = 0.0;
4408 DL2[0] = 0.0; DL2[1] = 2.0; DL2[2] = 0.0;
4409 DL3[0] = 0.0; DL3[1] = 0.0; DL3[2] = 2.0;
4410 DL4[0] = 2.0; DL4[1] = 2.0; DL4[2] = 0.0;
4411 DL5[0] = 0.0; DL5[1] = 2.0; DL5[2] = 2.0;
4413 for (i = 0; i < 10; i++)
4414 for (j = 0; j < 3; j++)
4421 for (j = 0; j < 3; j++)
4423 dshape(0,j) = DL0[j];
4424 dshape(4,j) = DL1[j];
4425 dshape(5,j) = DL2[j];
4426 dshape(6,j) = DL3[j];
4431 for (j = 0; j < 3; j++)
4433 dshape(4,j) = DL0[j];
4434 dshape(1,j) = DL1[j];
4435 dshape(7,j) = DL2[j];
4436 dshape(8,j) = DL3[j];
4441 for (j = 0; j < 3; j++)
4443 dshape(5,j) = DL0[j];
4444 dshape(7,j) = DL1[j];
4445 dshape(2,j) = DL2[j];
4446 dshape(9,j) = DL3[j];
4451 for (j = 0; j < 3; j++)
4453 dshape(6,j) = DL0[j];
4454 dshape(8,j) = DL1[j];
4455 dshape(9,j) = DL2[j];
4456 dshape(3,j) = DL3[j];
4459 else if ((L4 <= 1.0) && (L5 <= 1.0))
4461 for (j = 0; j < 3; j++)
4463 dshape(4,j) = - DL5[j];
4464 dshape(5,j) = DL2[j];
4465 dshape(6,j) = - DL4[j];
4466 dshape(8,j) = - DL0[j];
4469 else if ((L4 >= 1.0) && (L5 <= 1.0))
4471 for (j = 0; j < 3; j++)
4473 dshape(4,j) = - DL5[j];
4474 dshape(5,j) = - DL1[j];
4475 dshape(7,j) = DL4[j];
4476 dshape(8,j) = DL3[j];
4479 else if ((L4 <= 1.0) && (L5 >= 1.0))
4481 for (j = 0; j < 3; j++)
4483 dshape(5,j) = - DL3[j];
4484 dshape(6,j) = - DL4[j];
4485 dshape(8,j) = DL1[j];
4486 dshape(9,j) = DL5[j];
4489 else if ((L4 >= 1.0) && (L5 >= 1.0))
4491 for (j = 0; j < 3; j++)
4493 dshape(5,j) = DL0[j];
4494 dshape(7,j) = DL4[j];
4495 dshape(8,j) = - DL2[j];
4496 dshape(9,j) = DL5[j];
4529 double x = ip.
x, y = ip.
y;
4531 Lx = 2.0 * ( 1. - x );
4532 Ly = 2.0 * ( 1. - y );
4541 for (i = 0; i < 9; i++)
4546 if ((x <= 0.5) && (y <= 0.5))
4548 shape(0) = (Lx - 1.0) * (Ly - 1.0);
4549 shape(4) = (2.0 - Lx) * (Ly - 1.0);
4550 shape(8) = (2.0 - Lx) * (2.0 - Ly);
4551 shape(7) = (Lx - 1.0) * (2.0 - Ly);
4553 else if ((x >= 0.5) && (y <= 0.5))
4555 shape(4) = Lx * (Ly - 1.0);
4556 shape(1) = (1.0 - Lx) * (Ly - 1.0);
4557 shape(5) = (1.0 - Lx) * (2.0 - Ly);
4558 shape(8) = Lx * (2.0 - Ly);
4560 else if ((x >= 0.5) && (y >= 0.5))
4562 shape(8) = Lx * Ly ;
4563 shape(5) = (1.0 - Lx) * Ly ;
4564 shape(2) = (1.0 - Lx) * (1.0 - Ly);
4565 shape(6) = Lx * (1.0 - Ly);
4567 else if ((x <= 0.5) && (y >= 0.5))
4569 shape(7) = (Lx - 1.0) * Ly ;
4570 shape(8) = (2.0 - Lx) * Ly ;
4571 shape(6) = (2.0 - Lx) * (1.0 - Ly);
4572 shape(3) = (Lx - 1.0) * (1.0 - Ly);
4580 double x = ip.
x, y = ip.
y;
4582 Lx = 2.0 * ( 1. - x );
4583 Ly = 2.0 * ( 1. - y );
4585 for (i = 0; i < 9; i++)
4586 for (j = 0; j < 2; j++)
4591 if ((x <= 0.5) && (y <= 0.5))
4593 dshape(0,0) = 2.0 * (1.0 - Ly);
4594 dshape(0,1) = 2.0 * (1.0 - Lx);
4596 dshape(4,0) = 2.0 * (Ly - 1.0);
4597 dshape(4,1) = -2.0 * (2.0 - Lx);
4599 dshape(8,0) = 2.0 * (2.0 - Ly);
4600 dshape(8,1) = 2.0 * (2.0 - Lx);
4602 dshape(7,0) = -2.0 * (2.0 - Ly);
4603 dshape(7,0) = 2.0 * (Lx - 1.0);
4605 else if ((x >= 0.5) && (y <= 0.5))
4607 dshape(4,0) = -2.0 * (Ly - 1.0);
4608 dshape(4,1) = -2.0 * Lx;
4610 dshape(1,0) = 2.0 * (Ly - 1.0);
4611 dshape(1,1) = -2.0 * (1.0 - Lx);
4613 dshape(5,0) = 2.0 * (2.0 - Ly);
4614 dshape(5,1) = 2.0 * (1.0 - Lx);
4616 dshape(8,0) = -2.0 * (2.0 - Ly);
4617 dshape(8,1) = 2.0 * Lx;
4619 else if ((x >= 0.5) && (y >= 0.5))
4621 dshape(8,0) = -2.0 * Ly;
4622 dshape(8,1) = -2.0 * Lx;
4624 dshape(5,0) = 2.0 * Ly;
4625 dshape(5,1) = -2.0 * (1.0 - Lx);
4627 dshape(2,0) = 2.0 * (1.0 - Ly);
4628 dshape(2,1) = 2.0 * (1.0 - Lx);
4630 dshape(6,0) = -2.0 * (1.0 - Ly);
4631 dshape(6,1) = 2.0 * Lx;
4633 else if ((x <= 0.5) && (y >= 0.5))
4635 dshape(7,0) = -2.0 * Ly;
4636 dshape(7,1) = -2.0 * (Lx - 1.0);
4638 dshape(8,0) = 2.0 * Ly ;
4639 dshape(8,1) = -2.0 * (2.0 - Lx);
4641 dshape(6,0) = 2.0 * (1.0 - Ly);
4642 dshape(6,1) = 2.0 * (2.0 - Lx);
4644 dshape(3,0) = -2.0 * (1.0 - Ly);
4645 dshape(3,1) = 2.0 * (Lx - 1.0);
4656 I[ 0] = 0.0; J[ 0] = 0.0; K[ 0] = 0.0;
4657 I[ 1] = 1.0; J[ 1] = 0.0; K[ 1] = 0.0;
4658 I[ 2] = 1.0; J[ 2] = 1.0; K[ 2] = 0.0;
4659 I[ 3] = 0.0; J[ 3] = 1.0; K[ 3] = 0.0;
4660 I[ 4] = 0.0; J[ 4] = 0.0; K[ 4] = 1.0;
4661 I[ 5] = 1.0; J[ 5] = 0.0; K[ 5] = 1.0;
4662 I[ 6] = 1.0; J[ 6] = 1.0; K[ 6] = 1.0;
4663 I[ 7] = 0.0; J[ 7] = 1.0; K[ 7] = 1.0;
4665 I[ 8] = 0.5; J[ 8] = 0.0; K[ 8] = 0.0;
4666 I[ 9] = 1.0; J[ 9] = 0.5; K[ 9] = 0.0;
4667 I[10] = 0.5; J[10] = 1.0; K[10] = 0.0;
4668 I[11] = 0.0; J[11] = 0.5; K[11] = 0.0;
4669 I[12] = 0.5; J[12] = 0.0; K[12] = 1.0;
4670 I[13] = 1.0; J[13] = 0.5; K[13] = 1.0;
4671 I[14] = 0.5; J[14] = 1.0; K[14] = 1.0;
4672 I[15] = 0.0; J[15] = 0.5; K[15] = 1.0;
4673 I[16] = 0.0; J[16] = 0.0; K[16] = 0.5;
4674 I[17] = 1.0; J[17] = 0.0; K[17] = 0.5;
4675 I[18] = 1.0; J[18] = 1.0; K[18] = 0.5;
4676 I[19] = 0.0; J[19] = 1.0; K[19] = 0.5;
4678 I[20] = 0.5; J[20] = 0.5; K[20] = 0.0;
4679 I[21] = 0.5; J[21] = 0.0; K[21] = 0.5;
4680 I[22] = 1.0; J[22] = 0.5; K[22] = 0.5;
4681 I[23] = 0.5; J[23] = 1.0; K[23] = 0.5;
4682 I[24] = 0.0; J[24] = 0.5; K[24] = 0.5;
4683 I[25] = 0.5; J[25] = 0.5; K[25] = 1.0;
4685 I[26] = 0.5; J[26] = 0.5; K[26] = 0.5;
4687 for (
int n = 0; n < 27; n++)
4700 double x = ip.
x, y = ip.
y, z = ip.
z;
4702 for (i = 0; i < 27; i++)
4707 if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5))
4722 else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5))
4737 else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5))
4752 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5))
4767 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5))
4782 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5))
4797 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5))
4828 shape(N[0]) = Lx * Ly * Lz;
4829 shape(N[1]) = (1 - Lx) * Ly * Lz;
4830 shape(N[2]) = (1 - Lx) * (1 - Ly) * Lz;
4831 shape(N[3]) = Lx * (1 - Ly) * Lz;
4832 shape(N[4]) = Lx * Ly * (1 - Lz);
4833 shape(N[5]) = (1 - Lx) * Ly * (1 - Lz);
4834 shape(N[6]) = (1 - Lx) * (1 - Ly) * (1 - Lz);
4835 shape(N[7]) = Lx * (1 - Ly) * (1 - Lz);
4843 double x = ip.
x, y = ip.
y, z = ip.
z;
4845 for (i = 0; i < 27; i++)
4846 for (j = 0; j < 3; j++)
4851 if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5))
4866 else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5))
4881 else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5))
4896 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5))
4911 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5))
4926 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5))
4941 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5))
4972 dshape(N[0],0) = -2.0 * Ly * Lz ;
4973 dshape(N[0],1) = -2.0 * Lx * Lz ;
4974 dshape(N[0],2) = -2.0 * Lx * Ly ;
4976 dshape(N[1],0) = 2.0 * Ly * Lz ;
4977 dshape(N[1],1) = -2.0 * (1 - Lx) * Lz ;
4978 dshape(N[1],2) = -2.0 * (1 - Lx) * Ly ;
4980 dshape(N[2],0) = 2.0 * (1 - Ly) * Lz ;
4981 dshape(N[2],1) = 2.0 * (1 - Lx) * Lz ;
4982 dshape(N[2],2) = -2.0 * (1 - Lx) * (1 - Ly);
4984 dshape(N[3],0) = -2.0 * (1 - Ly) * Lz ;
4985 dshape(N[3],1) = 2.0 * Lx * Lz ;
4986 dshape(N[3],2) = -2.0 * Lx * (1 - Ly);
4988 dshape(N[4],0) = -2.0 * Ly * (1 - Lz);
4989 dshape(N[4],1) = -2.0 * Lx * (1 - Lz);
4990 dshape(N[4],2) = 2.0 * Lx * Ly ;
4992 dshape(N[5],0) = 2.0 * Ly * (1 - Lz);
4993 dshape(N[5],1) = -2.0 * (1 - Lx) * (1 - Lz);
4994 dshape(N[5],2) = 2.0 * (1 - Lx) * Ly ;
4996 dshape(N[6],0) = 2.0 * (1 - Ly) * (1 - Lz);
4997 dshape(N[6],1) = 2.0 * (1 - Lx) * (1 - Lz);
4998 dshape(N[6],2) = 2.0 * (1 - Lx) * (1 - Ly);
5000 dshape(N[7],0) = -2.0 * (1 - Ly) * (1 - Lz);
5001 dshape(N[7],1) = 2.0 * Lx * (1 - Lz);
5002 dshape(N[7],2) = 2.0 * Lx * (1 - Ly);
5062 double x = ip.
x, y = ip.
y, z = ip.
z;
5064 shape(0,0) = (1. - y) * (1. - z);
5068 shape(2,0) = y * (1. - z);
5072 shape(4,0) = z * (1. - y);
5081 shape(1,1) = x * (1. - z);
5085 shape(3,1) = (1. - x) * (1. - z);
5093 shape(7,1) = (1. - x) * z;
5098 shape(8,2) = (1. - x) * (1. - y);
5102 shape(9,2) = x * (1. - y);
5106 shape(10,2) = x * y;
5110 shape(11,2) = y * (1. - x);
5118 double x = ip.
x, y = ip.
y, z = ip.
z;
5120 curl_shape(0,0) = 0.;
5121 curl_shape(0,1) = y - 1.;
5122 curl_shape(0,2) = 1. - z;
5124 curl_shape(2,0) = 0.;
5125 curl_shape(2,1) = -y;
5126 curl_shape(2,2) = z - 1.;
5128 curl_shape(4,0) = 0;
5129 curl_shape(4,1) = 1. - y;
5130 curl_shape(4,2) = z;
5132 curl_shape(6,0) = 0.;
5133 curl_shape(6,1) = y;
5134 curl_shape(6,2) = -z;
5136 curl_shape(1,0) = x;
5137 curl_shape(1,1) = 0.;
5138 curl_shape(1,2) = 1. - z;
5140 curl_shape(3,0) = 1. - x;
5141 curl_shape(3,1) = 0.;
5142 curl_shape(3,2) = z - 1.;
5144 curl_shape(5,0) = -x;
5145 curl_shape(5,1) = 0.;
5146 curl_shape(5,2) = z;
5148 curl_shape(7,0) = x - 1.;
5149 curl_shape(7,1) = 0.;
5150 curl_shape(7,2) = -z;
5152 curl_shape(8,0) = x - 1.;
5153 curl_shape(8,1) = 1. - y;
5154 curl_shape(8,2) = 0.;
5156 curl_shape(9,0) = -x;
5157 curl_shape(9,1) = y - 1.;
5158 curl_shape(9,2) = 0;
5160 curl_shape(10,0) = x;
5161 curl_shape(10,1) = -y;
5162 curl_shape(10,2) = 0.;
5164 curl_shape(11,0) = 1. - x;
5165 curl_shape(11,1) = y;
5166 curl_shape(11,2) = 0.;
5169 const double Nedelec1HexFiniteElement::tk[12][3] =
5171 {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
5172 {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
5173 {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}
5180 #ifdef MFEM_THREAD_SAFE
5185 for (k = 0; k < 12; k++)
5188 for (j = 0; j < 12; j++)
5190 double d = (
vshape(j,0)*tk[k][0] +
vshape(j,1)*tk[k][1] +
5192 if (j == k) { d -= 1.0; }
5193 if (fabs(d) > 1.0e-12)
5195 mfem::err <<
"Nedelec1HexFiniteElement::GetLocalInterpolation (...)\n"
5196 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5204 ip.
x = ip.
y = ip.
z = 0.0;
5211 for (k = 0; k < 12; k++)
5214 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5217 vk[0] =
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2];
5218 vk[1] =
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2];
5219 vk[2] =
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2];
5220 for (j = 0; j < 12; j++)
5221 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5222 vshape(j,2)*vk[2])) < 1.0e-12)
5236 for (
int k = 0; k < 12; k++)
5244 vk[0] * (
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2] ) +
5245 vk[1] * (
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2] ) +
5246 vk[2] * (
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2] );
5283 double x = ip.
x, y = ip.
y, z = ip.
z;
5285 shape(0,0) = 1. - y - z;
5290 shape(1,1) = 1. - x - z;
5295 shape(2,2) = 1. - x - y;
5314 curl_shape(0,0) = 0.;
5315 curl_shape(0,1) = -2.;
5316 curl_shape(0,2) = 2.;
5318 curl_shape(1,0) = 2.;
5319 curl_shape(1,1) = 0.;
5320 curl_shape(1,2) = -2.;
5322 curl_shape(2,0) = -2.;
5323 curl_shape(2,1) = 2.;
5324 curl_shape(2,2) = 0.;
5326 curl_shape(3,0) = 0.;
5327 curl_shape(3,1) = 0.;
5328 curl_shape(3,2) = 2.;
5330 curl_shape(4,0) = 0.;
5331 curl_shape(4,1) = -2.;
5332 curl_shape(4,2) = 0.;
5334 curl_shape(5,0) = 2.;
5335 curl_shape(5,1) = 0.;
5336 curl_shape(5,2) = 0.;
5339 const double Nedelec1TetFiniteElement::tk[6][3] =
5340 {{1,0,0}, {0,1,0}, {0,0,1}, {-1,1,0}, {-1,0,1}, {0,-1,1}};
5346 #ifdef MFEM_THREAD_SAFE
5351 for (k = 0; k < 6; k++)
5354 for (j = 0; j < 6; j++)
5356 double d = (
vshape(j,0)*tk[k][0] +
vshape(j,1)*tk[k][1] +
5358 if (j == k) { d -= 1.0; }
5359 if (fabs(d) > 1.0e-12)
5361 mfem::err <<
"Nedelec1TetFiniteElement::GetLocalInterpolation (...)\n"
5362 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5370 ip.
x = ip.
y = ip.
z = 0.0;
5377 for (k = 0; k < 6; k++)
5380 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5383 vk[0] =
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2];
5384 vk[1] =
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2];
5385 vk[2] =
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2];
5386 for (j = 0; j < 6; j++)
5387 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5388 vshape(j,2)*vk[2])) < 1.0e-12)
5402 for (
int k = 0; k < 6; k++)
5410 vk[0] * (
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2] ) +
5411 vk[1] * (
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2] ) +
5412 vk[2] * (
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2] );
5449 double x = ip.
x, y = ip.
y, z = ip.
z;
5453 shape(0,2) = z - 1.;
5456 shape(1,1) = y - 1.;
5467 shape(4,0) = x - 1.;
5487 const double RT0HexFiniteElement::nk[6][3] =
5488 {{0,0,-1}, {0,-1,0}, {1,0,0}, {0,1,0}, {-1,0,0}, {0,0,1}};
5494 #ifdef MFEM_THREAD_SAFE
5500 for (k = 0; k < 6; k++)
5503 for (j = 0; j < 6; j++)
5505 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5507 if (j == k) { d -= 1.0; }
5508 if (fabs(d) > 1.0e-12)
5510 mfem::err <<
"RT0HexFiniteElement::GetLocalInterpolation (...)\n"
5511 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5519 ip.
x = ip.
y = ip.
z = 0.0;
5527 for (k = 0; k < 6; k++)
5530 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5533 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2];
5534 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2];
5535 vk[2] =
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2];
5536 for (j = 0; j < 6; j++)
5537 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5538 vshape(j,2)*vk[2])) < 1.0e-12)
5551 #ifdef MFEM_THREAD_SAFE
5555 for (
int k = 0; k < 6; k++)
5564 vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2] ) +
5565 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2] ) +
5566 vk[2] * (
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2] );
5695 double x = ip.
x, y = ip.
y, z = ip.
z;
5699 shape(2,2) = -(1. - 3.*z + 2.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5702 shape(3,2) = -(1. - 3.*z + 2.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5705 shape(0,2) = -(1. - 3.*z + 2.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5708 shape(1,2) = -(1. - 3.*z + 2.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5711 shape(4,1) = -(1. - 3.*y + 2.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5714 shape(5,1) = -(1. - 3.*y + 2.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5717 shape(6,1) = -(1. - 3.*y + 2.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5720 shape(7,1) = -(1. - 3.*y + 2.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5723 shape(8,0) = (-x + 2.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5726 shape(9,0) = (-x + 2.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5729 shape(10,0) = (-x + 2.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5732 shape(11,0) = (-x + 2.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5737 shape(13,1) = (-y + 2.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5740 shape(12,1) = (-y + 2.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5743 shape(15,1) = (-y + 2.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5746 shape(14,1) = (-y + 2.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5749 shape(17,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5752 shape(16,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5755 shape(19,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5758 shape(18,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5764 shape(20,2) = (-z + 2.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5767 shape(21,2) = (-z + 2.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5770 shape(22,2) = (-z + 2.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5773 shape(23,2) = (-z + 2.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5775 shape(24,0) = (4.*x - 4.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5778 shape(25,0) = (4.*x - 4.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5781 shape(26,0) = (4.*x - 4.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5784 shape(27,0) = (4.*x - 4.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5789 shape(28,1) = (4.*y - 4.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5792 shape(29,1) = (4.*y - 4.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5795 shape(30,1) = (4.*y - 4.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5798 shape(31,1) = (4.*y - 4.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5803 shape(32,2) = (4.*z - 4.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5806 shape(33,2) = (4.*z - 4.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5809 shape(34,2) = (4.*z - 4.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5812 shape(35,2) = (4.*z - 4.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5818 double x = ip.
x, y = ip.
y, z = ip.
z;
5820 divshape(2) = -(-3. + 4.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5821 divshape(3) = -(-3. + 4.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5822 divshape(0) = -(-3. + 4.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5823 divshape(1) = -(-3. + 4.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5825 divshape(4) = -(-3. + 4.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5826 divshape(5) = -(-3. + 4.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5827 divshape(6) = -(-3. + 4.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5828 divshape(7) = -(-3. + 4.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5830 divshape(8) = (-1. + 4.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5831 divshape(9) = (-1. + 4.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5832 divshape(10) = (-1. + 4.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5833 divshape(11) = (-1. + 4.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5835 divshape(13) = (-1. + 4.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5836 divshape(12) = (-1. + 4.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5837 divshape(15) = (-1. + 4.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5838 divshape(14) = (-1. + 4.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5840 divshape(17) = -(-3. + 4.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5841 divshape(16) = -(-3. + 4.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5842 divshape(19) = -(-3. + 4.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5843 divshape(18) = -(-3. + 4.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5845 divshape(20) = (-1. + 4.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5846 divshape(21) = (-1. + 4.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5847 divshape(22) = (-1. + 4.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5848 divshape(23) = (-1. + 4.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5850 divshape(24) = ( 4. - 8.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5851 divshape(25) = ( 4. - 8.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5852 divshape(26) = ( 4. - 8.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5853 divshape(27) = ( 4. - 8.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5855 divshape(28) = ( 4. - 8.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5856 divshape(29) = ( 4. - 8.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5857 divshape(30) = ( 4. - 8.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5858 divshape(31) = ( 4. - 8.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5860 divshape(32) = ( 4. - 8.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5861 divshape(33) = ( 4. - 8.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5862 divshape(34) = ( 4. - 8.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5863 divshape(35) = ( 4. - 8.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5866 const double RT1HexFiniteElement::nk[36][3] =
5868 {0, 0,-1}, {0, 0,-1}, {0, 0,-1}, {0, 0,-1},
5869 {0,-1, 0}, {0,-1, 0}, {0,-1, 0}, {0,-1, 0},
5870 {1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0},
5871 {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0},
5872 {-1,0, 0}, {-1,0, 0}, {-1,0, 0}, {-1,0, 0},
5873 {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1},
5874 {1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0},
5875 {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0},
5876 {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1}
5883 #ifdef MFEM_THREAD_SAFE
5889 for (k = 0; k < 36; k++)
5892 for (j = 0; j < 36; j++)
5894 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5896 if (j == k) { d -= 1.0; }
5897 if (fabs(d) > 1.0e-12)
5899 mfem::err <<
"RT0HexFiniteElement::GetLocalInterpolation (...)\n"
5900 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5908 ip.
x = ip.
y = ip.
z = 0.0;
5916 for (k = 0; k < 36; k++)
5919 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5922 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2];
5923 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2];
5924 vk[2] =
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2];
5925 for (j = 0; j < 36; j++)
5926 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5927 vshape(j,2)*vk[2])) < 1.0e-12)
5940 #ifdef MFEM_THREAD_SAFE
5944 for (
int k = 0; k < 36; k++)
5953 vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2] ) +
5954 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2] ) +
5955 vk[2] * (
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2] );
5983 double x2 = 2.0*ip.
x, y2 = 2.0*ip.
y, z2 = 2.0*ip.
z;
5989 shape(1,0) = x2 - 2.0;
5994 shape(2,1) = y2 - 2.0;
5999 shape(3,2) = z2 - 2.0;
6011 const double RT0TetFiniteElement::nk[4][3] =
6012 {{.5,.5,.5}, {-.5,0,0}, {0,-.5,0}, {0,0,-.5}};
6018 #ifdef MFEM_THREAD_SAFE
6024 for (k = 0; k < 4; k++)
6027 for (j = 0; j < 4; j++)
6029 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
6031 if (j == k) { d -= 1.0; }
6032 if (fabs(d) > 1.0e-12)
6034 mfem::err <<
"RT0TetFiniteElement::GetLocalInterpolation (...)\n"
6035 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
6043 ip.
x = ip.
y = ip.
z = 0.0;
6051 for (k = 0; k < 4; k++)
6054 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
6057 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2];
6058 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2];
6059 vk[2] =
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2];
6060 for (j = 0; j < 4; j++)
6061 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
6062 vshape(j,2)*vk[2])) < 1.0e-12)
6075 #ifdef MFEM_THREAD_SAFE
6079 for (
int k = 0; k < 4; k++)
6088 vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2] ) +
6089 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2] ) +
6090 vk[2] * (
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2] );
6125 double x = 2. * ip.
x - 1.;
6126 double y = 2. * ip.
y - 1.;
6127 double z = 2. * ip.
z - 1.;
6128 double f5 = x * x - y * y;
6129 double f6 = y * y - z * z;
6131 shape(0) = (1./6.) * (1. - 3. * z - f5 - 2. * f6);
6132 shape(1) = (1./6.) * (1. - 3. * y - f5 + f6);
6133 shape(2) = (1./6.) * (1. + 3. * x + 2. * f5 + f6);
6134 shape(3) = (1./6.) * (1. + 3. * y - f5 + f6);
6135 shape(4) = (1./6.) * (1. - 3. * x + 2. * f5 + f6);
6136 shape(5) = (1./6.) * (1. + 3. * z - f5 - 2. * f6);
6142 const double a = 2./3.;
6144 double xt = a * (1. - 2. * ip.
x);
6145 double yt = a * (1. - 2. * ip.
y);
6146 double zt = a * (1. - 2. * ip.
z);
6150 dshape(0,2) = -1. - 2. * zt;
6153 dshape(1,1) = -1. - 2. * yt;
6156 dshape(2,0) = 1. - 2. * xt;
6161 dshape(3,1) = 1. - 2. * yt;
6164 dshape(4,0) = -1. - 2. * xt;
6170 dshape(5,2) = 1. - 2. * zt;
6175 : x(p + 1), w(p + 1)
6181 for (
int i = 0; i <= p; i++)
6184 for (
int j = 0; j <= p; j++)
6197 for (
int i = 0; i <= p; i++)
6199 for (
int j = 0; j < i; j++)
6201 double xij = x(i) - x(j);
6206 for (
int i = 0; i <= p; i++)
6213 for (
int i = 0; i < p; i++)
6217 mfem_error(
"Poly_1D::Basis::Basis : nodes are not increasing!");
6233 int i, k, p = x.Size() - 1;
6243 for (k = 0; k < p; k++)
6245 if (y >= (x(k) + x(k+1))/2)
6251 for (i = k+1; i <= p; i++)
6258 l = lk * (y - x(k));
6260 for (i = 0; i < k; i++)
6262 u(i) = l * w(i) / (y - x(i));
6265 for (i++; i <= p; i++)
6267 u(i) = l * w(i) / (y - x(i));
6282 int i, k, p = x.Size() - 1;
6283 double l, lp, lk, sk, si;
6293 for (k = 0; k < p; k++)
6295 if (y >= (x(k) + x(k+1))/2)
6301 for (i = k+1; i <= p; i++)
6308 l = lk * (y - x(k));
6311 for (i = 0; i < k; i++)
6313 si = 1.0/(y - x(i));
6315 u(i) = l * si * w(i);
6318 for (i++; i <= p; i++)
6320 si = 1.0/(y - x(i));
6322 u(i) = l * si * w(i);
6326 for (i = 0; i < k; i++)
6328 d(i) = (lp * w(i) - u(i))/(y - x(i));
6331 for (i++; i <= p; i++)
6333 d(i) = (lp * w(i) - u(i))/(y - x(i));
6343 for (
int i = 0; i <= p; i++)
6345 binom(i,0) = binom(i,i) = 1;
6346 for (
int j = 1; j < i; j++)
6348 binom(i,j) = binom(i-1,j) + binom(i-1,j-1);
6357 for (
int i = 0; i <= p; i++)
6360 double s = sin(M_PI_2*(i + 0.5)/(p + 1));
6365 void Poly_1D::CalcMono(
const int p,
const double x,
double *u)
6369 for (
int n = 1; n <= p; n++)
6375 void Poly_1D::CalcMono(
const int p,
const double x,
double *u,
double *d)
6380 for (
int n = 1; n <= p; n++)
6397 const int *b =
Binom(p);
6400 for (i = 1; i < p; i++)
6407 for (i--; i > 0; i--)
6417 double *u,
double *d)
6427 const int *b =
Binom(p);
6428 const double xpy = x + y, ptx = p*x;
6431 for (i = 1; i < p; i++)
6433 d[i] = b[i]*z*(i*xpy - ptx);
6440 for (i--; i > 0; i--)
6461 const int *b =
Binom(p);
6462 const double xpy = x + y, ptx = p*x;
6465 for (i = 1; i < p; i++)
6467 d[i] = b[i]*z*(i*xpy - ptx);
6472 for (i--; i > 0; i--)
6481 void Poly_1D::CalcLegendre(
const int p,
const double x,
double *u)
6487 if (p == 0) {
return; }
6488 u[1] = z = 2.*x - 1.;
6489 for (
int n = 1; n < p; n++)
6491 u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
6495 void Poly_1D::CalcLegendre(
const int p,
const double x,
double *u,
double *d)
6504 if (p == 0) {
return; }
6505 u[1] = z = 2.*x - 1.;
6507 for (
int n = 1; n < p; n++)
6509 u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
6510 d[n+1] = (4*n + 2)*u[n] + d[n-1];
6514 void Poly_1D::CalcChebyshev(
const int p,
const double x,
double *u)
6521 if (p == 0) {
return; }
6522 u[1] = z = 2.*x - 1.;
6523 for (
int n = 1; n < p; n++)
6525 u[n+1] = 2*z*u[n] - u[n-1];
6529 void Poly_1D::CalcChebyshev(
const int p,
const double x,
double *u,
double *d)
6542 if (p == 0) {
return; }
6543 u[1] = z = 2.*x - 1.;
6545 for (
int n = 1; n < p; n++)
6547 u[n+1] = 2*z*u[n] - u[n-1];
6548 d[n+1] = (n + 1)*(z*d[n]/n + 2*u[n]);
6556 if (points_container.find(type) == points_container.end())
6561 if (pts.
Size() <= p)
6567 pts[p] =
new double[p + 1];
6577 if ( bases_container.find(type) == bases_container.end() )
6583 if (bases.
Size() <= p)
6587 if (bases[p] == NULL)
6596 for (std::map<
int,
Array<double*>*>::iterator it = points_container.begin();
6597 it != points_container.end() ; ++it)
6600 for (
int i = 0 ; i < pts.
Size() ; ++i )
6607 for (std::map<
int,
Array<Basis*>*>::iterator it = bases_container.begin();
6608 it != bases_container.end() ; ++it )
6611 for (
int i = 0 ; i < bases.
Size() ; ++i )
6625 pt_type(VerifyClosed(type)),
6631 #ifndef MFEM_THREAD_SAFE
6640 for (
int i = 1; i < p; i++)
6650 const int p =
Order;
6652 #ifdef MFEM_THREAD_SAFE
6656 basis1d.
Eval(ip.
x, shape_x);
6658 shape(0) = shape_x(0);
6659 shape(1) = shape_x(p);
6660 for (
int i = 1; i < p; i++)
6662 shape(i+1) = shape_x(i);
6669 const int p =
Order;
6671 #ifdef MFEM_THREAD_SAFE
6672 Vector shape_x(p+1), dshape_x(p+1);
6675 basis1d.
Eval(ip.
x, shape_x, dshape_x);
6677 dshape(0,0) = dshape_x(0);
6678 dshape(1,0) = dshape_x(p);
6679 for (
int i = 1; i < p; i++)
6681 dshape(i+1,0) = dshape_x(i);
6687 const int p =
Order;
6695 for (
int i = 1; i < p; i++)
6704 for (
int i = 1; i < p; i++)
6716 pt_type(VerifyClosed(type)),
6717 basis1d(
poly1d.ClosedBasis(p, pt_type)),
6718 dof_map((p + 1)*(p + 1))
6722 const int p1 = p + 1;
6724 #ifndef MFEM_THREAD_SAFE
6732 dof_map[0 + 0*p1] = 0;
6733 dof_map[p + 0*p1] = 1;
6734 dof_map[p + p*p1] = 2;
6735 dof_map[0 + p*p1] = 3;
6739 for (
int i = 1; i < p; i++)
6741 dof_map[i + 0*p1] = o++;
6743 for (
int i = 1; i < p; i++)
6745 dof_map[p + i*p1] = o++;
6747 for (
int i = 1; i < p; i++)
6749 dof_map[(p-i) + p*p1] = o++;
6751 for (
int i = 1; i < p; i++)
6753 dof_map[0 + (p-i)*p1] = o++;
6757 for (
int j = 1; j < p; j++)
6758 for (
int i = 1; i < p; i++)
6760 dof_map[i + j*p1] = o++;
6764 for (
int j = 0; j <= p; j++)
6766 for (
int i = 0; i <= p; i++)
6776 const int p =
Order;
6778 #ifdef MFEM_THREAD_SAFE
6779 Vector shape_x(p+1), shape_y(p+1);
6782 basis1d.
Eval(ip.
x, shape_x);
6783 basis1d.
Eval(ip.
y, shape_y);
6785 for (
int o = 0, j = 0; j <= p; j++)
6786 for (
int i = 0; i <= p; i++)
6788 shape(dof_map[o++]) = shape_x(i)*shape_y(j);
6795 const int p =
Order;
6797 #ifdef MFEM_THREAD_SAFE
6798 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
6801 basis1d.
Eval(ip.
x, shape_x, dshape_x);
6802 basis1d.
Eval(ip.
y, shape_y, dshape_y);
6804 for (
int o = 0, j = 0; j <= p; j++)
6806 for (
int i = 0; i <= p; i++)
6808 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j);
6809 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j); o++;
6816 const int p =
Order;
6819 #ifdef MFEM_THREAD_SAFE
6820 Vector shape_x(p+1), shape_y(p+1);
6823 for (
int i = 0; i <= p; i++)
6832 for (
int o = 0, j = 0; j <= p; j++)
6833 for (
int i = 0; i <= p; i++)
6835 dofs(dof_map[o++]) = shape_x(i)*shape_x(j);
6839 for (
int o = 0, j = 0; j <= p; j++)
6840 for (
int i = 0; i <= p; i++)
6842 dofs(dof_map[o++]) = shape_y(i)*shape_x(j);
6846 for (
int o = 0, j = 0; j <= p; j++)
6847 for (
int i = 0; i <= p; i++)
6849 dofs(dof_map[o++]) = shape_y(i)*shape_y(j);
6853 for (
int o = 0, j = 0; j <= p; j++)
6854 for (
int i = 0; i <= p; i++)
6856 dofs(dof_map[o++]) = shape_x(i)*shape_y(j);
6866 pt_type(VerifyClosed(type)),
6867 basis1d(
poly1d.ClosedBasis(p, pt_type)),
6868 dof_map((p + 1)*(p + 1)*(p + 1))
6872 const int p1 = p + 1;
6874 #ifndef MFEM_THREAD_SAFE
6884 dof_map[0 + (0 + 0*p1)*p1] = 0;
6885 dof_map[p + (0 + 0*p1)*p1] = 1;
6886 dof_map[p + (p + 0*p1)*p1] = 2;
6887 dof_map[0 + (p + 0*p1)*p1] = 3;
6888 dof_map[0 + (0 + p*p1)*p1] = 4;
6889 dof_map[p + (0 + p*p1)*p1] = 5;
6890 dof_map[p + (p + p*p1)*p1] = 6;
6891 dof_map[0 + (p + p*p1)*p1] = 7;
6895 for (
int i = 1; i < p; i++)
6897 dof_map[i + (0 + 0*p1)*p1] = o++;
6899 for (
int i = 1; i < p; i++)
6901 dof_map[p + (i + 0*p1)*p1] = o++;
6903 for (
int i = 1; i < p; i++)
6905 dof_map[i + (p + 0*p1)*p1] = o++;
6907 for (
int i = 1; i < p; i++)
6909 dof_map[0 + (i + 0*p1)*p1] = o++;
6911 for (
int i = 1; i < p; i++)
6913 dof_map[i + (0 + p*p1)*p1] = o++;
6915 for (
int i = 1; i < p; i++)
6917 dof_map[p + (i + p*p1)*p1] = o++;
6919 for (
int i = 1; i < p; i++)
6921 dof_map[i + (p + p*p1)*p1] = o++;
6923 for (
int i = 1; i < p; i++)
6925 dof_map[0 + (i + p*p1)*p1] = o++;
6927 for (
int i = 1; i < p; i++)
6929 dof_map[0 + (0 + i*p1)*p1] = o++;
6931 for (
int i = 1; i < p; i++)
6933 dof_map[p + (0 + i*p1)*p1] = o++;
6935 for (
int i = 1; i < p; i++)
6937 dof_map[p + (p + i*p1)*p1] = o++;
6939 for (
int i = 1; i < p; i++)
6941 dof_map[0 + (p + i*p1)*p1] = o++;
6945 for (
int j = 1; j < p; j++)
6946 for (
int i = 1; i < p; i++)
6948 dof_map[i + ((p-j) + 0*p1)*p1] = o++;
6950 for (
int j = 1; j < p; j++)
6951 for (
int i = 1; i < p; i++)
6953 dof_map[i + (0 + j*p1)*p1] = o++;
6955 for (
int j = 1; j < p; j++)
6956 for (
int i = 1; i < p; i++)
6958 dof_map[p + (i + j*p1)*p1] = o++;
6960 for (
int j = 1; j < p; j++)
6961 for (
int i = 1; i < p; i++)
6963 dof_map[(p-i) + (p + j*p1)*p1] = o++;
6965 for (
int j = 1; j < p; j++)
6966 for (
int i = 1; i < p; i++)
6968 dof_map[0 + ((p-i) + j*p1)*p1] = o++;
6970 for (
int j = 1; j < p; j++)
6971 for (
int i = 1; i < p; i++)
6973 dof_map[i + (j + p*p1)*p1] = o++;
6977 for (
int k = 1; k < p; k++)
6978 for (
int j = 1; j < p; j++)
6979 for (
int i = 1; i < p; i++)
6981 dof_map[i + (j + k*p1)*p1] = o++;
6985 for (
int k = 0; k <= p; k++)
6986 for (
int j = 0; j <= p; j++)
6987 for (
int i = 0; i <= p; i++)
6996 const int p =
Order;
6998 #ifdef MFEM_THREAD_SAFE
6999 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7002 basis1d.
Eval(ip.
x, shape_x);
7003 basis1d.
Eval(ip.
y, shape_y);
7004 basis1d.
Eval(ip.
z, shape_z);
7006 for (
int o = 0, k = 0; k <= p; k++)
7007 for (
int j = 0; j <= p; j++)
7008 for (
int i = 0; i <= p; i++)
7010 shape(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_z(k);
7017 const int p =
Order;
7019 #ifdef MFEM_THREAD_SAFE
7020 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7021 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
7024 basis1d.
Eval(ip.
x, shape_x, dshape_x);
7025 basis1d.
Eval(ip.
y, shape_y, dshape_y);
7026 basis1d.
Eval(ip.
z, shape_z, dshape_z);
7028 for (
int o = 0, k = 0; k <= p; k++)
7029 for (
int j = 0; j <= p; j++)
7030 for (
int i = 0; i <= p; i++)
7032 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j)* shape_z(k);
7033 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j)* shape_z(k);
7034 dshape(dof_map[o],2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
7040 const int p =
Order;
7043 #ifdef MFEM_THREAD_SAFE
7044 Vector shape_x(p+1), shape_y(p+1);
7047 for (
int i = 0; i <= p; i++)
7056 for (
int o = 0, k = 0; k <= p; k++)
7057 for (
int j = 0; j <= p; j++)
7058 for (
int i = 0; i <= p; i++)
7060 dofs(dof_map[o++]) = shape_x(i)*shape_x(j)*shape_x(k);
7064 for (
int o = 0, k = 0; k <= p; k++)
7065 for (
int j = 0; j <= p; j++)
7066 for (
int i = 0; i <= p; i++)
7068 dofs(dof_map[o++]) = shape_y(i)*shape_x(j)*shape_x(k);
7072 for (
int o = 0, k = 0; k <= p; k++)
7073 for (
int j = 0; j <= p; j++)
7074 for (
int i = 0; i <= p; i++)
7076 dofs(dof_map[o++]) = shape_y(i)*shape_y(j)*shape_x(k);
7080 for (
int o = 0, k = 0; k <= p; k++)
7081 for (
int j = 0; j <= p; j++)
7082 for (
int i = 0; i <= p; i++)
7084 dofs(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_x(k);
7088 for (
int o = 0, k = 0; k <= p; k++)
7089 for (
int j = 0; j <= p; j++)
7090 for (
int i = 0; i <= p; i++)
7092 dofs(dof_map[o++]) = shape_x(i)*shape_x(j)*shape_y(k);
7096 for (
int o = 0, k = 0; k <= p; k++)
7097 for (
int j = 0; j <= p; j++)
7098 for (
int i = 0; i <= p; i++)
7100 dofs(dof_map[o++]) = shape_y(i)*shape_x(j)*shape_y(k);
7104 for (
int o = 0, k = 0; k <= p; k++)
7105 for (
int j = 0; j <= p; j++)
7106 for (
int i = 0; i <= p; i++)
7108 dofs(dof_map[o++]) = shape_y(i)*shape_y(j)*shape_y(k);
7112 for (
int o = 0, k = 0; k <= p; k++)
7113 for (
int j = 0; j <= p; j++)
7114 for (
int i = 0; i <= p; i++)
7116 dofs(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_y(k);
7127 #ifndef MFEM_THREAD_SAFE
7138 for (
int i = 1; i < p; i++)
7148 const int p =
Order;
7150 #ifdef MFEM_THREAD_SAFE
7157 shape(0) = shape_x(0);
7158 shape(1) = shape_x(p);
7159 for (
int i = 1; i < p; i++)
7161 shape(i+1) = shape_x(i);
7168 const int p =
Order;
7170 #ifdef MFEM_THREAD_SAFE
7171 Vector shape_x(p+1), dshape_x(p+1);
7177 dshape(0,0) = dshape_x(0);
7178 dshape(1,0) = dshape_x(p);
7179 for (
int i = 1; i < p; i++)
7181 dshape(i+1,0) = dshape_x(i);
7195 dof_map((p + 1)*(p + 1))
7197 const int p1 = p + 1;
7199 #ifndef MFEM_THREAD_SAFE
7208 dof_map[0 + 0*p1] = 0;
7209 dof_map[p + 0*p1] = 1;
7210 dof_map[p + p*p1] = 2;
7211 dof_map[0 + p*p1] = 3;
7215 for (
int i = 1; i < p; i++)
7217 dof_map[i + 0*p1] = o++;
7219 for (
int i = 1; i < p; i++)
7221 dof_map[p + i*p1] = o++;
7223 for (
int i = 1; i < p; i++)
7225 dof_map[(p-i) + p*p1] = o++;
7227 for (
int i = 1; i < p; i++)
7229 dof_map[0 + (p-i)*p1] = o++;
7233 for (
int j = 1; j < p; j++)
7234 for (
int i = 1; i < p; i++)
7236 dof_map[i + j*p1] = o++;
7240 for (
int j = 0; j <= p; j++)
7241 for (
int i = 0; i <= p; i++)
7250 const int p =
Order;
7252 #ifdef MFEM_THREAD_SAFE
7253 Vector shape_x(p+1), shape_y(p+1);
7260 for (
int o = 0, j = 0; j <= p; j++)
7261 for (
int i = 0; i <= p; i++)
7263 shape(dof_map[o++]) = shape_x(i)*shape_y(j);
7270 const int p =
Order;
7272 #ifdef MFEM_THREAD_SAFE
7273 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
7280 for (
int o = 0, j = 0; j <= p; j++)
7281 for (
int i = 0; i <= p; i++)
7283 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j);
7284 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j); o++;
7298 dof_map((p + 1)*(p + 1)*(p + 1))
7300 const int p1 = p + 1;
7302 #ifndef MFEM_THREAD_SAFE
7313 dof_map[0 + (0 + 0*p1)*p1] = 0;
7314 dof_map[p + (0 + 0*p1)*p1] = 1;
7315 dof_map[p + (p + 0*p1)*p1] = 2;
7316 dof_map[0 + (p + 0*p1)*p1] = 3;
7317 dof_map[0 + (0 + p*p1)*p1] = 4;
7318 dof_map[p + (0 + p*p1)*p1] = 5;
7319 dof_map[p + (p + p*p1)*p1] = 6;
7320 dof_map[0 + (p + p*p1)*p1] = 7;
7324 for (
int i = 1; i < p; i++)
7326 dof_map[i + (0 + 0*p1)*p1] = o++;
7328 for (
int i = 1; i < p; i++)
7330 dof_map[p + (i + 0*p1)*p1] = o++;
7332 for (
int i = 1; i < p; i++)
7334 dof_map[i + (p + 0*p1)*p1] = o++;
7336 for (
int i = 1; i < p; i++)
7338 dof_map[0 + (i + 0*p1)*p1] = o++;
7340 for (
int i = 1; i < p; i++)
7342 dof_map[i + (0 + p*p1)*p1] = o++;
7344 for (
int i = 1; i < p; i++)
7346 dof_map[p + (i + p*p1)*p1] = o++;
7348 for (
int i = 1; i < p; i++)
7350 dof_map[i + (p + p*p1)*p1] = o++;
7352 for (
int i = 1; i < p; i++)
7354 dof_map[0 + (i + p*p1)*p1] = o++;
7356 for (
int i = 1; i < p; i++)
7358 dof_map[0 + (0 + i*p1)*p1] = o++;
7360 for (
int i = 1; i < p; i++)
7362 dof_map[p + (0 + i*p1)*p1] = o++;
7364 for (
int i = 1; i < p; i++)
7366 dof_map[p + (p + i*p1)*p1] = o++;
7368 for (
int i = 1; i < p; i++)
7370 dof_map[0 + (p + i*p1)*p1] = o++;
7374 for (
int j = 1; j < p; j++)
7375 for (
int i = 1; i < p; i++)
7377 dof_map[i + ((p-j) + 0*p1)*p1] = o++;
7379 for (
int j = 1; j < p; j++)
7380 for (
int i = 1; i < p; i++)
7382 dof_map[i + (0 + j*p1)*p1] = o++;
7384 for (
int j = 1; j < p; j++)
7385 for (
int i = 1; i < p; i++)
7387 dof_map[p + (i + j*p1)*p1] = o++;
7389 for (
int j = 1; j < p; j++)
7390 for (
int i = 1; i < p; i++)
7392 dof_map[(p-i) + (p + j*p1)*p1] = o++;
7394 for (
int j = 1; j < p; j++)
7395 for (
int i = 1; i < p; i++)
7397 dof_map[0 + ((p-i) + j*p1)*p1] = o++;
7399 for (
int j = 1; j < p; j++)
7400 for (
int i = 1; i < p; i++)
7402 dof_map[i + (j + p*p1)*p1] = o++;
7406 for (
int k = 1; k < p; k++)
7407 for (
int j = 1; j < p; j++)
7408 for (
int i = 1; i < p; i++)
7410 dof_map[i + (j + k*p1)*p1] = o++;
7414 for (
int k = 0; k <= p; k++)
7415 for (
int j = 0; j <= p; j++)
7416 for (
int i = 0; i <= p; i++)
7424 const int p =
Order;
7426 #ifdef MFEM_THREAD_SAFE
7427 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7434 for (
int o = 0, k = 0; k <= p; k++)
7435 for (
int j = 0; j <= p; j++)
7436 for (
int i = 0; i <= p; i++)
7438 shape(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_z(k);
7445 const int p =
Order;
7447 #ifdef MFEM_THREAD_SAFE
7448 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7449 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
7456 for (
int o = 0, k = 0; k <= p; k++)
7457 for (
int j = 0; j <= p; j++)
7458 for (
int i = 0; i <= p; i++)
7460 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j)* shape_z(k);
7461 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j)* shape_z(k);
7462 dshape(dof_map[o],2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
7479 #ifndef MFEM_THREAD_SAFE
7489 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
7499 for (
int i = 1; i < p; i++)
7503 for (
int i = 1; i < p; i++)
7507 for (
int i = 1; i < p; i++)
7513 for (
int j = 1; j < p; j++)
7514 for (
int i = 1; i + j < p; i++)
7516 const double w = cp[i] + cp[j] + cp[p-i-j];
7521 for (
int k = 0; k <
Dof; k++)
7529 for (
int j = 0; j <= p; j++)
7530 for (
int i = 0; i + j <= p; i++)
7532 T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
7543 const int p =
Order;
7545 #ifdef MFEM_THREAD_SAFE
7546 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(
Dof);
7553 for (
int o = 0, j = 0; j <= p; j++)
7554 for (
int i = 0; i + j <= p; i++)
7556 u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
7565 const int p =
Order;
7567 #ifdef MFEM_THREAD_SAFE
7568 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
7569 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
7577 for (
int o = 0, j = 0; j <= p; j++)
7578 for (
int i = 0; i + j <= p; i++)
7581 du(o,0) = ((dshape_x(i)* shape_l(k)) -
7582 ( shape_x(i)*dshape_l(k)))*shape_y(j);
7583 du(o,1) = ((dshape_y(j)* shape_l(k)) -
7584 ( shape_y(j)*dshape_l(k)))*shape_x(i);
7588 Ti.
Mult(du, dshape);
7598 #ifndef MFEM_THREAD_SAFE
7610 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7621 for (
int i = 1; i < p; i++)
7625 for (
int i = 1; i < p; i++)
7629 for (
int i = 1; i < p; i++)
7633 for (
int i = 1; i < p; i++)
7637 for (
int i = 1; i < p; i++)
7641 for (
int i = 1; i < p; i++)
7647 for (
int j = 1; j < p; j++)
7648 for (
int i = 1; i + j < p; i++)
7650 double w = cp[i] + cp[j] + cp[p-i-j];
7653 for (
int j = 1; j < p; j++)
7654 for (
int i = 1; i + j < p; i++)
7656 double w = cp[i] + cp[j] + cp[p-i-j];
7659 for (
int j = 1; j < p; j++)
7660 for (
int i = 1; i + j < p; i++)
7662 double w = cp[i] + cp[j] + cp[p-i-j];
7665 for (
int j = 1; j < p; j++)
7666 for (
int i = 1; i + j < p; i++)
7668 double w = cp[i] + cp[j] + cp[p-i-j];
7673 for (
int k = 1; k < p; k++)
7674 for (
int j = 1; j + k < p; j++)
7675 for (
int i = 1; i + j + k < p; i++)
7677 double w = cp[i] + cp[j] + cp[k] + cp[p-i-j-k];
7682 for (
int m = 0; m <
Dof; m++)
7691 for (
int k = 0; k <= p; k++)
7692 for (
int j = 0; j + k <= p; j++)
7693 for (
int i = 0; i + j + k <= p; i++)
7695 T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
7706 const int p =
Order;
7708 #ifdef MFEM_THREAD_SAFE
7709 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7718 for (
int o = 0, k = 0; k <= p; k++)
7719 for (
int j = 0; j + k <= p; j++)
7720 for (
int i = 0; i + j + k <= p; i++)
7722 u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
7731 const int p =
Order;
7733 #ifdef MFEM_THREAD_SAFE
7734 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7735 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
7744 for (
int o = 0, k = 0; k <= p; k++)
7745 for (
int j = 0; j + k <= p; j++)
7746 for (
int i = 0; i + j + k <= p; i++)
7748 int l = p - i - j - k;
7749 du(o,0) = ((dshape_x(i)* shape_l(l)) -
7750 ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
7751 du(o,1) = ((dshape_y(j)* shape_l(l)) -
7752 ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
7753 du(o,2) = ((dshape_z(k)* shape_l(l)) -
7754 ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
7758 Ti.
Mult(du, dshape);
7766 #ifndef MFEM_THREAD_SAFE
7776 Index(
int p) { p2p3 = 2*p + 3; }
7777 int operator()(
int i,
int j) {
return ((p2p3-j)*j)/2+i; }
7791 for (
int i = 1; i < p; i++)
7796 for (
int i = 1; i < p; i++)
7801 for (
int i = 1; i < p; i++)
7808 for (
int j = 1; j < p; j++)
7809 for (
int i = 1; i + j < p; i++)
7818 const int p,
const double l1,
const double l2,
double *shape)
7820 const double l3 = 1. - l1 - l2;
7830 for (
int o = 0, j = 0; j <= p; j++)
7834 for (
int i = 0; i <= p - j; i++)
7844 const int p,
const double l1,
const double l2,
7845 double *dshape_1d,
double *dshape)
7847 const int dof = ((p + 1)*(p + 2))/2;
7848 const double l3 = 1. - l1 - l2;
7852 for (
int o = 0, j = 0; j <= p; j++)
7856 for (
int i = 0; i <= p - j; i++)
7858 dshape[o++] = s*dshape_1d[i];
7863 for (
int i = 0; i <= p; i++)
7867 for (
int o = i, j = 0; j <= p - i; j++)
7869 dshape[dof + o] = s*dshape_1d[j];
7879 #ifdef MFEM_THREAD_SAFE
7883 for (
int i = 0; i <
Dof; i++)
7892 #ifdef MFEM_THREAD_SAFE
7897 for (
int d = 0; d < 2; d++)
7899 for (
int i = 0; i <
Dof; i++)
7911 #ifndef MFEM_THREAD_SAFE
7921 int tri(
int k) {
return (k*(k + 1))/2; }
7922 int tet(
int k) {
return (k*(k + 1)*(k + 2))/6; }
7923 Index(
int p_) { p = p_; dof = tet(p + 1); }
7924 int operator()(
int i,
int j,
int k)
7925 {
return dof - tet(p - k) - tri(p + 1 - k - j) + i; }
7941 for (
int i = 1; i < p; i++)
7946 for (
int i = 1; i < p; i++)
7951 for (
int i = 1; i < p; i++)
7956 for (
int i = 1; i < p; i++)
7961 for (
int i = 1; i < p; i++)
7966 for (
int i = 1; i < p; i++)
7973 for (
int j = 1; j < p; j++)
7974 for (
int i = 1; i + j < p; i++)
7979 for (
int j = 1; j < p; j++)
7980 for (
int i = 1; i + j < p; i++)
7985 for (
int j = 1; j < p; j++)
7986 for (
int i = 1; i + j < p; i++)
7991 for (
int j = 1; j < p; j++)
7992 for (
int i = 1; i + j < p; i++)
7999 for (
int k = 1; k < p; k++)
8000 for (
int j = 1; j + k < p; j++)
8001 for (
int i = 1; i + j + k < p; i++)
8010 const int p,
const double l1,
const double l2,
const double l3,
8013 const double l4 = 1. - l1 - l2 - l3;
8022 for (
int o = 0, k = 0; k <= p; k++)
8025 const double ek = bp[k]*l3k;
8027 for (
int j = 0; j <= p - k; j++)
8030 double ekj = ek*bpk[j]*l2j;
8031 for (
int i = 0; i <= p - k - j; i++)
8043 const int p,
const double l1,
const double l2,
const double l3,
8044 double *dshape_1d,
double *dshape)
8046 const int dof = ((p + 1)*(p + 2)*(p + 3))/6;
8047 const double l4 = 1. - l1 - l2 - l3;
8055 for (
int o = 0, k = 0; k <= p; k++)
8058 const double ek = bp[k]*l3k;
8060 for (
int j = 0; j <= p - k; j++)
8063 double ekj = ek*bpk[j]*l2j;
8064 for (
int i = 0; i <= p - k - j; i++)
8066 dshape[o++] = dshape_1d[i]*ekj;
8077 for (
int ok = 0, k = 0; k <= p; k++)
8080 const double ek = bp[k]*l3k;
8082 for (
int i = 0; i <= p - k; i++)
8085 double eki = ek*bpk[i]*l1i;
8087 for (
int j = 0; j <= p - k - i; j++)
8089 dshape[dof + o] = dshape_1d[j]*eki;
8095 ok += ((p - k + 2)*(p - k + 1))/2;
8102 for (
int j = 0; j <= p; j++)
8105 const double ej = bp[j]*l2j;
8107 for (
int i = 0; i <= p - j; i++)
8110 double eji = ej*bpj[i]*l1i;
8111 int m = ((p + 2)*(p + 1))/2;
8112 int n = ((p - j + 2)*(p - j + 1))/2;
8113 for (
int o = i, k = 0; k <= p - j - i; k++)
8118 dshape[2*dof + o - n] = dshape_1d[k]*eji;
8131 #ifdef MFEM_THREAD_SAFE
8135 for (
int i = 0; i <
Dof; i++)
8144 #ifdef MFEM_THREAD_SAFE
8149 for (
int d = 0; d < 3; d++)
8151 for (
int i = 0; i <
Dof; i++)
8161 type(VerifyOpen(type)),
8162 basis1d(
poly1d.OpenBasis(p, type))
8166 #ifndef MFEM_THREAD_SAFE
8171 for (
int i = 0; i <= p; i++)
8180 basis1d.
Eval(ip.
x, shape);
8186 #ifdef MFEM_THREAD_SAFE
8191 basis1d.
Eval(ip.
x, shape_x, dshape_x);
8196 const int p =
Order;
8202 for (
int i = 0; i <= p; i++)
8209 for (
int i = 0; i <= p; i++)
8221 #ifndef MFEM_THREAD_SAFE
8232 for (
int i = 0; i <= p; i++)
8248 #ifdef MFEM_THREAD_SAFE
8259 dofs[vertex*
Order] = 1.0;
8266 type(VerifyOpen(_type)),
8267 basis1d(
poly1d.OpenBasis(p, type))
8271 #ifndef MFEM_THREAD_SAFE
8278 for (
int o = 0, j = 0; j <= p; j++)
8279 for (
int i = 0; i <= p; i++)
8288 const int p =
Order;
8290 #ifdef MFEM_THREAD_SAFE
8291 Vector shape_x(p+1), shape_y(p+1);
8294 basis1d.
Eval(ip.
x, shape_x);
8295 basis1d.
Eval(ip.
y, shape_y);
8297 for (
int o = 0, j = 0; j <= p; j++)
8298 for (
int i = 0; i <= p; i++)
8300 shape(o++) = shape_x(i)*shape_y(j);
8307 const int p =
Order;
8309 #ifdef MFEM_THREAD_SAFE
8310 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
8313 basis1d.
Eval(ip.
x, shape_x, dshape_x);
8314 basis1d.
Eval(ip.
y, shape_y, dshape_y);
8316 for (
int o = 0, j = 0; j <= p; j++)
8317 for (
int i = 0; i <= p; i++)
8319 dshape(o,0) = dshape_x(i)* shape_y(j);
8320 dshape(o,1) = shape_x(i)*dshape_y(j); o++;
8326 const int p =
Order;
8329 #ifdef MFEM_THREAD_SAFE
8330 Vector shape_x(p+1), shape_y(p+1);
8333 for (
int i = 0; i <= p; i++)
8342 for (
int o = 0, j = 0; j <= p; j++)
8343 for (
int i = 0; i <= p; i++)
8345 dofs[o++] = shape_x(i)*shape_x(j);
8349 for (
int o = 0, j = 0; j <= p; j++)
8350 for (
int i = 0; i <= p; i++)
8352 dofs[o++] = shape_y(i)*shape_x(j);
8356 for (
int o = 0, j = 0; j <= p; j++)
8357 for (
int i = 0; i <= p; i++)
8359 dofs[o++] = shape_y(i)*shape_y(j);
8363 for (
int o = 0, j = 0; j <= p; j++)
8364 for (
int i = 0; i <= p; i++)
8366 dofs[o++] = shape_x(i)*shape_y(j);
8377 #ifndef MFEM_THREAD_SAFE
8390 for (
int o = 0, j = 0; j <= p; j++)
8391 for (
int i = 0; i <= p; i++)
8401 const int p =
Order;
8403 #ifdef MFEM_THREAD_SAFE
8404 Vector shape_x(p+1), shape_y(p+1);
8410 for (
int o = 0, j = 0; j <= p; j++)
8411 for (
int i = 0; i <= p; i++)
8413 shape(o++) = shape_x(i)*shape_y(j);
8420 const int p =
Order;
8422 #ifdef MFEM_THREAD_SAFE
8423 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
8429 for (
int o = 0, j = 0; j <= p; j++)
8430 for (
int i = 0; i <= p; i++)
8432 dshape(o,0) = dshape_x(i)* shape_y(j);
8433 dshape(o,1) = shape_x(i)*dshape_y(j); o++;
8439 const int p =
Order;
8444 case 0: dofs[0] = 1.0;
break;
8445 case 1: dofs[p] = 1.0;
break;
8446 case 2: dofs[p*(p + 2)] = 1.0;
break;
8447 case 3: dofs[p*(p + 1)] = 1.0;
break;
8455 type(VerifyOpen(_type)),
8456 basis1d(
poly1d.OpenBasis(p, type))
8460 #ifndef MFEM_THREAD_SAFE
8469 for (
int o = 0, k = 0; k <= p; k++)
8470 for (
int j = 0; j <= p; j++)
8471 for (
int i = 0; i <= p; i++)
8480 const int p =
Order;
8482 #ifdef MFEM_THREAD_SAFE
8483 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
8486 basis1d.
Eval(ip.
x, shape_x);
8487 basis1d.
Eval(ip.
y, shape_y);
8488 basis1d.
Eval(ip.
z, shape_z);
8490 for (
int o = 0, k = 0; k <= p; k++)
8491 for (
int j = 0; j <= p; j++)
8492 for (
int i = 0; i <= p; i++)
8494 shape(o++) = shape_x(i)*shape_y(j)*shape_z(k);
8501 const int p =
Order;
8503 #ifdef MFEM_THREAD_SAFE
8504 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
8505 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
8508 basis1d.
Eval(ip.
x, shape_x, dshape_x);
8509 basis1d.
Eval(ip.
y, shape_y, dshape_y);
8510 basis1d.
Eval(ip.
z, shape_z, dshape_z);
8512 for (
int o = 0, k = 0; k <= p; k++)
8513 for (
int j = 0; j <= p; j++)
8514 for (
int i = 0; i <= p; i++)
8516 dshape(o,0) = dshape_x(i)* shape_y(j)* shape_z(k);
8517 dshape(o,1) = shape_x(i)*dshape_y(j)* shape_z(k);
8518 dshape(o,2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
8524 const int p =
Order;
8527 #ifdef MFEM_THREAD_SAFE
8528 Vector shape_x(p+1), shape_y(p+1);
8531 for (
int i = 0; i <= p; i++)
8540 for (
int o = 0, k = 0; k <= p; k++)
8541 for (
int j = 0; j <= p; j++)
8542 for (
int i = 0; i <= p; i++)
8544 dofs[o++] = shape_x(i)*shape_x(j)*shape_x(k);
8548 for (
int o = 0, k = 0; k <= p; k++)
8549 for (
int j = 0; j <= p; j++)
8550 for (
int i = 0; i <= p; i++)
8552 dofs[o++] = shape_y(i)*shape_x(j)*shape_x(k);
8556 for (
int o = 0, k = 0; k <= p; k++)
8557 for (
int j = 0; j <= p; j++)
8558 for (
int i = 0; i <= p; i++)
8560 dofs[o++] = shape_y(i)*shape_y(j)*shape_x(k);
8564 for (
int o = 0, k = 0; k <= p; k++)
8565 for (
int j = 0; j <= p; j++)
8566 for (
int i = 0; i <= p; i++)
8568 dofs[o++] = shape_x(i)*shape_y(j)*shape_x(k);
8572 for (
int o = 0, k = 0; k <= p; k++)
8573 for (
int j = 0; j <= p; j++)
8574 for (
int i = 0; i <= p; i++)
8576 dofs[o++] = shape_x(i)*shape_x(j)*shape_y(k);
8580 for (
int o = 0, k = 0; k <= p; k++)
8581 for (
int j = 0; j <= p; j++)
8582 for (
int i = 0; i <= p; i++)
8584 dofs[o++] = shape_y(i)*shape_x(j)*shape_y(k);
8588 for (
int o = 0, k = 0; k <= p; k++)
8589 for (
int j = 0; j <= p; j++)
8590 for (
int i = 0; i <= p; i++)
8592 dofs[o++] = shape_y(i)*shape_y(j)*shape_y(k);
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++)
8600 dofs[o++] = shape_x(i)*shape_y(j)*shape_y(k);
8611 #ifndef MFEM_THREAD_SAFE
8626 for (
int o = 0, k = 0; k <= p; k++)
8627 for (
int j = 0; j <= p; j++)
8628 for (
int i = 0; i <= p; i++)
8638 const int p =
Order;
8640 #ifdef MFEM_THREAD_SAFE
8641 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
8648 for (
int o = 0, k = 0; k <= p; k++)
8649 for (
int j = 0; j <= p; j++)
8650 for (
int i = 0; i <= p; i++)
8652 shape(o++) = shape_x(i)*shape_y(j)*shape_z(k);
8659 const int p =
Order;
8661 #ifdef MFEM_THREAD_SAFE
8662 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
8663 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
8670 for (
int o = 0, k = 0; k <= p; k++)
8671 for (
int j = 0; j <= p; j++)
8672 for (
int i = 0; i <= p; i++)
8674 dshape(o,0) = dshape_x(i)* shape_y(j)* shape_z(k);
8675 dshape(o,1) = shape_x(i)*dshape_y(j)* shape_z(k);
8676 dshape(o,2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
8682 const int p =
Order;
8687 case 0: dofs[0] = 1.0;
break;
8688 case 1: dofs[p] = 1.0;
break;
8689 case 2: dofs[p*(p + 2)] = 1.0;
break;
8690 case 3: dofs[p*(p + 1)] = 1.0;
break;
8691 case 4: dofs[p*(p + 1)*(p + 1)] = 1.0;
break;
8692 case 5: dofs[p + p*(p + 1)*(p + 1)] = 1.0;
break;
8693 case 6: dofs[
Dof - 1] = 1.0;
break;
8694 case 7: dofs[
Dof - p - 1] = 1.0;
break;
8705 #ifndef MFEM_THREAD_SAFE
8715 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
8718 for (
int o = 0, j = 0; j <= p; j++)
8719 for (
int i = 0; i + j <= p; i++)
8721 double w = op[i] + op[j] + op[p-i-j];
8726 for (
int k = 0; k <
Dof; k++)
8733 for (
int o = 0, j = 0; j <= p; j++)
8734 for (
int i = 0; i + j <= p; i++)
8736 T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
8747 const int p =
Order;
8749 #ifdef MFEM_THREAD_SAFE
8750 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(
Dof);
8757 for (
int o = 0, j = 0; j <= p; j++)
8758 for (
int i = 0; i + j <= p; i++)
8760 u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
8769 const int p =
Order;
8771 #ifdef MFEM_THREAD_SAFE
8772 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
8773 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
8781 for (
int o = 0, j = 0; j <= p; j++)
8782 for (
int i = 0; i + j <= p; i++)
8785 du(o,0) = ((dshape_x(i)* shape_l(k)) -
8786 ( shape_x(i)*dshape_l(k)))*shape_y(j);
8787 du(o,1) = ((dshape_y(j)* shape_l(k)) -
8788 ( shape_y(j)*dshape_l(k)))*shape_x(i);
8792 Ti.
Mult(du, dshape);
8800 for (
int i = 0; i <
Dof; i++)
8803 dofs[i] = pow(1.0 - ip.
x - ip.
y,
Order);
8807 for (
int i = 0; i <
Dof; i++)
8810 dofs[i] = pow(ip.
x,
Order);
8814 for (
int i = 0; i <
Dof; i++)
8817 dofs[i] = pow(ip.
y,
Order);
8828 #ifndef MFEM_THREAD_SAFE
8838 for (
int o = 0, j = 0; j <= p; j++)
8839 for (
int i = 0; i + j <= p; i++)
8855 #ifdef MFEM_THREAD_SAFE
8868 case 0: dofs[0] = 1.0;
break;
8869 case 1: dofs[
Order] = 1.0;
break;
8870 case 2: dofs[
Dof-1] = 1.0;
break;
8881 #ifndef MFEM_THREAD_SAFE
8893 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8896 for (
int o = 0, k = 0; k <= p; k++)
8897 for (
int j = 0; j + k <= p; j++)
8898 for (
int i = 0; i + j + k <= p; i++)
8900 double w = op[i] + op[j] + op[k] + op[p-i-j-k];
8905 for (
int m = 0; m <
Dof; m++)
8913 for (
int o = 0, k = 0; k <= p; k++)
8914 for (
int j = 0; j + k <= p; j++)
8915 for (
int i = 0; i + j + k <= p; i++)
8917 T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
8928 const int p =
Order;
8930 #ifdef MFEM_THREAD_SAFE
8931 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8940 for (
int o = 0, k = 0; k <= p; k++)
8941 for (
int j = 0; j + k <= p; j++)
8942 for (
int i = 0; i + j + k <= p; i++)
8944 u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
8953 const int p =
Order;
8955 #ifdef MFEM_THREAD_SAFE
8956 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8957 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
8966 for (
int o = 0, k = 0; k <= p; k++)
8967 for (
int j = 0; j + k <= p; j++)
8968 for (
int i = 0; i + j + k <= p; i++)
8970 int l = p - i - j - k;
8971 du(o,0) = ((dshape_x(i)* shape_l(l)) -
8972 ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
8973 du(o,1) = ((dshape_y(j)* shape_l(l)) -
8974 ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
8975 du(o,2) = ((dshape_z(k)* shape_l(l)) -
8976 ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
8980 Ti.
Mult(du, dshape);
8988 for (
int i = 0; i <
Dof; i++)
8991 dofs[i] = pow(1.0 - ip.
x - ip.
y - ip.
z,
Order);
8995 for (
int i = 0; i <
Dof; i++)
8998 dofs[i] = pow(ip.
x,
Order);
9002 for (
int i = 0; i <
Dof; i++)
9005 dofs[i] = pow(ip.
y,
Order);
9008 for (
int i = 0; i <
Dof; i++)
9011 dofs[i] = pow(ip.
z,
Order);
9022 #ifndef MFEM_THREAD_SAFE
9032 for (
int o = 0, k = 0; k <= p; k++)
9033 for (
int j = 0; j + k <= p; j++)
9034 for (
int i = 0; i + j + k <= p; i++)
9051 #ifdef MFEM_THREAD_SAFE
9064 case 0: dofs[0] = 1.0;
break;
9065 case 1: dofs[
Order] = 1.0;
break;
9066 case 2: dofs[(
Order*(
Order+3))/2] = 1.0;
break;
9067 case 3: dofs[
Dof-1] = 1.0;
break;
9072 const double RT_QuadrilateralElement::nk[8] =
9073 { 0., -1., 1., 0., 0., 1., -1., 0. };
9080 cbasis1d(
poly1d.ClosedBasis(p + 1, VerifyClosed(cp_type))),
9081 obasis1d(
poly1d.OpenBasis(p, VerifyOpen(op_type))),
9082 dof_map(Dof), dof2nk(Dof)
9086 const int dof2 =
Dof/2;
9088 #ifndef MFEM_THREAD_SAFE
9099 for (
int i = 0; i <= p; i++)
9101 dof_map[1*dof2 + i + 0*(p + 1)] = o++;
9103 for (
int i = 0; i <= p; i++)
9105 dof_map[0*dof2 + (p + 1) + i*(p + 2)] = o++;
9107 for (
int i = 0; i <= p; i++)
9109 dof_map[1*dof2 + (p - i) + (p + 1)*(p + 1)] = o++;
9111 for (
int i = 0; i <= p; i++)
9113 dof_map[0*dof2 + 0 + (p - i)*(p + 2)] = o++;
9117 for (
int j = 0; j <= p; j++)
9118 for (
int i = 1; i <= p; i++)
9120 dof_map[0*dof2 + i + j*(p + 2)] = o++;
9122 for (
int j = 1; j <= p; j++)
9123 for (
int i = 0; i <= p; i++)
9125 dof_map[1*dof2 + i + j*(p + 1)] = o++;
9130 for (
int j = 0; j <= p; j++)
9131 for (
int i = 0; i <= p/2; i++)
9133 int idx = 0*dof2 + i + j*(p + 2);
9134 dof_map[idx] = -1 - dof_map[idx];
9137 for (
int j = p/2 + 1; j <= p; j++)
9139 int idx = 0*dof2 + (p/2 + 1) + j*(p + 2);
9140 dof_map[idx] = -1 - dof_map[idx];
9143 for (
int j = 0; j <= p/2; j++)
9144 for (
int i = 0; i <= p; i++)
9146 int idx = 1*dof2 + i + j*(p + 1);
9147 dof_map[idx] = -1 - dof_map[idx];
9150 for (
int i = 0; i <= p/2; i++)
9152 int idx = 1*dof2 + i + (p/2 + 1)*(p + 1);
9153 dof_map[idx] = -1 - dof_map[idx];
9157 for (
int j = 0; j <= p; j++)
9158 for (
int i = 0; i <= p + 1; i++)
9161 if ((idx = dof_map[o++]) < 0)
9172 for (
int j = 0; j <= p + 1; j++)
9173 for (
int i = 0; i <= p; i++)
9176 if ((idx = dof_map[o++]) < 0)
9192 const int pp1 =
Order;
9194 #ifdef MFEM_THREAD_SAFE
9195 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9198 cbasis1d.
Eval(ip.
x, shape_cx);
9199 obasis1d.
Eval(ip.
x, shape_ox);
9200 cbasis1d.
Eval(ip.
y, shape_cy);
9201 obasis1d.
Eval(ip.
y, shape_oy);
9204 for (
int j = 0; j < pp1; j++)
9205 for (
int i = 0; i <= pp1; i++)
9208 if ((idx = dof_map[o++]) < 0)
9210 idx = -1 - idx, s = -1;
9216 shape(idx,0) = s*shape_cx(i)*shape_oy(j);
9219 for (
int j = 0; j <= pp1; j++)
9220 for (
int i = 0; i < pp1; i++)
9223 if ((idx = dof_map[o++]) < 0)
9225 idx = -1 - idx, s = -1;
9232 shape(idx,1) = s*shape_ox(i)*shape_cy(j);
9239 const int pp1 =
Order;
9241 #ifdef MFEM_THREAD_SAFE
9242 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9243 Vector dshape_cx(pp1 + 1), dshape_cy(pp1 + 1);
9246 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
9247 obasis1d.
Eval(ip.
x, shape_ox);
9248 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
9249 obasis1d.
Eval(ip.
y, shape_oy);
9252 for (
int j = 0; j < pp1; j++)
9253 for (
int i = 0; i <= pp1; i++)
9256 if ((idx = dof_map[o++]) < 0)
9258 idx = -1 - idx, s = -1;
9264 divshape(idx) = s*dshape_cx(i)*shape_oy(j);
9266 for (
int j = 0; j <= pp1; j++)
9267 for (
int i = 0; i < pp1; i++)
9270 if ((idx = dof_map[o++]) < 0)
9272 idx = -1 - idx, s = -1;
9278 divshape(idx) = s*shape_ox(i)*dshape_cy(j);
9283 const double RT_HexahedronElement::nk[18] =
9284 { 0.,0.,-1., 0.,-1.,0., 1.,0.,0., 0.,1.,0., -1.,0.,0., 0.,0.,1. };
9291 cbasis1d(
poly1d.ClosedBasis(p + 1, VerifyClosed(cp_type))),
9292 obasis1d(
poly1d.OpenBasis(p, VerifyOpen(op_type))),
9293 dof_map(Dof), dof2nk(Dof)
9297 const int dof3 =
Dof/3;
9299 #ifndef MFEM_THREAD_SAFE
9313 for (
int j = 0; j <= p; j++)
9314 for (
int i = 0; i <= p; i++)
9316 dof_map[2*dof3 + i + ((p - j) + 0*(p + 1))*(p + 1)] = o++;
9318 for (
int j = 0; j <= p; j++)
9319 for (
int i = 0; i <= p; i++)
9321 dof_map[1*dof3 + i + (0 + j*(p + 2))*(p + 1)] = o++;
9323 for (
int j = 0; j <= p; j++)
9324 for (
int i = 0; i <= p; i++)
9326 dof_map[0*dof3 + (p + 1) + (i + j*(p + 1))*(p + 2)] = o++;
9328 for (
int j = 0; j <= p; j++)
9329 for (
int i = 0; i <= p; i++)
9331 dof_map[1*dof3 + (p - i) + ((p + 1) + j*(p + 2))*(p + 1)] = o++;
9333 for (
int j = 0; j <= p; j++)
9334 for (
int i = 0; i <= p; i++)
9336 dof_map[0*dof3 + 0 + ((p - i) + j*(p + 1))*(p + 2)] = o++;
9338 for (
int j = 0; j <= p; j++)
9339 for (
int i = 0; i <= p; i++)
9341 dof_map[2*dof3 + i + (j + (p + 1)*(p + 1))*(p + 1)] = o++;
9346 for (
int k = 0; k <= p; k++)
9347 for (
int j = 0; j <= p; j++)
9348 for (
int i = 1; i <= p; i++)
9350 dof_map[0*dof3 + i + (j + k*(p + 1))*(p + 2)] = o++;
9353 for (
int k = 0; k <= p; k++)
9354 for (
int j = 1; j <= p; j++)
9355 for (
int i = 0; i <= p; i++)
9357 dof_map[1*dof3 + i + (j + k*(p + 2))*(p + 1)] = o++;
9360 for (
int k = 1; k <= p; k++)
9361 for (
int j = 0; j <= p; j++)
9362 for (
int i = 0; i <= p; i++)
9364 dof_map[2*dof3 + i + (j + k*(p + 1))*(p + 1)] = o++;
9372 for (
int k = 0; k <= p; k++)
9373 for (
int j = 0; j <= p; j++)
9374 for (
int i = 0; i <= p/2; i++)
9376 int idx = 0*dof3 + i + (j + k*(p + 1))*(p + 2);
9377 dof_map[idx] = -1 - dof_map[idx];
9380 for (
int k = 0; k <= p; k++)
9381 for (
int j = 0; j <= p/2; j++)
9382 for (
int i = 0; i <= p; i++)
9384 int idx = 1*dof3 + i + (j + k*(p + 2))*(p + 1);
9385 dof_map[idx] = -1 - dof_map[idx];
9388 for (
int k = 0; k <= p/2; k++)
9389 for (
int j = 0; j <= p; j++)
9390 for (
int i = 0; i <= p; i++)
9392 int idx = 2*dof3 + i + (j + k*(p + 1))*(p + 1);
9393 dof_map[idx] = -1 - dof_map[idx];
9398 for (
int k = 0; k <= p; k++)
9399 for (
int j = 0; j <= p; j++)
9400 for (
int i = 0; i <= p + 1; i++)
9403 if ((idx = dof_map[o++]) < 0)
9415 for (
int k = 0; k <= p; k++)
9416 for (
int j = 0; j <= p + 1; j++)
9417 for (
int i = 0; i <= p; i++)
9420 if ((idx = dof_map[o++]) < 0)
9432 for (
int k = 0; k <= p + 1; k++)
9433 for (
int j = 0; j <= p; j++)
9434 for (
int i = 0; i <= p; i++)
9437 if ((idx = dof_map[o++]) < 0)
9453 const int pp1 =
Order;
9455 #ifdef MFEM_THREAD_SAFE
9456 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9457 Vector shape_cz(pp1 + 1), shape_oz(pp1);
9460 cbasis1d.
Eval(ip.
x, shape_cx);
9461 obasis1d.
Eval(ip.
x, shape_ox);
9462 cbasis1d.
Eval(ip.
y, shape_cy);
9463 obasis1d.
Eval(ip.
y, shape_oy);
9464 cbasis1d.
Eval(ip.
z, shape_cz);
9465 obasis1d.
Eval(ip.
z, shape_oz);
9469 for (
int k = 0; k < pp1; k++)
9470 for (
int j = 0; j < pp1; j++)
9471 for (
int i = 0; i <= pp1; i++)
9474 if ((idx = dof_map[o++]) < 0)
9476 idx = -1 - idx, s = -1;
9482 shape(idx,0) = s*shape_cx(i)*shape_oy(j)*shape_oz(k);
9487 for (
int k = 0; k < pp1; k++)
9488 for (
int j = 0; j <= pp1; j++)
9489 for (
int i = 0; i < pp1; i++)
9492 if ((idx = dof_map[o++]) < 0)
9494 idx = -1 - idx, s = -1;
9501 shape(idx,1) = s*shape_ox(i)*shape_cy(j)*shape_oz(k);
9505 for (
int k = 0; k <= pp1; k++)
9506 for (
int j = 0; j < pp1; j++)
9507 for (
int i = 0; i < pp1; i++)
9510 if ((idx = dof_map[o++]) < 0)
9512 idx = -1 - idx, s = -1;
9520 shape(idx,2) = s*shape_ox(i)*shape_oy(j)*shape_cz(k);
9527 const int pp1 =
Order;
9529 #ifdef MFEM_THREAD_SAFE
9530 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9531 Vector shape_cz(pp1 + 1), shape_oz(pp1);
9532 Vector dshape_cx(pp1 + 1), dshape_cy(pp1 + 1), dshape_cz(pp1 + 1);
9535 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
9536 obasis1d.
Eval(ip.
x, shape_ox);
9537 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
9538 obasis1d.
Eval(ip.
y, shape_oy);
9539 cbasis1d.
Eval(ip.
z, shape_cz, dshape_cz);
9540 obasis1d.
Eval(ip.
z, shape_oz);
9544 for (
int k = 0; k < pp1; k++)
9545 for (
int j = 0; j < pp1; j++)
9546 for (
int i = 0; i <= pp1; i++)
9549 if ((idx = dof_map[o++]) < 0)
9551 idx = -1 - idx, s = -1;
9557 divshape(idx) = s*dshape_cx(i)*shape_oy(j)*shape_oz(k);
9560 for (
int k = 0; k < pp1; k++)
9561 for (
int j = 0; j <= pp1; j++)
9562 for (
int i = 0; i < pp1; i++)
9565 if ((idx = dof_map[o++]) < 0)
9567 idx = -1 - idx, s = -1;
9573 divshape(idx) = s*shape_ox(i)*dshape_cy(j)*shape_oz(k);
9576 for (
int k = 0; k <= pp1; k++)
9577 for (
int j = 0; j < pp1; j++)
9578 for (
int i = 0; i < pp1; i++)
9581 if ((idx = dof_map[o++]) < 0)
9583 idx = -1 - idx, s = -1;
9589 divshape(idx) = s*shape_ox(i)*shape_oy(j)*dshape_cz(k);
9594 const double RT_TriangleElement::nk[6] =
9595 { 0., -1., 1., 1., -1., 0. };
9597 const double RT_TriangleElement::c = 1./3.;
9607 #ifndef MFEM_THREAD_SAFE
9617 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
9622 for (
int i = 0; i <= p; i++)
9627 for (
int i = 0; i <= p; i++)
9632 for (
int i = 0; i <= p; i++)
9639 for (
int j = 0; j < p; j++)
9640 for (
int i = 0; i + j < p; i++)
9642 double w = iop[i] + iop[j] + iop[p-1-i-j];
9650 for (
int k = 0; k <
Dof; k++)
9656 const double *n_k = nk + 2*dof2nk[k];
9659 for (
int j = 0; j <= p; j++)
9660 for (
int i = 0; i + j <= p; i++)
9662 double s = shape_x(i)*shape_y(j)*shape_l(p-i-j);
9663 T(o++, k) = s*n_k[0];
9664 T(o++, k) = s*n_k[1];
9666 for (
int i = 0; i <= p; i++)
9668 double s = shape_x(i)*shape_y(p-i);
9669 T(o++, k) = s*((ip.
x - c)*n_k[0] + (ip.
y - c)*n_k[1]);
9680 const int p =
Order - 1;
9682 #ifdef MFEM_THREAD_SAFE
9683 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
9692 for (
int j = 0; j <= p; j++)
9693 for (
int i = 0; i + j <= p; i++)
9695 double s = shape_x(i)*shape_y(j)*shape_l(p-i-j);
9696 u(o,0) = s; u(o,1) = 0; o++;
9697 u(o,0) = 0; u(o,1) = s; o++;
9699 for (
int i = 0; i <= p; i++)
9701 double s = shape_x(i)*shape_y(p-i);
9702 u(o,0) = (ip.
x - c)*s;
9703 u(o,1) = (ip.
y - c)*s;
9713 const int p =
Order - 1;
9715 #ifdef MFEM_THREAD_SAFE
9716 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
9717 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
9726 for (
int j = 0; j <= p; j++)
9727 for (
int i = 0; i + j <= p; i++)
9730 divu(o++) = (dshape_x(i)*shape_l(k) -
9731 shape_x(i)*dshape_l(k))*shape_y(j);
9732 divu(o++) = (dshape_y(j)*shape_l(k) -
9733 shape_y(j)*dshape_l(k))*shape_x(i);
9735 for (
int i = 0; i <= p; i++)
9738 divu(o++) = ((shape_x(i) + (ip.
x - c)*dshape_x(i))*shape_y(j) +
9739 (shape_y(j) + (ip.
y - c)*dshape_y(j))*shape_x(i));
9742 Ti.
Mult(divu, divshape);
9746 const double RT_TetrahedronElement::nk[12] =
9747 { 1,1,1, -1,0,0, 0,-1,0, 0,0,-1 };
9750 const double RT_TetrahedronElement::c = 1./4.;
9760 #ifndef MFEM_THREAD_SAFE
9772 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
9778 for (
int j = 0; j <= p; j++)
9779 for (
int i = 0; i + j <= p; i++)
9781 double w = bop[i] + bop[j] + bop[p-i-j];
9785 for (
int j = 0; j <= p; j++)
9786 for (
int i = 0; i + j <= p; i++)
9788 double w = bop[i] + bop[j] + bop[p-i-j];
9792 for (
int j = 0; j <= p; j++)
9793 for (
int i = 0; i + j <= p; i++)
9795 double w = bop[i] + bop[j] + bop[p-i-j];
9799 for (
int j = 0; j <= p; j++)
9800 for (
int i = 0; i + j <= p; i++)
9802 double w = bop[i] + bop[j] + bop[p-i-j];
9808 for (
int k = 0; k < p; k++)
9809 for (
int j = 0; j + k < p; j++)
9810 for (
int i = 0; i + j + k < p; i++)
9812 double w = iop[i] + iop[j] + iop[k] + iop[p-1-i-j-k];
9822 for (
int m = 0; m <
Dof; m++)
9829 const double *nm = nk + 3*dof2nk[m];
9832 for (
int k = 0; k <= p; k++)
9833 for (
int j = 0; j + k <= p; j++)
9834 for (
int i = 0; i + j + k <= p; i++)
9836 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
9837 T(o++, m) = s * nm[0];
9838 T(o++, m) = s * nm[1];
9839 T(o++, m) = s * nm[2];
9841 for (
int j = 0; j <= p; j++)
9842 for (
int i = 0; i + j <= p; i++)
9844 double s = shape_x(i)*shape_y(j)*shape_z(p-i-j);
9845 T(o++, m) = s*((ip.
x - c)*nm[0] + (ip.
y - c)*nm[1] +
9857 const int p =
Order - 1;
9859 #ifdef MFEM_THREAD_SAFE
9860 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
9870 for (
int k = 0; k <= p; k++)
9871 for (
int j = 0; j + k <= p; j++)
9872 for (
int i = 0; i + j + k <= p; i++)
9874 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
9875 u(o,0) = s; u(o,1) = 0; u(o,2) = 0; o++;
9876 u(o,0) = 0; u(o,1) = s; u(o,2) = 0; o++;
9877 u(o,0) = 0; u(o,1) = 0; u(o,2) = s; o++;
9879 for (
int j = 0; j <= p; j++)
9880 for (
int i = 0; i + j <= p; i++)
9882 double s = shape_x(i)*shape_y(j)*shape_z(p-i-j);
9883 u(o,0) = (ip.
x - c)*s; u(o,1) = (ip.
y - c)*s; u(o,2) = (ip.
z - c)*s;
9893 const int p =
Order - 1;
9895 #ifdef MFEM_THREAD_SAFE
9896 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
9897 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
9907 for (
int k = 0; k <= p; k++)
9908 for (
int j = 0; j + k <= p; j++)
9909 for (
int i = 0; i + j + k <= p; i++)
9911 int l = p - i - j - k;
9912 divu(o++) = (dshape_x(i)*shape_l(l) -
9913 shape_x(i)*dshape_l(l))*shape_y(j)*shape_z(k);
9914 divu(o++) = (dshape_y(j)*shape_l(l) -
9915 shape_y(j)*dshape_l(l))*shape_x(i)*shape_z(k);
9916 divu(o++) = (dshape_z(k)*shape_l(l) -
9917 shape_z(k)*dshape_l(l))*shape_x(i)*shape_y(j);
9919 for (
int j = 0; j <= p; j++)
9920 for (
int i = 0; i + j <= p; i++)
9924 (shape_x(i) + (ip.
x - c)*dshape_x(i))*shape_y(j)*shape_z(k) +
9925 (shape_y(j) + (ip.
y - c)*dshape_y(j))*shape_x(i)*shape_z(k) +
9926 (shape_z(k) + (ip.
z - c)*dshape_z(k))*shape_x(i)*shape_y(j);
9929 Ti.
Mult(divu, divshape);
9933 const double ND_HexahedronElement::tk[18] =
9934 { 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,0.,0., 0.,-1.,0., 0.,0.,-1. };
9937 const int cp_type,
const int op_type)
9940 cbasis1d(
poly1d.ClosedBasis(p, VerifyClosed(cp_type))),
9941 obasis1d(
poly1d.OpenBasis(p - 1, VerifyOpen(op_type))),
9942 dof_map(Dof), dof2tk(Dof)
9946 const int dof3 =
Dof/3;
9948 #ifndef MFEM_THREAD_SAFE
9962 for (
int i = 0; i < p; i++)
9964 dof_map[0*dof3 + i + (0 + 0*(p + 1))*p] = o++;
9966 for (
int i = 0; i < p; i++)
9968 dof_map[1*dof3 + p + (i + 0*p)*(p + 1)] = o++;
9970 for (
int i = 0; i < p; i++)
9972 dof_map[0*dof3 + i + (p + 0*(p + 1))*p] = o++;
9974 for (
int i = 0; i < p; i++)
9976 dof_map[1*dof3 + 0 + (i + 0*p)*(p + 1)] = o++;
9978 for (
int i = 0; i < p; i++)
9980 dof_map[0*dof3 + i + (0 + p*(p + 1))*p] = o++;
9982 for (
int i = 0; i < p; i++)
9984 dof_map[1*dof3 + p + (i + p*p)*(p + 1)] = o++;
9986 for (
int i = 0; i < p; i++)
9988 dof_map[0*dof3 + i + (p + p*(p + 1))*p] = o++;
9990 for (
int i = 0; i < p; i++)
9992 dof_map[1*dof3 + 0 + (i + p*p)*(p + 1)] = o++;
9994 for (
int i = 0; i < p; i++)
9996 dof_map[2*dof3 + 0 + (0 + i*(p + 1))*(p + 1)] = o++;
9998 for (
int i = 0; i < p; i++)
10000 dof_map[2*dof3 + p + (0 + i*(p + 1))*(p + 1)] = o++;
10002 for (
int i = 0; i < p; i++)
10004 dof_map[2*dof3 + p + (p + i*(p + 1))*(p + 1)] = o++;
10006 for (
int i = 0; i < p; i++)
10008 dof_map[2*dof3 + 0 + (p + i*(p + 1))*(p + 1)] = o++;
10013 for (
int j = 1; j < p; j++)
10014 for (
int i = 0; i < p; i++)
10016 dof_map[0*dof3 + i + ((p - j) + 0*(p + 1))*p] = o++;
10018 for (
int j = 0; j < p; j++)
10019 for (
int i = 1; i < p; i++)
10021 dof_map[1*dof3 + i + ((p - 1 - j) + 0*p)*(p + 1)] = -1 - (o++);
10024 for (
int k = 1; k < p; k++)
10025 for (
int i = 0; i < p; i++)
10027 dof_map[0*dof3 + i + (0 + k*(p + 1))*p] = o++;
10029 for (
int k = 0; k < p; k++)
10030 for (
int i = 1; i < p; i++ )
10032 dof_map[2*dof3 + i + (0 + k*(p + 1))*(p + 1)] = o++;
10035 for (
int k = 1; k < p; k++)
10036 for (
int j = 0; j < p; j++)
10038 dof_map[1*dof3 + p + (j + k*p)*(p + 1)] = o++;
10040 for (
int k = 0; k < p; k++)
10041 for (
int j = 1; j < p; j++)
10043 dof_map[2*dof3 + p + (j + k*(p + 1))*(p + 1)] = o++;
10046 for (
int k = 1; k < p; k++)
10047 for (
int i = 0; i < p; i++)
10049 dof_map[0*dof3 + (p - 1 - i) + (p + k*(p + 1))*p] = -1 - (o++);
10051 for (
int k = 0; k < p; k++)
10052 for (
int i = 1; i < p; i++)
10054 dof_map[2*dof3 + (p - i) + (p + k*(p + 1))*(p + 1)] = o++;
10057 for (
int k = 1; k < p; k++)
10058 for (
int j = 0; j < p; j++)
10060 dof_map[1*dof3 + 0 + ((p - 1 - j) + k*p)*(p + 1)] = -1 - (o++);
10062 for (
int k = 0; k < p; k++)
10063 for (
int j = 1; j < p; j++)
10065 dof_map[2*dof3 + 0 + ((p - j) + k*(p + 1))*(p + 1)] = o++;
10068 for (
int j = 1; j < p; j++)
10069 for (
int i = 0; i < p; i++)
10071 dof_map[0*dof3 + i + (j + p*(p + 1))*p] = o++;
10073 for (
int j = 0; j < p; j++)
10074 for (
int i = 1; i < p; i++)
10076 dof_map[1*dof3 + i + (j + p*p)*(p + 1)] = o++;
10081 for (
int k = 1; k < p; k++)
10082 for (
int j = 1; j < p; j++)
10083 for (
int i = 0; i < p; i++)
10085 dof_map[0*dof3 + i + (j + k*(p + 1))*p] = o++;
10088 for (
int k = 1; k < p; k++)
10089 for (
int j = 0; j < p; j++)
10090 for (
int i = 1; i < p; i++)
10092 dof_map[1*dof3 + i + (j + k*p)*(p + 1)] = o++;
10095 for (
int k = 0; k < p; k++)
10096 for (
int j = 1; j < p; j++)
10097 for (
int i = 1; i < p; i++)
10099 dof_map[2*dof3 + i + (j + k*(p + 1))*(p + 1)] = o++;
10105 for (
int k = 0; k <= p; k++)
10106 for (
int j = 0; j <= p; j++)
10107 for (
int i = 0; i < p; i++)
10110 if ((idx = dof_map[o++]) < 0)
10112 dof2tk[idx = -1 - idx] = 3;
10121 for (
int k = 0; k <= p; k++)
10122 for (
int j = 0; j < p; j++)
10123 for (
int i = 0; i <= p; i++)
10126 if ((idx = dof_map[o++]) < 0)
10128 dof2tk[idx = -1 - idx] = 4;
10137 for (
int k = 0; k < p; k++)
10138 for (
int j = 0; j <= p; j++)
10139 for (
int i = 0; i <= p; i++)
10142 if ((idx = dof_map[o++]) < 0)
10144 dof2tk[idx = -1 - idx] = 5;
10157 const int p =
Order;
10159 #ifdef MFEM_THREAD_SAFE
10160 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10161 Vector shape_cz(p + 1), shape_oz(p);
10164 cbasis1d.
Eval(ip.
x, shape_cx);
10165 obasis1d.
Eval(ip.
x, shape_ox);
10166 cbasis1d.
Eval(ip.
y, shape_cy);
10167 obasis1d.
Eval(ip.
y, shape_oy);
10168 cbasis1d.
Eval(ip.
z, shape_cz);
10169 obasis1d.
Eval(ip.
z, shape_oz);
10173 for (
int k = 0; k <= p; k++)
10174 for (
int j = 0; j <= p; j++)
10175 for (
int i = 0; i < p; i++)
10178 if ((idx = dof_map[o++]) < 0)
10180 idx = -1 - idx, s = -1;
10186 shape(idx,0) = s*shape_ox(i)*shape_cy(j)*shape_cz(k);
10191 for (
int k = 0; k <= p; k++)
10192 for (
int j = 0; j < p; j++)
10193 for (
int i = 0; i <= p; i++)
10196 if ((idx = dof_map[o++]) < 0)
10198 idx = -1 - idx, s = -1;
10205 shape(idx,1) = s*shape_cx(i)*shape_oy(j)*shape_cz(k);
10209 for (
int k = 0; k < p; k++)
10210 for (
int j = 0; j <= p; j++)
10211 for (
int i = 0; i <= p; i++)
10214 if ((idx = dof_map[o++]) < 0)
10216 idx = -1 - idx, s = -1;
10224 shape(idx,2) = s*shape_cx(i)*shape_cy(j)*shape_oz(k);
10231 const int p =
Order;
10233 #ifdef MFEM_THREAD_SAFE
10234 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10235 Vector shape_cz(p + 1), shape_oz(p);
10236 Vector dshape_cx(p + 1), dshape_cy(p + 1), dshape_cz(p + 1);
10239 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
10240 obasis1d.
Eval(ip.
x, shape_ox);
10241 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
10242 obasis1d.
Eval(ip.
y, shape_oy);
10243 cbasis1d.
Eval(ip.
z, shape_cz, dshape_cz);
10244 obasis1d.
Eval(ip.
z, shape_oz);
10248 for (
int k = 0; k <= p; k++)
10249 for (
int j = 0; j <= p; j++)
10250 for (
int i = 0; i < p; i++)
10253 if ((idx = dof_map[o++]) < 0)
10255 idx = -1 - idx, s = -1;
10261 curl_shape(idx,0) = 0.;
10262 curl_shape(idx,1) = s*shape_ox(i)* shape_cy(j)*dshape_cz(k);
10263 curl_shape(idx,2) = -s*shape_ox(i)*dshape_cy(j)* shape_cz(k);
10266 for (
int k = 0; k <= p; k++)
10267 for (
int j = 0; j < p; j++)
10268 for (
int i = 0; i <= p; i++)
10271 if ((idx = dof_map[o++]) < 0)
10273 idx = -1 - idx, s = -1;
10279 curl_shape(idx,0) = -s* shape_cx(i)*shape_oy(j)*dshape_cz(k);
10280 curl_shape(idx,1) = 0.;
10281 curl_shape(idx,2) = s*dshape_cx(i)*shape_oy(j)* shape_cz(k);
10284 for (
int k = 0; k < p; k++)
10285 for (
int j = 0; j <= p; j++)
10286 for (
int i = 0; i <= p; i++)
10289 if ((idx = dof_map[o++]) < 0)
10291 idx = -1 - idx, s = -1;
10297 curl_shape(idx,0) = s* shape_cx(i)*dshape_cy(j)*shape_oz(k);
10298 curl_shape(idx,1) = -s*dshape_cx(i)* shape_cy(j)*shape_oz(k);
10299 curl_shape(idx,2) = 0.;
10304 const double ND_QuadrilateralElement::tk[8] =
10305 { 1.,0., 0.,1., -1.,0., 0.,-1. };
10312 cbasis1d(
poly1d.ClosedBasis(p, VerifyClosed(cp_type))),
10313 obasis1d(
poly1d.OpenBasis(p - 1, VerifyOpen(op_type))),
10314 dof_map(Dof), dof2tk(Dof)
10318 const int dof2 =
Dof/2;
10320 #ifndef MFEM_THREAD_SAFE
10331 for (
int i = 0; i < p; i++)
10333 dof_map[0*dof2 + i + 0*p] = o++;
10335 for (
int j = 0; j < p; j++)
10337 dof_map[1*dof2 + p + j*(p + 1)] = o++;
10339 for (
int i = 0; i < p; i++)
10341 dof_map[0*dof2 + (p - 1 - i) + p*p] = -1 - (o++);
10343 for (
int j = 0; j < p; j++)
10345 dof_map[1*dof2 + 0 + (p - 1 - j)*(p + 1)] = -1 - (o++);
10350 for (
int j = 1; j < p; j++)
10351 for (
int i = 0; i < p; i++)
10353 dof_map[0*dof2 + i + j*p] = o++;
10356 for (
int j = 0; j < p; j++)
10357 for (
int i = 1; i < p; i++)
10359 dof_map[1*dof2 + i + j*(p + 1)] = o++;
10365 for (
int j = 0; j <= p; j++)
10366 for (
int i = 0; i < p; i++)
10369 if ((idx = dof_map[o++]) < 0)
10371 dof2tk[idx = -1 - idx] = 2;
10380 for (
int j = 0; j < p; j++)
10381 for (
int i = 0; i <= p; i++)
10384 if ((idx = dof_map[o++]) < 0)
10386 dof2tk[idx = -1 - idx] = 3;
10399 const int p =
Order;
10401 #ifdef MFEM_THREAD_SAFE
10402 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10405 cbasis1d.
Eval(ip.
x, shape_cx);
10406 obasis1d.
Eval(ip.
x, shape_ox);
10407 cbasis1d.
Eval(ip.
y, shape_cy);
10408 obasis1d.
Eval(ip.
y, shape_oy);
10412 for (
int j = 0; j <= p; j++)
10413 for (
int i = 0; i < p; i++)
10416 if ((idx = dof_map[o++]) < 0)
10418 idx = -1 - idx, s = -1;
10424 shape(idx,0) = s*shape_ox(i)*shape_cy(j);
10428 for (
int j = 0; j < p; j++)
10429 for (
int i = 0; i <= p; i++)
10432 if ((idx = dof_map[o++]) < 0)
10434 idx = -1 - idx, s = -1;
10441 shape(idx,1) = s*shape_cx(i)*shape_oy(j);
10448 const int p =
Order;
10450 #ifdef MFEM_THREAD_SAFE
10451 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10452 Vector dshape_cx(p + 1), dshape_cy(p + 1);
10455 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
10456 obasis1d.
Eval(ip.
x, shape_ox);
10457 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
10458 obasis1d.
Eval(ip.
y, shape_oy);
10462 for (
int j = 0; j <= p; j++)
10463 for (
int i = 0; i < p; i++)
10466 if ((idx = dof_map[o++]) < 0)
10468 idx = -1 - idx, s = -1;
10474 curl_shape(idx,0) = -s*shape_ox(i)*dshape_cy(j);
10477 for (
int j = 0; j < p; j++)
10478 for (
int i = 0; i <= p; i++)
10481 if ((idx = dof_map[o++]) < 0)
10483 idx = -1 - idx, s = -1;
10489 curl_shape(idx,0) = s*dshape_cx(i)*shape_oy(j);
10494 const double ND_TetrahedronElement::tk[18] =
10495 { 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,1.,0., -1.,0.,1., 0.,-1.,1. };
10497 const double ND_TetrahedronElement::c = 1./4.;
10507 const int pm1 = p - 1, pm2 = p - 2, pm3 = p - 3;
10509 #ifndef MFEM_THREAD_SAFE
10520 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
10525 for (
int i = 0; i < p; i++)
10530 for (
int i = 0; i < p; i++)
10535 for (
int i = 0; i < p; i++)
10540 for (
int i = 0; i < p; i++)
10545 for (
int i = 0; i < p; i++)
10550 for (
int i = 0; i < p; i++)
10557 for (
int j = 0; j <= pm2; j++)
10558 for (
int i = 0; i + j <= pm2; i++)
10560 double w = fop[i] + fop[j] + fop[pm2-i-j];
10566 for (
int j = 0; j <= pm2; j++)
10567 for (
int i = 0; i + j <= pm2; i++)
10569 double w = fop[i] + fop[j] + fop[pm2-i-j];
10575 for (
int j = 0; j <= pm2; j++)
10576 for (
int i = 0; i + j <= pm2; i++)
10578 double w = fop[i] + fop[j] + fop[pm2-i-j];
10584 for (
int j = 0; j <= pm2; j++)
10585 for (
int i = 0; i + j <= pm2; i++)
10587 double w = fop[i] + fop[j] + fop[pm2-i-j];
10595 for (
int k = 0; k <= pm3; k++)
10596 for (
int j = 0; j + k <= pm3; j++)
10597 for (
int i = 0; i + j + k <= pm3; i++)
10599 double w = iop[i] + iop[j] + iop[k] + iop[pm3-i-j-k];
10609 for (
int m = 0; m <
Dof; m++)
10612 const double *tm = tk + 3*dof2tk[m];
10620 for (
int k = 0; k <= pm1; k++)
10621 for (
int j = 0; j + k <= pm1; j++)
10622 for (
int i = 0; i + j + k <= pm1; i++)
10624 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
10625 T(o++, m) = s * tm[0];
10626 T(o++, m) = s * tm[1];
10627 T(o++, m) = s * tm[2];
10629 for (
int k = 0; k <= pm1; k++)
10630 for (
int j = 0; j + k <= pm1; j++)
10632 double s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
10633 T(o++, m) = s*((ip.
y - c)*tm[0] - (ip.
x - c)*tm[1]);
10634 T(o++, m) = s*((ip.
z - c)*tm[0] - (ip.
x - c)*tm[2]);
10636 for (
int k = 0; k <= pm1; k++)
10639 shape_y(pm1-k)*shape_z(k)*((ip.
z - c)*tm[1] - (ip.
y - c)*tm[2]);
10650 const int pm1 =
Order - 1;
10652 #ifdef MFEM_THREAD_SAFE
10653 const int p =
Order;
10654 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
10664 for (
int k = 0; k <= pm1; k++)
10665 for (
int j = 0; j + k <= pm1; j++)
10666 for (
int i = 0; i + j + k <= pm1; i++)
10668 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
10669 u(n,0) = s; u(n,1) = 0.; u(n,2) = 0.; n++;
10670 u(n,0) = 0.; u(n,1) = s; u(n,2) = 0.; n++;
10671 u(n,0) = 0.; u(n,1) = 0.; u(n,2) = s; n++;
10673 for (
int k = 0; k <= pm1; k++)
10674 for (
int j = 0; j + k <= pm1; j++)
10676 double s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
10677 u(n,0) = s*(ip.
y - c); u(n,1) = -s*(ip.
x - c); u(n,2) = 0.; n++;
10678 u(n,0) = s*(ip.
z - c); u(n,1) = 0.; u(n,2) = -s*(ip.
x - c); n++;
10680 for (
int k = 0; k <= pm1; k++)
10682 double s = shape_y(pm1-k)*shape_z(k);
10683 u(n,0) = 0.; u(n,1) = s*(ip.
z - c); u(n,2) = -s*(ip.
y - c); n++;
10692 const int pm1 =
Order - 1;
10694 #ifdef MFEM_THREAD_SAFE
10695 const int p =
Order;
10696 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
10697 Vector dshape_x(p), dshape_y(p), dshape_z(p), dshape_l(p);
10707 for (
int k = 0; k <= pm1; k++)
10708 for (
int j = 0; j + k <= pm1; j++)
10709 for (
int i = 0; i + j + k <= pm1; i++)
10712 const double dx = (dshape_x(i)*shape_l(l) -
10713 shape_x(i)*dshape_l(l))*shape_y(j)*shape_z(k);
10714 const double dy = (dshape_y(j)*shape_l(l) -
10715 shape_y(j)*dshape_l(l))*shape_x(i)*shape_z(k);
10716 const double dz = (dshape_z(k)*shape_l(l) -
10717 shape_z(k)*dshape_l(l))*shape_x(i)*shape_y(j);
10719 u(n,0) = 0.; u(n,1) = dz; u(n,2) = -dy; n++;
10720 u(n,0) = -dz; u(n,1) = 0.; u(n,2) = dx; n++;
10721 u(n,0) = dy; u(n,1) = -dx; u(n,2) = 0.; n++;
10723 for (
int k = 0; k <= pm1; k++)
10724 for (
int j = 0; j + k <= pm1; j++)
10726 int i = pm1 - j - k;
10729 u(n,0) = shape_x(i)*(ip.
x - c)*shape_y(j)*dshape_z(k);
10730 u(n,1) = shape_x(i)*shape_y(j)*(ip.
y - c)*dshape_z(k);
10732 -((dshape_x(i)*(ip.
x - c) + shape_x(i))*shape_y(j)*shape_z(k) +
10733 (dshape_y(j)*(ip.
y - c) + shape_y(j))*shape_x(i)*shape_z(k));
10736 u(n,0) = -shape_x(i)*(ip.
x - c)*dshape_y(j)*shape_z(k);
10737 u(n,1) = (shape_x(i)*shape_y(j)*(dshape_z(k)*(ip.
z - c) + shape_z(k)) +
10738 (dshape_x(i)*(ip.
x - c) + shape_x(i))*shape_y(j)*shape_z(k));
10739 u(n,2) = -shape_x(i)*dshape_y(j)*shape_z(k)*(ip.
z - c);
10742 for (
int k = 0; k <= pm1; k++)
10746 u(n,0) = -((dshape_y(j)*(ip.
y - c) + shape_y(j))*shape_z(k) +
10747 shape_y(j)*(dshape_z(k)*(ip.
z - c) + shape_z(k)));
10752 Ti.
Mult(u, curl_shape);
10756 const double ND_TriangleElement::tk[8] =
10757 { 1.,0., -1.,1., 0.,-1., 0.,1. };
10759 const double ND_TriangleElement::c = 1./3.;
10769 const int pm1 = p - 1, pm2 = p - 2;
10771 #ifndef MFEM_THREAD_SAFE
10781 Vector shape_x(p), shape_y(p), shape_l(p);
10786 for (
int i = 0; i < p; i++)
10791 for (
int i = 0; i < p; i++)
10796 for (
int i = 0; i < p; i++)
10803 for (
int j = 0; j <= pm2; j++)
10804 for (
int i = 0; i + j <= pm2; i++)
10806 double w = iop[i] + iop[j] + iop[pm2-i-j];
10814 for (
int m = 0; m <
Dof; m++)
10817 const double *tm = tk + 2*dof2tk[m];
10824 for (
int j = 0; j <= pm1; j++)
10825 for (
int i = 0; i + j <= pm1; i++)
10827 double s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
10828 T(n++, m) = s * tm[0];
10829 T(n++, m) = s * tm[1];
10831 for (
int j = 0; j <= pm1; j++)
10834 shape_x(pm1-j)*shape_y(j)*((ip.
y - c)*tm[0] - (ip.
x - c)*tm[1]);
10845 const int pm1 =
Order - 1;
10847 #ifdef MFEM_THREAD_SAFE
10848 const int p =
Order;
10849 Vector shape_x(p), shape_y(p), shape_l(p);
10858 for (
int j = 0; j <= pm1; j++)
10859 for (
int i = 0; i + j <= pm1; i++)
10861 double s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
10862 u(n,0) = s; u(n,1) = 0; n++;
10863 u(n,0) = 0; u(n,1) = s; n++;
10865 for (
int j = 0; j <= pm1; j++)
10867 double s = shape_x(pm1-j)*shape_y(j);
10868 u(n,0) = s*(ip.
y - c);
10869 u(n,1) = -s*(ip.
x - c);
10879 const int pm1 =
Order - 1;
10881 #ifdef MFEM_THREAD_SAFE
10882 const int p =
Order;
10883 Vector shape_x(p), shape_y(p), shape_l(p);
10884 Vector dshape_x(p), dshape_y(p), dshape_l(p);
10893 for (
int j = 0; j <= pm1; j++)
10894 for (
int i = 0; i + j <= pm1; i++)
10897 const double dx = (dshape_x(i)*shape_l(l) -
10898 shape_x(i)*dshape_l(l)) * shape_y(j);
10899 const double dy = (dshape_y(j)*shape_l(l) -
10900 shape_y(j)*dshape_l(l)) * shape_x(i);
10906 for (
int j = 0; j <= pm1; j++)
10910 curlu(n++) = -((dshape_x(i)*(ip.
x - c) + shape_x(i)) * shape_y(j) +
10911 (dshape_y(j)*(ip.
y - c) + shape_y(j)) * shape_x(i));
10915 Ti.
Mult(curlu, curl2d);
10919 const double ND_SegmentElement::tk[1] = { 1. };
10924 obasis1d(
poly1d.OpenBasis(p - 1, VerifyOpen(op_type))),
10930 for (
int i = 0; i < p; i++)
10949 kv[0]->CalcShape(shape,
ijk[0], ip.
x);
10952 for (
int i = 0; i <=
Order; i++)
10954 sum += (shape(i) *=
weights(i));
10966 kv[0]->CalcDShape(grad,
ijk[0], ip.
x);
10968 double sum = 0.0, dsum = 0.0;
10969 for (
int i = 0; i <=
Order; i++)
10972 dsum += ( grad(i) *=
weights(i));
10976 add(sum, grad, -dsum*sum*sum,
shape_x, grad);
10986 for (
int o = 0, j = 0; j <=
Order; j++)
10988 const double sy =
shape_y(j);
10989 for (
int i = 0; i <=
Order; i++, o++)
11001 double sum, dsum[2];
11009 sum = dsum[0] = dsum[1] = 0.0;
11010 for (
int o = 0, j = 0; j <=
Order; j++)
11013 for (
int i = 0; i <=
Order; i++, o++)
11023 dsum[0] *= sum*sum;
11024 dsum[1] *= sum*sum;
11026 for (
int o = 0; o <
Dof; o++)
11028 dshape(o,0) = dshape(o,0)*sum -
u(o)*dsum[0];
11029 dshape(o,1) = dshape(o,1)*sum -
u(o)*dsum[1];
11041 for (
int o = 0, k = 0; k <=
Order; k++)
11043 const double sz =
shape_z(k);
11044 for (
int j = 0; j <=
Order; j++)
11046 const double sy_sz =
shape_y(j)*sz;
11047 for (
int i = 0; i <=
Order; i++, o++)
11060 double sum, dsum[3];
11070 sum = dsum[0] = dsum[1] = dsum[2] = 0.0;
11071 for (
int o = 0, k = 0; k <=
Order; k++)
11074 for (
int j = 0; j <=
Order; j++)
11076 const double sy_sz =
shape_y(j)* sz;
11077 const double dsy_sz =
dshape_y(j)* sz;
11078 const double sy_dsz =
shape_y(j)*dsz;
11079 for (
int i = 0; i <=
Order; i++, o++)
11091 dsum[0] *= sum*sum;
11092 dsum[1] *= sum*sum;
11093 dsum[2] *= sum*sum;
11095 for (
int o = 0; o <
Dof; o++)
11097 dshape(o,0) = dshape(o,0)*sum -
u(o)*dsum[0];
11098 dshape(o,1) = dshape(o,1)*sum -
u(o)*dsum[1];
11099 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()
int Size() const
Logical size of the array.
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
L2_SegmentElement(const int p, const int type=Quadrature1D::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...
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 double CalcDelta(const int p, const double x)
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
RT1TriangleFiniteElement()
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.
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...
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)
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
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
const double * OpenPoints(const int p, const int type=Quadrature1D::GaussLegendre)
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...
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
L2_TriangleElement(const int p, const int type=Quadrature1D::GaussLegendre)
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
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...
void LocalInterpolation_ND(const double *tk, const Array< int > &d2t, ElementTransformation &Trans, DenseMatrix &I) const
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &h) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
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)
BiQuadPos2DFiniteElement()
int GetOrder() const
Returns the order of the finite element.
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...
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
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 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
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
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
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
static void CalcShape(const int p, const double x, const double y, const double z, double *shape)
H1_SegmentElement(const int p, const int type=Quadrature1D::GaussLobatto)
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
H1_TetrahedronElement(const int p, const int type=Quadrature1D::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...
Linear2DFiniteElement()
Construct a linear FE on triangle.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
void add(const Vector &v1, const Vector &v2, Vector &v)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
L2_QuadrilateralElement(const int p, const int _type=Quadrature1D::GaussLegendre)
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
virtual void 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)
RT_QuadrilateralElement(const int p, const int cp_type=Quadrature1D::GaussLobatto, const int op_type=Quadrature1D::GaussLegendre)
const double * GetPoints(const int p, const int type)
Get the coordinates of the points of the given Quadrature1D type.
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...
Implements CalcDivShape methods.
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)
H1_HexahedronElement(const int p, const int type=Quadrature1D::GaussLobatto)
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
evaluate shape function - constant 1
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &h) const
static void CalcBinomTerms(const int p, const double x, const double y, double *u)
Compute the terms in the expansion of the binomial (x + y)^p.
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &h) const
H1Pos_TriangleElement(const int p)
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
static int VerifyOpen(int pt_type)
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
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...
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
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...
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
void LocalInterpolation_RT(const double *nk, const Array< int > &d2n, ElementTransformation &Trans, DenseMatrix &I) const
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
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...
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...
For scalar fields; preserves point values.
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
For scalar fields; preserves volume integrals.
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void ProjectDelta(int vertex, Vector &dofs) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
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...
Implements CalcCurlShape methods.
void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
VVt += a * v v^t.
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...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void 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...
void ProjectGrad_ND(const double *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
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
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
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
~LagrangeHexFiniteElement()
void CalcVShape_ND(ElementTransformation &Trans, DenseMatrix &shape) const
No derivatives implemented.
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 mfem_error(const char *msg)
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
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 poiner 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...
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
ND_HexahedronElement(const int p, const int cp_type=Quadrature1D::GaussLobatto, const int op_type=Quadrature1D::GaussLegendre)
Quad1DFiniteElement()
Construct a quadratic FE on interval.
RT_HexahedronElement(const int p, const int cp_type=Quadrature1D::GaussLobatto, const int op_type=Quadrature1D::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 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
Basis & ClosedBasis(const int p, const int type=Quadrature1D::GaussLobatto)
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)
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
Basis & GetBasis(const int p, const int type)
Get a Poly_1D::Basis object of the given degree and Quadrature1D type.
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...
virtual void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const
static int VerifyClosed(int pt_type)
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &elmat)
virtual void ProjectDelta(int vertex, Vector &dofs) const
RefinedLinear2DFiniteElement()
Construct a quadratic FE on triangle.
void NodalLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const NodalFiniteElement &fine_fe) const
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
void ProjectCurl_2D(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
RT_TetrahedronElement(const int p)
ND_QuadrilateralElement(const int p, const int cp_type=Quadrature1D::GaussLobatto, const int op_type=Quadrature1D::GaussLegendre)
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
L2_TetrahedronElement(const int p, const int type=Quadrature1D::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...
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.
const double * ClosedPoints(const int p, const int type=Quadrature1D::GaussLobatto)
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 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.
H1_QuadrilateralElement(const int p, const int type=Quadrature1D::GaussLobatto)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
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...
void PositiveLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const PositiveFiniteElement &fine_fe) const
GaussBiLinear2DFiniteElement()
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...
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
H1_TriangleElement(const int p, const int type=Quadrature1D::GaussLobatto)
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...
Basis(const int p, const double *nodes, const int _mode=1)
static void CalcDBinomTerms(const int p, const double x, const double y, double *d)
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
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.
ND_SegmentElement(const int p, const int op_type=Quadrature1D::GaussLegendre)
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...
BiLinear2DFiniteElement()
Construct a bilinear FE on quadrilateral.
L2_HexahedronElement(const int p, const int _type=Quadrature1D::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...
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