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 !");
129 mfem_error(
"FiniteElement::ProjectDelta(...) is not implemented for "
136 mfem_error(
"FiniteElement::Project(...) (fe version) is not implemented "
137 "for this element!");
144 mfem_error(
"FiniteElement::ProjectGrad(...) is not implemented for "
152 mfem_error(
"FiniteElement::ProjectCurl(...) is not implemented for "
160 mfem_error(
"FiniteElement::ProjectDiv(...) is not implemented for "
178 #ifdef MFEM_THREAD_SAFE
193 #ifdef MFEM_THREAD_SAFE
199 for (
int i = 0; i < fine_fe.
Dof; i++)
204 for (
int j = 0; j <
Dof; j++)
205 if (fabs (I (i,j) = c_shape (j)) < 1.0e-12)
227 for (
int i = 0; i <
Dof; i++)
230 for (
int j = 0; j < fe.
GetDof(); j++)
232 curl(i,j) = curl_shape(j,0);
240 for (
int i = 0; i <
Dof; i++)
246 dofs(i) = coeff.
Eval (Trans, ip);
249 dofs(i) *= Trans.
Weight();
257 MFEM_ASSERT(vc.
GetVDim() <= 3,
"");
262 for (
int i = 0; i <
Dof; i++)
266 vc.
Eval (x, Trans, ip);
271 for (
int j = 0; j < x.Size(); j++)
273 dofs(Dof*j+i) = v[j];
288 for (
int k = 0; k <
Dof; k++)
291 for (
int j = 0; j < shape.Size(); j++)
293 I(k,j) = (fabs(shape(j)) < 1e-12) ? 0.0 : shape(j);
302 for (
int k = 0; k <
Dof; k++)
313 I(k+d*Dof,j) =
vshape(j,d);
329 for (
int k = 0; k <
Dof; k++)
335 Mult(dshape, Jinv, grad_k);
340 for (
int j = 0; j < grad_k.Height(); j++)
341 for (
int d = 0; d <
Dim; d++)
343 grad(k+d*Dof,j) = grad_k(j,d);
356 for (
int k = 0; k <
Dof; k++)
364 for (
int j = 0; j < div_shape.Size(); j++)
366 div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j)/detJ;
371 for (
int j = 0; j < div_shape.Size(); j++)
373 div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j);
391 Vector fine_shape(fs), coarse_shape(cs);
392 DenseMatrix fine_mass(fs), fine_coarse_mass(fs, cs);
409 fine_mass_inv.
Mult(fine_coarse_mass, I);
422 for (
int i = 0; i <
Dof; i++)
426 dofs(i) = coeff.
Eval(Trans, ip);
452 pos_mass_inv.
Mult(mixed_mass, I);
457 void VectorFiniteElement::CalcShape (
460 mfem_error (
"Error: Cannot use scalar CalcShape(...) function with\n"
461 " VectorFiniteElements!");
464 void VectorFiniteElement::CalcDShape (
465 const IntegrationPoint &ip, DenseMatrix &dshape )
const
467 mfem_error (
"Error: Cannot use scalar CalcDShape(...) function with\n"
468 " VectorFiniteElements!");
500 MFEM_ABORT(
"Invalid dimension, Dim = " <<
Dim);
504 MFEM_ABORT(
"Invalid MapType = " <<
MapType);
512 #ifdef MFEM_THREAD_SAFE
517 shape *= (1.0 / Trans.
Weight());
524 #ifdef MFEM_THREAD_SAFE
537 MFEM_ASSERT(vc.
GetVDim() == sdim,
"");
539 const bool square_J = (
Dim == sdim);
541 for (
int k = 0; k <
Dof; k++)
547 if (!square_J) { dofs(k) /= Trans.
Weight(); }
562 for (
int k = 0; k <
Dof; k++)
571 double w = 1.0/Trans.
Weight();
572 for (
int d = 0; d <
Dim; d++)
578 for (
int j = 0; j < shape.Size(); j++)
585 for (
int d = 0; d < sdim; d++)
587 I(k,j+d*shape.Size()) = s*vk[d];
594 mfem_error(
"VectorFiniteElement::Project_RT (fe version)");
604 mfem_error(
"VectorFiniteElement::ProjectGrad_RT works only in 2D!");
612 for (
int k = 0; k <
Dof; k++)
615 tk[0] = nk[d2n[k]*
Dim+1];
616 tk[1] = -nk[d2n[k]*
Dim];
617 dshape.Mult(tk, grad_k);
618 for (
int j = 0; j < grad_k.Size(); j++)
620 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
629 #ifdef MFEM_THREAD_SAFE
642 for (
int k = 0; k <
Dof; k++)
649 J *= 1.0 / Trans.
Weight();
656 for (
int j = 0; j < curl_k.Size(); j++)
658 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
671 for (
int k = 0; k <
Dof; k++)
674 curl_shape.Mult(nk + d2n[k]*
Dim, curl_k);
675 for (
int j = 0; j < curl_k.Size(); j++)
677 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
689 for (
int k = 0; k <
Dof; k++)
710 for (
int k = 0; k <
Dof; k++)
719 double w = 1.0/Trans.
Weight();
720 for (
int d = 0; d < sdim; d++)
726 for (
int j = 0; j < shape.Size(); j++)
733 for (
int d = 0; d < sdim; d++)
735 I(k, j + d*shape.Size()) = s*vk[d];
742 mfem_error(
"VectorFiniteElement::Project_ND (fe version)");
756 for (
int k = 0; k <
Dof; k++)
759 dshape.Mult(tk + d2t[k]*
Dim, grad_k);
760 for (
int j = 0; j < grad_k.Size(); j++)
762 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
774 #ifdef MFEM_THREAD_SAFE
781 for (
int k = 0; k <
Dof; k++)
789 for (
int j = 0; j <
Dof; j++)
792 for (
int i = 0; i <
Dim; i++)
794 Ikj +=
vshape(j, i) * vk[i];
796 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
808 #ifdef MFEM_THREAD_SAFE
815 for (
int k = 0; k <
Dof; k++)
823 for (
int j = 0; j <
Dof; j++)
826 for (
int i = 0; i <
Dim; i++)
828 Ikj +=
vshape(j, i) * vk[i];
830 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
864 shape(0) = 1. - ip.
x;
889 shape(0) = 1. - ip.
x - ip.
y;
897 dshape(0,0) = -1.; dshape(0,1) = -1.;
898 dshape(1,0) = 1.; dshape(1,1) = 0.;
899 dshape(2,0) = 0.; dshape(2,1) = 1.;
918 shape(0) = (1. - ip.
x) * (1. - ip.
y) ;
919 shape(1) = ip.
x * (1. - ip.
y) ;
920 shape(2) = ip.
x * ip.
y ;
921 shape(3) = (1. - ip.
x) * ip.
y ;
927 dshape(0,0) = -1. + ip.
y; dshape(0,1) = -1. + ip.
x ;
928 dshape(1,0) = 1. - ip.
y; dshape(1,1) = -ip.
x ;
929 dshape(2,0) = ip.
y ; dshape(2,1) = ip.
x ;
930 dshape(3,0) = -ip.
y ; dshape(3,1) = 1. - ip.
x ;
936 h(0,0) = 0.; h(0,1) = 1.; h(0,2) = 0.;
937 h(1,0) = 0.; h(1,1) = -1.; h(1,2) = 0.;
938 h(2,0) = 0.; h(2,1) = 1.; h(2,2) = 0.;
939 h(3,0) = 0.; h(3,1) = -1.; h(3,2) = 0.;
957 const double x = ip.
x, y = ip.
y;
959 shape(0) = 5./3. - 2. * (x + y);
960 shape(1) = 2. * (x - 1./6.);
961 shape(2) = 2. * (y - 1./6.);
967 dshape(0,0) = -2.; dshape(0,1) = -2.;
968 dshape(1,0) = 2.; dshape(1,1) = 0.;
969 dshape(2,0) = 0.; dshape(2,1) = 2.;
974 dofs(vertex) = 2./3.;
975 dofs((vertex+1)%3) = 1./6.;
976 dofs((vertex+2)%3) = 1./6.;
981 const double GaussBiLinear2DFiniteElement::p[] =
982 { 0.2113248654051871177454256, 0.7886751345948128822545744 };
1000 const double x = ip.
x, y = ip.
y;
1002 shape(0) = 3. * (p[1] - x) * (p[1] - y);
1003 shape(1) = 3. * (x - p[0]) * (p[1] - y);
1004 shape(2) = 3. * (x - p[0]) * (y - p[0]);
1005 shape(3) = 3. * (p[1] - x) * (y - p[0]);
1011 const double x = ip.
x, y = ip.
y;
1013 dshape(0,0) = 3. * (y - p[1]); dshape(0,1) = 3. * (x - p[1]);
1014 dshape(1,0) = 3. * (p[1] - y); dshape(1,1) = 3. * (p[0] - x);
1015 dshape(2,0) = 3. * (y - p[0]); dshape(2,1) = 3. * (x - p[0]);
1016 dshape(3,0) = 3. * (p[0] - y); dshape(3,1) = 3. * (p[1] - x);
1022 dofs(vertex) = p[1]*p[1];
1023 dofs((vertex+1)%4) = p[0]*p[1];
1024 dofs((vertex+2)%4) = p[0]*p[0];
1025 dofs((vertex+3)%4) = p[0]*p[1];
1046 shape(0) = 1. - ip.
x - ip.
y;
1054 dshape(0,0) = -1.; dshape(0,1) = -1.;
1055 dshape(1,0) = 1.; dshape(1,1) = 0.;
1056 dshape(2,0) = 0.; dshape(2,1) = 1.;
1072 double l1 = 1.0 - x, l2 = x, l3 = 2. * x - 1.;
1074 shape(0) = l1 * (-l3);
1076 shape(2) = 4. * l1 * l2;
1084 dshape(0,0) = 4. * x - 3.;
1085 dshape(1,0) = 4. * x - 1.;
1086 dshape(2,0) = 4. - 8. * x;
1101 const double x = ip.
x, x1 = 1. - x;
1105 shape(2) = 2. * x * x1;
1111 const double x = ip.
x;
1113 dshape(0,0) = 2. * x - 2.;
1114 dshape(1,0) = 2. * x;
1115 dshape(2,0) = 2. - 4. * x;
1138 double x = ip.
x, y = ip.
y;
1139 double l1 = 1.-x-y, l2 = x, l3 = y;
1141 shape(0) = l1 * (2. * l1 - 1.);
1142 shape(1) = l2 * (2. * l2 - 1.);
1143 shape(2) = l3 * (2. * l3 - 1.);
1144 shape(3) = 4. * l1 * l2;
1145 shape(4) = 4. * l2 * l3;
1146 shape(5) = 4. * l3 * l1;
1152 double x = ip.
x, y = ip.
y;
1155 dshape(0,1) = 4. * (x + y) - 3.;
1157 dshape(1,0) = 4. * x - 1.;
1161 dshape(2,1) = 4. * y - 1.;
1163 dshape(3,0) = -4. * (2. * x + y - 1.);
1164 dshape(3,1) = -4. * x;
1166 dshape(4,0) = 4. * y;
1167 dshape(4,1) = 4. * x;
1169 dshape(5,0) = -4. * y;
1170 dshape(5,1) = -4. * (x + 2. * y - 1.);
1210 case 0: dofs(3) = 0.25; dofs(5) = 0.25;
break;
1211 case 1: dofs(3) = 0.25; dofs(4) = 0.25;
break;
1212 case 2: dofs(4) = 0.25; dofs(5) = 0.25;
break;
1218 const double GaussQuad2DFiniteElement::p[] =
1219 { 0.0915762135097707434595714634022015, 0.445948490915964886318329253883051 };
1237 for (
int i = 0; i < 6; i++)
1254 const double x = ip.
x, y = ip.
y;
1268 const double x = ip.
x, y = ip.
y;
1269 D(0,0) = 0.; D(0,1) = 0.;
1270 D(1,0) = 1.; D(1,1) = 0.;
1271 D(2,0) = 0.; D(2,1) = 1.;
1272 D(3,0) = 2. * x; D(3,1) = 0.;
1273 D(4,0) = y; D(4,1) = x;
1274 D(5,0) = 0.; D(5,1) = 2. * y;
1306 double x = ip.
x, y = ip.
y;
1307 double l1x, l2x, l3x, l1y, l2y, l3y;
1309 l1x = (x - 1.) * (2. * x - 1);
1310 l2x = 4. * x * (1. - x);
1311 l3x = x * (2. * x - 1.);
1312 l1y = (y - 1.) * (2. * y - 1);
1313 l2y = 4. * y * (1. - y);
1314 l3y = y * (2. * y - 1.);
1316 shape(0) = l1x * l1y;
1317 shape(4) = l2x * l1y;
1318 shape(1) = l3x * l1y;
1319 shape(7) = l1x * l2y;
1320 shape(8) = l2x * l2y;
1321 shape(5) = l3x * l2y;
1322 shape(3) = l1x * l3y;
1323 shape(6) = l2x * l3y;
1324 shape(2) = l3x * l3y;
1330 double x = ip.
x, y = ip.
y;
1331 double l1x, l2x, l3x, l1y, l2y, l3y;
1332 double d1x, d2x, d3x, d1y, d2y, d3y;
1334 l1x = (x - 1.) * (2. * x - 1);
1335 l2x = 4. * x * (1. - x);
1336 l3x = x * (2. * x - 1.);
1337 l1y = (y - 1.) * (2. * y - 1);
1338 l2y = 4. * y * (1. - y);
1339 l3y = y * (2. * y - 1.);
1348 dshape(0,0) = d1x * l1y;
1349 dshape(0,1) = l1x * d1y;
1351 dshape(4,0) = d2x * l1y;
1352 dshape(4,1) = l2x * d1y;
1354 dshape(1,0) = d3x * l1y;
1355 dshape(1,1) = l3x * d1y;
1357 dshape(7,0) = d1x * l2y;
1358 dshape(7,1) = l1x * d2y;
1360 dshape(8,0) = d2x * l2y;
1361 dshape(8,1) = l2x * d2y;
1363 dshape(5,0) = d3x * l2y;
1364 dshape(5,1) = l3x * d2y;
1366 dshape(3,0) = d1x * l3y;
1367 dshape(3,1) = l1x * d3y;
1369 dshape(6,0) = d2x * l3y;
1370 dshape(6,1) = l2x * d3y;
1372 dshape(2,0) = d3x * l3y;
1373 dshape(2,1) = l3x * d3y;
1385 case 0: dofs(4) = 0.25; dofs(7) = 0.25;
break;
1386 case 1: dofs(4) = 0.25; dofs(5) = 0.25;
break;
1387 case 2: dofs(5) = 0.25; dofs(6) = 0.25;
break;
1388 case 3: dofs(6) = 0.25; dofs(7) = 0.25;
break;
1420 double x = ip.
x, y = ip.
y;
1421 double l1x, l2x, l3x, l1y, l2y, l3y;
1423 l1x = (1. - x) * (1. - x);
1424 l2x = 2. * x * (1. - x);
1426 l1y = (1. - y) * (1. - y);
1427 l2y = 2. * y * (1. - y);
1430 shape(0) = l1x * l1y;
1431 shape(4) = l2x * l1y;
1432 shape(1) = l3x * l1y;
1433 shape(7) = l1x * l2y;
1434 shape(8) = l2x * l2y;
1435 shape(5) = l3x * l2y;
1436 shape(3) = l1x * l3y;
1437 shape(6) = l2x * l3y;
1438 shape(2) = l3x * l3y;
1444 double x = ip.
x, y = ip.
y;
1445 double l1x, l2x, l3x, l1y, l2y, l3y;
1446 double d1x, d2x, d3x, d1y, d2y, d3y;
1448 l1x = (1. - x) * (1. - x);
1449 l2x = 2. * x * (1. - x);
1451 l1y = (1. - y) * (1. - y);
1452 l2y = 2. * y * (1. - y);
1462 dshape(0,0) = d1x * l1y;
1463 dshape(0,1) = l1x * d1y;
1465 dshape(4,0) = d2x * l1y;
1466 dshape(4,1) = l2x * d1y;
1468 dshape(1,0) = d3x * l1y;
1469 dshape(1,1) = l3x * d1y;
1471 dshape(7,0) = d1x * l2y;
1472 dshape(7,1) = l1x * d2y;
1474 dshape(8,0) = d2x * l2y;
1475 dshape(8,1) = l2x * d2y;
1477 dshape(5,0) = d3x * l2y;
1478 dshape(5,1) = l3x * d2y;
1480 dshape(3,0) = d1x * l3y;
1481 dshape(3,1) = l1x * d3y;
1483 dshape(6,0) = d2x * l3y;
1484 dshape(6,1) = l2x * d3y;
1486 dshape(2,0) = d3x * l3y;
1487 dshape(2,1) = l3x * d3y;
1495 Vector xx(&tr_ip.
x, 2), shape(s, 9);
1497 for (
int i = 0; i < 9; i++)
1501 for (
int j = 0; j < 9; j++)
1502 if (fabs(I(i,j) = s[j]) < 1.0e-12)
1507 for (
int i = 0; i < 9; i++)
1509 double *d = &I(0,i);
1510 d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
1511 d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
1512 d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
1513 d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
1514 d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
1515 0.25 * (d[0] + d[1] + d[2] + d[3]);
1524 for (
int i = 0; i < 9; i++)
1528 d[i] = coeff.
Eval(Trans, ip);
1530 d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
1531 d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
1532 d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
1533 d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
1534 d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
1535 0.25 * (d[0] + d[1] + d[2] + d[3]);
1545 for (
int i = 0; i < 9; i++)
1549 vc.
Eval (x, Trans, ip);
1550 for (
int j = 0; j < x.Size(); j++)
1555 for (
int j = 0; j < x.Size(); j++)
1557 double *d = &dofs(9*j);
1559 d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
1560 d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
1561 d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
1562 d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
1563 d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
1564 0.25 * (d[0] + d[1] + d[2] + d[3]);
1572 const double p1 = 0.5*(1.-sqrt(3./5.));
1597 const double a = sqrt(5./3.);
1598 const double p1 = 0.5*(1.-sqrt(3./5.));
1600 double x = a*(ip.
x-p1), y = a*(ip.
y-p1);
1601 double l1x, l2x, l3x, l1y, l2y, l3y;
1603 l1x = (x - 1.) * (2. * x - 1);
1604 l2x = 4. * x * (1. - x);
1605 l3x = x * (2. * x - 1.);
1606 l1y = (y - 1.) * (2. * y - 1);
1607 l2y = 4. * y * (1. - y);
1608 l3y = y * (2. * y - 1.);
1610 shape(0) = l1x * l1y;
1611 shape(4) = l2x * l1y;
1612 shape(1) = l3x * l1y;
1613 shape(7) = l1x * l2y;
1614 shape(8) = l2x * l2y;
1615 shape(5) = l3x * l2y;
1616 shape(3) = l1x * l3y;
1617 shape(6) = l2x * l3y;
1618 shape(2) = l3x * l3y;
1624 const double a = sqrt(5./3.);
1625 const double p1 = 0.5*(1.-sqrt(3./5.));
1627 double x = a*(ip.
x-p1), y = a*(ip.
y-p1);
1628 double l1x, l2x, l3x, l1y, l2y, l3y;
1629 double d1x, d2x, d3x, d1y, d2y, d3y;
1631 l1x = (x - 1.) * (2. * x - 1);
1632 l2x = 4. * x * (1. - x);
1633 l3x = x * (2. * x - 1.);
1634 l1y = (y - 1.) * (2. * y - 1);
1635 l2y = 4. * y * (1. - y);
1636 l3y = y * (2. * y - 1.);
1638 d1x = a * (4. * x - 3.);
1639 d2x = a * (4. - 8. * x);
1640 d3x = a * (4. * x - 1.);
1641 d1y = a * (4. * y - 3.);
1642 d2y = a * (4. - 8. * y);
1643 d3y = a * (4. * y - 1.);
1645 dshape(0,0) = d1x * l1y;
1646 dshape(0,1) = l1x * d1y;
1648 dshape(4,0) = d2x * l1y;
1649 dshape(4,1) = l2x * d1y;
1651 dshape(1,0) = d3x * l1y;
1652 dshape(1,1) = l3x * d1y;
1654 dshape(7,0) = d1x * l2y;
1655 dshape(7,1) = l1x * d2y;
1657 dshape(8,0) = d2x * l2y;
1658 dshape(8,1) = l2x * d2y;
1660 dshape(5,0) = d3x * l2y;
1661 dshape(5,1) = l3x * d2y;
1663 dshape(3,0) = d1x * l3y;
1664 dshape(3,1) = l1x * d3y;
1666 dshape(6,0) = d2x * l3y;
1667 dshape(6,1) = l2x * d3y;
1669 dshape(2,0) = d3x * l3y;
1670 dshape(2,1) = l3x * d3y;
1713 double x = ip.
x, y = ip.
y;
1715 double w1x, w2x, w3x, w1y, w2y, w3y;
1716 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
1718 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
1719 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
1721 l0x = (- 4.5) * w1x * w2x * w3x;
1722 l1x = ( 13.5) * x * w2x * w3x;
1723 l2x = (-13.5) * x * w1x * w3x;
1724 l3x = ( 4.5) * x * w1x * w2x;
1726 l0y = (- 4.5) * w1y * w2y * w3y;
1727 l1y = ( 13.5) * y * w2y * w3y;
1728 l2y = (-13.5) * y * w1y * w3y;
1729 l3y = ( 4.5) * y * w1y * w2y;
1731 shape(0) = l0x * l0y;
1732 shape(1) = l3x * l0y;
1733 shape(2) = l3x * l3y;
1734 shape(3) = l0x * l3y;
1735 shape(4) = l1x * l0y;
1736 shape(5) = l2x * l0y;
1737 shape(6) = l3x * l1y;
1738 shape(7) = l3x * l2y;
1739 shape(8) = l2x * l3y;
1740 shape(9) = l1x * l3y;
1741 shape(10) = l0x * l2y;
1742 shape(11) = l0x * l1y;
1743 shape(12) = l1x * l1y;
1744 shape(13) = l2x * l1y;
1745 shape(14) = l1x * l2y;
1746 shape(15) = l2x * l2y;
1752 double x = ip.
x, y = ip.
y;
1754 double w1x, w2x, w3x, w1y, w2y, w3y;
1755 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
1756 double d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
1758 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
1759 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
1761 l0x = (- 4.5) * w1x * w2x * w3x;
1762 l1x = ( 13.5) * x * w2x * w3x;
1763 l2x = (-13.5) * x * w1x * w3x;
1764 l3x = ( 4.5) * x * w1x * w2x;
1766 l0y = (- 4.5) * w1y * w2y * w3y;
1767 l1y = ( 13.5) * y * w2y * w3y;
1768 l2y = (-13.5) * y * w1y * w3y;
1769 l3y = ( 4.5) * y * w1y * w2y;
1771 d0x = -5.5 + ( 18. - 13.5 * x) * x;
1772 d1x = 9. + (-45. + 40.5 * x) * x;
1773 d2x = -4.5 + ( 36. - 40.5 * x) * x;
1774 d3x = 1. + (- 9. + 13.5 * x) * x;
1776 d0y = -5.5 + ( 18. - 13.5 * y) * y;
1777 d1y = 9. + (-45. + 40.5 * y) * y;
1778 d2y = -4.5 + ( 36. - 40.5 * y) * y;
1779 d3y = 1. + (- 9. + 13.5 * y) * y;
1781 dshape( 0,0) = d0x * l0y; dshape( 0,1) = l0x * d0y;
1782 dshape( 1,0) = d3x * l0y; dshape( 1,1) = l3x * d0y;
1783 dshape( 2,0) = d3x * l3y; dshape( 2,1) = l3x * d3y;
1784 dshape( 3,0) = d0x * l3y; dshape( 3,1) = l0x * d3y;
1785 dshape( 4,0) = d1x * l0y; dshape( 4,1) = l1x * d0y;
1786 dshape( 5,0) = d2x * l0y; dshape( 5,1) = l2x * d0y;
1787 dshape( 6,0) = d3x * l1y; dshape( 6,1) = l3x * d1y;
1788 dshape( 7,0) = d3x * l2y; dshape( 7,1) = l3x * d2y;
1789 dshape( 8,0) = d2x * l3y; dshape( 8,1) = l2x * d3y;
1790 dshape( 9,0) = d1x * l3y; dshape( 9,1) = l1x * d3y;
1791 dshape(10,0) = d0x * l2y; dshape(10,1) = l0x * d2y;
1792 dshape(11,0) = d0x * l1y; dshape(11,1) = l0x * d1y;
1793 dshape(12,0) = d1x * l1y; dshape(12,1) = l1x * d1y;
1794 dshape(13,0) = d2x * l1y; dshape(13,1) = l2x * d1y;
1795 dshape(14,0) = d1x * l2y; dshape(14,1) = l1x * d2y;
1796 dshape(15,0) = d2x * l2y; dshape(15,1) = l2x * d2y;
1802 double x = ip.
x, y = ip.
y;
1804 double w1x, w2x, w3x, w1y, w2y, w3y;
1805 double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
1806 double d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
1807 double h0x, h1x, h2x, h3x, h0y, h1y, h2y, h3y;
1809 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
1810 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
1812 l0x = (- 4.5) * w1x * w2x * w3x;
1813 l1x = ( 13.5) * x * w2x * w3x;
1814 l2x = (-13.5) * x * w1x * w3x;
1815 l3x = ( 4.5) * x * w1x * w2x;
1817 l0y = (- 4.5) * w1y * w2y * w3y;
1818 l1y = ( 13.5) * y * w2y * w3y;
1819 l2y = (-13.5) * y * w1y * w3y;
1820 l3y = ( 4.5) * y * w1y * w2y;
1822 d0x = -5.5 + ( 18. - 13.5 * x) * x;
1823 d1x = 9. + (-45. + 40.5 * x) * x;
1824 d2x = -4.5 + ( 36. - 40.5 * x) * x;
1825 d3x = 1. + (- 9. + 13.5 * x) * x;
1827 d0y = -5.5 + ( 18. - 13.5 * y) * y;
1828 d1y = 9. + (-45. + 40.5 * y) * y;
1829 d2y = -4.5 + ( 36. - 40.5 * y) * y;
1830 d3y = 1. + (- 9. + 13.5 * y) * y;
1832 h0x = -27. * x + 18.;
1833 h1x = 81. * x - 45.;
1834 h2x = -81. * x + 36.;
1837 h0y = -27. * y + 18.;
1838 h1y = 81. * y - 45.;
1839 h2y = -81. * y + 36.;
1842 h( 0,0) = h0x * l0y; h( 0,1) = d0x * d0y; h( 0,2) = l0x * h0y;
1843 h( 1,0) = h3x * l0y; h( 1,1) = d3x * d0y; h( 1,2) = l3x * h0y;
1844 h( 2,0) = h3x * l3y; h( 2,1) = d3x * d3y; h( 2,2) = l3x * h3y;
1845 h( 3,0) = h0x * l3y; h( 3,1) = d0x * d3y; h( 3,2) = l0x * h3y;
1846 h( 4,0) = h1x * l0y; h( 4,1) = d1x * d0y; h( 4,2) = l1x * h0y;
1847 h( 5,0) = h2x * l0y; h( 5,1) = d2x * d0y; h( 5,2) = l2x * h0y;
1848 h( 6,0) = h3x * l1y; h( 6,1) = d3x * d1y; h( 6,2) = l3x * h1y;
1849 h( 7,0) = h3x * l2y; h( 7,1) = d3x * d2y; h( 7,2) = l3x * h2y;
1850 h( 8,0) = h2x * l3y; h( 8,1) = d2x * d3y; h( 8,2) = l2x * h3y;
1851 h( 9,0) = h1x * l3y; h( 9,1) = d1x * d3y; h( 9,2) = l1x * h3y;
1852 h(10,0) = h0x * l2y; h(10,1) = d0x * d2y; h(10,2) = l0x * h2y;
1853 h(11,0) = h0x * l1y; h(11,1) = d0x * d1y; h(11,2) = l0x * h1y;
1854 h(12,0) = h1x * l1y; h(12,1) = d1x * d1y; h(12,2) = l1x * h1y;
1855 h(13,0) = h2x * l1y; h(13,1) = d2x * d1y; h(13,2) = l2x * h1y;
1856 h(14,0) = h1x * l2y; h(14,1) = d1x * d2y; h(14,2) = l1x * h2y;
1857 h(15,0) = h2x * l2y; h(15,1) = d2x * d2y; h(15,2) = l2x * h2y;
1876 l3 = (0.33333333333333333333-x),
1877 l4 = (0.66666666666666666667-x);
1879 shape(0) = 4.5 * l2 * l3 * l4;
1880 shape(1) = 4.5 * l1 * l3 * l4;
1881 shape(2) = 13.5 * l1 * l2 * l4;
1882 shape(3) = -13.5 * l1 * l2 * l3;
1890 dshape(0,0) = -5.5 + x * (18. - 13.5 * x);
1891 dshape(1,0) = 1. - x * (9. - 13.5 * x);
1892 dshape(2,0) = 9. - x * (45. - 40.5 * x);
1893 dshape(3,0) = -4.5 + x * (36. - 40.5 * x);
1925 double x = ip.
x, y = ip.
y;
1926 double l1 = (-1. + x + y),
1930 shape(0) = -0.5*l1*(3.*l1 + 1.)*(3.*l1 + 2.);
1931 shape(1) = 0.5*x*(lx - 1.)*lx;
1932 shape(2) = 0.5*y*(-1. + ly)*ly;
1933 shape(3) = 4.5*x*l1*(3.*l1 + 1.);
1934 shape(4) = -4.5*x*lx*l1;
1935 shape(5) = 4.5*x*lx*y;
1936 shape(6) = 4.5*x*y*ly;
1937 shape(7) = -4.5*y*l1*ly;
1938 shape(8) = 4.5*y*l1*(1. + 3.*l1);
1939 shape(9) = -27.*x*y*l1;
1945 double x = ip.
x, y = ip.
y;
1947 dshape(0,0) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
1948 dshape(1,0) = 1. + 4.5*x*(-2. + 3.*x);
1950 dshape(3,0) = 4.5*(2. + 9.*x*x - 5.*y + 3.*y*y + 2.*x*(-5. + 6.*y));
1951 dshape(4,0) = -4.5*(1. - 1.*y + x*(-8. + 9.*x + 6.*y));
1952 dshape(5,0) = 4.5*(-1. + 6.*x)*y;
1953 dshape(6,0) = 4.5*y*(-1. + 3.*y);
1954 dshape(7,0) = 4.5*(1. - 3.*y)*y;
1955 dshape(8,0) = 4.5*y*(-5. + 6.*x + 6.*y);
1956 dshape(9,0) = -27.*y*(-1. + 2.*x + y);
1958 dshape(0,1) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
1960 dshape(2,1) = 1. + 4.5*y*(-2. + 3.*y);
1961 dshape(3,1) = 4.5*x*(-5. + 6.*x + 6.*y);
1962 dshape(4,1) = 4.5*(1. - 3.*x)*x;
1963 dshape(5,1) = 4.5*x*(-1. + 3.*x);
1964 dshape(6,1) = 4.5*x*(-1. + 6.*y);
1965 dshape(7,1) = -4.5*(1. + x*(-1. + 6.*y) + y*(-8. + 9.*y));
1966 dshape(8,1) = 4.5*(2. + 3.*x*x + y*(-10. + 9.*y) + x*(-5. + 12.*y));
1967 dshape(9,1) = -27.*x*(-1. + x + 2.*y);
1973 double x = ip.
x, y = ip.
y;
1975 h(0,0) = 18.-27.*(x+y);
1976 h(0,1) = 18.-27.*(x+y);
1977 h(0,2) = 18.-27.*(x+y);
1987 h(3,0) = -45.+81.*x+54.*y;
1988 h(3,1) = -22.5+54.*x+27.*y;
1991 h(4,0) = 36.-81.*x-27.*y;
1996 h(5,1) = -4.5+27.*x;
2000 h(6,1) = -4.5+27.*y;
2005 h(7,2) = 36.-27.*x-81.*y;
2008 h(8,1) = -22.5+27.*x+54.*y;
2009 h(8,2) = -45.+54.*x+81.*y;
2012 h(9,1) = 27.-54.*(x+y);
2085 double x = ip.
x, y = ip.
y, z = ip.
z;
2087 shape(0) = -((-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z)*
2088 (-1 + 3*x + 3*y + 3*z))/2.;
2089 shape(4) = (9*x*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
2090 shape(5) = (-9*x*(-1 + 3*x)*(-1 + x + y + z))/2.;
2091 shape(1) = (x*(2 + 9*(-1 + x)*x))/2.;
2092 shape(6) = (9*y*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
2093 shape(19) = -27*x*y*(-1 + x + y + z);
2094 shape(10) = (9*x*(-1 + 3*x)*y)/2.;
2095 shape(7) = (-9*y*(-1 + 3*y)*(-1 + x + y + z))/2.;
2096 shape(11) = (9*x*y*(-1 + 3*y))/2.;
2097 shape(2) = (y*(2 + 9*(-1 + y)*y))/2.;
2098 shape(8) = (9*z*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
2099 shape(18) = -27*x*z*(-1 + x + y + z);
2100 shape(12) = (9*x*(-1 + 3*x)*z)/2.;
2101 shape(17) = -27*y*z*(-1 + x + y + z);
2102 shape(16) = 27*x*y*z;
2103 shape(14) = (9*y*(-1 + 3*y)*z)/2.;
2104 shape(9) = (-9*z*(-1 + x + y + z)*(-1 + 3*z))/2.;
2105 shape(13) = (9*x*z*(-1 + 3*z))/2.;
2106 shape(15) = (9*y*z*(-1 + 3*z))/2.;
2107 shape(3) = (z*(2 + 9*(-1 + z)*z))/2.;
2113 double x = ip.
x, y = ip.
y, z = ip.
z;
2115 dshape(0,0) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
2116 x*(-4 + 6*y + 6*z)))/2.;
2117 dshape(0,1) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
2118 x*(-4 + 6*y + 6*z)))/2.;
2119 dshape(0,2) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
2120 x*(-4 + 6*y + 6*z)))/2.;
2121 dshape(4,0) = (9*(9*pow(x,2) + (-1 + y + z)*(-2 + 3*y + 3*z) +
2122 2*x*(-5 + 6*y + 6*z)))/2.;
2123 dshape(4,1) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
2124 dshape(4,2) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
2125 dshape(5,0) = (-9*(1 - y - z + x*(-8 + 9*x + 6*y + 6*z)))/2.;
2126 dshape(5,1) = (9*(1 - 3*x)*x)/2.;
2127 dshape(5,2) = (9*(1 - 3*x)*x)/2.;
2128 dshape(1,0) = 1 + (9*x*(-2 + 3*x))/2.;
2131 dshape(6,0) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
2132 dshape(6,1) = (9*(2 + 3*pow(x,2) - 10*y - 5*z + 3*(y + z)*(3*y + z) +
2133 x*(-5 + 12*y + 6*z)))/2.;
2134 dshape(6,2) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
2135 dshape(19,0) = -27*y*(-1 + 2*x + y + z);
2136 dshape(19,1) = -27*x*(-1 + x + 2*y + z);
2137 dshape(19,2) = -27*x*y;
2138 dshape(10,0) = (9*(-1 + 6*x)*y)/2.;
2139 dshape(10,1) = (9*x*(-1 + 3*x))/2.;
2141 dshape(7,0) = (9*(1 - 3*y)*y)/2.;
2142 dshape(7,1) = (-9*(1 + x*(-1 + 6*y) - z + y*(-8 + 9*y + 6*z)))/2.;
2143 dshape(7,2) = (9*(1 - 3*y)*y)/2.;
2144 dshape(11,0) = (9*y*(-1 + 3*y))/2.;
2145 dshape(11,1) = (9*x*(-1 + 6*y))/2.;
2148 dshape(2,1) = 1 + (9*y*(-2 + 3*y))/2.;
2150 dshape(8,0) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
2151 dshape(8,1) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
2152 dshape(8,2) = (9*(2 + 3*pow(x,2) - 5*y - 10*z + 3*(y + z)*(y + 3*z) +
2153 x*(-5 + 6*y + 12*z)))/2.;
2154 dshape(18,0) = -27*z*(-1 + 2*x + y + z);
2155 dshape(18,1) = -27*x*z;
2156 dshape(18,2) = -27*x*(-1 + x + y + 2*z);
2157 dshape(12,0) = (9*(-1 + 6*x)*z)/2.;
2159 dshape(12,2) = (9*x*(-1 + 3*x))/2.;
2160 dshape(17,0) = -27*y*z;
2161 dshape(17,1) = -27*z*(-1 + x + 2*y + z);
2162 dshape(17,2) = -27*y*(-1 + x + y + 2*z);
2163 dshape(16,0) = 27*y*z;
2164 dshape(16,1) = 27*x*z;
2165 dshape(16,2) = 27*x*y;
2167 dshape(14,1) = (9*(-1 + 6*y)*z)/2.;
2168 dshape(14,2) = (9*y*(-1 + 3*y))/2.;
2169 dshape(9,0) = (9*(1 - 3*z)*z)/2.;
2170 dshape(9,1) = (9*(1 - 3*z)*z)/2.;
2171 dshape(9,2) = (9*(-1 + x + y + 8*z - 6*(x + y)*z - 9*pow(z,2)))/2.;
2172 dshape(13,0) = (9*z*(-1 + 3*z))/2.;
2174 dshape(13,2) = (9*x*(-1 + 6*z))/2.;
2176 dshape(15,1) = (9*z*(-1 + 3*z))/2.;
2177 dshape(15,2) = (9*y*(-1 + 6*z))/2.;
2180 dshape(3,2) = 1 + (9*z*(-2 + 3*z))/2.;
2246 shape(0) = 1. - ip.
x - ip.
y - ip.
z;
2255 if (dshape.
Height() == 4)
2257 double *A = &dshape(0,0);
2258 A[0] = -1.; A[4] = -1.; A[8] = -1.;
2259 A[1] = 1.; A[5] = 0.; A[9] = 0.;
2260 A[2] = 0.; A[6] = 1.; A[10] = 0.;
2261 A[3] = 0.; A[7] = 0.; A[11] = 1.;
2265 dshape(0,0) = -1.; dshape(0,1) = -1.; dshape(0,2) = -1.;
2266 dshape(1,0) = 1.; dshape(1,1) = 0.; dshape(1,2) = 0.;
2267 dshape(2,0) = 0.; dshape(2,1) = 1.; dshape(2,2) = 0.;
2268 dshape(3,0) = 0.; dshape(3,1) = 0.; dshape(3,2) = 1.;
2275 static int face_dofs[4][3] = {{1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {0, 1, 2}};
2278 *dofs = face_dofs[face];
2320 double L0, L1, L2, L3;
2322 L0 = 1. - ip.
x - ip.
y - ip.
z;
2327 shape(0) = L0 * ( 2.0 * L0 - 1.0 );
2328 shape(1) = L1 * ( 2.0 * L1 - 1.0 );
2329 shape(2) = L2 * ( 2.0 * L2 - 1.0 );
2330 shape(3) = L3 * ( 2.0 * L3 - 1.0 );
2331 shape(4) = 4.0 * L0 * L1;
2332 shape(5) = 4.0 * L0 * L2;
2333 shape(6) = 4.0 * L0 * L3;
2334 shape(7) = 4.0 * L1 * L2;
2335 shape(8) = 4.0 * L1 * L3;
2336 shape(9) = 4.0 * L2 * L3;
2347 L0 = 1.0 - x - y - z;
2349 dshape(0,0) = dshape(0,1) = dshape(0,2) = 1.0 - 4.0 * L0;
2350 dshape(1,0) = -1.0 + 4.0 * x; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
2351 dshape(2,0) = 0.0; dshape(2,1) = -1.0 + 4.0 * y; dshape(2,2) = 0.0;
2352 dshape(3,0) = dshape(3,1) = 0.0; dshape(3,2) = -1.0 + 4.0 * z;
2353 dshape(4,0) = 4.0 * (L0 - x); dshape(4,1) = dshape(4,2) = -4.0 * x;
2354 dshape(5,0) = dshape(5,2) = -4.0 * y; dshape(5,1) = 4.0 * (L0 - y);
2355 dshape(6,0) = dshape(6,1) = -4.0 * z; dshape(6,2) = 4.0 * (L0 - z);
2356 dshape(7,0) = 4.0 * y; dshape(7,1) = 4.0 * x; dshape(7,2) = 0.0;
2357 dshape(8,0) = 4.0 * z; dshape(8,1) = 0.0; dshape(8,2) = 4.0 * x;
2358 dshape(9,0) = 0.0; dshape(9,1) = 4.0 * z; dshape(9,2) = 4.0 * y;
2400 double x = ip.
x, y = ip.
y, z = ip.
z;
2401 double ox = 1.-x, oy = 1.-y, oz = 1.-z;
2403 shape(0) = ox * oy * oz;
2404 shape(1) = x * oy * oz;
2405 shape(2) = x * y * oz;
2406 shape(3) = ox * y * oz;
2407 shape(4) = ox * oy * z;
2408 shape(5) = x * oy * z;
2409 shape(6) = x * y * z;
2410 shape(7) = ox * y * z;
2416 double x = ip.
x, y = ip.
y, z = ip.
z;
2417 double ox = 1.-x, oy = 1.-y, oz = 1.-z;
2419 dshape(0,0) = - oy * oz;
2420 dshape(0,1) = - ox * oz;
2421 dshape(0,2) = - ox * oy;
2423 dshape(1,0) = oy * oz;
2424 dshape(1,1) = - x * oz;
2425 dshape(1,2) = - x * oy;
2427 dshape(2,0) = y * oz;
2428 dshape(2,1) = x * oz;
2429 dshape(2,2) = - x * y;
2431 dshape(3,0) = - y * oz;
2432 dshape(3,1) = ox * oz;
2433 dshape(3,2) = - ox * y;
2435 dshape(4,0) = - oy * z;
2436 dshape(4,1) = - ox * z;
2437 dshape(4,2) = ox * oy;
2439 dshape(5,0) = oy * z;
2440 dshape(5,1) = - x * z;
2441 dshape(5,2) = x * oy;
2443 dshape(6,0) = y * z;
2444 dshape(6,1) = x * z;
2445 dshape(6,2) = x * y;
2447 dshape(7,0) = - y * z;
2448 dshape(7,1) = ox * z;
2449 dshape(7,2) = ox * y;
2484 shape(0) = 1.0 - 2.0 * ip.
y;
2485 shape(1) = -1.0 + 2.0 * ( ip.
x + ip.
y );
2486 shape(2) = 1.0 - 2.0 * ip.
x;
2492 dshape(0,0) = 0.0; dshape(0,1) = -2.0;
2493 dshape(1,0) = 2.0; dshape(1,1) = 2.0;
2494 dshape(2,0) = -2.0; dshape(2,1) = 0.0;
2515 const double l1 = ip.
x+ip.
y-0.5, l2 = 1.-l1, l3 = ip.
x-ip.
y+0.5, l4 = 1.-l3;
2526 const double x2 = 2.*ip.
x, y2 = 2.*ip.
y;
2528 dshape(0,0) = 1. - x2; dshape(0,1) = -2. + y2;
2529 dshape(1,0) = x2; dshape(1,1) = 1. - y2;
2530 dshape(2,0) = 1. - x2; dshape(2,1) = y2;
2531 dshape(3,0) = -2. + x2; dshape(3,1) = 1. - y2;
2549 double x = ip.
x, y = ip.
y;
2552 shape(0,1) = y - 1.;
2555 shape(2,0) = x - 1.;
2567 const double RT0TriangleFiniteElement::nk[3][2] =
2568 { {0, -1}, {1, 1}, {-1, 0} };
2574 #ifdef MFEM_THREAD_SAFE
2580 for (k = 0; k < 3; k++)
2583 for (j = 0; j < 3; j++)
2586 if (j == k) { d -= 1.0; }
2587 if (fabs(d) > 1.0e-12)
2589 cerr <<
"RT0TriangleFiniteElement::GetLocalInterpolation (...)\n"
2590 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2606 for (k = 0; k < 3; k++)
2609 ip.
x = vk[0]; ip.
y = vk[1];
2612 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2613 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2614 for (j = 0; j < 3; j++)
2615 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2628 #ifdef MFEM_THREAD_SAFE
2632 for (
int k = 0; k < 3; k++)
2640 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2641 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
2661 double x = ip.
x, y = ip.
y;
2664 shape(0,1) = y - 1.;
2669 shape(3,0) = x - 1.;
2682 const double RT0QuadFiniteElement::nk[4][2] =
2683 { {0, -1}, {1, 0}, {0, 1}, {-1, 0} };
2689 #ifdef MFEM_THREAD_SAFE
2695 for (k = 0; k < 4; k++)
2698 for (j = 0; j < 4; j++)
2701 if (j == k) { d -= 1.0; }
2702 if (fabs(d) > 1.0e-12)
2704 cerr <<
"RT0QuadFiniteElement::GetLocalInterpolation (...)\n"
2705 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2721 for (k = 0; k < 4; k++)
2724 ip.
x = vk[0]; ip.
y = vk[1];
2727 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2728 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2729 for (j = 0; j < 4; j++)
2730 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2743 #ifdef MFEM_THREAD_SAFE
2747 for (
int k = 0; k < 4; k++)
2755 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2756 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
2784 double x = ip.
x, y = ip.
y;
2786 shape(0,0) = -2 * x * (-1 + x + 2 * y);
2787 shape(0,1) = -2 * (-1 + y) * (-1 + x + 2 * y);
2788 shape(1,0) = 2 * x * (x - y);
2789 shape(1,1) = 2 * (x - y) * (-1 + y);
2790 shape(2,0) = 2 * x * (-1 + 2 * x + y);
2791 shape(2,1) = 2 * y * (-1 + 2 * x + y);
2792 shape(3,0) = 2 * x * (-1 + x + 2 * y);
2793 shape(3,1) = 2 * y * (-1 + x + 2 * y);
2794 shape(4,0) = -2 * (-1 + x) * (x - y);
2795 shape(4,1) = 2 * y * (-x + y);
2796 shape(5,0) = -2 * (-1 + x) * (-1 + 2 * x + y);
2797 shape(5,1) = -2 * y * (-1 + 2 * x + y);
2798 shape(6,0) = -3 * x * (-2 + 2 * x + y);
2799 shape(6,1) = -3 * y * (-1 + 2 * x + y);
2800 shape(7,0) = -3 * x * (-1 + x + 2 * y);
2801 shape(7,1) = -3 * y * (-2 + x + 2 * y);
2807 double x = ip.
x, y = ip.
y;
2809 divshape(0) = -2 * (-4 + 3 * x + 6 * y);
2810 divshape(1) = 2 + 6 * x - 6 * y;
2811 divshape(2) = -4 + 12 * x + 6 * y;
2812 divshape(3) = -4 + 6 * x + 12 * y;
2813 divshape(4) = 2 - 6 * x + 6 * y;
2814 divshape(5) = -2 * (-4 + 6 * x + 3 * y);
2815 divshape(6) = -9 * (-1 + 2 * x + y);
2816 divshape(7) = -9 * (-1 + x + 2 * y);
2819 const double RT1TriangleFiniteElement::nk[8][2] =
2831 #ifdef MFEM_THREAD_SAFE
2837 for (k = 0; k < 8; k++)
2840 for (j = 0; j < 8; j++)
2843 if (j == k) { d -= 1.0; }
2844 if (fabs(d) > 1.0e-12)
2846 cerr <<
"RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
2847 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
2863 for (k = 0; k < 8; k++)
2866 ip.
x = vk[0]; ip.
y = vk[1];
2869 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
2870 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
2871 for (j = 0; j < 8; j++)
2872 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
2884 #ifdef MFEM_THREAD_SAFE
2888 for (
int k = 0; k < 8; k++)
2896 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
2897 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
2940 double x = ip.
x, y = ip.
y;
2944 shape(0,1) = -( 1. - 3.*y + 2.*y*y)*( 2. - 3.*x);
2946 shape(1,1) = -( 1. - 3.*y + 2.*y*y)*(-1. + 3.*x);
2948 shape(2,0) = (-x + 2.*x*x)*( 2. - 3.*y);
2950 shape(3,0) = (-x + 2.*x*x)*(-1. + 3.*y);
2954 shape(4,1) = (-y + 2.*y*y)*(-1. + 3.*x);
2956 shape(5,1) = (-y + 2.*y*y)*( 2. - 3.*x);
2958 shape(6,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y);
2960 shape(7,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y);
2963 shape(8,0) = (4.*x - 4.*x*x)*( 2. - 3.*y);
2965 shape(9,0) = (4.*x - 4.*x*x)*(-1. + 3.*y);
2969 shape(10,1) = (4.*y - 4.*y*y)*( 2. - 3.*x);
2971 shape(11,1) = (4.*y - 4.*y*y)*(-1. + 3.*x);
2977 double x = ip.
x, y = ip.
y;
2979 divshape(0) = -(-3. + 4.*y)*( 2. - 3.*x);
2980 divshape(1) = -(-3. + 4.*y)*(-1. + 3.*x);
2981 divshape(2) = (-1. + 4.*x)*( 2. - 3.*y);
2982 divshape(3) = (-1. + 4.*x)*(-1. + 3.*y);
2983 divshape(4) = (-1. + 4.*y)*(-1. + 3.*x);
2984 divshape(5) = (-1. + 4.*y)*( 2. - 3.*x);
2985 divshape(6) = -(-3. + 4.*x)*(-1. + 3.*y);
2986 divshape(7) = -(-3. + 4.*x)*( 2. - 3.*y);
2987 divshape(8) = ( 4. - 8.*x)*( 2. - 3.*y);
2988 divshape(9) = ( 4. - 8.*x)*(-1. + 3.*y);
2989 divshape(10) = ( 4. - 8.*y)*( 2. - 3.*x);
2990 divshape(11) = ( 4. - 8.*y)*(-1. + 3.*x);
2993 const double RT1QuadFiniteElement::nk[12][2] =
3013 #ifdef MFEM_THREAD_SAFE
3019 for (k = 0; k < 12; k++)
3022 for (j = 0; j < 12; j++)
3025 if (j == k) { d -= 1.0; }
3026 if (fabs(d) > 1.0e-12)
3028 cerr <<
"RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
3029 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
3045 for (k = 0; k < 12; k++)
3048 ip.
x = vk[0]; ip.
y = vk[1];
3051 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
3052 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
3053 for (j = 0; j < 12; j++)
3054 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
3066 #ifdef MFEM_THREAD_SAFE
3070 for (
int k = 0; k < 12; k++)
3078 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
3079 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
3083 const double RT2TriangleFiniteElement::M[15][15] =
3086 0, -5.3237900077244501311, 5.3237900077244501311, 16.647580015448900262,
3087 0, 24.442740046346700787, -16.647580015448900262, -12.,
3088 -19.118950038622250656, -47.237900077244501311, 0, -34.414110069520051180,
3089 12., 30.590320061795601049, 15.295160030897800524
3092 0, 1.5, -1.5, -15., 0, 2.625, 15., 15., -4.125, 30., 0, -14.625, -15.,
3096 0, -0.67620999227554986889, 0.67620999227554986889, 7.3524199845510997378,
3097 0, -3.4427400463467007866, -7.3524199845510997378, -12.,
3098 4.1189500386222506555, -0.76209992275549868892, 0, 7.4141100695200511800,
3099 12., -6.5903200617956010489, -3.2951600308978005244
3102 0, 0, 1.5, 0, 0, 1.5, -11.471370023173350393, 0, 2.4713700231733503933,
3103 -11.471370023173350393, 0, 2.4713700231733503933, 15.295160030897800524,
3104 0, -3.2951600308978005244
3107 0, 0, 4.875, 0, 0, 4.875, -16.875, 0, -16.875, -16.875, 0, -16.875, 10.5,
3111 0, 0, 1.5, 0, 0, 1.5, 2.4713700231733503933, 0, -11.471370023173350393,
3112 2.4713700231733503933, 0, -11.471370023173350393, -3.2951600308978005244,
3113 0, 15.295160030897800524
3116 -0.67620999227554986889, 0, -3.4427400463467007866, 0,
3117 7.3524199845510997378, 0.67620999227554986889, 7.4141100695200511800, 0,
3118 -0.76209992275549868892, 4.1189500386222506555, -12.,
3119 -7.3524199845510997378, -3.2951600308978005244, -6.5903200617956010489,
3123 1.5, 0, 2.625, 0, -15., -1.5, -14.625, 0, 30., -4.125, 15., 15., 10.5,
3127 -5.3237900077244501311, 0, 24.442740046346700787, 0, 16.647580015448900262,
3128 5.3237900077244501311, -34.414110069520051180, 0, -47.237900077244501311,
3129 -19.118950038622250656, -12., -16.647580015448900262, 15.295160030897800524,
3130 30.590320061795601049, 12.
3132 { 0, 0, 18., 0, 0, 6., -42., 0, -30., -26., 0, -14., 24., 32., 8.},
3133 { 0, 0, 6., 0, 0, 18., -14., 0, -26., -30., 0, -42., 8., 32., 24.},
3134 { 0, 0, -6., 0, 0, -4., 30., 0, 4., 22., 0, 4., -24., -16., 0},
3135 { 0, 0, -4., 0, 0, -8., 20., 0, 8., 36., 0, 8., -16., -32., 0},
3136 { 0, 0, -8., 0, 0, -4., 8., 0, 36., 8., 0, 20., 0, -32., -16.},
3137 { 0, 0, -4., 0, 0, -6., 4., 0, 22., 4., 0, 30., 0, -16., -24.}
3143 const double p = 0.11270166537925831148;
3180 double x = ip.
x, y = ip.
y;
3182 double Bx[15] = {1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y, 0., x*x*x,
3185 double By[15] = {0., 1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y,
3189 for (
int i = 0; i < 15; i++)
3191 double cx = 0.0, cy = 0.0;
3192 for (
int j = 0; j < 15; j++)
3194 cx += M[i][j] * Bx[j];
3195 cy += M[i][j] * By[j];
3205 double x = ip.
x, y = ip.
y;
3207 double DivB[15] = {0., 0., 1., 0., 0., 1., 2.*x, 0., y, x, 0., 2.*y,
3208 4.*x*x, 4.*x*y, 4.*y*y
3211 for (
int i = 0; i < 15; i++)
3214 for (
int j = 0; j < 15; j++)
3216 div += M[i][j] * DivB[j];
3222 const double RT2QuadFiniteElement::pt[4] = {0.,1./3.,2./3.,1.};
3224 const double RT2QuadFiniteElement::dpt[3] = {0.25,0.5,0.75};
3266 double x = ip.
x, y = ip.
y;
3268 double ax0 = pt[0] - x;
3269 double ax1 = pt[1] - x;
3270 double ax2 = pt[2] - x;
3271 double ax3 = pt[3] - x;
3273 double by0 = dpt[0] - y;
3274 double by1 = dpt[1] - y;
3275 double by2 = dpt[2] - y;
3277 double ay0 = pt[0] - y;
3278 double ay1 = pt[1] - y;
3279 double ay2 = pt[2] - y;
3280 double ay3 = pt[3] - y;
3282 double bx0 = dpt[0] - x;
3283 double bx1 = dpt[1] - x;
3284 double bx2 = dpt[2] - x;
3286 double A01 = pt[0] - pt[1];
3287 double A02 = pt[0] - pt[2];
3288 double A12 = pt[1] - pt[2];
3289 double A03 = pt[0] - pt[3];
3290 double A13 = pt[1] - pt[3];
3291 double A23 = pt[2] - pt[3];
3293 double B01 = dpt[0] - dpt[1];
3294 double B02 = dpt[0] - dpt[2];
3295 double B12 = dpt[1] - dpt[2];
3297 double tx0 = (bx1*bx2)/(B01*B02);
3298 double tx1 = -(bx0*bx2)/(B01*B12);
3299 double tx2 = (bx0*bx1)/(B02*B12);
3301 double ty0 = (by1*by2)/(B01*B02);
3302 double ty1 = -(by0*by2)/(B01*B12);
3303 double ty2 = (by0*by1)/(B02*B12);
3307 shape(0, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx0;
3309 shape(1, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx1;
3311 shape(2, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx2;
3313 shape(3, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty0;
3315 shape(4, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty1;
3317 shape(5, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty2;
3321 shape(6, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx2;
3323 shape(7, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx1;
3325 shape(8, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx0;
3327 shape(9, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty2;
3329 shape(10, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty1;
3331 shape(11, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty0;
3334 shape(12, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty0;
3336 shape(13, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty1;
3338 shape(14, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty2;
3341 shape(15, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty0;
3343 shape(16, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty1;
3345 shape(17, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty2;
3349 shape(18, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx0;
3351 shape(19, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx1;
3353 shape(20, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx2;
3356 shape(21, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx0;
3358 shape(22, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx1;
3360 shape(23, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx2;
3366 double x = ip.
x, y = ip.
y;
3368 double a01 = pt[0]*pt[1];
3369 double a02 = pt[0]*pt[2];
3370 double a12 = pt[1]*pt[2];
3371 double a03 = pt[0]*pt[3];
3372 double a13 = pt[1]*pt[3];
3373 double a23 = pt[2]*pt[3];
3375 double bx0 = dpt[0] - x;
3376 double bx1 = dpt[1] - x;
3377 double bx2 = dpt[2] - x;
3379 double by0 = dpt[0] - y;
3380 double by1 = dpt[1] - y;
3381 double by2 = dpt[2] - y;
3383 double A01 = pt[0] - pt[1];
3384 double A02 = pt[0] - pt[2];
3385 double A12 = pt[1] - pt[2];
3386 double A03 = pt[0] - pt[3];
3387 double A13 = pt[1] - pt[3];
3388 double A23 = pt[2] - pt[3];
3390 double A012 = pt[0] + pt[1] + pt[2];
3391 double A013 = pt[0] + pt[1] + pt[3];
3392 double A023 = pt[0] + pt[2] + pt[3];
3393 double A123 = pt[1] + pt[2] + pt[3];
3395 double B01 = dpt[0] - dpt[1];
3396 double B02 = dpt[0] - dpt[2];
3397 double B12 = dpt[1] - dpt[2];
3399 double tx0 = (bx1*bx2)/(B01*B02);
3400 double tx1 = -(bx0*bx2)/(B01*B12);
3401 double tx2 = (bx0*bx1)/(B02*B12);
3403 double ty0 = (by1*by2)/(B01*B02);
3404 double ty1 = -(by0*by2)/(B01*B12);
3405 double ty2 = (by0*by1)/(B02*B12);
3408 divshape(0) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx0;
3409 divshape(1) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx1;
3410 divshape(2) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx2;
3412 divshape(3) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty0;
3413 divshape(4) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty1;
3414 divshape(5) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty2;
3416 divshape(6) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx2;
3417 divshape(7) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx1;
3418 divshape(8) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx0;
3420 divshape(9) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty2;
3421 divshape(10) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty1;
3422 divshape(11) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty0;
3424 divshape(12) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty0;
3425 divshape(13) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty1;
3426 divshape(14) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty2;
3428 divshape(15) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty0;
3429 divshape(16) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty1;
3430 divshape(17) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty2;
3432 divshape(18) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx0;
3433 divshape(19) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx1;
3434 divshape(20) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx2;
3436 divshape(21) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx0;
3437 divshape(22) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx1;
3438 divshape(23) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx2;
3441 const double RT2QuadFiniteElement::nk[24][2] =
3444 {0,-1}, {0,-1}, {0,-1},
3446 {1, 0}, {1, 0}, {1, 0},
3448 {0, 1}, {0, 1}, {0, 1},
3450 {-1,0}, {-1,0}, {-1,0},
3452 {1, 0}, {1, 0}, {1, 0},
3454 {1, 0}, {1, 0}, {1, 0},
3456 {0, 1}, {0, 1}, {0, 1},
3458 {0, 1}, {0, 1}, {0, 1}
3465 #ifdef MFEM_THREAD_SAFE
3471 for (k = 0; k < 24; k++)
3474 for (j = 0; j < 24; j++)
3477 if (j == k) { d -= 1.0; }
3478 if (fabs(d) > 1.0e-12)
3480 cerr <<
"RT2QuadFiniteElement::GetLocalInterpolation (...)\n"
3481 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
3497 for (k = 0; k < 24; k++)
3500 ip.
x = vk[0]; ip.
y = vk[1];
3503 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1];
3504 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1];
3505 for (j = 0; j < 24; j++)
3506 if (fabs (I(k,j) =
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]) < 1.0e-12)
3518 #ifdef MFEM_THREAD_SAFE
3522 for (
int k = 0; k < 24; k++)
3530 dofs(k) = (vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1] ) +
3531 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1] ));
3547 shape(0) = 2. - 3. * x;
3548 shape(1) = 3. * x - 1.;
3562 const double p = 0.11270166537925831148;
3572 const double p = 0.11270166537925831148;
3573 const double w = 1./((1-2*p)*(1-2*p));
3576 shape(0) = (2*x-1)*(x-1+p)*w;
3577 shape(1) = 4*(x-1+p)*(p-x)*w;
3578 shape(2) = (2*x-1)*(x-p)*w;
3584 const double p = 0.11270166537925831148;
3585 const double w = 1./((1-2*p)*(1-2*p));
3588 dshape(0,0) = (-3+4*x+2*p)*w;
3589 dshape(1,0) = (4-8*x)*w;
3590 dshape(2,0) = (-1+4*x-2*p)*w;
3601 for (i = 1; i < m; i++)
3607 #ifndef MFEM_THREAD_SAFE
3612 for (i = 1; i <= m; i++)
3614 rwk(i) = rwk(i-1) * ( (double)(m) / (double)(i) );
3616 for (i = 0; i < m/2+1; i++)
3618 rwk(m-i) = ( rwk(i) *= rwk(m-i) );
3620 for (i = m-1; i >= 0; i -= 2)
3629 double w, wk, x = ip.
x;
3632 #ifdef MFEM_THREAD_SAFE
3636 k = (int) floor ( m * x + 0.5 );
3638 for (i = 0; i <= m; i++)
3641 wk *= ( rxxk(i) = x - (double)(i) / m );
3643 w = wk * ( rxxk(k) = x - (double)(k) / m );
3647 shape(0) = w * rwk(0) / rxxk(0);
3651 shape(0) = wk * rwk(0);
3655 shape(1) = w * rwk(m) / rxxk(m);
3659 shape(1) = wk * rwk(k);
3661 for (i = 1; i < m; i++)
3664 shape(i+1) = w * rwk(i) / rxxk(i);
3668 shape(k+1) = wk * rwk(k);
3675 double s, srx, w, wk, x = ip.
x;
3678 #ifdef MFEM_THREAD_SAFE
3682 k = (int) floor ( m * x + 0.5 );
3684 for (i = 0; i <= m; i++)
3687 wk *= ( rxxk(i) = x - (double)(i) / m );
3689 w = wk * ( rxxk(k) = x - (double)(k) / m );
3691 for (i = 0; i <= m; i++)
3693 rxxk(i) = 1.0 / rxxk(i);
3696 for (i = 0; i <= m; i++)
3705 dshape(0,0) = (s - w * rxxk(0)) * rwk(0) * rxxk(0);
3709 dshape(0,0) = wk * srx * rwk(0);
3713 dshape(1,0) = (s - w * rxxk(m)) * rwk(m) * rxxk(m);
3717 dshape(1,0) = wk * srx * rwk(k);
3719 for (i = 1; i < m; i++)
3722 dshape(i+1,0) = (s - w * rxxk(i)) * rwk(i) * rxxk(i);
3726 dshape(k+1,0) = wk * srx * rwk(k);
3755 double L0, L1, L2, L3;
3757 L1 = ip.
x; L2 = ip.
y; L3 = ip.
z; L0 = 1.0 - L1 - L2 - L3;
3758 shape(0) = 1.0 - 3.0 * L0;
3759 shape(1) = 1.0 - 3.0 * L1;
3760 shape(2) = 1.0 - 3.0 * L2;
3761 shape(3) = 1.0 - 3.0 * L3;
3767 dshape(0,0) = 3.0; dshape(0,1) = 3.0; dshape(0,2) = 3.0;
3768 dshape(1,0) = -3.0; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
3769 dshape(2,0) = 0.0; dshape(2,1) = -3.0; dshape(2,2) = 0.0;
3770 dshape(3,0) = 0.0; dshape(3,1) = 0.0; dshape(3,2) = -3.0;
3791 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3812 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3826 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3827 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3828 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3829 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3830 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3831 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3832 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3833 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3835 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3836 I[ 9] = 1; J[ 9] = 2; K[ 9] = 0;
3837 I[10] = 2; J[10] = 1; K[10] = 0;
3838 I[11] = 0; J[11] = 2; K[11] = 0;
3839 I[12] = 2; J[12] = 0; K[12] = 1;
3840 I[13] = 1; J[13] = 2; K[13] = 1;
3841 I[14] = 2; J[14] = 1; K[14] = 1;
3842 I[15] = 0; J[15] = 2; K[15] = 1;
3843 I[16] = 0; J[16] = 0; K[16] = 2;
3844 I[17] = 1; J[17] = 0; K[17] = 2;
3845 I[18] = 1; J[18] = 1; K[18] = 2;
3846 I[19] = 0; J[19] = 1; K[19] = 2;
3848 I[20] = 2; J[20] = 2; K[20] = 0;
3849 I[21] = 2; J[21] = 0; K[21] = 2;
3850 I[22] = 1; J[22] = 2; K[22] = 2;
3851 I[23] = 2; J[23] = 1; K[23] = 2;
3852 I[24] = 0; J[24] = 2; K[24] = 2;
3853 I[25] = 2; J[25] = 2; K[25] = 1;
3855 I[26] = 2; J[26] = 2; K[26] = 2;
3857 else if (degree == 3)
3863 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3864 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3865 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3866 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3867 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3868 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3869 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3870 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3872 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3873 I[ 9] = 3; J[ 9] = 0; K[ 9] = 0;
3874 I[10] = 1; J[10] = 2; K[10] = 0;
3875 I[11] = 1; J[11] = 3; K[11] = 0;
3876 I[12] = 2; J[12] = 1; K[12] = 0;
3877 I[13] = 3; J[13] = 1; K[13] = 0;
3878 I[14] = 0; J[14] = 2; K[14] = 0;
3879 I[15] = 0; J[15] = 3; K[15] = 0;
3880 I[16] = 2; J[16] = 0; K[16] = 1;
3881 I[17] = 3; J[17] = 0; K[17] = 1;
3882 I[18] = 1; J[18] = 2; K[18] = 1;
3883 I[19] = 1; J[19] = 3; K[19] = 1;
3884 I[20] = 2; J[20] = 1; K[20] = 1;
3885 I[21] = 3; J[21] = 1; K[21] = 1;
3886 I[22] = 0; J[22] = 2; K[22] = 1;
3887 I[23] = 0; J[23] = 3; K[23] = 1;
3888 I[24] = 0; J[24] = 0; K[24] = 2;
3889 I[25] = 0; J[25] = 0; K[25] = 3;
3890 I[26] = 1; J[26] = 0; K[26] = 2;
3891 I[27] = 1; J[27] = 0; K[27] = 3;
3892 I[28] = 1; J[28] = 1; K[28] = 2;
3893 I[29] = 1; J[29] = 1; K[29] = 3;
3894 I[30] = 0; J[30] = 1; K[30] = 2;
3895 I[31] = 0; J[31] = 1; K[31] = 3;
3897 I[32] = 2; J[32] = 3; K[32] = 0;
3898 I[33] = 3; J[33] = 3; K[33] = 0;
3899 I[34] = 2; J[34] = 2; K[34] = 0;
3900 I[35] = 3; J[35] = 2; K[35] = 0;
3901 I[36] = 2; J[36] = 0; K[36] = 2;
3902 I[37] = 3; J[37] = 0; K[37] = 2;
3903 I[38] = 2; J[38] = 0; K[38] = 3;
3904 I[39] = 3; J[39] = 0; K[39] = 3;
3905 I[40] = 1; J[40] = 2; K[40] = 2;
3906 I[41] = 1; J[41] = 3; K[41] = 2;
3907 I[42] = 1; J[42] = 2; K[42] = 3;
3908 I[43] = 1; J[43] = 3; K[43] = 3;
3909 I[44] = 3; J[44] = 1; K[44] = 2;
3910 I[45] = 2; J[45] = 1; K[45] = 2;
3911 I[46] = 3; J[46] = 1; K[46] = 3;
3912 I[47] = 2; J[47] = 1; K[47] = 3;
3913 I[48] = 0; J[48] = 3; K[48] = 2;
3914 I[49] = 0; J[49] = 2; K[49] = 2;
3915 I[50] = 0; J[50] = 3; K[50] = 3;
3916 I[51] = 0; J[51] = 2; K[51] = 3;
3917 I[52] = 2; J[52] = 2; K[52] = 1;
3918 I[53] = 3; J[53] = 2; K[53] = 1;
3919 I[54] = 2; J[54] = 3; K[54] = 1;
3920 I[55] = 3; J[55] = 3; K[55] = 1;
3922 I[56] = 2; J[56] = 2; K[56] = 2;
3923 I[57] = 3; J[57] = 2; K[57] = 2;
3924 I[58] = 3; J[58] = 3; K[58] = 2;
3925 I[59] = 2; J[59] = 3; K[59] = 2;
3926 I[60] = 2; J[60] = 2; K[60] = 3;
3927 I[61] = 3; J[61] = 2; K[61] = 3;
3928 I[62] = 3; J[62] = 3; K[62] = 3;
3929 I[63] = 2; J[63] = 3; K[63] = 3;
3933 mfem_error (
"LagrangeHexFiniteElement::LagrangeHexFiniteElement");
3937 dof1d = fe1d ->
GetDof();
3939 #ifndef MFEM_THREAD_SAFE
3949 for (
int n = 0; n <
Dof; n++)
3964 #ifdef MFEM_THREAD_SAFE
3965 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
3972 for (
int n = 0; n <
Dof; n++)
3974 shape(n) = shape1dx(I[n]) * shape1dy(J[n]) * shape1dz(K[n]);
3985 #ifdef MFEM_THREAD_SAFE
3986 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
3987 DenseMatrix dshape1dx(dof1d,1), dshape1dy(dof1d,1), dshape1dz(dof1d,1);
3998 for (
int n = 0; n <
Dof; n++)
4000 dshape(n,0) = dshape1dx(I[n],0) * shape1dy(J[n]) * shape1dz(K[n]);
4001 dshape(n,1) = shape1dx(I[n]) * dshape1dy(J[n],0) * shape1dz(K[n]);
4002 dshape(n,2) = shape1dx(I[n]) * shape1dy(J[n]) * dshape1dz(K[n],0);
4031 shape(0) = 1.0 - 2.0 * x;
4038 shape(1) = 2.0 * x - 1.0;
4039 shape(2) = 2.0 - 2.0 * x;
4050 dshape(0,0) = - 2.0;
4058 dshape(2,0) = - 2.0;
4085 L0 = 2.0 * ( 1. - ip.
x - ip.
y );
4086 L1 = 2.0 * ( ip.
x );
4087 L2 = 2.0 * ( ip.
y );
4096 for (i = 0; i < 6; i++)
4103 shape(0) = L0 - 1.0;
4110 shape(1) = L1 - 1.0;
4117 shape(2) = L2 - 1.0;
4121 shape(3) = 1.0 - L2;
4122 shape(4) = 1.0 - L0;
4123 shape(5) = 1.0 - L1;
4133 L0 = 2.0 * ( 1. - ip.
x - ip.
y );
4134 L1 = 2.0 * ( ip.
x );
4135 L2 = 2.0 * ( ip.
y );
4137 double DL0[2], DL1[2], DL2[2];
4138 DL0[0] = -2.0; DL0[1] = -2.0;
4139 DL1[0] = 2.0; DL1[1] = 0.0;
4140 DL2[0] = 0.0; DL2[1] = 2.0;
4142 for (i = 0; i < 6; i++)
4143 for (j = 0; j < 2; j++)
4150 for (j = 0; j < 2; j++)
4152 dshape(0,j) = DL0[j];
4153 dshape(3,j) = DL1[j];
4154 dshape(5,j) = DL2[j];
4159 for (j = 0; j < 2; j++)
4161 dshape(3,j) = DL0[j];
4162 dshape(1,j) = DL1[j];
4163 dshape(4,j) = DL2[j];
4168 for (j = 0; j < 2; j++)
4170 dshape(5,j) = DL0[j];
4171 dshape(4,j) = DL1[j];
4172 dshape(2,j) = DL2[j];
4177 for (j = 0; j < 2; j++)
4179 dshape(3,j) = - DL2[j];
4180 dshape(4,j) = - DL0[j];
4181 dshape(5,j) = - DL1[j];
4226 double L0, L1, L2, L3, L4, L5;
4227 L0 = 2.0 * ( 1. - ip.
x - ip.
y - ip.
z );
4228 L1 = 2.0 * ( ip.
x );
4229 L2 = 2.0 * ( ip.
y );
4230 L3 = 2.0 * ( ip.
z );
4231 L4 = 2.0 * ( ip.
x + ip.
y );
4232 L5 = 2.0 * ( ip.
y + ip.
z );
4245 for (i = 0; i < 10; i++)
4252 shape(0) = L0 - 1.0;
4260 shape(1) = L1 - 1.0;
4268 shape(2) = L2 - 1.0;
4276 shape(3) = L3 - 1.0;
4278 else if ((L4 <= 1.0) && (L5 <= 1.0))
4280 shape(4) = 1.0 - L5;
4282 shape(6) = 1.0 - L4;
4283 shape(8) = 1.0 - L0;
4285 else if ((L4 >= 1.0) && (L5 <= 1.0))
4287 shape(4) = 1.0 - L5;
4288 shape(5) = 1.0 - L1;
4289 shape(7) = L4 - 1.0;
4292 else if ((L4 <= 1.0) && (L5 >= 1.0))
4294 shape(5) = 1.0 - L3;
4295 shape(6) = 1.0 - L4;
4297 shape(9) = L5 - 1.0;
4299 else if ((L4 >= 1.0) && (L5 >= 1.0))
4302 shape(7) = L4 - 1.0;
4303 shape(8) = 1.0 - L2;
4304 shape(9) = L5 - 1.0;
4313 double L0, L1, L2, L3, L4, L5;
4314 L0 = 2.0 * ( 1. - ip.
x - ip.
y - ip.
z );
4315 L1 = 2.0 * ( ip.
x );
4316 L2 = 2.0 * ( ip.
y );
4317 L3 = 2.0 * ( ip.
z );
4318 L4 = 2.0 * ( ip.
x + ip.
y );
4319 L5 = 2.0 * ( ip.
y + ip.
z );
4321 double DL0[3], DL1[3], DL2[3], DL3[3], DL4[3], DL5[3];
4322 DL0[0] = -2.0; DL0[1] = -2.0; DL0[2] = -2.0;
4323 DL1[0] = 2.0; DL1[1] = 0.0; DL1[2] = 0.0;
4324 DL2[0] = 0.0; DL2[1] = 2.0; DL2[2] = 0.0;
4325 DL3[0] = 0.0; DL3[1] = 0.0; DL3[2] = 2.0;
4326 DL4[0] = 2.0; DL4[1] = 2.0; DL4[2] = 0.0;
4327 DL5[0] = 0.0; DL5[1] = 2.0; DL5[2] = 2.0;
4329 for (i = 0; i < 10; i++)
4330 for (j = 0; j < 3; j++)
4337 for (j = 0; j < 3; j++)
4339 dshape(0,j) = DL0[j];
4340 dshape(4,j) = DL1[j];
4341 dshape(5,j) = DL2[j];
4342 dshape(6,j) = DL3[j];
4347 for (j = 0; j < 3; j++)
4349 dshape(4,j) = DL0[j];
4350 dshape(1,j) = DL1[j];
4351 dshape(7,j) = DL2[j];
4352 dshape(8,j) = DL3[j];
4357 for (j = 0; j < 3; j++)
4359 dshape(5,j) = DL0[j];
4360 dshape(7,j) = DL1[j];
4361 dshape(2,j) = DL2[j];
4362 dshape(9,j) = DL3[j];
4367 for (j = 0; j < 3; j++)
4369 dshape(6,j) = DL0[j];
4370 dshape(8,j) = DL1[j];
4371 dshape(9,j) = DL2[j];
4372 dshape(3,j) = DL3[j];
4375 else if ((L4 <= 1.0) && (L5 <= 1.0))
4377 for (j = 0; j < 3; j++)
4379 dshape(4,j) = - DL5[j];
4380 dshape(5,j) = DL2[j];
4381 dshape(6,j) = - DL4[j];
4382 dshape(8,j) = - DL0[j];
4385 else if ((L4 >= 1.0) && (L5 <= 1.0))
4387 for (j = 0; j < 3; j++)
4389 dshape(4,j) = - DL5[j];
4390 dshape(5,j) = - DL1[j];
4391 dshape(7,j) = DL4[j];
4392 dshape(8,j) = DL3[j];
4395 else if ((L4 <= 1.0) && (L5 >= 1.0))
4397 for (j = 0; j < 3; j++)
4399 dshape(5,j) = - DL3[j];
4400 dshape(6,j) = - DL4[j];
4401 dshape(8,j) = DL1[j];
4402 dshape(9,j) = DL5[j];
4405 else if ((L4 >= 1.0) && (L5 >= 1.0))
4407 for (j = 0; j < 3; j++)
4409 dshape(5,j) = DL0[j];
4410 dshape(7,j) = DL4[j];
4411 dshape(8,j) = - DL2[j];
4412 dshape(9,j) = DL5[j];
4445 double x = ip.
x, y = ip.
y;
4447 Lx = 2.0 * ( 1. - x );
4448 Ly = 2.0 * ( 1. - y );
4457 for (i = 0; i < 9; i++)
4462 if ((x <= 0.5) && (y <= 0.5))
4464 shape(0) = (Lx - 1.0) * (Ly - 1.0);
4465 shape(4) = (2.0 - Lx) * (Ly - 1.0);
4466 shape(8) = (2.0 - Lx) * (2.0 - Ly);
4467 shape(7) = (Lx - 1.0) * (2.0 - Ly);
4469 else if ((x >= 0.5) && (y <= 0.5))
4471 shape(4) = Lx * (Ly - 1.0);
4472 shape(1) = (1.0 - Lx) * (Ly - 1.0);
4473 shape(5) = (1.0 - Lx) * (2.0 - Ly);
4474 shape(8) = Lx * (2.0 - Ly);
4476 else if ((x >= 0.5) && (y >= 0.5))
4478 shape(8) = Lx * Ly ;
4479 shape(5) = (1.0 - Lx) * Ly ;
4480 shape(2) = (1.0 - Lx) * (1.0 - Ly);
4481 shape(6) = Lx * (1.0 - Ly);
4483 else if ((x <= 0.5) && (y >= 0.5))
4485 shape(7) = (Lx - 1.0) * Ly ;
4486 shape(8) = (2.0 - Lx) * Ly ;
4487 shape(6) = (2.0 - Lx) * (1.0 - Ly);
4488 shape(3) = (Lx - 1.0) * (1.0 - Ly);
4496 double x = ip.
x, y = ip.
y;
4498 Lx = 2.0 * ( 1. - x );
4499 Ly = 2.0 * ( 1. - y );
4501 for (i = 0; i < 9; i++)
4502 for (j = 0; j < 2; j++)
4507 if ((x <= 0.5) && (y <= 0.5))
4509 dshape(0,0) = 2.0 * (1.0 - Ly);
4510 dshape(0,1) = 2.0 * (1.0 - Lx);
4512 dshape(4,0) = 2.0 * (Ly - 1.0);
4513 dshape(4,1) = -2.0 * (2.0 - Lx);
4515 dshape(8,0) = 2.0 * (2.0 - Ly);
4516 dshape(8,1) = 2.0 * (2.0 - Lx);
4518 dshape(7,0) = -2.0 * (2.0 - Ly);
4519 dshape(7,0) = 2.0 * (Lx - 1.0);
4521 else if ((x >= 0.5) && (y <= 0.5))
4523 dshape(4,0) = -2.0 * (Ly - 1.0);
4524 dshape(4,1) = -2.0 * Lx;
4526 dshape(1,0) = 2.0 * (Ly - 1.0);
4527 dshape(1,1) = -2.0 * (1.0 - Lx);
4529 dshape(5,0) = 2.0 * (2.0 - Ly);
4530 dshape(5,1) = 2.0 * (1.0 - Lx);
4532 dshape(8,0) = -2.0 * (2.0 - Ly);
4533 dshape(8,1) = 2.0 * Lx;
4535 else if ((x >= 0.5) && (y >= 0.5))
4537 dshape(8,0) = -2.0 * Ly;
4538 dshape(8,1) = -2.0 * Lx;
4540 dshape(5,0) = 2.0 * Ly;
4541 dshape(5,1) = -2.0 * (1.0 - Lx);
4543 dshape(2,0) = 2.0 * (1.0 - Ly);
4544 dshape(2,1) = 2.0 * (1.0 - Lx);
4546 dshape(6,0) = -2.0 * (1.0 - Ly);
4547 dshape(6,1) = 2.0 * Lx;
4549 else if ((x <= 0.5) && (y >= 0.5))
4551 dshape(7,0) = -2.0 * Ly;
4552 dshape(7,1) = -2.0 * (Lx - 1.0);
4554 dshape(8,0) = 2.0 * Ly ;
4555 dshape(8,1) = -2.0 * (2.0 - Lx);
4557 dshape(6,0) = 2.0 * (1.0 - Ly);
4558 dshape(6,1) = 2.0 * (2.0 - Lx);
4560 dshape(3,0) = -2.0 * (1.0 - Ly);
4561 dshape(3,1) = 2.0 * (Lx - 1.0);
4572 I[ 0] = 0.0; J[ 0] = 0.0; K[ 0] = 0.0;
4573 I[ 1] = 1.0; J[ 1] = 0.0; K[ 1] = 0.0;
4574 I[ 2] = 1.0; J[ 2] = 1.0; K[ 2] = 0.0;
4575 I[ 3] = 0.0; J[ 3] = 1.0; K[ 3] = 0.0;
4576 I[ 4] = 0.0; J[ 4] = 0.0; K[ 4] = 1.0;
4577 I[ 5] = 1.0; J[ 5] = 0.0; K[ 5] = 1.0;
4578 I[ 6] = 1.0; J[ 6] = 1.0; K[ 6] = 1.0;
4579 I[ 7] = 0.0; J[ 7] = 1.0; K[ 7] = 1.0;
4581 I[ 8] = 0.5; J[ 8] = 0.0; K[ 8] = 0.0;
4582 I[ 9] = 1.0; J[ 9] = 0.5; K[ 9] = 0.0;
4583 I[10] = 0.5; J[10] = 1.0; K[10] = 0.0;
4584 I[11] = 0.0; J[11] = 0.5; K[11] = 0.0;
4585 I[12] = 0.5; J[12] = 0.0; K[12] = 1.0;
4586 I[13] = 1.0; J[13] = 0.5; K[13] = 1.0;
4587 I[14] = 0.5; J[14] = 1.0; K[14] = 1.0;
4588 I[15] = 0.0; J[15] = 0.5; K[15] = 1.0;
4589 I[16] = 0.0; J[16] = 0.0; K[16] = 0.5;
4590 I[17] = 1.0; J[17] = 0.0; K[17] = 0.5;
4591 I[18] = 1.0; J[18] = 1.0; K[18] = 0.5;
4592 I[19] = 0.0; J[19] = 1.0; K[19] = 0.5;
4594 I[20] = 0.5; J[20] = 0.5; K[20] = 0.0;
4595 I[21] = 0.5; J[21] = 0.0; K[21] = 0.5;
4596 I[22] = 1.0; J[22] = 0.5; K[22] = 0.5;
4597 I[23] = 0.5; J[23] = 1.0; K[23] = 0.5;
4598 I[24] = 0.0; J[24] = 0.5; K[24] = 0.5;
4599 I[25] = 0.5; J[25] = 0.5; K[25] = 1.0;
4601 I[26] = 0.5; J[26] = 0.5; K[26] = 0.5;
4603 for (
int n = 0; n < 27; n++)
4616 double x = ip.
x, y = ip.
y, z = ip.
z;
4618 for (i = 0; i < 27; i++)
4623 if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5))
4638 else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5))
4653 else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5))
4668 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5))
4683 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5))
4698 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5))
4713 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5))
4744 shape(N[0]) = Lx * Ly * Lz;
4745 shape(N[1]) = (1 - Lx) * Ly * Lz;
4746 shape(N[2]) = (1 - Lx) * (1 - Ly) * Lz;
4747 shape(N[3]) = Lx * (1 - Ly) * Lz;
4748 shape(N[4]) = Lx * Ly * (1 - Lz);
4749 shape(N[5]) = (1 - Lx) * Ly * (1 - Lz);
4750 shape(N[6]) = (1 - Lx) * (1 - Ly) * (1 - Lz);
4751 shape(N[7]) = Lx * (1 - Ly) * (1 - Lz);
4759 double x = ip.
x, y = ip.
y, z = ip.
z;
4761 for (i = 0; i < 27; i++)
4762 for (j = 0; j < 3; j++)
4767 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))
4812 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5))
4827 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5))
4842 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5))
4857 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5))
4888 dshape(N[0],0) = -2.0 * Ly * Lz ;
4889 dshape(N[0],1) = -2.0 * Lx * Lz ;
4890 dshape(N[0],2) = -2.0 * Lx * Ly ;
4892 dshape(N[1],0) = 2.0 * Ly * Lz ;
4893 dshape(N[1],1) = -2.0 * (1 - Lx) * Lz ;
4894 dshape(N[1],2) = -2.0 * (1 - Lx) * Ly ;
4896 dshape(N[2],0) = 2.0 * (1 - Ly) * Lz ;
4897 dshape(N[2],1) = 2.0 * (1 - Lx) * Lz ;
4898 dshape(N[2],2) = -2.0 * (1 - Lx) * (1 - Ly);
4900 dshape(N[3],0) = -2.0 * (1 - Ly) * Lz ;
4901 dshape(N[3],1) = 2.0 * Lx * Lz ;
4902 dshape(N[3],2) = -2.0 * Lx * (1 - Ly);
4904 dshape(N[4],0) = -2.0 * Ly * (1 - Lz);
4905 dshape(N[4],1) = -2.0 * Lx * (1 - Lz);
4906 dshape(N[4],2) = 2.0 * Lx * Ly ;
4908 dshape(N[5],0) = 2.0 * Ly * (1 - Lz);
4909 dshape(N[5],1) = -2.0 * (1 - Lx) * (1 - Lz);
4910 dshape(N[5],2) = 2.0 * (1 - Lx) * Ly ;
4912 dshape(N[6],0) = 2.0 * (1 - Ly) * (1 - Lz);
4913 dshape(N[6],1) = 2.0 * (1 - Lx) * (1 - Lz);
4914 dshape(N[6],2) = 2.0 * (1 - Lx) * (1 - Ly);
4916 dshape(N[7],0) = -2.0 * (1 - Ly) * (1 - Lz);
4917 dshape(N[7],1) = 2.0 * Lx * (1 - Lz);
4918 dshape(N[7],2) = 2.0 * Lx * (1 - Ly);
4978 double x = ip.
x, y = ip.
y, z = ip.
z;
4980 shape(0,0) = (1. - y) * (1. - z);
4984 shape(2,0) = y * (1. - z);
4988 shape(4,0) = z * (1. - y);
4997 shape(1,1) = x * (1. - z);
5001 shape(3,1) = (1. - x) * (1. - z);
5009 shape(7,1) = (1. - x) * z;
5014 shape(8,2) = (1. - x) * (1. - y);
5018 shape(9,2) = x * (1. - y);
5022 shape(10,2) = x * y;
5026 shape(11,2) = y * (1. - x);
5034 double x = ip.
x, y = ip.
y, z = ip.
z;
5036 curl_shape(0,0) = 0.;
5037 curl_shape(0,1) = y - 1.;
5038 curl_shape(0,2) = 1. - z;
5040 curl_shape(2,0) = 0.;
5041 curl_shape(2,1) = -y;
5042 curl_shape(2,2) = z - 1.;
5044 curl_shape(4,0) = 0;
5045 curl_shape(4,1) = 1. - y;
5046 curl_shape(4,2) = z;
5048 curl_shape(6,0) = 0.;
5049 curl_shape(6,1) = y;
5050 curl_shape(6,2) = -z;
5052 curl_shape(1,0) = x;
5053 curl_shape(1,1) = 0.;
5054 curl_shape(1,2) = 1. - z;
5056 curl_shape(3,0) = 1. - x;
5057 curl_shape(3,1) = 0.;
5058 curl_shape(3,2) = z - 1.;
5060 curl_shape(5,0) = -x;
5061 curl_shape(5,1) = 0.;
5062 curl_shape(5,2) = z;
5064 curl_shape(7,0) = x - 1.;
5065 curl_shape(7,1) = 0.;
5066 curl_shape(7,2) = -z;
5068 curl_shape(8,0) = x - 1.;
5069 curl_shape(8,1) = 1. - y;
5070 curl_shape(8,2) = 0.;
5072 curl_shape(9,0) = -x;
5073 curl_shape(9,1) = y - 1.;
5074 curl_shape(9,2) = 0;
5076 curl_shape(10,0) = x;
5077 curl_shape(10,1) = -y;
5078 curl_shape(10,2) = 0.;
5080 curl_shape(11,0) = 1. - x;
5081 curl_shape(11,1) = y;
5082 curl_shape(11,2) = 0.;
5085 const double Nedelec1HexFiniteElement::tk[12][3] =
5087 {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
5088 {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
5089 {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}
5096 #ifdef MFEM_THREAD_SAFE
5101 for (k = 0; k < 12; k++)
5104 for (j = 0; j < 12; j++)
5106 double d = (
vshape(j,0)*tk[k][0] +
vshape(j,1)*tk[k][1] +
5108 if (j == k) { d -= 1.0; }
5109 if (fabs(d) > 1.0e-12)
5111 cerr <<
"Nedelec1HexFiniteElement::GetLocalInterpolation (...)\n"
5112 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5120 ip.
x = ip.
y = ip.
z = 0.0;
5127 for (k = 0; k < 12; k++)
5130 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5133 vk[0] =
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2];
5134 vk[1] =
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2];
5135 vk[2] =
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2];
5136 for (j = 0; j < 12; j++)
5137 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5138 vshape(j,2)*vk[2])) < 1.0e-12)
5152 for (
int k = 0; k < 12; k++)
5160 vk[0] * (
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2] ) +
5161 vk[1] * (
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2] ) +
5162 vk[2] * (
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2] );
5199 double x = ip.
x, y = ip.
y, z = ip.
z;
5201 shape(0,0) = 1. - y - z;
5206 shape(1,1) = 1. - x - z;
5211 shape(2,2) = 1. - x - y;
5230 curl_shape(0,0) = 0.;
5231 curl_shape(0,1) = -2.;
5232 curl_shape(0,2) = 2.;
5234 curl_shape(1,0) = 2.;
5235 curl_shape(1,1) = 0.;
5236 curl_shape(1,2) = -2.;
5238 curl_shape(2,0) = -2.;
5239 curl_shape(2,1) = 2.;
5240 curl_shape(2,2) = 0.;
5242 curl_shape(3,0) = 0.;
5243 curl_shape(3,1) = 0.;
5244 curl_shape(3,2) = 2.;
5246 curl_shape(4,0) = 0.;
5247 curl_shape(4,1) = -2.;
5248 curl_shape(4,2) = 0.;
5250 curl_shape(5,0) = 2.;
5251 curl_shape(5,1) = 0.;
5252 curl_shape(5,2) = 0.;
5255 const double Nedelec1TetFiniteElement::tk[6][3] =
5256 {{1,0,0}, {0,1,0}, {0,0,1}, {-1,1,0}, {-1,0,1}, {0,-1,1}};
5262 #ifdef MFEM_THREAD_SAFE
5267 for (k = 0; k < 6; k++)
5270 for (j = 0; j < 6; j++)
5272 double d = (
vshape(j,0)*tk[k][0] +
vshape(j,1)*tk[k][1] +
5274 if (j == k) { d -= 1.0; }
5275 if (fabs(d) > 1.0e-12)
5277 cerr <<
"Nedelec1TetFiniteElement::GetLocalInterpolation (...)\n"
5278 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5286 ip.
x = ip.
y = ip.
z = 0.0;
5293 for (k = 0; k < 6; k++)
5296 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5299 vk[0] =
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2];
5300 vk[1] =
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2];
5301 vk[2] =
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2];
5302 for (j = 0; j < 6; j++)
5303 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5304 vshape(j,2)*vk[2])) < 1.0e-12)
5318 for (
int k = 0; k < 6; k++)
5326 vk[0] * (
J(0,0)*tk[k][0]+
J(0,1)*tk[k][1]+
J(0,2)*tk[k][2] ) +
5327 vk[1] * (
J(1,0)*tk[k][0]+
J(1,1)*tk[k][1]+
J(1,2)*tk[k][2] ) +
5328 vk[2] * (
J(2,0)*tk[k][0]+
J(2,1)*tk[k][1]+
J(2,2)*tk[k][2] );
5365 double x = ip.
x, y = ip.
y, z = ip.
z;
5369 shape(0,2) = z - 1.;
5372 shape(1,1) = y - 1.;
5383 shape(4,0) = x - 1.;
5403 const double RT0HexFiniteElement::nk[6][3] =
5404 {{0,0,-1}, {0,-1,0}, {1,0,0}, {0,1,0}, {-1,0,0}, {0,0,1}};
5410 #ifdef MFEM_THREAD_SAFE
5416 for (k = 0; k < 6; k++)
5419 for (j = 0; j < 6; j++)
5421 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5423 if (j == k) { d -= 1.0; }
5424 if (fabs(d) > 1.0e-12)
5426 cerr <<
"RT0HexFiniteElement::GetLocalInterpolation (...)\n"
5427 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5435 ip.
x = ip.
y = ip.
z = 0.0;
5443 for (k = 0; k < 6; k++)
5446 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5449 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2];
5450 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2];
5451 vk[2] =
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2];
5452 for (j = 0; j < 6; j++)
5453 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5454 vshape(j,2)*vk[2])) < 1.0e-12)
5467 #ifdef MFEM_THREAD_SAFE
5471 for (
int k = 0; k < 6; k++)
5480 vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2] ) +
5481 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2] ) +
5482 vk[2] * (
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2] );
5611 double x = ip.
x, y = ip.
y, z = ip.
z;
5615 shape(2,2) = -(1. - 3.*z + 2.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5618 shape(3,2) = -(1. - 3.*z + 2.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5621 shape(0,2) = -(1. - 3.*z + 2.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5624 shape(1,2) = -(1. - 3.*z + 2.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5627 shape(4,1) = -(1. - 3.*y + 2.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5630 shape(5,1) = -(1. - 3.*y + 2.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5633 shape(6,1) = -(1. - 3.*y + 2.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5636 shape(7,1) = -(1. - 3.*y + 2.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5639 shape(8,0) = (-x + 2.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5642 shape(9,0) = (-x + 2.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5645 shape(10,0) = (-x + 2.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5648 shape(11,0) = (-x + 2.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5653 shape(13,1) = (-y + 2.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5656 shape(12,1) = (-y + 2.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5659 shape(15,1) = (-y + 2.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5662 shape(14,1) = (-y + 2.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5665 shape(17,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5668 shape(16,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5671 shape(19,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5674 shape(18,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5680 shape(20,2) = (-z + 2.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5683 shape(21,2) = (-z + 2.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5686 shape(22,2) = (-z + 2.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5689 shape(23,2) = (-z + 2.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5691 shape(24,0) = (4.*x - 4.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5694 shape(25,0) = (4.*x - 4.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5697 shape(26,0) = (4.*x - 4.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5700 shape(27,0) = (4.*x - 4.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5705 shape(28,1) = (4.*y - 4.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5708 shape(29,1) = (4.*y - 4.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5711 shape(30,1) = (4.*y - 4.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5714 shape(31,1) = (4.*y - 4.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5719 shape(32,2) = (4.*z - 4.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5722 shape(33,2) = (4.*z - 4.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5725 shape(34,2) = (4.*z - 4.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5728 shape(35,2) = (4.*z - 4.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5734 double x = ip.
x, y = ip.
y, z = ip.
z;
5736 divshape(2) = -(-3. + 4.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5737 divshape(3) = -(-3. + 4.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5738 divshape(0) = -(-3. + 4.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5739 divshape(1) = -(-3. + 4.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5741 divshape(4) = -(-3. + 4.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5742 divshape(5) = -(-3. + 4.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5743 divshape(6) = -(-3. + 4.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5744 divshape(7) = -(-3. + 4.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5746 divshape(8) = (-1. + 4.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5747 divshape(9) = (-1. + 4.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5748 divshape(10) = (-1. + 4.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5749 divshape(11) = (-1. + 4.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5751 divshape(13) = (-1. + 4.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5752 divshape(12) = (-1. + 4.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5753 divshape(15) = (-1. + 4.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5754 divshape(14) = (-1. + 4.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5756 divshape(17) = -(-3. + 4.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5757 divshape(16) = -(-3. + 4.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5758 divshape(19) = -(-3. + 4.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5759 divshape(18) = -(-3. + 4.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5761 divshape(20) = (-1. + 4.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5762 divshape(21) = (-1. + 4.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5763 divshape(22) = (-1. + 4.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5764 divshape(23) = (-1. + 4.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5766 divshape(24) = ( 4. - 8.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5767 divshape(25) = ( 4. - 8.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5768 divshape(26) = ( 4. - 8.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5769 divshape(27) = ( 4. - 8.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5771 divshape(28) = ( 4. - 8.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5772 divshape(29) = ( 4. - 8.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5773 divshape(30) = ( 4. - 8.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5774 divshape(31) = ( 4. - 8.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5776 divshape(32) = ( 4. - 8.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5777 divshape(33) = ( 4. - 8.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5778 divshape(34) = ( 4. - 8.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5779 divshape(35) = ( 4. - 8.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5782 const double RT1HexFiniteElement::nk[36][3] =
5784 {0, 0,-1}, {0, 0,-1}, {0, 0,-1}, {0, 0,-1},
5785 {0,-1, 0}, {0,-1, 0}, {0,-1, 0}, {0,-1, 0},
5786 {1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0},
5787 {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0},
5788 {-1,0, 0}, {-1,0, 0}, {-1,0, 0}, {-1,0, 0},
5789 {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1},
5790 {1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0},
5791 {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0},
5792 {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1}
5799 #ifdef MFEM_THREAD_SAFE
5805 for (k = 0; k < 36; k++)
5808 for (j = 0; j < 36; j++)
5810 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5812 if (j == k) { d -= 1.0; }
5813 if (fabs(d) > 1.0e-12)
5815 cerr <<
"RT0HexFiniteElement::GetLocalInterpolation (...)\n"
5816 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5824 ip.
x = ip.
y = ip.
z = 0.0;
5832 for (k = 0; k < 36; k++)
5835 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5838 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2];
5839 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2];
5840 vk[2] =
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2];
5841 for (j = 0; j < 36; j++)
5842 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5843 vshape(j,2)*vk[2])) < 1.0e-12)
5856 #ifdef MFEM_THREAD_SAFE
5860 for (
int k = 0; k < 36; k++)
5869 vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2] ) +
5870 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2] ) +
5871 vk[2] * (
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2] );
5899 double x2 = 2.0*ip.
x, y2 = 2.0*ip.
y, z2 = 2.0*ip.
z;
5905 shape(1,0) = x2 - 2.0;
5910 shape(2,1) = y2 - 2.0;
5915 shape(3,2) = z2 - 2.0;
5927 const double RT0TetFiniteElement::nk[4][3] =
5928 {{.5,.5,.5}, {-.5,0,0}, {0,-.5,0}, {0,0,-.5}};
5934 #ifdef MFEM_THREAD_SAFE
5940 for (k = 0; k < 4; k++)
5943 for (j = 0; j < 4; j++)
5945 double d = (
vshape(j,0)*nk[k][0] +
vshape(j,1)*nk[k][1] +
5947 if (j == k) { d -= 1.0; }
5948 if (fabs(d) > 1.0e-12)
5950 cerr <<
"RT0TetFiniteElement::GetLocalInterpolation (...)\n"
5951 " k = " << k <<
", j = " << j <<
", d = " << d << endl;
5959 ip.
x = ip.
y = ip.
z = 0.0;
5967 for (k = 0; k < 4; k++)
5970 ip.
x = vk[0]; ip.
y = vk[1]; ip.
z = vk[2];
5973 vk[0] =
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2];
5974 vk[1] =
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2];
5975 vk[2] =
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2];
5976 for (j = 0; j < 4; j++)
5977 if (fabs (I(k,j) = (
vshape(j,0)*vk[0]+
vshape(j,1)*vk[1]+
5978 vshape(j,2)*vk[2])) < 1.0e-12)
5991 #ifdef MFEM_THREAD_SAFE
5995 for (
int k = 0; k < 4; k++)
6004 vk[0] * (
Jinv(0,0)*nk[k][0]+
Jinv(0,1)*nk[k][1]+
Jinv(0,2)*nk[k][2] ) +
6005 vk[1] * (
Jinv(1,0)*nk[k][0]+
Jinv(1,1)*nk[k][1]+
Jinv(1,2)*nk[k][2] ) +
6006 vk[2] * (
Jinv(2,0)*nk[k][0]+
Jinv(2,1)*nk[k][1]+
Jinv(2,2)*nk[k][2] );
6041 double x = 2. * ip.
x - 1.;
6042 double y = 2. * ip.
y - 1.;
6043 double z = 2. * ip.
z - 1.;
6044 double f5 = x * x - y * y;
6045 double f6 = y * y - z * z;
6047 shape(0) = (1./6.) * (1. - 3. * z - f5 - 2. * f6);
6048 shape(1) = (1./6.) * (1. - 3. * y - f5 + f6);
6049 shape(2) = (1./6.) * (1. + 3. * x + 2. * f5 + f6);
6050 shape(3) = (1./6.) * (1. + 3. * y - f5 + f6);
6051 shape(4) = (1./6.) * (1. - 3. * x + 2. * f5 + f6);
6052 shape(5) = (1./6.) * (1. + 3. * z - f5 - 2. * f6);
6058 const double a = 2./3.;
6060 double xt = a * (1. - 2. * ip.
x);
6061 double yt = a * (1. - 2. * ip.
y);
6062 double zt = a * (1. - 2. * ip.
z);
6066 dshape(0,2) = -1. - 2. * zt;
6069 dshape(1,1) = -1. - 2. * yt;
6072 dshape(2,0) = 1. - 2. * xt;
6077 dshape(3,1) = 1. - 2. * yt;
6080 dshape(4,0) = -1. - 2. * xt;
6086 dshape(5,2) = 1. - 2. * zt;
6091 : x(p + 1), w(p + 1)
6097 for (
int i = 0; i <= p; i++)
6100 for (
int j = 0; j <= p; j++)
6113 for (
int i = 0; i <= p; i++)
6115 for (
int j = 0; j < i; j++)
6117 double xij = x(i) - x(j);
6122 for (
int i = 0; i <= p; i++)
6129 for (
int i = 0; i < p; i++)
6133 mfem_error(
"Poly_1D::Basis::Basis : nodes are not increasing!");
6149 int i, k, p = x.Size() - 1;
6159 for (k = 0; k < p; k++)
6161 if (y >= (x(k) + x(k+1))/2)
6167 for (i = k+1; i <= p; i++)
6174 l = lk * (y - x(k));
6176 for (i = 0; i < k; i++)
6178 u(i) = l * w(i) / (y - x(i));
6181 for (i++; i <= p; i++)
6183 u(i) = l * w(i) / (y - x(i));
6198 int i, k, p = x.Size() - 1;
6199 double l, lp, lk, sk, si;
6209 for (k = 0; k < p; k++)
6211 if (y >= (x(k) + x(k+1))/2)
6217 for (i = k+1; i <= p; i++)
6224 l = lk * (y - x(k));
6227 for (i = 0; i < k; i++)
6229 si = 1.0/(y - x(i));
6231 u(i) = l * si * w(i);
6234 for (i++; i <= p; i++)
6236 si = 1.0/(y - x(i));
6238 u(i) = l * si * w(i);
6242 for (i = 0; i < k; i++)
6244 d(i) = (lp * w(i) - u(i))/(y - x(i));
6247 for (i++; i <= p; i++)
6249 d(i) = (lp * w(i) - u(i))/(y - x(i));
6259 for (
int i = 0; i <= p; i++)
6261 binom(i,0) = binom(i,i) = 1;
6262 for (
int j = 1; j < i; j++)
6264 binom(i,j) = binom(i-1,j) + binom(i-1,j-1);
6273 for (
int i = 0; i <= p; i++)
6276 double s = sin(M_PI_2*(i + 0.5)/(p + 1));
6281 void Poly_1D::CalcMono(
const int p,
const double x,
double *u)
6285 for (
int n = 1; n <= p; n++)
6291 void Poly_1D::CalcMono(
const int p,
const double x,
double *u,
double *d)
6296 for (
int n = 1; n <= p; n++)
6313 const int *b =
Binom(p);
6316 for (i = 1; i < p; i++)
6323 for (i--; i > 0; i--)
6333 double *u,
double *d)
6343 const int *b =
Binom(p);
6344 const double xpy = x + y, ptx = p*x;
6347 for (i = 1; i < p; i++)
6349 d[i] = b[i]*z*(i*xpy - ptx);
6356 for (i--; i > 0; i--)
6377 const int *b =
Binom(p);
6378 const double xpy = x + y, ptx = p*x;
6381 for (i = 1; i < p; i++)
6383 d[i] = b[i]*z*(i*xpy - ptx);
6388 for (i--; i > 0; i--)
6397 void Poly_1D::CalcLegendre(
const int p,
const double x,
double *u)
6403 if (p == 0) {
return; }
6404 u[1] = z = 2.*x - 1.;
6405 for (
int n = 1; n < p; n++)
6407 u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
6411 void Poly_1D::CalcLegendre(
const int p,
const double x,
double *u,
double *d)
6420 if (p == 0) {
return; }
6421 u[1] = z = 2.*x - 1.;
6423 for (
int n = 1; n < p; n++)
6425 u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
6426 d[n+1] = (4*n + 2)*u[n] + d[n-1];
6430 void Poly_1D::CalcChebyshev(
const int p,
const double x,
double *u)
6437 if (p == 0) {
return; }
6438 u[1] = z = 2.*x - 1.;
6439 for (
int n = 1; n < p; n++)
6441 u[n+1] = 2*z*u[n] - u[n-1];
6445 void Poly_1D::CalcChebyshev(
const int p,
const double x,
double *u,
double *d)
6458 if (p == 0) {
return; }
6459 u[1] = z = 2.*x - 1.;
6461 for (
int n = 1; n < p; n++)
6463 u[n+1] = 2*z*u[n] - u[n-1];
6464 d[n+1] = (n + 1)*(z*d[n]/n + 2*u[n]);
6472 if (points_container.find(type) == points_container.end())
6477 if (pts.
Size() <= p)
6483 pts[p] =
new double[p + 1];
6493 if ( bases_container.find(type) == bases_container.end() )
6499 if (bases.
Size() <= p)
6503 if (bases[p] == NULL)
6512 for (std::map<
int,
Array<double*>*>::iterator it = points_container.begin();
6513 it != points_container.end() ; ++it)
6516 for (
int i = 0 ; i < pts.
Size() ; ++i )
6523 for (std::map<
int,
Array<Basis*>*>::iterator it = bases_container.begin();
6524 it != bases_container.end() ; ++it )
6527 for (
int i = 0 ; i < bases.
Size() ; ++i )
6541 pt_type(VerifyClosed(type)),
6547 #ifndef MFEM_THREAD_SAFE
6556 for (
int i = 1; i < p; i++)
6566 const int p =
Order;
6568 #ifdef MFEM_THREAD_SAFE
6572 basis1d.
Eval(ip.
x, shape_x);
6574 shape(0) = shape_x(0);
6575 shape(1) = shape_x(p);
6576 for (
int i = 1; i < p; i++)
6578 shape(i+1) = shape_x(i);
6585 const int p =
Order;
6587 #ifdef MFEM_THREAD_SAFE
6588 Vector shape_x(p+1), dshape_x(p+1);
6591 basis1d.
Eval(ip.
x, shape_x, dshape_x);
6593 dshape(0,0) = dshape_x(0);
6594 dshape(1,0) = dshape_x(p);
6595 for (
int i = 1; i < p; i++)
6597 dshape(i+1,0) = dshape_x(i);
6603 const int p =
Order;
6611 for (
int i = 1; i < p; i++)
6620 for (
int i = 1; i < p; i++)
6632 pt_type(VerifyClosed(type)),
6633 basis1d(
poly1d.ClosedBasis(p, pt_type)),
6634 dof_map((p + 1)*(p + 1))
6638 const int p1 = p + 1;
6640 #ifndef MFEM_THREAD_SAFE
6648 dof_map[0 + 0*p1] = 0;
6649 dof_map[p + 0*p1] = 1;
6650 dof_map[p + p*p1] = 2;
6651 dof_map[0 + p*p1] = 3;
6655 for (
int i = 1; i < p; i++)
6657 dof_map[i + 0*p1] = o++;
6659 for (
int i = 1; i < p; i++)
6661 dof_map[p + i*p1] = o++;
6663 for (
int i = 1; i < p; i++)
6665 dof_map[(p-i) + p*p1] = o++;
6667 for (
int i = 1; i < p; i++)
6669 dof_map[0 + (p-i)*p1] = o++;
6673 for (
int j = 1; j < p; j++)
6674 for (
int i = 1; i < p; i++)
6676 dof_map[i + j*p1] = o++;
6680 for (
int j = 0; j <= p; j++)
6682 for (
int i = 0; i <= p; i++)
6692 const int p =
Order;
6694 #ifdef MFEM_THREAD_SAFE
6695 Vector shape_x(p+1), shape_y(p+1);
6698 basis1d.
Eval(ip.
x, shape_x);
6699 basis1d.
Eval(ip.
y, shape_y);
6701 for (
int o = 0, j = 0; j <= p; j++)
6702 for (
int i = 0; i <= p; i++)
6704 shape(dof_map[o++]) = shape_x(i)*shape_y(j);
6711 const int p =
Order;
6713 #ifdef MFEM_THREAD_SAFE
6714 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
6717 basis1d.
Eval(ip.
x, shape_x, dshape_x);
6718 basis1d.
Eval(ip.
y, shape_y, dshape_y);
6720 for (
int o = 0, j = 0; j <= p; j++)
6722 for (
int i = 0; i <= p; i++)
6724 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j);
6725 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j); o++;
6732 const int p =
Order;
6735 #ifdef MFEM_THREAD_SAFE
6736 Vector shape_x(p+1), shape_y(p+1);
6739 for (
int i = 0; i <= p; i++)
6748 for (
int o = 0, j = 0; j <= p; j++)
6749 for (
int i = 0; i <= p; i++)
6751 dofs(dof_map[o++]) = shape_x(i)*shape_x(j);
6755 for (
int o = 0, j = 0; j <= p; j++)
6756 for (
int i = 0; i <= p; i++)
6758 dofs(dof_map[o++]) = shape_y(i)*shape_x(j);
6762 for (
int o = 0, j = 0; j <= p; j++)
6763 for (
int i = 0; i <= p; i++)
6765 dofs(dof_map[o++]) = shape_y(i)*shape_y(j);
6769 for (
int o = 0, j = 0; j <= p; j++)
6770 for (
int i = 0; i <= p; i++)
6772 dofs(dof_map[o++]) = shape_x(i)*shape_y(j);
6782 pt_type(VerifyClosed(type)),
6783 basis1d(
poly1d.ClosedBasis(p, pt_type)),
6784 dof_map((p + 1)*(p + 1)*(p + 1))
6788 const int p1 = p + 1;
6790 #ifndef MFEM_THREAD_SAFE
6800 dof_map[0 + (0 + 0*p1)*p1] = 0;
6801 dof_map[p + (0 + 0*p1)*p1] = 1;
6802 dof_map[p + (p + 0*p1)*p1] = 2;
6803 dof_map[0 + (p + 0*p1)*p1] = 3;
6804 dof_map[0 + (0 + p*p1)*p1] = 4;
6805 dof_map[p + (0 + p*p1)*p1] = 5;
6806 dof_map[p + (p + p*p1)*p1] = 6;
6807 dof_map[0 + (p + p*p1)*p1] = 7;
6811 for (
int i = 1; i < p; i++)
6813 dof_map[i + (0 + 0*p1)*p1] = o++;
6815 for (
int i = 1; i < p; i++)
6817 dof_map[p + (i + 0*p1)*p1] = o++;
6819 for (
int i = 1; i < p; i++)
6821 dof_map[i + (p + 0*p1)*p1] = o++;
6823 for (
int i = 1; i < p; i++)
6825 dof_map[0 + (i + 0*p1)*p1] = o++;
6827 for (
int i = 1; i < p; i++)
6829 dof_map[i + (0 + p*p1)*p1] = o++;
6831 for (
int i = 1; i < p; i++)
6833 dof_map[p + (i + p*p1)*p1] = o++;
6835 for (
int i = 1; i < p; i++)
6837 dof_map[i + (p + p*p1)*p1] = o++;
6839 for (
int i = 1; i < p; i++)
6841 dof_map[0 + (i + p*p1)*p1] = o++;
6843 for (
int i = 1; i < p; i++)
6845 dof_map[0 + (0 + i*p1)*p1] = o++;
6847 for (
int i = 1; i < p; i++)
6849 dof_map[p + (0 + i*p1)*p1] = o++;
6851 for (
int i = 1; i < p; i++)
6853 dof_map[p + (p + i*p1)*p1] = o++;
6855 for (
int i = 1; i < p; i++)
6857 dof_map[0 + (p + i*p1)*p1] = o++;
6861 for (
int j = 1; j < p; j++)
6862 for (
int i = 1; i < p; i++)
6864 dof_map[i + ((p-j) + 0*p1)*p1] = o++;
6866 for (
int j = 1; j < p; j++)
6867 for (
int i = 1; i < p; i++)
6869 dof_map[i + (0 + j*p1)*p1] = o++;
6871 for (
int j = 1; j < p; j++)
6872 for (
int i = 1; i < p; i++)
6874 dof_map[p + (i + j*p1)*p1] = o++;
6876 for (
int j = 1; j < p; j++)
6877 for (
int i = 1; i < p; i++)
6879 dof_map[(p-i) + (p + j*p1)*p1] = o++;
6881 for (
int j = 1; j < p; j++)
6882 for (
int i = 1; i < p; i++)
6884 dof_map[0 + ((p-i) + j*p1)*p1] = o++;
6886 for (
int j = 1; j < p; j++)
6887 for (
int i = 1; i < p; i++)
6889 dof_map[i + (j + p*p1)*p1] = o++;
6893 for (
int k = 1; k < p; k++)
6894 for (
int j = 1; j < p; j++)
6895 for (
int i = 1; i < p; i++)
6897 dof_map[i + (j + k*p1)*p1] = o++;
6901 for (
int k = 0; k <= p; k++)
6902 for (
int j = 0; j <= p; j++)
6903 for (
int i = 0; i <= p; i++)
6912 const int p =
Order;
6914 #ifdef MFEM_THREAD_SAFE
6915 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
6918 basis1d.
Eval(ip.
x, shape_x);
6919 basis1d.
Eval(ip.
y, shape_y);
6920 basis1d.
Eval(ip.
z, shape_z);
6922 for (
int o = 0, k = 0; k <= p; k++)
6923 for (
int j = 0; j <= p; j++)
6924 for (
int i = 0; i <= p; i++)
6926 shape(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_z(k);
6933 const int p =
Order;
6935 #ifdef MFEM_THREAD_SAFE
6936 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
6937 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
6940 basis1d.
Eval(ip.
x, shape_x, dshape_x);
6941 basis1d.
Eval(ip.
y, shape_y, dshape_y);
6942 basis1d.
Eval(ip.
z, shape_z, dshape_z);
6944 for (
int o = 0, k = 0; k <= p; k++)
6945 for (
int j = 0; j <= p; j++)
6946 for (
int i = 0; i <= p; i++)
6948 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j)* shape_z(k);
6949 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j)* shape_z(k);
6950 dshape(dof_map[o],2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
6956 const int p =
Order;
6959 #ifdef MFEM_THREAD_SAFE
6960 Vector shape_x(p+1), shape_y(p+1);
6963 for (
int i = 0; i <= p; i++)
6972 for (
int o = 0, k = 0; k <= p; k++)
6973 for (
int j = 0; j <= p; j++)
6974 for (
int i = 0; i <= p; i++)
6976 dofs(dof_map[o++]) = shape_x(i)*shape_x(j)*shape_x(k);
6980 for (
int o = 0, k = 0; k <= p; k++)
6981 for (
int j = 0; j <= p; j++)
6982 for (
int i = 0; i <= p; i++)
6984 dofs(dof_map[o++]) = shape_y(i)*shape_x(j)*shape_x(k);
6988 for (
int o = 0, k = 0; k <= p; k++)
6989 for (
int j = 0; j <= p; j++)
6990 for (
int i = 0; i <= p; i++)
6992 dofs(dof_map[o++]) = shape_y(i)*shape_y(j)*shape_x(k);
6996 for (
int o = 0, k = 0; k <= p; k++)
6997 for (
int j = 0; j <= p; j++)
6998 for (
int i = 0; i <= p; i++)
7000 dofs(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_x(k);
7004 for (
int o = 0, k = 0; k <= p; k++)
7005 for (
int j = 0; j <= p; j++)
7006 for (
int i = 0; i <= p; i++)
7008 dofs(dof_map[o++]) = shape_x(i)*shape_x(j)*shape_y(k);
7012 for (
int o = 0, k = 0; k <= p; k++)
7013 for (
int j = 0; j <= p; j++)
7014 for (
int i = 0; i <= p; i++)
7016 dofs(dof_map[o++]) = shape_y(i)*shape_x(j)*shape_y(k);
7020 for (
int o = 0, k = 0; k <= p; k++)
7021 for (
int j = 0; j <= p; j++)
7022 for (
int i = 0; i <= p; i++)
7024 dofs(dof_map[o++]) = shape_y(i)*shape_y(j)*shape_y(k);
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 dofs(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_y(k);
7043 #ifndef MFEM_THREAD_SAFE
7054 for (
int i = 1; i < p; i++)
7064 const int p =
Order;
7066 #ifdef MFEM_THREAD_SAFE
7073 shape(0) = shape_x(0);
7074 shape(1) = shape_x(p);
7075 for (
int i = 1; i < p; i++)
7077 shape(i+1) = shape_x(i);
7084 const int p =
Order;
7086 #ifdef MFEM_THREAD_SAFE
7087 Vector shape_x(p+1), dshape_x(p+1);
7093 dshape(0,0) = dshape_x(0);
7094 dshape(1,0) = dshape_x(p);
7095 for (
int i = 1; i < p; i++)
7097 dshape(i+1,0) = dshape_x(i);
7111 dof_map((p + 1)*(p + 1))
7113 const int p1 = p + 1;
7115 #ifndef MFEM_THREAD_SAFE
7124 dof_map[0 + 0*p1] = 0;
7125 dof_map[p + 0*p1] = 1;
7126 dof_map[p + p*p1] = 2;
7127 dof_map[0 + p*p1] = 3;
7131 for (
int i = 1; i < p; i++)
7133 dof_map[i + 0*p1] = o++;
7135 for (
int i = 1; i < p; i++)
7137 dof_map[p + i*p1] = o++;
7139 for (
int i = 1; i < p; i++)
7141 dof_map[(p-i) + p*p1] = o++;
7143 for (
int i = 1; i < p; i++)
7145 dof_map[0 + (p-i)*p1] = o++;
7149 for (
int j = 1; j < p; j++)
7150 for (
int i = 1; i < p; i++)
7152 dof_map[i + j*p1] = o++;
7156 for (
int j = 0; j <= p; j++)
7157 for (
int i = 0; i <= p; i++)
7166 const int p =
Order;
7168 #ifdef MFEM_THREAD_SAFE
7169 Vector shape_x(p+1), shape_y(p+1);
7176 for (
int o = 0, j = 0; j <= p; j++)
7177 for (
int i = 0; i <= p; i++)
7179 shape(dof_map[o++]) = shape_x(i)*shape_y(j);
7186 const int p =
Order;
7188 #ifdef MFEM_THREAD_SAFE
7189 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
7196 for (
int o = 0, j = 0; j <= p; j++)
7197 for (
int i = 0; i <= p; i++)
7199 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j);
7200 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j); o++;
7214 dof_map((p + 1)*(p + 1)*(p + 1))
7216 const int p1 = p + 1;
7218 #ifndef MFEM_THREAD_SAFE
7229 dof_map[0 + (0 + 0*p1)*p1] = 0;
7230 dof_map[p + (0 + 0*p1)*p1] = 1;
7231 dof_map[p + (p + 0*p1)*p1] = 2;
7232 dof_map[0 + (p + 0*p1)*p1] = 3;
7233 dof_map[0 + (0 + p*p1)*p1] = 4;
7234 dof_map[p + (0 + p*p1)*p1] = 5;
7235 dof_map[p + (p + p*p1)*p1] = 6;
7236 dof_map[0 + (p + p*p1)*p1] = 7;
7240 for (
int i = 1; i < p; i++)
7242 dof_map[i + (0 + 0*p1)*p1] = o++;
7244 for (
int i = 1; i < p; i++)
7246 dof_map[p + (i + 0*p1)*p1] = o++;
7248 for (
int i = 1; i < p; i++)
7250 dof_map[i + (p + 0*p1)*p1] = o++;
7252 for (
int i = 1; i < p; i++)
7254 dof_map[0 + (i + 0*p1)*p1] = o++;
7256 for (
int i = 1; i < p; i++)
7258 dof_map[i + (0 + p*p1)*p1] = o++;
7260 for (
int i = 1; i < p; i++)
7262 dof_map[p + (i + p*p1)*p1] = o++;
7264 for (
int i = 1; i < p; i++)
7266 dof_map[i + (p + p*p1)*p1] = o++;
7268 for (
int i = 1; i < p; i++)
7270 dof_map[0 + (i + p*p1)*p1] = o++;
7272 for (
int i = 1; i < p; i++)
7274 dof_map[0 + (0 + i*p1)*p1] = o++;
7276 for (
int i = 1; i < p; i++)
7278 dof_map[p + (0 + i*p1)*p1] = o++;
7280 for (
int i = 1; i < p; i++)
7282 dof_map[p + (p + i*p1)*p1] = o++;
7284 for (
int i = 1; i < p; i++)
7286 dof_map[0 + (p + i*p1)*p1] = o++;
7290 for (
int j = 1; j < p; j++)
7291 for (
int i = 1; i < p; i++)
7293 dof_map[i + ((p-j) + 0*p1)*p1] = o++;
7295 for (
int j = 1; j < p; j++)
7296 for (
int i = 1; i < p; i++)
7298 dof_map[i + (0 + j*p1)*p1] = o++;
7300 for (
int j = 1; j < p; j++)
7301 for (
int i = 1; i < p; i++)
7303 dof_map[p + (i + j*p1)*p1] = o++;
7305 for (
int j = 1; j < p; j++)
7306 for (
int i = 1; i < p; i++)
7308 dof_map[(p-i) + (p + j*p1)*p1] = o++;
7310 for (
int j = 1; j < p; j++)
7311 for (
int i = 1; i < p; i++)
7313 dof_map[0 + ((p-i) + j*p1)*p1] = o++;
7315 for (
int j = 1; j < p; j++)
7316 for (
int i = 1; i < p; i++)
7318 dof_map[i + (j + p*p1)*p1] = o++;
7322 for (
int k = 1; k < p; k++)
7323 for (
int j = 1; j < p; j++)
7324 for (
int i = 1; i < p; i++)
7326 dof_map[i + (j + k*p1)*p1] = o++;
7330 for (
int k = 0; k <= p; k++)
7331 for (
int j = 0; j <= p; j++)
7332 for (
int i = 0; i <= p; i++)
7340 const int p =
Order;
7342 #ifdef MFEM_THREAD_SAFE
7343 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7350 for (
int o = 0, k = 0; k <= p; k++)
7351 for (
int j = 0; j <= p; j++)
7352 for (
int i = 0; i <= p; i++)
7354 shape(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_z(k);
7361 const int p =
Order;
7363 #ifdef MFEM_THREAD_SAFE
7364 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
7365 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
7372 for (
int o = 0, k = 0; k <= p; k++)
7373 for (
int j = 0; j <= p; j++)
7374 for (
int i = 0; i <= p; i++)
7376 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j)* shape_z(k);
7377 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j)* shape_z(k);
7378 dshape(dof_map[o],2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
7395 #ifndef MFEM_THREAD_SAFE
7405 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
7415 for (
int i = 1; i < p; i++)
7419 for (
int i = 1; i < p; i++)
7423 for (
int i = 1; i < p; i++)
7429 for (
int j = 1; j < p; j++)
7430 for (
int i = 1; i + j < p; i++)
7432 const double w = cp[i] + cp[j] + cp[p-i-j];
7437 for (
int k = 0; k <
Dof; k++)
7445 for (
int j = 0; j <= p; j++)
7446 for (
int i = 0; i + j <= p; i++)
7448 T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
7459 const int p =
Order;
7461 #ifdef MFEM_THREAD_SAFE
7462 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(
Dof);
7469 for (
int o = 0, j = 0; j <= p; j++)
7470 for (
int i = 0; i + j <= p; i++)
7472 u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
7481 const int p =
Order;
7483 #ifdef MFEM_THREAD_SAFE
7484 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
7485 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
7493 for (
int o = 0, j = 0; j <= p; j++)
7494 for (
int i = 0; i + j <= p; i++)
7497 du(o,0) = ((dshape_x(i)* shape_l(k)) -
7498 ( shape_x(i)*dshape_l(k)))*shape_y(j);
7499 du(o,1) = ((dshape_y(j)* shape_l(k)) -
7500 ( shape_y(j)*dshape_l(k)))*shape_x(i);
7504 Ti.
Mult(du, dshape);
7514 #ifndef MFEM_THREAD_SAFE
7526 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7537 for (
int i = 1; i < p; i++)
7541 for (
int i = 1; i < p; i++)
7545 for (
int i = 1; i < p; i++)
7549 for (
int i = 1; i < p; i++)
7553 for (
int i = 1; i < p; i++)
7557 for (
int i = 1; i < p; i++)
7563 for (
int j = 1; j < p; j++)
7564 for (
int i = 1; i + j < p; i++)
7566 double w = cp[i] + cp[j] + cp[p-i-j];
7569 for (
int j = 1; j < p; j++)
7570 for (
int i = 1; i + j < p; i++)
7572 double w = cp[i] + cp[j] + cp[p-i-j];
7575 for (
int j = 1; j < p; j++)
7576 for (
int i = 1; i + j < p; i++)
7578 double w = cp[i] + cp[j] + cp[p-i-j];
7581 for (
int j = 1; j < p; j++)
7582 for (
int i = 1; i + j < p; i++)
7584 double w = cp[i] + cp[j] + cp[p-i-j];
7589 for (
int k = 1; k < p; k++)
7590 for (
int j = 1; j + k < p; j++)
7591 for (
int i = 1; i + j + k < p; i++)
7593 double w = cp[i] + cp[j] + cp[k] + cp[p-i-j-k];
7598 for (
int m = 0; m <
Dof; m++)
7607 for (
int k = 0; k <= p; k++)
7608 for (
int j = 0; j + k <= p; j++)
7609 for (
int i = 0; i + j + k <= p; i++)
7611 T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
7622 const int p =
Order;
7624 #ifdef MFEM_THREAD_SAFE
7625 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7634 for (
int o = 0, k = 0; k <= p; k++)
7635 for (
int j = 0; j + k <= p; j++)
7636 for (
int i = 0; i + j + k <= p; i++)
7638 u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
7647 const int p =
Order;
7649 #ifdef MFEM_THREAD_SAFE
7650 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
7651 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
7660 for (
int o = 0, k = 0; k <= p; k++)
7661 for (
int j = 0; j + k <= p; j++)
7662 for (
int i = 0; i + j + k <= p; i++)
7664 int l = p - i - j - k;
7665 du(o,0) = ((dshape_x(i)* shape_l(l)) -
7666 ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
7667 du(o,1) = ((dshape_y(j)* shape_l(l)) -
7668 ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
7669 du(o,2) = ((dshape_z(k)* shape_l(l)) -
7670 ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
7674 Ti.
Mult(du, dshape);
7682 #ifndef MFEM_THREAD_SAFE
7692 Index(
int p) { p2p3 = 2*p + 3; }
7693 int operator()(
int i,
int j) {
return ((p2p3-j)*j)/2+i; }
7707 for (
int i = 1; i < p; i++)
7712 for (
int i = 1; i < p; i++)
7717 for (
int i = 1; i < p; i++)
7724 for (
int j = 1; j < p; j++)
7725 for (
int i = 1; i + j < p; i++)
7734 const int p,
const double l1,
const double l2,
double *shape)
7736 const double l3 = 1. - l1 - l2;
7746 for (
int o = 0, j = 0; j <= p; j++)
7750 for (
int i = 0; i <= p - j; i++)
7760 const int p,
const double l1,
const double l2,
7761 double *dshape_1d,
double *dshape)
7763 const int dof = ((p + 1)*(p + 2))/2;
7764 const double l3 = 1. - l1 - l2;
7768 for (
int o = 0, j = 0; j <= p; j++)
7772 for (
int i = 0; i <= p - j; i++)
7774 dshape[o++] = s*dshape_1d[i];
7779 for (
int i = 0; i <= p; i++)
7783 for (
int o = i, j = 0; j <= p - i; j++)
7785 dshape[dof + o] = s*dshape_1d[j];
7795 #ifdef MFEM_THREAD_SAFE
7799 for (
int i = 0; i <
Dof; i++)
7808 #ifdef MFEM_THREAD_SAFE
7813 for (
int d = 0; d < 2; d++)
7815 for (
int i = 0; i <
Dof; i++)
7827 #ifndef MFEM_THREAD_SAFE
7837 int tri(
int k) {
return (k*(k + 1))/2; }
7838 int tet(
int k) {
return (k*(k + 1)*(k + 2))/6; }
7839 Index(
int p_) { p = p_; dof = tet(p + 1); }
7840 int operator()(
int i,
int j,
int k)
7841 {
return dof - tet(p - k) - tri(p + 1 - k - j) + i; }
7857 for (
int i = 1; i < p; i++)
7862 for (
int i = 1; i < p; i++)
7867 for (
int i = 1; i < p; i++)
7872 for (
int i = 1; i < p; i++)
7877 for (
int i = 1; i < p; i++)
7882 for (
int i = 1; i < p; i++)
7889 for (
int j = 1; j < p; j++)
7890 for (
int i = 1; i + j < p; i++)
7895 for (
int j = 1; j < p; j++)
7896 for (
int i = 1; i + j < p; i++)
7901 for (
int j = 1; j < p; j++)
7902 for (
int i = 1; i + j < p; i++)
7907 for (
int j = 1; j < p; j++)
7908 for (
int i = 1; i + j < p; i++)
7915 for (
int k = 1; k < p; k++)
7916 for (
int j = 1; j + k < p; j++)
7917 for (
int i = 1; i + j + k < p; i++)
7926 const int p,
const double l1,
const double l2,
const double l3,
7929 const double l4 = 1. - l1 - l2 - l3;
7938 for (
int o = 0, k = 0; k <= p; k++)
7941 const double ek = bp[k]*l3k;
7943 for (
int j = 0; j <= p - k; j++)
7946 double ekj = ek*bpk[j]*l2j;
7947 for (
int i = 0; i <= p - k - j; i++)
7959 const int p,
const double l1,
const double l2,
const double l3,
7960 double *dshape_1d,
double *dshape)
7962 const int dof = ((p + 1)*(p + 2)*(p + 3))/6;
7963 const double l4 = 1. - l1 - l2 - l3;
7971 for (
int o = 0, k = 0; k <= p; k++)
7974 const double ek = bp[k]*l3k;
7976 for (
int j = 0; j <= p - k; j++)
7979 double ekj = ek*bpk[j]*l2j;
7980 for (
int i = 0; i <= p - k - j; i++)
7982 dshape[o++] = dshape_1d[i]*ekj;
7993 for (
int ok = 0, k = 0; k <= p; k++)
7996 const double ek = bp[k]*l3k;
7998 for (
int i = 0; i <= p - k; i++)
8001 double eki = ek*bpk[i]*l1i;
8003 for (
int j = 0; j <= p - k - i; j++)
8005 dshape[dof + o] = dshape_1d[j]*eki;
8011 ok += ((p - k + 2)*(p - k + 1))/2;
8018 for (
int j = 0; j <= p; j++)
8021 const double ej = bp[j]*l2j;
8023 for (
int i = 0; i <= p - j; i++)
8026 double eji = ej*bpj[i]*l1i;
8027 int m = ((p + 2)*(p + 1))/2;
8028 int n = ((p - j + 2)*(p - j + 1))/2;
8029 for (
int o = i, k = 0; k <= p - j - i; k++)
8034 dshape[2*dof + o - n] = dshape_1d[k]*eji;
8047 #ifdef MFEM_THREAD_SAFE
8051 for (
int i = 0; i <
Dof; i++)
8060 #ifdef MFEM_THREAD_SAFE
8065 for (
int d = 0; d < 3; d++)
8067 for (
int i = 0; i <
Dof; i++)
8077 type(VerifyOpen(type)),
8078 basis1d(
poly1d.OpenBasis(p, type))
8082 #ifndef MFEM_THREAD_SAFE
8087 for (
int i = 0; i <= p; i++)
8096 basis1d.
Eval(ip.
x, shape);
8102 #ifdef MFEM_THREAD_SAFE
8107 basis1d.
Eval(ip.
x, shape_x, dshape_x);
8112 const int p =
Order;
8118 for (
int i = 0; i <= p; i++)
8125 for (
int i = 0; i <= p; i++)
8137 #ifndef MFEM_THREAD_SAFE
8148 for (
int i = 0; i <= p; i++)
8164 #ifdef MFEM_THREAD_SAFE
8175 dofs[vertex*
Order] = 1.0;
8182 type(VerifyOpen(_type)),
8183 basis1d(
poly1d.OpenBasis(p, type))
8187 #ifndef MFEM_THREAD_SAFE
8194 for (
int o = 0, j = 0; j <= p; j++)
8195 for (
int i = 0; i <= p; i++)
8204 const int p =
Order;
8206 #ifdef MFEM_THREAD_SAFE
8207 Vector shape_x(p+1), shape_y(p+1);
8210 basis1d.
Eval(ip.
x, shape_x);
8211 basis1d.
Eval(ip.
y, shape_y);
8213 for (
int o = 0, j = 0; j <= p; j++)
8214 for (
int i = 0; i <= p; i++)
8216 shape(o++) = shape_x(i)*shape_y(j);
8223 const int p =
Order;
8225 #ifdef MFEM_THREAD_SAFE
8226 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
8229 basis1d.
Eval(ip.
x, shape_x, dshape_x);
8230 basis1d.
Eval(ip.
y, shape_y, dshape_y);
8232 for (
int o = 0, j = 0; j <= p; j++)
8233 for (
int i = 0; i <= p; i++)
8235 dshape(o,0) = dshape_x(i)* shape_y(j);
8236 dshape(o,1) = shape_x(i)*dshape_y(j); o++;
8242 const int p =
Order;
8245 #ifdef MFEM_THREAD_SAFE
8246 Vector shape_x(p+1), shape_y(p+1);
8249 for (
int i = 0; i <= p; i++)
8258 for (
int o = 0, j = 0; j <= p; j++)
8259 for (
int i = 0; i <= p; i++)
8261 dofs[o++] = shape_x(i)*shape_x(j);
8265 for (
int o = 0, j = 0; j <= p; j++)
8266 for (
int i = 0; i <= p; i++)
8268 dofs[o++] = shape_y(i)*shape_x(j);
8272 for (
int o = 0, j = 0; j <= p; j++)
8273 for (
int i = 0; i <= p; i++)
8275 dofs[o++] = shape_y(i)*shape_y(j);
8279 for (
int o = 0, j = 0; j <= p; j++)
8280 for (
int i = 0; i <= p; i++)
8282 dofs[o++] = shape_x(i)*shape_y(j);
8293 #ifndef MFEM_THREAD_SAFE
8306 for (
int o = 0, j = 0; j <= p; j++)
8307 for (
int i = 0; i <= p; i++)
8317 const int p =
Order;
8319 #ifdef MFEM_THREAD_SAFE
8320 Vector shape_x(p+1), shape_y(p+1);
8326 for (
int o = 0, j = 0; j <= p; j++)
8327 for (
int i = 0; i <= p; i++)
8329 shape(o++) = shape_x(i)*shape_y(j);
8336 const int p =
Order;
8338 #ifdef MFEM_THREAD_SAFE
8339 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
8345 for (
int o = 0, j = 0; j <= p; j++)
8346 for (
int i = 0; i <= p; i++)
8348 dshape(o,0) = dshape_x(i)* shape_y(j);
8349 dshape(o,1) = shape_x(i)*dshape_y(j); o++;
8355 const int p =
Order;
8360 case 0: dofs[0] = 1.0;
break;
8361 case 1: dofs[p] = 1.0;
break;
8362 case 2: dofs[p*(p + 2)] = 1.0;
break;
8363 case 3: dofs[p*(p + 1)] = 1.0;
break;
8371 type(VerifyOpen(_type)),
8372 basis1d(
poly1d.OpenBasis(p, type))
8376 #ifndef MFEM_THREAD_SAFE
8385 for (
int o = 0, k = 0; k <= p; k++)
8386 for (
int j = 0; j <= p; j++)
8387 for (
int i = 0; i <= p; i++)
8396 const int p =
Order;
8398 #ifdef MFEM_THREAD_SAFE
8399 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
8402 basis1d.
Eval(ip.
x, shape_x);
8403 basis1d.
Eval(ip.
y, shape_y);
8404 basis1d.
Eval(ip.
z, shape_z);
8406 for (
int o = 0, k = 0; k <= p; k++)
8407 for (
int j = 0; j <= p; j++)
8408 for (
int i = 0; i <= p; i++)
8410 shape(o++) = shape_x(i)*shape_y(j)*shape_z(k);
8417 const int p =
Order;
8419 #ifdef MFEM_THREAD_SAFE
8420 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
8421 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
8424 basis1d.
Eval(ip.
x, shape_x, dshape_x);
8425 basis1d.
Eval(ip.
y, shape_y, dshape_y);
8426 basis1d.
Eval(ip.
z, shape_z, dshape_z);
8428 for (
int o = 0, k = 0; k <= p; k++)
8429 for (
int j = 0; j <= p; j++)
8430 for (
int i = 0; i <= p; i++)
8432 dshape(o,0) = dshape_x(i)* shape_y(j)* shape_z(k);
8433 dshape(o,1) = shape_x(i)*dshape_y(j)* shape_z(k);
8434 dshape(o,2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
8440 const int p =
Order;
8443 #ifdef MFEM_THREAD_SAFE
8444 Vector shape_x(p+1), shape_y(p+1);
8447 for (
int i = 0; i <= p; i++)
8456 for (
int o = 0, k = 0; k <= p; k++)
8457 for (
int j = 0; j <= p; j++)
8458 for (
int i = 0; i <= p; i++)
8460 dofs[o++] = shape_x(i)*shape_x(j)*shape_x(k);
8464 for (
int o = 0, k = 0; k <= p; k++)
8465 for (
int j = 0; j <= p; j++)
8466 for (
int i = 0; i <= p; i++)
8468 dofs[o++] = shape_y(i)*shape_x(j)*shape_x(k);
8472 for (
int o = 0, k = 0; k <= p; k++)
8473 for (
int j = 0; j <= p; j++)
8474 for (
int i = 0; i <= p; i++)
8476 dofs[o++] = shape_y(i)*shape_y(j)*shape_x(k);
8480 for (
int o = 0, k = 0; k <= p; k++)
8481 for (
int j = 0; j <= p; j++)
8482 for (
int i = 0; i <= p; i++)
8484 dofs[o++] = shape_x(i)*shape_y(j)*shape_x(k);
8488 for (
int o = 0, k = 0; k <= p; k++)
8489 for (
int j = 0; j <= p; j++)
8490 for (
int i = 0; i <= p; i++)
8492 dofs[o++] = shape_x(i)*shape_x(j)*shape_y(k);
8496 for (
int o = 0, k = 0; k <= p; k++)
8497 for (
int j = 0; j <= p; j++)
8498 for (
int i = 0; i <= p; i++)
8500 dofs[o++] = shape_y(i)*shape_x(j)*shape_y(k);
8504 for (
int o = 0, k = 0; k <= p; k++)
8505 for (
int j = 0; j <= p; j++)
8506 for (
int i = 0; i <= p; i++)
8508 dofs[o++] = shape_y(i)*shape_y(j)*shape_y(k);
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 dofs[o++] = shape_x(i)*shape_y(j)*shape_y(k);
8527 #ifndef MFEM_THREAD_SAFE
8542 for (
int o = 0, k = 0; k <= p; k++)
8543 for (
int j = 0; j <= p; j++)
8544 for (
int i = 0; i <= p; i++)
8554 const int p =
Order;
8556 #ifdef MFEM_THREAD_SAFE
8557 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
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 shape(o++) = shape_x(i)*shape_y(j)*shape_z(k);
8575 const int p =
Order;
8577 #ifdef MFEM_THREAD_SAFE
8578 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
8579 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
8586 for (
int o = 0, k = 0; k <= p; k++)
8587 for (
int j = 0; j <= p; j++)
8588 for (
int i = 0; i <= p; i++)
8590 dshape(o,0) = dshape_x(i)* shape_y(j)* shape_z(k);
8591 dshape(o,1) = shape_x(i)*dshape_y(j)* shape_z(k);
8592 dshape(o,2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
8598 const int p =
Order;
8603 case 0: dofs[0] = 1.0;
break;
8604 case 1: dofs[p] = 1.0;
break;
8605 case 2: dofs[p*(p + 2)] = 1.0;
break;
8606 case 3: dofs[p*(p + 1)] = 1.0;
break;
8607 case 4: dofs[p*(p + 1)*(p + 1)] = 1.0;
break;
8608 case 5: dofs[p + p*(p + 1)*(p + 1)] = 1.0;
break;
8609 case 6: dofs[
Dof - 1] = 1.0;
break;
8610 case 7: dofs[
Dof - p - 1] = 1.0;
break;
8621 #ifndef MFEM_THREAD_SAFE
8631 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
8634 for (
int o = 0, j = 0; j <= p; j++)
8635 for (
int i = 0; i + j <= p; i++)
8637 double w = op[i] + op[j] + op[p-i-j];
8642 for (
int k = 0; k <
Dof; k++)
8649 for (
int o = 0, j = 0; j <= p; j++)
8650 for (
int i = 0; i + j <= p; i++)
8652 T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
8663 const int p =
Order;
8665 #ifdef MFEM_THREAD_SAFE
8666 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(
Dof);
8673 for (
int o = 0, j = 0; j <= p; j++)
8674 for (
int i = 0; i + j <= p; i++)
8676 u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
8685 const int p =
Order;
8687 #ifdef MFEM_THREAD_SAFE
8688 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
8689 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
8697 for (
int o = 0, j = 0; j <= p; j++)
8698 for (
int i = 0; i + j <= p; i++)
8701 du(o,0) = ((dshape_x(i)* shape_l(k)) -
8702 ( shape_x(i)*dshape_l(k)))*shape_y(j);
8703 du(o,1) = ((dshape_y(j)* shape_l(k)) -
8704 ( shape_y(j)*dshape_l(k)))*shape_x(i);
8708 Ti.
Mult(du, dshape);
8716 for (
int i = 0; i <
Dof; i++)
8719 dofs[i] = pow(1.0 - ip.
x - ip.
y,
Order);
8723 for (
int i = 0; i <
Dof; i++)
8726 dofs[i] = pow(ip.
x,
Order);
8730 for (
int i = 0; i <
Dof; i++)
8733 dofs[i] = pow(ip.
y,
Order);
8744 #ifndef MFEM_THREAD_SAFE
8754 for (
int o = 0, j = 0; j <= p; j++)
8755 for (
int i = 0; i + j <= p; i++)
8771 #ifdef MFEM_THREAD_SAFE
8784 case 0: dofs[0] = 1.0;
break;
8785 case 1: dofs[
Order] = 1.0;
break;
8786 case 2: dofs[
Dof-1] = 1.0;
break;
8797 #ifndef MFEM_THREAD_SAFE
8809 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8812 for (
int o = 0, k = 0; k <= p; k++)
8813 for (
int j = 0; j + k <= p; j++)
8814 for (
int i = 0; i + j + k <= p; i++)
8816 double w = op[i] + op[j] + op[k] + op[p-i-j-k];
8821 for (
int m = 0; m <
Dof; m++)
8829 for (
int o = 0, k = 0; k <= p; k++)
8830 for (
int j = 0; j + k <= p; j++)
8831 for (
int i = 0; i + j + k <= p; i++)
8833 T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
8844 const int p =
Order;
8846 #ifdef MFEM_THREAD_SAFE
8847 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8856 for (
int o = 0, k = 0; k <= p; k++)
8857 for (
int j = 0; j + k <= p; j++)
8858 for (
int i = 0; i + j + k <= p; i++)
8860 u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
8869 const int p =
Order;
8871 #ifdef MFEM_THREAD_SAFE
8872 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
8873 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
8882 for (
int o = 0, k = 0; k <= p; k++)
8883 for (
int j = 0; j + k <= p; j++)
8884 for (
int i = 0; i + j + k <= p; i++)
8886 int l = p - i - j - k;
8887 du(o,0) = ((dshape_x(i)* shape_l(l)) -
8888 ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
8889 du(o,1) = ((dshape_y(j)* shape_l(l)) -
8890 ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
8891 du(o,2) = ((dshape_z(k)* shape_l(l)) -
8892 ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
8896 Ti.
Mult(du, dshape);
8904 for (
int i = 0; i <
Dof; i++)
8907 dofs[i] = pow(1.0 - ip.
x - ip.
y - ip.
z,
Order);
8911 for (
int i = 0; i <
Dof; i++)
8914 dofs[i] = pow(ip.
x,
Order);
8918 for (
int i = 0; i <
Dof; i++)
8921 dofs[i] = pow(ip.
y,
Order);
8924 for (
int i = 0; i <
Dof; i++)
8927 dofs[i] = pow(ip.
z,
Order);
8938 #ifndef MFEM_THREAD_SAFE
8948 for (
int o = 0, k = 0; k <= p; k++)
8949 for (
int j = 0; j + k <= p; j++)
8950 for (
int i = 0; i + j + k <= p; i++)
8967 #ifdef MFEM_THREAD_SAFE
8980 case 0: dofs[0] = 1.0;
break;
8981 case 1: dofs[
Order] = 1.0;
break;
8982 case 2: dofs[(
Order*(
Order+3))/2] = 1.0;
break;
8983 case 3: dofs[
Dof-1] = 1.0;
break;
8988 const double RT_QuadrilateralElement::nk[8] =
8989 { 0., -1., 1., 0., 0., 1., -1., 0. };
8996 cbasis1d(
poly1d.ClosedBasis(p + 1, VerifyClosed(cp_type))),
8997 obasis1d(
poly1d.OpenBasis(p, VerifyOpen(op_type))),
8998 dof_map(Dof), dof2nk(Dof)
9002 const int dof2 =
Dof/2;
9004 #ifndef MFEM_THREAD_SAFE
9015 for (
int i = 0; i <= p; i++)
9017 dof_map[1*dof2 + i + 0*(p + 1)] = o++;
9019 for (
int i = 0; i <= p; i++)
9021 dof_map[0*dof2 + (p + 1) + i*(p + 2)] = o++;
9023 for (
int i = 0; i <= p; i++)
9025 dof_map[1*dof2 + (p - i) + (p + 1)*(p + 1)] = o++;
9027 for (
int i = 0; i <= p; i++)
9029 dof_map[0*dof2 + 0 + (p - i)*(p + 2)] = o++;
9033 for (
int j = 0; j <= p; j++)
9034 for (
int i = 1; i <= p; i++)
9036 dof_map[0*dof2 + i + j*(p + 2)] = o++;
9038 for (
int j = 1; j <= p; j++)
9039 for (
int i = 0; i <= p; i++)
9041 dof_map[1*dof2 + i + j*(p + 1)] = o++;
9046 for (
int j = 0; j <= p; j++)
9047 for (
int i = 0; i <= p/2; i++)
9049 int idx = 0*dof2 + i + j*(p + 2);
9050 dof_map[idx] = -1 - dof_map[idx];
9053 for (
int j = p/2 + 1; j <= p; j++)
9055 int idx = 0*dof2 + (p/2 + 1) + j*(p + 2);
9056 dof_map[idx] = -1 - dof_map[idx];
9059 for (
int j = 0; j <= p/2; j++)
9060 for (
int i = 0; i <= p; i++)
9062 int idx = 1*dof2 + i + j*(p + 1);
9063 dof_map[idx] = -1 - dof_map[idx];
9066 for (
int i = 0; i <= p/2; i++)
9068 int idx = 1*dof2 + i + (p/2 + 1)*(p + 1);
9069 dof_map[idx] = -1 - dof_map[idx];
9073 for (
int j = 0; j <= p; j++)
9074 for (
int i = 0; i <= p + 1; i++)
9077 if ((idx = dof_map[o++]) < 0)
9088 for (
int j = 0; j <= p + 1; j++)
9089 for (
int i = 0; i <= p; i++)
9092 if ((idx = dof_map[o++]) < 0)
9108 const int pp1 =
Order;
9110 #ifdef MFEM_THREAD_SAFE
9111 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9114 cbasis1d.
Eval(ip.
x, shape_cx);
9115 obasis1d.
Eval(ip.
x, shape_ox);
9116 cbasis1d.
Eval(ip.
y, shape_cy);
9117 obasis1d.
Eval(ip.
y, shape_oy);
9120 for (
int j = 0; j < pp1; j++)
9121 for (
int i = 0; i <= pp1; i++)
9124 if ((idx = dof_map[o++]) < 0)
9126 idx = -1 - idx, s = -1;
9132 shape(idx,0) = s*shape_cx(i)*shape_oy(j);
9135 for (
int j = 0; j <= pp1; j++)
9136 for (
int i = 0; i < pp1; i++)
9139 if ((idx = dof_map[o++]) < 0)
9141 idx = -1 - idx, s = -1;
9148 shape(idx,1) = s*shape_ox(i)*shape_cy(j);
9155 const int pp1 =
Order;
9157 #ifdef MFEM_THREAD_SAFE
9158 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9159 Vector dshape_cx(pp1 + 1), dshape_cy(pp1 + 1);
9162 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
9163 obasis1d.
Eval(ip.
x, shape_ox);
9164 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
9165 obasis1d.
Eval(ip.
y, shape_oy);
9168 for (
int j = 0; j < pp1; j++)
9169 for (
int i = 0; i <= pp1; i++)
9172 if ((idx = dof_map[o++]) < 0)
9174 idx = -1 - idx, s = -1;
9180 divshape(idx) = s*dshape_cx(i)*shape_oy(j);
9182 for (
int j = 0; j <= pp1; j++)
9183 for (
int i = 0; i < pp1; i++)
9186 if ((idx = dof_map[o++]) < 0)
9188 idx = -1 - idx, s = -1;
9194 divshape(idx) = s*shape_ox(i)*dshape_cy(j);
9199 const double RT_HexahedronElement::nk[18] =
9200 { 0.,0.,-1., 0.,-1.,0., 1.,0.,0., 0.,1.,0., -1.,0.,0., 0.,0.,1. };
9207 cbasis1d(
poly1d.ClosedBasis(p + 1, VerifyClosed(cp_type))),
9208 obasis1d(
poly1d.OpenBasis(p, VerifyOpen(op_type))),
9209 dof_map(Dof), dof2nk(Dof)
9213 const int dof3 =
Dof/3;
9215 #ifndef MFEM_THREAD_SAFE
9229 for (
int j = 0; j <= p; j++)
9230 for (
int i = 0; i <= p; i++)
9232 dof_map[2*dof3 + i + ((p - j) + 0*(p + 1))*(p + 1)] = o++;
9234 for (
int j = 0; j <= p; j++)
9235 for (
int i = 0; i <= p; i++)
9237 dof_map[1*dof3 + i + (0 + j*(p + 2))*(p + 1)] = o++;
9239 for (
int j = 0; j <= p; j++)
9240 for (
int i = 0; i <= p; i++)
9242 dof_map[0*dof3 + (p + 1) + (i + j*(p + 1))*(p + 2)] = o++;
9244 for (
int j = 0; j <= p; j++)
9245 for (
int i = 0; i <= p; i++)
9247 dof_map[1*dof3 + (p - i) + ((p + 1) + j*(p + 2))*(p + 1)] = o++;
9249 for (
int j = 0; j <= p; j++)
9250 for (
int i = 0; i <= p; i++)
9252 dof_map[0*dof3 + 0 + ((p - i) + j*(p + 1))*(p + 2)] = o++;
9254 for (
int j = 0; j <= p; j++)
9255 for (
int i = 0; i <= p; i++)
9257 dof_map[2*dof3 + i + (j + (p + 1)*(p + 1))*(p + 1)] = o++;
9262 for (
int k = 0; k <= p; k++)
9263 for (
int j = 0; j <= p; j++)
9264 for (
int i = 1; i <= p; i++)
9266 dof_map[0*dof3 + i + (j + k*(p + 1))*(p + 2)] = o++;
9269 for (
int k = 0; k <= p; k++)
9270 for (
int j = 1; j <= p; j++)
9271 for (
int i = 0; i <= p; i++)
9273 dof_map[1*dof3 + i + (j + k*(p + 2))*(p + 1)] = o++;
9276 for (
int k = 1; k <= p; k++)
9277 for (
int j = 0; j <= p; j++)
9278 for (
int i = 0; i <= p; i++)
9280 dof_map[2*dof3 + i + (j + k*(p + 1))*(p + 1)] = o++;
9288 for (
int k = 0; k <= p; k++)
9289 for (
int j = 0; j <= p; j++)
9290 for (
int i = 0; i <= p/2; i++)
9292 int idx = 0*dof3 + i + (j + k*(p + 1))*(p + 2);
9293 dof_map[idx] = -1 - dof_map[idx];
9296 for (
int k = 0; k <= p; k++)
9297 for (
int j = 0; j <= p/2; j++)
9298 for (
int i = 0; i <= p; i++)
9300 int idx = 1*dof3 + i + (j + k*(p + 2))*(p + 1);
9301 dof_map[idx] = -1 - dof_map[idx];
9304 for (
int k = 0; k <= p/2; k++)
9305 for (
int j = 0; j <= p; j++)
9306 for (
int i = 0; i <= p; i++)
9308 int idx = 2*dof3 + i + (j + k*(p + 1))*(p + 1);
9309 dof_map[idx] = -1 - dof_map[idx];
9314 for (
int k = 0; k <= p; k++)
9315 for (
int j = 0; j <= p; j++)
9316 for (
int i = 0; i <= p + 1; i++)
9319 if ((idx = dof_map[o++]) < 0)
9331 for (
int k = 0; k <= p; k++)
9332 for (
int j = 0; j <= p + 1; j++)
9333 for (
int i = 0; i <= p; i++)
9336 if ((idx = dof_map[o++]) < 0)
9348 for (
int k = 0; k <= p + 1; k++)
9349 for (
int j = 0; j <= p; j++)
9350 for (
int i = 0; i <= p; i++)
9353 if ((idx = dof_map[o++]) < 0)
9369 const int pp1 =
Order;
9371 #ifdef MFEM_THREAD_SAFE
9372 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9373 Vector shape_cz(pp1 + 1), shape_oz(pp1);
9376 cbasis1d.
Eval(ip.
x, shape_cx);
9377 obasis1d.
Eval(ip.
x, shape_ox);
9378 cbasis1d.
Eval(ip.
y, shape_cy);
9379 obasis1d.
Eval(ip.
y, shape_oy);
9380 cbasis1d.
Eval(ip.
z, shape_cz);
9381 obasis1d.
Eval(ip.
z, shape_oz);
9385 for (
int k = 0; k < pp1; k++)
9386 for (
int j = 0; j < pp1; j++)
9387 for (
int i = 0; i <= pp1; i++)
9390 if ((idx = dof_map[o++]) < 0)
9392 idx = -1 - idx, s = -1;
9398 shape(idx,0) = s*shape_cx(i)*shape_oy(j)*shape_oz(k);
9403 for (
int k = 0; k < pp1; k++)
9404 for (
int j = 0; j <= pp1; j++)
9405 for (
int i = 0; i < pp1; i++)
9408 if ((idx = dof_map[o++]) < 0)
9410 idx = -1 - idx, s = -1;
9417 shape(idx,1) = s*shape_ox(i)*shape_cy(j)*shape_oz(k);
9421 for (
int k = 0; k <= pp1; k++)
9422 for (
int j = 0; j < pp1; j++)
9423 for (
int i = 0; i < pp1; i++)
9426 if ((idx = dof_map[o++]) < 0)
9428 idx = -1 - idx, s = -1;
9436 shape(idx,2) = s*shape_ox(i)*shape_oy(j)*shape_cz(k);
9443 const int pp1 =
Order;
9445 #ifdef MFEM_THREAD_SAFE
9446 Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
9447 Vector shape_cz(pp1 + 1), shape_oz(pp1);
9448 Vector dshape_cx(pp1 + 1), dshape_cy(pp1 + 1), dshape_cz(pp1 + 1);
9451 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
9452 obasis1d.
Eval(ip.
x, shape_ox);
9453 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
9454 obasis1d.
Eval(ip.
y, shape_oy);
9455 cbasis1d.
Eval(ip.
z, shape_cz, dshape_cz);
9456 obasis1d.
Eval(ip.
z, shape_oz);
9460 for (
int k = 0; k < pp1; k++)
9461 for (
int j = 0; j < pp1; j++)
9462 for (
int i = 0; i <= pp1; i++)
9465 if ((idx = dof_map[o++]) < 0)
9467 idx = -1 - idx, s = -1;
9473 divshape(idx) = s*dshape_cx(i)*shape_oy(j)*shape_oz(k);
9476 for (
int k = 0; k < pp1; k++)
9477 for (
int j = 0; j <= pp1; j++)
9478 for (
int i = 0; i < pp1; i++)
9481 if ((idx = dof_map[o++]) < 0)
9483 idx = -1 - idx, s = -1;
9489 divshape(idx) = s*shape_ox(i)*dshape_cy(j)*shape_oz(k);
9492 for (
int k = 0; k <= pp1; k++)
9493 for (
int j = 0; j < pp1; j++)
9494 for (
int i = 0; i < pp1; i++)
9497 if ((idx = dof_map[o++]) < 0)
9499 idx = -1 - idx, s = -1;
9505 divshape(idx) = s*shape_ox(i)*shape_oy(j)*dshape_cz(k);
9510 const double RT_TriangleElement::nk[6] =
9511 { 0., -1., 1., 1., -1., 0. };
9513 const double RT_TriangleElement::c = 1./3.;
9523 #ifndef MFEM_THREAD_SAFE
9533 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
9538 for (
int i = 0; i <= p; i++)
9543 for (
int i = 0; i <= p; i++)
9548 for (
int i = 0; i <= p; i++)
9555 for (
int j = 0; j < p; j++)
9556 for (
int i = 0; i + j < p; i++)
9558 double w = iop[i] + iop[j] + iop[p-1-i-j];
9566 for (
int k = 0; k <
Dof; k++)
9572 const double *n_k = nk + 2*dof2nk[k];
9575 for (
int j = 0; j <= p; j++)
9576 for (
int i = 0; i + j <= p; i++)
9578 double s = shape_x(i)*shape_y(j)*shape_l(p-i-j);
9579 T(o++, k) = s*n_k[0];
9580 T(o++, k) = s*n_k[1];
9582 for (
int i = 0; i <= p; i++)
9584 double s = shape_x(i)*shape_y(p-i);
9585 T(o++, k) = s*((ip.
x - c)*n_k[0] + (ip.
y - c)*n_k[1]);
9596 const int p =
Order - 1;
9598 #ifdef MFEM_THREAD_SAFE
9599 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
9608 for (
int j = 0; j <= p; j++)
9609 for (
int i = 0; i + j <= p; i++)
9611 double s = shape_x(i)*shape_y(j)*shape_l(p-i-j);
9612 u(o,0) = s; u(o,1) = 0; o++;
9613 u(o,0) = 0; u(o,1) = s; o++;
9615 for (
int i = 0; i <= p; i++)
9617 double s = shape_x(i)*shape_y(p-i);
9618 u(o,0) = (ip.
x - c)*s;
9619 u(o,1) = (ip.
y - c)*s;
9629 const int p =
Order - 1;
9631 #ifdef MFEM_THREAD_SAFE
9632 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
9633 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
9642 for (
int j = 0; j <= p; j++)
9643 for (
int i = 0; i + j <= p; i++)
9646 divu(o++) = (dshape_x(i)*shape_l(k) -
9647 shape_x(i)*dshape_l(k))*shape_y(j);
9648 divu(o++) = (dshape_y(j)*shape_l(k) -
9649 shape_y(j)*dshape_l(k))*shape_x(i);
9651 for (
int i = 0; i <= p; i++)
9654 divu(o++) = ((shape_x(i) + (ip.
x - c)*dshape_x(i))*shape_y(j) +
9655 (shape_y(j) + (ip.
y - c)*dshape_y(j))*shape_x(i));
9658 Ti.
Mult(divu, divshape);
9662 const double RT_TetrahedronElement::nk[12] =
9663 { 1,1,1, -1,0,0, 0,-1,0, 0,0,-1 };
9666 const double RT_TetrahedronElement::c = 1./4.;
9676 #ifndef MFEM_THREAD_SAFE
9688 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
9694 for (
int j = 0; j <= p; j++)
9695 for (
int i = 0; i + j <= p; i++)
9697 double w = bop[i] + bop[j] + bop[p-i-j];
9701 for (
int j = 0; j <= p; j++)
9702 for (
int i = 0; i + j <= p; i++)
9704 double w = bop[i] + bop[j] + bop[p-i-j];
9708 for (
int j = 0; j <= p; j++)
9709 for (
int i = 0; i + j <= p; i++)
9711 double w = bop[i] + bop[j] + bop[p-i-j];
9715 for (
int j = 0; j <= p; j++)
9716 for (
int i = 0; i + j <= p; i++)
9718 double w = bop[i] + bop[j] + bop[p-i-j];
9724 for (
int k = 0; k < p; k++)
9725 for (
int j = 0; j + k < p; j++)
9726 for (
int i = 0; i + j + k < p; i++)
9728 double w = iop[i] + iop[j] + iop[k] + iop[p-1-i-j-k];
9738 for (
int m = 0; m <
Dof; m++)
9745 const double *nm = nk + 3*dof2nk[m];
9748 for (
int k = 0; k <= p; k++)
9749 for (
int j = 0; j + k <= p; j++)
9750 for (
int i = 0; i + j + k <= p; i++)
9752 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
9753 T(o++, m) = s * nm[0];
9754 T(o++, m) = s * nm[1];
9755 T(o++, m) = s * nm[2];
9757 for (
int j = 0; j <= p; j++)
9758 for (
int i = 0; i + j <= p; i++)
9760 double s = shape_x(i)*shape_y(j)*shape_z(p-i-j);
9761 T(o++, m) = s*((ip.
x - c)*nm[0] + (ip.
y - c)*nm[1] +
9773 const int p =
Order - 1;
9775 #ifdef MFEM_THREAD_SAFE
9776 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
9786 for (
int k = 0; k <= p; k++)
9787 for (
int j = 0; j + k <= p; j++)
9788 for (
int i = 0; i + j + k <= p; i++)
9790 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
9791 u(o,0) = s; u(o,1) = 0; u(o,2) = 0; o++;
9792 u(o,0) = 0; u(o,1) = s; u(o,2) = 0; o++;
9793 u(o,0) = 0; u(o,1) = 0; u(o,2) = s; o++;
9795 for (
int j = 0; j <= p; j++)
9796 for (
int i = 0; i + j <= p; i++)
9798 double s = shape_x(i)*shape_y(j)*shape_z(p-i-j);
9799 u(o,0) = (ip.
x - c)*s; u(o,1) = (ip.
y - c)*s; u(o,2) = (ip.
z - c)*s;
9809 const int p =
Order - 1;
9811 #ifdef MFEM_THREAD_SAFE
9812 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
9813 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
9823 for (
int k = 0; k <= p; k++)
9824 for (
int j = 0; j + k <= p; j++)
9825 for (
int i = 0; i + j + k <= p; i++)
9827 int l = p - i - j - k;
9828 divu(o++) = (dshape_x(i)*shape_l(l) -
9829 shape_x(i)*dshape_l(l))*shape_y(j)*shape_z(k);
9830 divu(o++) = (dshape_y(j)*shape_l(l) -
9831 shape_y(j)*dshape_l(l))*shape_x(i)*shape_z(k);
9832 divu(o++) = (dshape_z(k)*shape_l(l) -
9833 shape_z(k)*dshape_l(l))*shape_x(i)*shape_y(j);
9835 for (
int j = 0; j <= p; j++)
9836 for (
int i = 0; i + j <= p; i++)
9840 (shape_x(i) + (ip.
x - c)*dshape_x(i))*shape_y(j)*shape_z(k) +
9841 (shape_y(j) + (ip.
y - c)*dshape_y(j))*shape_x(i)*shape_z(k) +
9842 (shape_z(k) + (ip.
z - c)*dshape_z(k))*shape_x(i)*shape_y(j);
9845 Ti.
Mult(divu, divshape);
9849 const double ND_HexahedronElement::tk[18] =
9850 { 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,0.,0., 0.,-1.,0., 0.,0.,-1. };
9853 const int cp_type,
const int op_type)
9856 cbasis1d(
poly1d.ClosedBasis(p, VerifyClosed(cp_type))),
9857 obasis1d(
poly1d.OpenBasis(p - 1, VerifyOpen(op_type))),
9858 dof_map(Dof), dof2tk(Dof)
9862 const int dof3 =
Dof/3;
9864 #ifndef MFEM_THREAD_SAFE
9878 for (
int i = 0; i < p; i++)
9880 dof_map[0*dof3 + i + (0 + 0*(p + 1))*p] = o++;
9882 for (
int i = 0; i < p; i++)
9884 dof_map[1*dof3 + p + (i + 0*p)*(p + 1)] = o++;
9886 for (
int i = 0; i < p; i++)
9888 dof_map[0*dof3 + i + (p + 0*(p + 1))*p] = o++;
9890 for (
int i = 0; i < p; i++)
9892 dof_map[1*dof3 + 0 + (i + 0*p)*(p + 1)] = o++;
9894 for (
int i = 0; i < p; i++)
9896 dof_map[0*dof3 + i + (0 + p*(p + 1))*p] = o++;
9898 for (
int i = 0; i < p; i++)
9900 dof_map[1*dof3 + p + (i + p*p)*(p + 1)] = o++;
9902 for (
int i = 0; i < p; i++)
9904 dof_map[0*dof3 + i + (p + p*(p + 1))*p] = o++;
9906 for (
int i = 0; i < p; i++)
9908 dof_map[1*dof3 + 0 + (i + p*p)*(p + 1)] = o++;
9910 for (
int i = 0; i < p; i++)
9912 dof_map[2*dof3 + 0 + (0 + i*(p + 1))*(p + 1)] = o++;
9914 for (
int i = 0; i < p; i++)
9916 dof_map[2*dof3 + p + (0 + i*(p + 1))*(p + 1)] = o++;
9918 for (
int i = 0; i < p; i++)
9920 dof_map[2*dof3 + p + (p + i*(p + 1))*(p + 1)] = o++;
9922 for (
int i = 0; i < p; i++)
9924 dof_map[2*dof3 + 0 + (p + i*(p + 1))*(p + 1)] = o++;
9929 for (
int j = 1; j < p; j++)
9930 for (
int i = 0; i < p; i++)
9932 dof_map[0*dof3 + i + ((p - j) + 0*(p + 1))*p] = o++;
9934 for (
int j = 0; j < p; j++)
9935 for (
int i = 1; i < p; i++)
9937 dof_map[1*dof3 + i + ((p - 1 - j) + 0*p)*(p + 1)] = -1 - (o++);
9940 for (
int k = 1; k < p; k++)
9941 for (
int i = 0; i < p; i++)
9943 dof_map[0*dof3 + i + (0 + k*(p + 1))*p] = o++;
9945 for (
int k = 0; k < p; k++)
9946 for (
int i = 1; i < p; i++ )
9948 dof_map[2*dof3 + i + (0 + k*(p + 1))*(p + 1)] = o++;
9951 for (
int k = 1; k < p; k++)
9952 for (
int j = 0; j < p; j++)
9954 dof_map[1*dof3 + p + (j + k*p)*(p + 1)] = o++;
9956 for (
int k = 0; k < p; k++)
9957 for (
int j = 1; j < p; j++)
9959 dof_map[2*dof3 + p + (j + k*(p + 1))*(p + 1)] = o++;
9962 for (
int k = 1; k < p; k++)
9963 for (
int i = 0; i < p; i++)
9965 dof_map[0*dof3 + (p - 1 - i) + (p + k*(p + 1))*p] = -1 - (o++);
9967 for (
int k = 0; k < p; k++)
9968 for (
int i = 1; i < p; i++)
9970 dof_map[2*dof3 + (p - i) + (p + k*(p + 1))*(p + 1)] = o++;
9973 for (
int k = 1; k < p; k++)
9974 for (
int j = 0; j < p; j++)
9976 dof_map[1*dof3 + 0 + ((p - 1 - j) + k*p)*(p + 1)] = -1 - (o++);
9978 for (
int k = 0; k < p; k++)
9979 for (
int j = 1; j < p; j++)
9981 dof_map[2*dof3 + 0 + ((p - j) + k*(p + 1))*(p + 1)] = o++;
9984 for (
int j = 1; j < p; j++)
9985 for (
int i = 0; i < p; i++)
9987 dof_map[0*dof3 + i + (j + p*(p + 1))*p] = o++;
9989 for (
int j = 0; j < p; j++)
9990 for (
int i = 1; i < p; i++)
9992 dof_map[1*dof3 + i + (j + p*p)*(p + 1)] = o++;
9997 for (
int k = 1; k < p; k++)
9998 for (
int j = 1; j < p; j++)
9999 for (
int i = 0; i < p; i++)
10001 dof_map[0*dof3 + i + (j + k*(p + 1))*p] = o++;
10004 for (
int k = 1; k < p; k++)
10005 for (
int j = 0; j < p; j++)
10006 for (
int i = 1; i < p; i++)
10008 dof_map[1*dof3 + i + (j + k*p)*(p + 1)] = o++;
10011 for (
int k = 0; k < p; k++)
10012 for (
int j = 1; j < p; j++)
10013 for (
int i = 1; i < p; i++)
10015 dof_map[2*dof3 + i + (j + k*(p + 1))*(p + 1)] = o++;
10021 for (
int k = 0; k <= p; k++)
10022 for (
int j = 0; j <= p; j++)
10023 for (
int i = 0; i < p; i++)
10026 if ((idx = dof_map[o++]) < 0)
10028 dof2tk[idx = -1 - idx] = 3;
10037 for (
int k = 0; k <= p; k++)
10038 for (
int j = 0; j < p; j++)
10039 for (
int i = 0; i <= p; i++)
10042 if ((idx = dof_map[o++]) < 0)
10044 dof2tk[idx = -1 - idx] = 4;
10053 for (
int k = 0; k < p; k++)
10054 for (
int j = 0; j <= p; j++)
10055 for (
int i = 0; i <= p; i++)
10058 if ((idx = dof_map[o++]) < 0)
10060 dof2tk[idx = -1 - idx] = 5;
10073 const int p =
Order;
10075 #ifdef MFEM_THREAD_SAFE
10076 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10077 Vector shape_cz(p + 1), shape_oz(p);
10080 cbasis1d.
Eval(ip.
x, shape_cx);
10081 obasis1d.
Eval(ip.
x, shape_ox);
10082 cbasis1d.
Eval(ip.
y, shape_cy);
10083 obasis1d.
Eval(ip.
y, shape_oy);
10084 cbasis1d.
Eval(ip.
z, shape_cz);
10085 obasis1d.
Eval(ip.
z, shape_oz);
10089 for (
int k = 0; k <= p; k++)
10090 for (
int j = 0; j <= p; j++)
10091 for (
int i = 0; i < p; i++)
10094 if ((idx = dof_map[o++]) < 0)
10096 idx = -1 - idx, s = -1;
10102 shape(idx,0) = s*shape_ox(i)*shape_cy(j)*shape_cz(k);
10107 for (
int k = 0; k <= p; k++)
10108 for (
int j = 0; j < p; j++)
10109 for (
int i = 0; i <= p; i++)
10112 if ((idx = dof_map[o++]) < 0)
10114 idx = -1 - idx, s = -1;
10121 shape(idx,1) = s*shape_cx(i)*shape_oy(j)*shape_cz(k);
10125 for (
int k = 0; k < p; k++)
10126 for (
int j = 0; j <= p; j++)
10127 for (
int i = 0; i <= p; i++)
10130 if ((idx = dof_map[o++]) < 0)
10132 idx = -1 - idx, s = -1;
10140 shape(idx,2) = s*shape_cx(i)*shape_cy(j)*shape_oz(k);
10147 const int p =
Order;
10149 #ifdef MFEM_THREAD_SAFE
10150 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10151 Vector shape_cz(p + 1), shape_oz(p);
10152 Vector dshape_cx(p + 1), dshape_cy(p + 1), dshape_cz(p + 1);
10155 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
10156 obasis1d.
Eval(ip.
x, shape_ox);
10157 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
10158 obasis1d.
Eval(ip.
y, shape_oy);
10159 cbasis1d.
Eval(ip.
z, shape_cz, dshape_cz);
10160 obasis1d.
Eval(ip.
z, shape_oz);
10164 for (
int k = 0; k <= p; k++)
10165 for (
int j = 0; j <= p; j++)
10166 for (
int i = 0; i < p; i++)
10169 if ((idx = dof_map[o++]) < 0)
10171 idx = -1 - idx, s = -1;
10177 curl_shape(idx,0) = 0.;
10178 curl_shape(idx,1) = s*shape_ox(i)* shape_cy(j)*dshape_cz(k);
10179 curl_shape(idx,2) = -s*shape_ox(i)*dshape_cy(j)* shape_cz(k);
10182 for (
int k = 0; k <= p; k++)
10183 for (
int j = 0; j < p; j++)
10184 for (
int i = 0; i <= p; i++)
10187 if ((idx = dof_map[o++]) < 0)
10189 idx = -1 - idx, s = -1;
10195 curl_shape(idx,0) = -s* shape_cx(i)*shape_oy(j)*dshape_cz(k);
10196 curl_shape(idx,1) = 0.;
10197 curl_shape(idx,2) = s*dshape_cx(i)*shape_oy(j)* shape_cz(k);
10200 for (
int k = 0; k < p; k++)
10201 for (
int j = 0; j <= p; j++)
10202 for (
int i = 0; i <= p; i++)
10205 if ((idx = dof_map[o++]) < 0)
10207 idx = -1 - idx, s = -1;
10213 curl_shape(idx,0) = s* shape_cx(i)*dshape_cy(j)*shape_oz(k);
10214 curl_shape(idx,1) = -s*dshape_cx(i)* shape_cy(j)*shape_oz(k);
10215 curl_shape(idx,2) = 0.;
10220 const double ND_QuadrilateralElement::tk[8] =
10221 { 1.,0., 0.,1., -1.,0., 0.,-1. };
10228 cbasis1d(
poly1d.ClosedBasis(p, VerifyClosed(cp_type))),
10229 obasis1d(
poly1d.OpenBasis(p - 1, VerifyOpen(op_type))),
10230 dof_map(Dof), dof2tk(Dof)
10234 const int dof2 =
Dof/2;
10236 #ifndef MFEM_THREAD_SAFE
10247 for (
int i = 0; i < p; i++)
10249 dof_map[0*dof2 + i + 0*p] = o++;
10251 for (
int j = 0; j < p; j++)
10253 dof_map[1*dof2 + p + j*(p + 1)] = o++;
10255 for (
int i = 0; i < p; i++)
10257 dof_map[0*dof2 + (p - 1 - i) + p*p] = -1 - (o++);
10259 for (
int j = 0; j < p; j++)
10261 dof_map[1*dof2 + 0 + (p - 1 - j)*(p + 1)] = -1 - (o++);
10266 for (
int j = 1; j < p; j++)
10267 for (
int i = 0; i < p; i++)
10269 dof_map[0*dof2 + i + j*p] = o++;
10272 for (
int j = 0; j < p; j++)
10273 for (
int i = 1; i < p; i++)
10275 dof_map[1*dof2 + i + j*(p + 1)] = o++;
10281 for (
int j = 0; j <= p; j++)
10282 for (
int i = 0; i < p; i++)
10285 if ((idx = dof_map[o++]) < 0)
10287 dof2tk[idx = -1 - idx] = 2;
10296 for (
int j = 0; j < p; j++)
10297 for (
int i = 0; i <= p; i++)
10300 if ((idx = dof_map[o++]) < 0)
10302 dof2tk[idx = -1 - idx] = 3;
10315 const int p =
Order;
10317 #ifdef MFEM_THREAD_SAFE
10318 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10321 cbasis1d.
Eval(ip.
x, shape_cx);
10322 obasis1d.
Eval(ip.
x, shape_ox);
10323 cbasis1d.
Eval(ip.
y, shape_cy);
10324 obasis1d.
Eval(ip.
y, shape_oy);
10328 for (
int j = 0; j <= p; j++)
10329 for (
int i = 0; i < p; i++)
10332 if ((idx = dof_map[o++]) < 0)
10334 idx = -1 - idx, s = -1;
10340 shape(idx,0) = s*shape_ox(i)*shape_cy(j);
10344 for (
int j = 0; j < p; j++)
10345 for (
int i = 0; i <= p; i++)
10348 if ((idx = dof_map[o++]) < 0)
10350 idx = -1 - idx, s = -1;
10357 shape(idx,1) = s*shape_cx(i)*shape_oy(j);
10364 const int p =
Order;
10366 #ifdef MFEM_THREAD_SAFE
10367 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
10368 Vector dshape_cx(p + 1), dshape_cy(p + 1);
10371 cbasis1d.
Eval(ip.
x, shape_cx, dshape_cx);
10372 obasis1d.
Eval(ip.
x, shape_ox);
10373 cbasis1d.
Eval(ip.
y, shape_cy, dshape_cy);
10374 obasis1d.
Eval(ip.
y, shape_oy);
10378 for (
int j = 0; j <= p; j++)
10379 for (
int i = 0; i < p; i++)
10382 if ((idx = dof_map[o++]) < 0)
10384 idx = -1 - idx, s = -1;
10390 curl_shape(idx,0) = -s*shape_ox(i)*dshape_cy(j);
10393 for (
int j = 0; j < p; j++)
10394 for (
int i = 0; i <= p; i++)
10397 if ((idx = dof_map[o++]) < 0)
10399 idx = -1 - idx, s = -1;
10405 curl_shape(idx,0) = s*dshape_cx(i)*shape_oy(j);
10410 const double ND_TetrahedronElement::tk[18] =
10411 { 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,1.,0., -1.,0.,1., 0.,-1.,1. };
10413 const double ND_TetrahedronElement::c = 1./4.;
10423 const int pm1 = p - 1, pm2 = p - 2, pm3 = p - 3;
10425 #ifndef MFEM_THREAD_SAFE
10436 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
10441 for (
int i = 0; i < p; i++)
10446 for (
int i = 0; i < p; i++)
10451 for (
int i = 0; i < p; i++)
10456 for (
int i = 0; i < p; i++)
10461 for (
int i = 0; i < p; i++)
10466 for (
int i = 0; i < p; i++)
10473 for (
int j = 0; j <= pm2; j++)
10474 for (
int i = 0; i + j <= pm2; i++)
10476 double w = fop[i] + fop[j] + fop[pm2-i-j];
10482 for (
int j = 0; j <= pm2; j++)
10483 for (
int i = 0; i + j <= pm2; i++)
10485 double w = fop[i] + fop[j] + fop[pm2-i-j];
10491 for (
int j = 0; j <= pm2; j++)
10492 for (
int i = 0; i + j <= pm2; i++)
10494 double w = fop[i] + fop[j] + fop[pm2-i-j];
10500 for (
int j = 0; j <= pm2; j++)
10501 for (
int i = 0; i + j <= pm2; i++)
10503 double w = fop[i] + fop[j] + fop[pm2-i-j];
10511 for (
int k = 0; k <= pm3; k++)
10512 for (
int j = 0; j + k <= pm3; j++)
10513 for (
int i = 0; i + j + k <= pm3; i++)
10515 double w = iop[i] + iop[j] + iop[k] + iop[pm3-i-j-k];
10525 for (
int m = 0; m <
Dof; m++)
10528 const double *tm = tk + 3*dof2tk[m];
10536 for (
int k = 0; k <= pm1; k++)
10537 for (
int j = 0; j + k <= pm1; j++)
10538 for (
int i = 0; i + j + k <= pm1; i++)
10540 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
10541 T(o++, m) = s * tm[0];
10542 T(o++, m) = s * tm[1];
10543 T(o++, m) = s * tm[2];
10545 for (
int k = 0; k <= pm1; k++)
10546 for (
int j = 0; j + k <= pm1; j++)
10548 double s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
10549 T(o++, m) = s*((ip.
y - c)*tm[0] - (ip.
x - c)*tm[1]);
10550 T(o++, m) = s*((ip.
z - c)*tm[0] - (ip.
x - c)*tm[2]);
10552 for (
int k = 0; k <= pm1; k++)
10555 shape_y(pm1-k)*shape_z(k)*((ip.
z - c)*tm[1] - (ip.
y - c)*tm[2]);
10566 const int pm1 =
Order - 1;
10568 #ifdef MFEM_THREAD_SAFE
10569 const int p =
Order;
10570 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
10580 for (
int k = 0; k <= pm1; k++)
10581 for (
int j = 0; j + k <= pm1; j++)
10582 for (
int i = 0; i + j + k <= pm1; i++)
10584 double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
10585 u(n,0) = s; u(n,1) = 0.; u(n,2) = 0.; n++;
10586 u(n,0) = 0.; u(n,1) = s; u(n,2) = 0.; n++;
10587 u(n,0) = 0.; u(n,1) = 0.; u(n,2) = s; n++;
10589 for (
int k = 0; k <= pm1; k++)
10590 for (
int j = 0; j + k <= pm1; j++)
10592 double s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
10593 u(n,0) = s*(ip.
y - c); u(n,1) = -s*(ip.
x - c); u(n,2) = 0.; n++;
10594 u(n,0) = s*(ip.
z - c); u(n,1) = 0.; u(n,2) = -s*(ip.
x - c); n++;
10596 for (
int k = 0; k <= pm1; k++)
10598 double s = shape_y(pm1-k)*shape_z(k);
10599 u(n,0) = 0.; u(n,1) = s*(ip.
z - c); u(n,2) = -s*(ip.
y - c); n++;
10608 const int pm1 =
Order - 1;
10610 #ifdef MFEM_THREAD_SAFE
10611 const int p =
Order;
10612 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
10613 Vector dshape_x(p), dshape_y(p), dshape_z(p), dshape_l(p);
10623 for (
int k = 0; k <= pm1; k++)
10624 for (
int j = 0; j + k <= pm1; j++)
10625 for (
int i = 0; i + j + k <= pm1; i++)
10628 const double dx = (dshape_x(i)*shape_l(l) -
10629 shape_x(i)*dshape_l(l))*shape_y(j)*shape_z(k);
10630 const double dy = (dshape_y(j)*shape_l(l) -
10631 shape_y(j)*dshape_l(l))*shape_x(i)*shape_z(k);
10632 const double dz = (dshape_z(k)*shape_l(l) -
10633 shape_z(k)*dshape_l(l))*shape_x(i)*shape_y(j);
10635 u(n,0) = 0.; u(n,1) = dz; u(n,2) = -dy; n++;
10636 u(n,0) = -dz; u(n,1) = 0.; u(n,2) = dx; n++;
10637 u(n,0) = dy; u(n,1) = -dx; u(n,2) = 0.; n++;
10639 for (
int k = 0; k <= pm1; k++)
10640 for (
int j = 0; j + k <= pm1; j++)
10642 int i = pm1 - j - k;
10645 u(n,0) = shape_x(i)*(ip.
x - c)*shape_y(j)*dshape_z(k);
10646 u(n,1) = shape_x(i)*shape_y(j)*(ip.
y - c)*dshape_z(k);
10648 -((dshape_x(i)*(ip.
x - c) + shape_x(i))*shape_y(j)*shape_z(k) +
10649 (dshape_y(j)*(ip.
y - c) + shape_y(j))*shape_x(i)*shape_z(k));
10652 u(n,0) = -shape_x(i)*(ip.
x - c)*dshape_y(j)*shape_z(k);
10653 u(n,1) = (shape_x(i)*shape_y(j)*(dshape_z(k)*(ip.
z - c) + shape_z(k)) +
10654 (dshape_x(i)*(ip.
x - c) + shape_x(i))*shape_y(j)*shape_z(k));
10655 u(n,2) = -shape_x(i)*dshape_y(j)*shape_z(k)*(ip.
z - c);
10658 for (
int k = 0; k <= pm1; k++)
10662 u(n,0) = -((dshape_y(j)*(ip.
y - c) + shape_y(j))*shape_z(k) +
10663 shape_y(j)*(dshape_z(k)*(ip.
z - c) + shape_z(k)));
10668 Ti.
Mult(u, curl_shape);
10672 const double ND_TriangleElement::tk[8] =
10673 { 1.,0., -1.,1., 0.,-1., 0.,1. };
10675 const double ND_TriangleElement::c = 1./3.;
10685 const int pm1 = p - 1, pm2 = p - 2;
10687 #ifndef MFEM_THREAD_SAFE
10697 Vector shape_x(p), shape_y(p), shape_l(p);
10702 for (
int i = 0; i < p; i++)
10707 for (
int i = 0; i < p; i++)
10712 for (
int i = 0; i < p; i++)
10719 for (
int j = 0; j <= pm2; j++)
10720 for (
int i = 0; i + j <= pm2; i++)
10722 double w = iop[i] + iop[j] + iop[pm2-i-j];
10730 for (
int m = 0; m <
Dof; m++)
10733 const double *tm = tk + 2*dof2tk[m];
10740 for (
int j = 0; j <= pm1; j++)
10741 for (
int i = 0; i + j <= pm1; i++)
10743 double s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
10744 T(n++, m) = s * tm[0];
10745 T(n++, m) = s * tm[1];
10747 for (
int j = 0; j <= pm1; j++)
10750 shape_x(pm1-j)*shape_y(j)*((ip.
y - c)*tm[0] - (ip.
x - c)*tm[1]);
10761 const int pm1 =
Order - 1;
10763 #ifdef MFEM_THREAD_SAFE
10764 const int p =
Order;
10765 Vector shape_x(p), shape_y(p), shape_l(p);
10774 for (
int j = 0; j <= pm1; j++)
10775 for (
int i = 0; i + j <= pm1; i++)
10777 double s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
10778 u(n,0) = s; u(n,1) = 0; n++;
10779 u(n,0) = 0; u(n,1) = s; n++;
10781 for (
int j = 0; j <= pm1; j++)
10783 double s = shape_x(pm1-j)*shape_y(j);
10784 u(n,0) = s*(ip.
y - c);
10785 u(n,1) = -s*(ip.
x - c);
10795 const int pm1 =
Order - 1;
10797 #ifdef MFEM_THREAD_SAFE
10798 const int p =
Order;
10799 Vector shape_x(p), shape_y(p), shape_l(p);
10800 Vector dshape_x(p), dshape_y(p), dshape_l(p);
10809 for (
int j = 0; j <= pm1; j++)
10810 for (
int i = 0; i + j <= pm1; i++)
10813 const double dx = (dshape_x(i)*shape_l(l) -
10814 shape_x(i)*dshape_l(l)) * shape_y(j);
10815 const double dy = (dshape_y(j)*shape_l(l) -
10816 shape_y(j)*dshape_l(l)) * shape_x(i);
10822 for (
int j = 0; j <= pm1; j++)
10826 curlu(n++) = -((dshape_x(i)*(ip.
x - c) + shape_x(i)) * shape_y(j) +
10827 (dshape_y(j)*(ip.
y - c) + shape_y(j)) * shape_x(i));
10831 Ti.
Mult(curlu, curl2d);
10835 const double ND_SegmentElement::tk[1] = { 1. };
10840 obasis1d(
poly1d.OpenBasis(p - 1, VerifyOpen(op_type))),
10846 for (
int i = 0; i < p; i++)
10865 kv[0]->CalcShape(shape,
ijk[0], ip.
x);
10868 for (
int i = 0; i <=
Order; i++)
10870 sum += (shape(i) *=
weights(i));
10882 kv[0]->CalcDShape(grad,
ijk[0], ip.
x);
10884 double sum = 0.0, dsum = 0.0;
10885 for (
int i = 0; i <=
Order; i++)
10888 dsum += ( grad(i) *=
weights(i));
10892 add(sum, grad, -dsum*sum*sum,
shape_x, grad);
10902 for (
int o = 0, j = 0; j <=
Order; j++)
10904 const double sy =
shape_y(j);
10905 for (
int i = 0; i <=
Order; i++, o++)
10917 double sum, dsum[2];
10925 sum = dsum[0] = dsum[1] = 0.0;
10926 for (
int o = 0, j = 0; j <=
Order; j++)
10929 for (
int i = 0; i <=
Order; i++, o++)
10939 dsum[0] *= sum*sum;
10940 dsum[1] *= sum*sum;
10942 for (
int o = 0; o <
Dof; o++)
10944 dshape(o,0) = dshape(o,0)*sum -
u(o)*dsum[0];
10945 dshape(o,1) = dshape(o,1)*sum -
u(o)*dsum[1];
10957 for (
int o = 0, k = 0; k <=
Order; k++)
10959 const double sz =
shape_z(k);
10960 for (
int j = 0; j <=
Order; j++)
10962 const double sy_sz =
shape_y(j)*sz;
10963 for (
int i = 0; i <=
Order; i++, o++)
10976 double sum, dsum[3];
10986 sum = dsum[0] = dsum[1] = dsum[2] = 0.0;
10987 for (
int o = 0, k = 0; k <=
Order; k++)
10990 for (
int j = 0; j <=
Order; j++)
10992 const double sy_sz =
shape_y(j)* sz;
10993 const double dsy_sz =
dshape_y(j)* sz;
10994 const double sy_dsz =
shape_y(j)*dsz;
10995 for (
int i = 0; i <=
Order; i++, o++)
11007 dsum[0] *= sum*sum;
11008 dsum[1] *= sum*sum;
11009 dsum[2] *= sum*sum;
11011 for (
int o = 0; o <
Dof; o++)
11013 dshape(o,0) = dshape(o,0)*sum -
u(o)*dsum[0];
11014 dshape(o,1) = dshape(o,1)*sum -
u(o)*dsum[1];
11015 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
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)
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
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()
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)
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...
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
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