35#ifndef MFEM_THREAD_SAFE
43 MFEM_ABORT(
"method is not implemented for this class");
49 MFEM_ABORT(
"method is not implemented for this class");
55 MFEM_ABORT(
"method is not implemented for this class");
62 div_shape *= (1.0 / Trans.
Weight());
68 MFEM_ABORT(
"method is not implemented for this class");
78#ifdef MFEM_THREAD_SAFE
83 curl_shape *= (1.0 / Trans.
Weight());
89 curl_shape *= (1.0 / Trans.
Weight());
92 MFEM_ABORT(
"Invalid dimension, Dim = " <<
dim);
98 MFEM_ABORT(
"method is not overloaded");
104 MFEM_ABORT(
"method is not overloaded");
110 MFEM_ABORT(
"method is not overloaded");
116 MFEM_ABORT(
"method is not overloaded");
123 MFEM_ABORT(
"method is not overloaded");
129 MFEM_ABORT(
"method is not overloaded");
135 MFEM_ABORT(
"method is not overloaded");
141 mfem_error(
"FiniteElement::ProjectFromNodes() (vector) is not overloaded!");
147 MFEM_ABORT(
"method is not overloaded");
152 MFEM_ABORT(
"method is not implemented for this element");
158 MFEM_ABORT(
"method is not implemented for this element");
165 MFEM_ABORT(
"method is not implemented for this element");
172 MFEM_ABORT(
"method is not implemented for this element");
179 MFEM_ABORT(
"method is not implemented for this element");
196#ifdef MFEM_THREAD_SAFE
216 int size = (
dim*(
dim+1))/2;
222 for (
int nd = 0; nd <
dof; nd++)
224 Laplacian[nd] = hess(nd,0) + hess(nd,4) + hess(nd,5);
229 for (
int nd = 0; nd <
dof; nd++)
231 Laplacian[nd] = hess(nd,0) + hess(nd,2);
236 for (
int nd = 0; nd <
dof; nd++)
238 Laplacian[nd] = hess(nd,0);
248 int size = (
dim*(
dim+1))/2;
259 scale[1] = 2*Gij(0,1);
260 scale[2] = 2*Gij(0,2);
262 scale[3] = 2*Gij(1,2);
270 scale[1] = 2*Gij(0,1);
278 for (
int nd = 0; nd <
dof; nd++)
281 for (
int ii = 0; ii < size; ii++)
283 Laplacian[nd] += hess(nd,ii)*scale[ii];
323 int size = (
dim*(
dim+1))/2;
341 for (
int i = 0; i <
dim; i++)
343 for (
int j = 0; j <
dim; j++)
345 for (
int k = 0; k <
dim; k++)
347 for (
int l = 0; l <
dim; l++)
349 lhm(map[i*
dim+j],map[k*
dim+l]) += invJ(i,k)*invJ(j,l);
357 for (
int i = 0; i <
dim*
dim; i++) { mult[map[i]]++; }
362 Mult(hess, lhm, Hessian);
371#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
372 #pragma omp critical (DofToQuad)
378 if (d2q->
IntRule != &ir || d2q->
mode != mode) { d2q =
nullptr; }
382#ifdef MFEM_THREAD_SAFE
401 for (
int i = 0; i < nqpt; i++)
405 for (
int j = 0; j <
dof; j++)
407 d2q->
B[i+nqpt*j] = d2q->
Bt[j+
dof*i] = shape(j);
417 for (
int i = 0; i < nqpt; i++)
421 for (
int d = 0; d <
dim; d++)
423 for (
int j = 0; j <
dof; j++)
425 d2q->
B[i+nqpt*(d+
dim*j)] =
443 for (
int i = 0; i < nqpt; i++)
447 for (
int d = 0; d <
dim; d++)
449 for (
int j = 0; j <
dof; j++)
451 d2q->
G[i+nqpt*(d+
dim*j)] =
465 for (
int i = 0; i < nqpt; i++)
467 const IntegrationPoint &ip = ir.
IntPoint(i);
469 for (
int j = 0; j <
dof; j++)
471 d2q->
G[i+nqpt*j] = d2q->
Gt[j+
dof*i] = divshape(j);
482 for (
int i = 0; i < nqpt; i++)
484 const IntegrationPoint &ip = ir.
IntPoint(i);
486 for (
int d = 0; d <
cdim; d++)
488 for (
int j = 0; j <
dof; j++)
490 d2q->
G[i+nqpt*(d+
cdim*j)] =
491 d2q->
Gt[j+
dof*(i+nqpt*d)] = curlshape(j, d);
510 MFEM_ABORT(
"method is not implemented for this element");
530#ifdef MFEM_THREAD_SAFE
540 for (
int i = 0; i < fine_fe.
dof; i++)
545 for (
int j = 0; j <
dof; j++)
547 if (fabs(I(i,j) = shape(j)) < 1.0e-12)
573 Vector fine_shape(fs), coarse_shape(cs);
574 DenseMatrix fine_mass(fs), fine_coarse_mass(fs, cs);
592 fine_mass_inv.
Mult(fine_coarse_mass, I);
612 Vector fine_shape(fs), coarse_shape(cs);
613 DenseMatrix coarse_mass(cs), coarse_fine_mass(cs, fs);
641 coarse_mass_inv.
Mult(coarse_fine_mass, R);
647 R *= 1.0 / Trans.
Weight();
651void NodalFiniteElement::CreateLexicographicFullMap(
const IntegrationRule &ir)
663 for (
int i = 0; i < nqpt; i++)
665 for (
int d = 0; d < b_dim; d++)
667 for (
int j = 0; j <
dof; j++)
669 const double val = d2q.
B[i + nqpt*(d+b_dim*
lex_ordering[j])];
670 d2q_new->B[i+nqpt*(d+b_dim*j)] = val;
671 d2q_new->Bt[j+
dof*(i+nqpt*d)] = val;
676 const int g_dim = [
this]()
687 for (
int i = 0; i < nqpt; i++)
689 for (
int d = 0; d < g_dim; d++)
691 for (
int j = 0; j <
dof; j++)
693 const double val = d2q.
G[i + nqpt*(d+g_dim*
lex_ordering[j])];
694 d2q_new->G[i+nqpt*(d+g_dim*j)] = val;
695 d2q_new->Gt[j+
dof*(i+nqpt*d)] = val;
710 if (d2q.
IntRule == &ir && d2q.
mode == mode) {
return d2q; }
719 CreateLexicographicFullMap(ir);
731 for (
int i = 0; i <
dof; i++)
741 for (
int j = 0; j < fe.
GetDof(); j++)
743 curl(i,j) = w * curl_shape(j,0);
754 trans.Transform(p0, x);
761 trans.InverseJacobian().Mult(v, x);
770#ifdef MFEM_THREAD_SAFE
779 for (
int j = 0; j <
dof; j++)
799 for (
int i = 0; i <
dof; i++)
805 dofs(i) = coeff.
Eval(Trans, ip);
808 dofs(i) *= Trans.
Weight();
819 for (
int i = 0; i <
dof; i++)
823 vc.
Eval (x, Trans, ip);
828 for (
int j = 0; j < x.
Size(); j++)
830 dofs(
dof*j+i) = x(j);
842 for (
int k = 0; k <
dof; k++)
847 for (
int r = 0; r < MQ.
Height(); r++)
849 for (
int d = 0; d < MQ.
Width(); d++)
851 dofs(k+
dof*(d+MQ.
Width()*r)) = MQ(r,d);
867 for (
int k = 0; k <
dof; k++)
870 for (
int j = 0; j < shape.
Size(); j++)
872 I(k,j) = (fabs(shape(j)) < 1e-12) ? 0.0 : shape(j);
878 for (
int k = 0; k <
dof; k++)
886 for (
int j = 0; j < shape.
Size(); j++)
888 I(k,j) = (fabs(shape(j)) < 1e-12) ? 0.0 : shape(j);
899 for (
int k = 0; k <
dof; k++)
926 for (
int k = 0; k <
dof; k++)
932 Mult(dshape, Jinv, grad_k);
937 for (
int j = 0; j < grad_k.Height(); j++)
938 for (
int d = 0; d <
dim; d++)
940 grad(k+d*
dof,j) = grad_k(j,d);
953 for (
int k = 0; k <
dof; k++)
961 for (
int j = 0; j < div_shape.
Size(); j++)
963 div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j)/detJ;
968 for (
int j = 0; j < div_shape.
Size(); j++)
970 div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j);
980 MFEM_ASSERT(dofs.
Size() == ncomp *
dof,
"Wrong input size.");
983 for (
int i = 0; i <
dof; i++)
985 for (
int c = 0; c < ncomp; c++)
994 int Do,
int O,
int M,
int F)
1008void VectorFiniteElement::CalcShape(
1011 mfem_error(
"Error: Cannot use scalar CalcShape(...) function with\n"
1012 " VectorFiniteElements!");
1015void VectorFiniteElement::CalcDShape(
1016 const IntegrationPoint &ip, DenseMatrix &dshape)
const
1018 mfem_error(
"Error: Cannot use scalar CalcDShape(...) function with\n"
1019 " VectorFiniteElements!");
1051 MFEM_ABORT(
"Invalid dimension, Dim = " <<
dim);
1055 MFEM_ABORT(
"Invalid MapType = " <<
map_type);
1063#ifdef MFEM_THREAD_SAFE
1068 shape *= (1.0 / Trans.
Weight());
1075#ifdef MFEM_THREAD_SAFE
1088 MFEM_ASSERT(vc.
GetVDim() == sdim,
"");
1090 const bool square_J = (
dim == sdim);
1092 for (
int k = 0; k <
dof; k++)
1098 if (!square_J) { dofs(k) /= Trans.
Weight(); }
1107 const bool square_J = (
dim == sdim);
1109 for (
int k = 0; k <
dof; k++)
1114 &vc[k*sdim], nk + d2n[k]*
dim);
1115 if (!square_J) { dofs(k) /= Trans.
Weight(); }
1126 MFEM_ASSERT(mc.
GetWidth() == sdim,
"");
1127 const bool square_J = (
dim == sdim);
1132 for (
int k = 0; k <
dof; k++)
1138 if (!square_J) { nk_phys /= T.
Weight(); }
1139 MQ.
Mult(nk_phys, dofs_k);
1140 for (
int r = 0; r < MQ.
Height(); r++)
1142 dofs(k+
dof*r) = dofs_k(r);
1158 for (
int k = 0; k <
dof; k++)
1170 for (
int d = 0; d <
dim; d++)
1176 for (
int j = 0; j < shape.
Size(); j++)
1179 if (fabs(s) < 1e-12)
1185 for (
int d = 0; d < sdim; d++)
1187 I(k,j+d*shape.
Size()) = s*vk[d];
1198 const bool square_J = (
dim == sdim);
1201 for (
int k = 0; k <
dof; k++)
1213 if (!square_J) { vshapenk /= Trans.
Weight(); }
1214 for (
int j=0; j<vshapenk.
Size(); j++)
1216 I(k,j) = vshapenk(j);
1228 mfem_error(
"VectorFiniteElement::ProjectGrad_RT works only in 2D!");
1236 for (
int k = 0; k <
dof; k++)
1239 tk[0] = nk[d2n[k]*
dim+1];
1240 tk[1] = -nk[d2n[k]*
dim];
1241 dshape.
Mult(tk, grad_k);
1242 for (
int j = 0; j < grad_k.
Size(); j++)
1244 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
1253#ifdef MFEM_THREAD_SAFE
1266 for (
int k = 0; k <
dof; k++)
1280 for (
int j = 0; j < curl_k.
Size(); j++)
1282 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
1295 for (
int k = 0; k <
dof; k++)
1298 curl_shape.
Mult(nk + d2n[k]*
dim, curl_k);
1299 for (
int j = 0; j < curl_k.
Size(); j++)
1301 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
1313 for (
int k = 0; k <
dof; k++)
1327 for (
int k = 0; k <
dof; k++)
1342 MFEM_ASSERT(mc.
GetWidth() == sdim,
"");
1347 for (
int k = 0; k <
dof; k++)
1353 MQ.
Mult(tk_phys, dofs_k);
1354 for (
int r = 0; r < MQ.
Height(); r++)
1356 dofs(k+
dof*r) = dofs_k(r);
1372 for (
int k = 0; k <
dof; k++)
1384 for (
int d = 0; d < sdim; d++)
1390 for (
int j = 0; j < shape.
Size(); j++)
1393 if (fabs(s) < 1e-12)
1399 for (
int d = 0; d < sdim; d++)
1401 I(k, j + d*shape.
Size()) = s*vk[d];
1414 for (
int k = 0; k <
dof; k++)
1426 for (
int j=0; j<vshapetk.
Size(); j++)
1428 I(k, j) = vshapetk(j);
1444 for (
int k = 0; k <
dof; k++)
1447 dshape.
Mult(tk + d2t[k]*
dim, grad_k);
1448 for (
int j = 0; j < grad_k.
Size(); j++)
1450 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
1465 DenseMatrix fine_mass(fs), fine_coarse_mass(fs, cs);
1466 const int ir_order =
1482 for (
int k=0; k<fs; ++k)
1484 for (
int j=0; j<cs; ++j)
1487 for (
int d1=0; d1<
dim; ++d1)
1489 for (
int d2=0; d2<
dim; ++d2)
1491 Mkj += w*fine_shape(k,d1)*adjJ(d2,d1)*coarse_shape(j,d2);
1494 fine_coarse_mass(k,j) += (fabs(Mkj) < 1e-12) ? 0.0 : Mkj;
1499 fine_mass_inv.
Mult(fine_coarse_mass, I);
1513#ifdef MFEM_THREAD_SAFE
1523 for (
int k = 0; k <
dof; k++)
1534 for (
int i = 0; i <
dim; i++)
1536 Ikj +=
vshape(j, i) * vk[i];
1538 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
1553 DenseMatrix fine_mass(fs), fine_coarse_mass(fs, cs);
1554 const int ir_order =
1569 for (
int k=0; k<fs; ++k)
1571 for (
int j=0; j<cs; ++j)
1574 for (
int d1=0; d1<
dim; ++d1)
1576 for (
int d2=0; d2<
dim; ++d2)
1578 Mkj += ip.
weight*fine_shape(k,d1)*J(d1,d2)*coarse_shape(j,d2);
1581 fine_coarse_mass(k,j) += (fabs(Mkj) < 1e-12) ? 0.0 : Mkj;
1586 fine_mass_inv.
Mult(fine_coarse_mass, I);
1598#ifdef MFEM_THREAD_SAFE
1608 for (
int k = 0; k <
dof; k++)
1619 for (
int i = 0; i <
dim; i++)
1621 Ikj +=
vshape(j, i) * vk[i];
1623 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
1636#ifdef MFEM_THREAD_SAFE
1643 for (
int j = 0; j <
dof; j++)
1652 for (
int k = 0; k <
dof; k++)
1655 for (
int d = 0; d <
dim; d++)
1657 R_jk +=
vshape(k,d)*pt_data[d];
1679#ifdef MFEM_THREAD_SAFE
1685 for (
int j = 0; j <
dof; j++)
1692 Jinv.
Mult(tk+
dim*d2t[j], pt_data);
1693 for (
int k = 0; k <
dof; k++)
1696 for (
int d = 0; d <
dim; d++)
1698 R_jk +=
vshape(k,d)*pt_data[d];
1714 : etype(etype), auxiliary_basis(NULL), scale_integrated(false)
1723 for (
int i = 0; i <=
p; i++)
1737 for (
int i = 0; i <=
p; i++)
1739 for (
int j = 0; j < i; j++)
1741 real_t xij = x(i) - x(j);
1746 for (
int i = 0; i <=
p; i++)
1753 for (
int i = 0; i <
p; i++)
1757 mfem_error(
"Poly_1D::Basis::Basis : nodes are not increasing!");
1767 auxiliary_basis =
new Basis(
1789 int i, k,
p = x.Size() - 1;
1799 for (k = 0; k <
p; k++)
1801 if (y >= (x(k) + x(k+1))/2)
1807 for (i = k+1; i <=
p; i++)
1814 l = lk * (y - x(k));
1816 for (i = 0; i < k; i++)
1818 u(i) = l * w(i) / (y - x(i));
1821 for (i++; i <=
p; i++)
1823 u(i) = l * w(i) / (y - x(i));
1831 auxiliary_basis->Eval(y, u_aux, d_aux);
1832 EvalIntegrated(d_aux,
u);
1851 int i, k,
p = x.Size() - 1;
1852 real_t l, lp, lk, sk, si;
1862 for (k = 0; k <
p; k++)
1864 if (y >= (x(k) + x(k+1))/2)
1870 for (i = k+1; i <=
p; i++)
1877 l = lk * (y - x(k));
1880 for (i = 0; i < k; i++)
1882 si = 1.0/(y - x(i));
1884 u(i) = l * si * w(i);
1887 for (i++; i <=
p; i++)
1889 si = 1.0/(y - x(i));
1891 u(i) = l * si * w(i);
1895 for (i = 0; i < k; i++)
1897 d(i) = (lp * w(i) -
u(i))/(y - x(i));
1900 for (i++; i <=
p; i++)
1902 d(i) = (lp * w(i) -
u(i))/(y - x(i));
1910 auxiliary_basis->Eval(y, u_aux, d_aux, d2_aux);
1911 EvalIntegrated(d_aux,
u);
1912 EvalIntegrated(d2_aux,d);
1922 "Basis::Eval with second order derivatives not implemented for"
1923 " etype = " << etype);
1936 int i, k,
p = x.Size() - 1;
1937 real_t l, lp, lp2, lk, sk, si, sk2;
1948 for (k = 0; k <
p; k++)
1950 if (y >= (x(k) + x(k+1))/2)
1956 for (i = k+1; i <=
p; i++)
1963 l = lk * (y - x(k));
1967 for (i = 0; i < k; i++)
1969 si = 1.0/(y - x(i));
1972 u(i) = l * si * w(i);
1975 for (i++; i <=
p; i++)
1977 si = 1.0/(y - x(i));
1980 u(i) = l * si * w(i);
1983 lp2 = lp * sk + l * sk2 + sk * lk;
1985 for (i = 0; i < k; i++)
1987 d(i) = (lp * w(i) -
u(i))/(y - x(i));
1988 d2(i) = (lp2 * w(i) - 2 * d(i))/(y - x(i));
1991 d2(k) = sk2 *
u(k) + sk * d(k);
1992 for (i++; i <=
p; i++)
1994 d(i) = (lp * w(i) -
u(i))/(y - x(i));
1995 d2(i) = (lp2 * w(i) - 2 * d(i))/(y - x(i));
2003 MFEM_ABORT(
"Integrated basis must be evaluated with EvalIntegrated");
2012 "EvalIntegrated is only valid for Integrated basis type");
2013 int p = d_aux_.
Size() - 1;
2017 for (
int j=1; j<
p; ++j)
2019 u[j] =
u[j-1] - d_aux_[j];
2024 if (scale_integrated)
2026 Vector &aux_nodes = auxiliary_basis->x;
2027 for (
int j=0; j<aux_nodes.
Size()-1; ++j)
2029 u[j] *= aux_nodes[j+1] - aux_nodes[j];
2036 scale_integrated = scale_integrated_;
2041 delete auxiliary_basis;
2049 for (
int i = 0; i <=
p; i++)
2051 binom(i,0) = binom(i,i) = 1;
2052 for (
int j = 1; j < i; j++)
2054 binom(i,j) = binom(i-1,j) + binom(i-1,j-1);
2063 for (
int i = 0; i <=
p; i++)
2066 real_t s = sin(M_PI_2*(i + 0.5)/(
p + 1));
2075 for (
int n = 1; n <=
p; n++)
2086 for (
int n = 1; n <=
p; n++)
2106 for (i = 1; i <
p; i++)
2113 for (i--; i > 0; i--)
2134 const real_t xpy = x + y, ptx =
p*x;
2137 for (i = 1; i <
p; i++)
2139 d[i] =
b[i]*z*(i*xpy - ptx);
2146 for (i--; i > 0; i--)
2168 const real_t xpy = x + y, ptx =
p*x;
2171 for (i = 1; i <
p; i++)
2173 d[i] =
b[i]*z*(i*xpy - ptx);
2178 for (i--; i > 0; i--)
2200 for (i = 1; i <
p; i++)
2207 for (i--; i > 0; i--)
2229 for (i = 1; i <
p; i++)
2236 for (i--; i > 0; i--)
2238 u[i] *= (
p - i) * z;
2251 if (
p == 0) {
return; }
2252 u[1] = z = 2.*x - 1.;
2253 for (
int n = 1; n <
p; n++)
2255 u[n+1] = ((2*n + 1)*z*
u[n] - n*
u[n-1])/(n + 1);
2268 if (
p == 0) {
return; }
2269 u[1] = z = 2.*x - 1.;
2271 for (
int n = 1; n <
p; n++)
2273 u[n+1] = ((2*n + 1)*z*
u[n] - n*
u[n-1])/(n + 1);
2274 d[n+1] = (4*n + 2)*
u[n] + d[n-1];
2278void Poly_1D::CalcChebyshev(
const int p,
const real_t x,
real_t *
u)
2285 if (
p == 0) {
return; }
2286 u[1] = z = 2.*x - 1.;
2287 for (
int n = 1; n <
p; n++)
2289 u[n+1] = 2*z*
u[n] -
u[n-1];
2306 if (
p == 0) {
return; }
2307 u[1] = z = 2.*x - 1.;
2309 for (
int n = 1; n <
p; n++)
2311 u[n+1] = 2*z*
u[n] -
u[n-1];
2312 d[n+1] = (n + 1)*(z*d[n]/n + 2*
u[n]);
2332 if (
p == 0) {
return; }
2333 u[1] = z = 2.*x - 1.;
2336 for (
int n = 1; n <
p; n++)
2338 u[n+1] = 2*z*
u[n] -
u[n-1];
2339 d[n+1] = (n + 1)*(z*d[n]/n + 2*
u[n]);
2340 dd[n+1] = (n + 1)*(2.*(n + 1)*d[n] + z*dd[n])/n;
2351#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
2352 #pragma omp critical (Poly1DGetPoints)
2355 std::pair<int, int> key(btype,
p);
2356 auto it = points_container.find(key);
2357 if (it == points_container.end())
2359 it = points_container.emplace(key,
new Array<real_t>(
p + 1, h_mt)).first;
2360 val = it->second.get();
2366 val = it->second.get();
2377#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
2378 #pragma omp critical (Poly1DGetBasis)
2381 std::pair<int, int> key(btype,
p);
2382 auto it = bases_container.find(key);
2383 if (it == bases_container.end())
2389 it = bases_container
2393 val = it->second.get();
2413 for (
int i = 1; i <
p; i++)
2421 const int p1 =
p + 1;
2432 for (
int i = 1; i <
p; i++)
2436 for (
int i = 1; i <
p; i++)
2440 for (
int i = 1; i <
p; i++)
2444 for (
int i = 1; i <
p; i++)
2450 for (
int j = 1; j <
p; j++)
2452 for (
int i = 1; i <
p; i++)
2461 const int p1 =
p + 1;
2465 dof_map[0 + (0 + 0*p1)*p1] = 0;
2477 for (
int i = 1; i <
p; i++)
2479 dof_map[i + (0 + 0*p1)*p1] = o++;
2481 for (
int i = 1; i <
p; i++)
2485 for (
int i = 1; i <
p; i++)
2489 for (
int i = 1; i <
p; i++)
2491 dof_map[0 + (i + 0*p1)*p1] = o++;
2493 for (
int i = 1; i <
p; i++)
2497 for (
int i = 1; i <
p; i++)
2501 for (
int i = 1; i <
p; i++)
2505 for (
int i = 1; i <
p; i++)
2509 for (
int i = 1; i <
p; i++)
2511 dof_map[0 + (0 + i*p1)*p1] = o++;
2513 for (
int i = 1; i <
p; i++)
2517 for (
int i = 1; i <
p; i++)
2521 for (
int i = 1; i <
p; i++)
2527 for (
int j = 1; j <
p; j++)
2529 for (
int i = 1; i <
p; i++)
2531 dof_map[i + ((
p-j) + 0*p1)*p1] = o++;
2534 for (
int j = 1; j <
p; j++)
2536 for (
int i = 1; i <
p; i++)
2538 dof_map[i + (0 + j*p1)*p1] = o++;
2541 for (
int j = 1; j <
p; j++)
2543 for (
int i = 1; i <
p; i++)
2548 for (
int j = 1; j <
p; j++)
2550 for (
int i = 1; i <
p; i++)
2555 for (
int j = 1; j <
p; j++)
2557 for (
int i = 1; i <
p; i++)
2559 dof_map[0 + ((
p-i) + j*p1)*p1] = o++;
2562 for (
int j = 1; j <
p; j++)
2564 for (
int i = 1; i <
p; i++)
2571 for (
int k = 1; k <
p; k++)
2573 for (
int j = 1; j <
p; j++)
2575 for (
int i = 1; i <
p; i++)
2577 dof_map[i + (j + k*p1)*p1] = o++;
2584 MFEM_ABORT(
"invalid dimension: " << dims);
2595 MFEM_ABORT(
"invalid DofMapType: " << dmtype);
2607#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
2608 #pragma omp critical (DofToQuad)
2611 for (
int i = 0; i < dof2quad_array.
Size(); i++)
2613 d2q = dof2quad_array[i];
2614 if (d2q->
IntRule != &ir || d2q->
mode != mode) { d2q =
nullptr; }
2630 Vector val(ndof), grad(ndof);
2631 for (
int i = 0; i < nqpt; i++)
2636 for (
int j = 0; j < ndof; j++)
2638 d2q->
B[i+nqpt*j] = d2q->
Bt[j+ndof*i] = val(j);
2639 d2q->
G[i+nqpt*j] = d2q->
Gt[j+ndof*i] = grad(j);
2642 dof2quad_array.
Append(d2q);
2689 internal::GetTensorFaceMap(
dim,
order, face_id, face_map);
2702 obasis1d(
poly1d.GetBasis(
p - 1, VerifyOpen(obtype)))
2704 MFEM_VERIFY(dims > 1,
"Constructor for VectorTensorFiniteElement with both "
2705 "open and closed bases is not valid for 1D elements.");
2717 obasis1d(
poly1d.GetBasis(
p, VerifyOpen(obtype)))
2719 MFEM_VERIFY(dims == 1,
"Constructor for VectorTensorFiniteElement without "
2720 "closed basis is only valid for 1D elements.");
2725 for (
int i = 0; i < dof2quad_array_open.Size(); i++)
2727 delete dof2quad_array_open[i];
void SetSize(int m, int n)
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
int Size() const
Return the logical size of the array.
int Append(const T &el)
Append element 'el' to array, resize if necessary.
T * HostWrite()
Shortcut for mfem::Write(a.GetMemory(), a.Size(), false).
static int GetQuadrature1D(int b_type)
Get the corresponding Quadrature1D constant, when that makes sense; otherwise return Quadrature1D::In...
static int Check(int b_type)
If the input does not represents a valid BasisType, abort with an error; otherwise return the input.
@ GaussLobatto
Closed type.
@ Positive
Bernstein polynomials.
@ IntegratedGLL
Integrated GLL indicator functions.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
void Factor()
Factor the current DenseMatrix, *a.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication with the inverse of dense matrix.
Data type dense matrix using column-major storage.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
void Threshold(real_t eps)
Replace small entries, abs(a_ij) <= eps, with zero.
void SetRow(int r, const real_t *row)
void GetColumnReference(int c, Vector &col)
real_t InnerProduct(const real_t *x, const real_t *y) const
Compute y^t A x.
real_t * Data() const
Returns the matrix data array.
real_t * GetData() const
Returns the matrix data array.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
void Invert()
Replaces the current matrix with its inverse.
void InvRightScaling(const Vector &s)
InvRightScaling: this = this * diag(1./s);.
void GetColumn(int c, Vector &col) const
real_t FNorm2() const
Compute the square of the Frobenius norm of the matrix.
Structure representing the matrices/tensors needed to evaluate (in reference space) the values,...
Mode mode
Describes the contents of the B, Bt, G, and Gt arrays, see Mode.
const IntegrationRule * IntRule
IntegrationRule that defines the quadrature points at which the basis functions of the FE are evaluat...
Array< real_t > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Mode
Type of data stored in the arrays B, Bt, G, and Gt.
@ FULL
Full multidimensional representation which does not use tensor product structure. The ordering of the...
@ LEXICOGRAPHIC_FULL
Full multidimensional representation which does not use tensor product structure. The ordering of the...
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Array< real_t > B
Basis functions evaluated at quadrature points.
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number.
Array< real_t > Gt
Transpose of G.
const class FiniteElement * FE
The FiniteElement that created and owns this object.
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Array< real_t > Bt
Transpose of B.
Abstract class for all finite elements.
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...
int dof
Number of degrees of freedom.
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &Hessian) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
virtual void ProjectMatrixCoefficient(MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Given a matrix coefficient and a transformation, compute an approximation ("projection") in the local...
virtual ~FiniteElement()
Deconstruct the FiniteElement.
virtual void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const
Compute the discrete divergence matrix from the given FiniteElement onto 'this' FiniteElement....
virtual void ProjectFromNodes(Vector &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector of values at the finite element nodes and a transformation, compute its projection (ap...
virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const
Get the dofs associated with the given face. *dofs is set to an internal array of the local dofc on t...
virtual const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
int GetRangeDim() const
Returns the vector dimension for vector-valued finite elements, which is also the dimension of the in...
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
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...
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 ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Compute the discrete gradient matrix from the given FiniteElement onto 'this' FiniteElement....
virtual void GetFaceMap(const int face_id, Array< int > &face_map) const
Return the mapping from lexicographic face DOFs to lexicographic element DOFs for the given local fac...
int GetDim() const
Returns the reference space dimension for the finite element.
int vdim
Vector dimension of vector-valued basis functions.
void CalcPhysHessian(ElementTransformation &Trans, DenseMatrix &Hessian) const
Evaluate the Hessian of all shape functions of a scalar finite element in physical space at the given...
FiniteElement(int D, Geometry::Type G, int Do, int O, int F=FunctionSpace::Pk)
Construct FiniteElement with given.
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
int orders[Geometry::MaxDim]
Anisotropic orders.
virtual void GetTransferMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const
Return interpolation matrix, I, which maps dofs from a coarse element, fe, to the fine dofs on this f...
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
virtual void GetLocalRestriction(ElementTransformation &Trans, DenseMatrix &R) const
Return a local restriction matrix R (Dof x Dof) mapping fine dofs to coarse dofs.
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...
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
int cdim
Dimension of curl for vector-valued basis functions.
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
@ DIV
Implements CalcDivShape methods.
@ NONE
No derivatives implemented.
@ CURL
Implements CalcCurlShape methods.
@ GRAD
Implements CalcDShape methods.
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Geometry::Type geom_type
Geometry::Type of the reference element.
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...
Array< DofToQuad * > dof2quad_array
Container for all DofToQuad objects created by the FiniteElement.
void CalcPhysLaplacian(ElementTransformation &Trans, Vector &Laplacian) const
Evaluate the Laplacian of all shape functions of a scalar finite element in physical space at the giv...
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...
virtual void ProjectCurl(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Compute the discrete curl matrix from the given FiniteElement onto 'this' FiniteElement....
virtual 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...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
int order
Order/degree of the shape functions.
void CalcPhysLinLaplacian(ElementTransformation &Trans, Vector &Laplacian) const
Evaluate the Laplacian of all shape functions of a scalar finite element in physical space at the giv...
int dim
Dimension of reference space.
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 ...
Describes the function space on each element.
const IntegrationPoint & GetCenter(int GeomType) const
Return the center of the given Geometry::Type, GeomType.
static bool CheckPoint(int GeomType, const IntegrationPoint &ip)
Check if the given point is inside the given reference element.
Class for integration point with weight.
void Get(real_t *p, const int dim) const
void Set(const real_t *p, const int dim)
void Set3(const real_t x1, const real_t x2, const real_t x3)
Class for an integration rule - an Array of IntegrationPoint.
int GetNPoints() const
Returns the number of the points in the integration rule.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Base class for Matrix Coefficients that optionally depend on time and space.
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the matrix coefficient in the element described by T at the point ip, storing the result in ...
int GetWidth() const
Get the width of the matrix.
int GetHeight() const
Get the height of the matrix.
Class for standard nodal finite elements.
void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const override
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
void ProjectCurl_2D(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
void ReorderLexToNative(int ncomp, Vector &dofs) const
void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const override
Compute the discrete divergence matrix from the given FiniteElement onto 'this' FiniteElement....
void GetLocalRestriction(ElementTransformation &Trans, DenseMatrix &R) const override
Return a local restriction matrix R (Dof x Dof) mapping fine dofs to coarse dofs.
void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const override
Compute the discrete gradient matrix from the given FiniteElement onto 'this' FiniteElement....
void ProjectMatrixCoefficient(MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const override
Given a matrix coefficient and a transformation, compute an approximation ("projection") in the local...
Array< int > lex_ordering
const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const override
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
void GetFaceMap(const int face_id, Array< int > &face_map) const override
Return the mapping from lexicographic face DOFs to lexicographic element DOFs for the given local fac...
NodalTensorFiniteElement(const int dims, const int p, const int btype, const DofMapType dmtype)
const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const override
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
void SetMapType(const int map_type_) override
Set the FiniteElement::MapType of the element to either VALUE or INTEGRAL. Also sets the FiniteElemen...
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Class for evaluating 1D nodal, positive (Bernstein), or integrated (Gerritsma) bases.
void Eval(const real_t x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
bool IsIntegratedType() const
Returns true if the basis is "integrated", false otherwise.
void ScaleIntegrated(bool scale_integrated_)
Set whether the "integrated" basis should be scaled by the subcell sizes. Has no effect for non-integ...
void EvalIntegrated(const Vector &d, Vector &i) const
Evaluate the "integrated" basis type using pre-computed closed basis derivatives.
Basis(const int p, const real_t *nodes, EvalType etype=Barycentric)
Create a nodal or positive (Bernstein) basis of degree p.
static void CalcDBinomTerms(const int p, const real_t x, const real_t y, real_t *d)
Compute the derivatives (w.r.t. x) of the terms in the expansion of the binomial (x + y)^p assuming t...
const real_t * GetPoints(const int p, const int btype, bool on_device=false)
Get the coordinates of the points of the given BasisType, btype.
static const int * Binom(const int p)
Get a pointer to an array containing the binomial coefficients "pchoose k" for k=0,...
Basis & GetBasis(const int p, const int btype)
Get a Poly_1D::Basis object of the given degree and BasisType, btype.
const Array< real_t > * GetPointsArray(const int p, const int btype)
Get the coordinates of the points of the given BasisType, btype.
static void CalcDyBinomTerms(const int p, const real_t x, const real_t y, real_t *d)
Compute the derivatives (w.r.t. y) of the terms in the expansion of the binomial (x + y)^p....
static void CalcLegendre(const int p, const real_t x, real_t *u)
static void CalcBernstein(const int p, const real_t x, real_t *u)
Compute the values of the Bernstein basis functions of order p at coordinate x and store the results ...
EvalType
One-dimensional basis evaluation type.
@ ChangeOfBasis
Use change of basis, O(p^2) Evals.
@ Integrated
Integrated indicator functions (cf. Gerritsma)
@ Positive
Fast evaluation of Bernstein polynomials.
@ Barycentric
Use barycentric Lagrangian interpolation, O(p) Evals.
static void ChebyshevPoints(const int p, real_t *x)
Compute the points for the Chebyshev polynomials of order p and place them in the already allocated x...
static void CalcBasis(const int p, const real_t x, real_t *u)
Evaluate the values of a hierarchical 1D basis at point x hierarchical = k-th basis function is degre...
static void CalcBinomTerms(const int p, const real_t x, const real_t y, real_t *u)
Compute the p terms in the expansion of the binomial (x + y)^p and store them in the already allocate...
static void CalcDxBinomTerms(const int p, const real_t x, const real_t y, real_t *d)
Compute the derivatives (w.r.t. x) of the terms in the expansion of the binomial (x + y)^p....
static void GivePolyPoints(const int np, real_t *pts, const int type)
Class for finite elements with basis functions that return scalar values.
void NodalLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const ScalarFiniteElement &fine_fe) const
Get the matrix I that defines nodal interpolation between this element and the refined element fine_f...
virtual void SetMapType(int M)
Set the FiniteElement::MapType of the element to either VALUE or INTEGRAL. Also sets the FiniteElemen...
void ScalarLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const ScalarFiniteElement &fine_fe) const
Get matrix I "Interpolation" defined through local L2-projection in the space defined by the fine_fe.
void ScalarLocalL2Restriction(ElementTransformation &Trans, DenseMatrix &R, const ScalarFiniteElement &coarse_fe) const
Get restriction matrix R defined through local L2-projection in the space defined by the coarse_fe.
static const DofToQuad & GetTensorDofToQuad(const FiniteElement &fe, const IntegrationRule &ir, DofToQuad::Mode mode, const Poly_1D::Basis &basis, bool closed, Array< DofToQuad * > &dof2quad_array)
TensorBasisElement(const int dims, const int p, const int btype, const DofMapType dmtype)
Base class for vector Coefficients that optionally depend on time and space.
int GetVDim()
Returns dimension of the vector.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
Intermediate class for finite elements whose basis functions return vector values.
void ProjectGrad_RT(const real_t *nk, const Array< int > &d2n, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
void LocalRestriction_RT(const real_t *nk, const Array< int > &d2n, ElementTransformation &Trans, DenseMatrix &R) const
void Project_ND(const real_t *tk, const Array< int > &d2t, VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Project a vector coefficient onto the ND basis functions.
void Project_RT(const real_t *nk, const Array< int > &d2n, VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Project a vector coefficient onto the RT basis functions.
void LocalInterpolation_RT(const VectorFiniteElement &cfe, const real_t *nk, const Array< int > &d2n, ElementTransformation &Trans, DenseMatrix &I) const
void ProjectMatrixCoefficient_ND(const real_t *tk, const Array< int > &d2t, MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Project the rows of the matrix coefficient in an ND space.
void LocalRestriction_ND(const real_t *tk, const Array< int > &d2t, ElementTransformation &Trans, DenseMatrix &R) const
void LocalL2Projection_RT(const VectorFiniteElement &cfe, ElementTransformation &Trans, DenseMatrix &I) const
void ProjectGrad_ND(const real_t *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
void CalcVShape_ND(ElementTransformation &Trans, DenseMatrix &shape) const
void CalcVShape_RT(ElementTransformation &Trans, DenseMatrix &shape) const
VectorFiniteElement(int D, Geometry::Type G, int Do, int O, int M, int F=FunctionSpace::Pk)
void LocalInterpolation_ND(const VectorFiniteElement &cfe, const real_t *tk, const Array< int > &d2t, ElementTransformation &Trans, DenseMatrix &I) const
void ProjectCurl_RT(const real_t *nk, const Array< int > &d2n, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
void LocalL2Projection_ND(const VectorFiniteElement &cfe, ElementTransformation &Trans, DenseMatrix &I) const
void ProjectCurl_ND(const real_t *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
void ProjectMatrixCoefficient_RT(const real_t *nk, const Array< int > &d2n, MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Project the rows of the matrix coefficient in an RT space.
virtual ~VectorTensorFiniteElement()
VectorTensorFiniteElement(const int dims, const int d, const int p, const int cbtype, const int obtype, const int M, const DofMapType dmtype)
void SetDataAndSize(real_t *d, int s)
Set the Vector data and size.
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
void trans(const Vector &u, Vector &x)
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
real_t u(const Vector &xvec)
void mfem_error(const char *msg)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
Multiply a matrix A with the transpose of a matrix B: A*Bt.
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
void AddMult_a_VWt(const real_t a, const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += a * v w^t.
void AddMult_a_VVt(const real_t a, const Vector &v, DenseMatrix &VVt)
VVt += a * v v^t.
void AddMult_a_AAt(real_t a, const DenseMatrix &A, DenseMatrix &AAt)
AAt += a * A * A^t.
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
void InvertLinearTrans(ElementTransformation &trans, const IntegrationPoint &pt, Vector &x)
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
real_t p(const Vector &x, real_t t)
std::array< int, NCMesh::MaxFaceNodes > nodes