phynexis-deps/mfem
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Using max of space dim and new vector dim to fix support for ND and RT on surfaces

Browse Source
This commit is contained in:
Stowell, Mark L
2022-03-04 13:55:38 -08:00
parent 8c775bcb25
commit fa8be935c9
5 changed files with 21 additions and 13 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fem/fe/fe_base.cpp
+7 -5
View File
@@ -752,7 +752,8 @@ void NodalFiniteElement::Project(
}
else
{
DenseMatrix vshape(fe.GetDof(), fe.GetVDim());
DenseMatrix vshape(fe.GetDof(), std::max(Trans.GetSpaceDim(),
fe.GetVDim()));
I.SetSize(vshape.Width()*dof, fe.GetDof());
for (int k = 0; k < dof; k++)
@@ -1036,7 +1037,7 @@ void VectorFiniteElement::Project_RT(
{
int sdim = Trans.GetSpaceDim();
double vk[Geometry::MaxDim];
DenseMatrix vshape(fe.GetDof(), fe.GetVDim());
DenseMatrix vshape(fe.GetDof(), sdim);
Vector vshapenk(fe.GetDof());
const bool square_J = (dim == sdim);
@@ -1221,7 +1222,7 @@ void VectorFiniteElement::Project_ND(
Trans.SetIntPoint(&ip);
// Transform ND edge tengents from reference to physical space
// vk = J tk
Trans.Jacobian().Mult(tk + d2t[k]*vdim, vk);
Trans.Jacobian().Mult(tk + d2t[k]*dim, vk);
if (fe.GetMapType() == INTEGRAL)
{
double w = 1.0/Trans.Weight();
@@ -1249,8 +1250,9 @@ void VectorFiniteElement::Project_ND(
}
else
{
int sdim = Trans.GetSpaceDim();
double vk[Geometry::MaxDim];
DenseMatrix vshape(fe.GetDof(), fe.GetVDim());
DenseMatrix vshape(fe.GetDof(), sdim);
Vector vshapetk(fe.GetDof());
I.SetSize(dof, fe.GetDof());
@@ -1287,7 +1289,7 @@ void VectorFiniteElement::ProjectGrad_ND(
for (int k = 0; k < dof; k++)
{
fe.CalcDShape(Nodes.IntPoint(k), dshape);
dshape.Mult(tk + d2t[k]*vdim, grad_k);
dshape.Mult(tk + d2t[k]*dim, grad_k);
for (int j = 0; j < grad_k.Size(); j++)
{
grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);