Commit fd5eb756 authored by Dave Moxey's avatar Dave Moxey

Merge branch 'fix/linearise_mapping' into 'master'

fix/linearise

This MR fixes the linearise threshold routine which was being far too over aggressive in the linearisation of prism elements.
This is because it was measuring linear deformation when considering the quality, this has now been fixed.

See merge request !717
parents 200fd1d0 f7866e7c
......@@ -54,7 +54,7 @@ ProcessLinear::ProcessLinear(MeshSharedPtr m) : ProcessModule(m)
m_config["invalid"] =
ConfigOption(false, "0", "remove curve nodes if element is invalid.");
m_config["prismonly"] =
ConfigOption(false, "", "only acts on prims");
ConfigOption(true, "0", "only acts on prims");
m_config["extract"] =
ConfigOption(false, "", "dump a mesh of the extracted elements");
}
......@@ -259,36 +259,61 @@ bool ProcessLinear::Invalid(ElementSharedPtr el, NekDouble thr)
return true;
}
vector<NodeSharedPtr> ns = el->GetVertexList();
ElmtConfig c = el->GetConf();
c.m_order = 1;
c.m_faceNodes = false;
c.m_volumeNodes = false;
c.m_reorient = false;
ElementSharedPtr elL = GetElementFactory().CreateInstance(
c.m_e, c, ns, el->GetTagList());
SpatialDomains::GeometrySharedPtr geomL = elL->GetGeom(m_mesh->m_spaceDim);
SpatialDomains::GeomFactorsSharedPtr gfacL = geomL->GetGeomFactors();
LibUtilities::PointsKeyVector p = geom->GetPointsKeys();
SpatialDomains::DerivStorage deriv = gfac->GetDeriv(p);
SpatialDomains::DerivStorage derivL = gfacL->GetDeriv(p);
const int pts = deriv[0][0].num_elements();
Array<OneD,NekDouble> jc(pts);
Array<OneD,NekDouble> jcL(pts);
for (int k = 0; k < pts; ++k)
{
DNekMat jac(m_mesh->m_expDim, m_mesh->m_expDim, 0.0, eFULL);
DNekMat jacL(m_mesh->m_expDim, m_mesh->m_expDim, 0.0, eFULL);
for (int l = 0; l < m_mesh->m_expDim; ++l)
{
for (int j = 0; j < m_mesh->m_expDim; ++j)
{
jac(j,l) = deriv[l][j][k];
jacL(j,l) = derivL[l][j][k];
}
}
if(m_mesh->m_expDim == 2)
{
jc[k] = jac(0,0) * jac(1,1) - jac(0,1)*jac(1,0);
jcL[k] = jacL(0,0) * jacL(1,1) - jacL(0,1)*jacL(1,0);
}
else if(m_mesh->m_expDim == 3)
{
jc[k] = jac(0,0) * (jac(1,1)*jac(2,2) - jac(2,1)*jac(1,2)) -
jac(0,1) * (jac(1,0)*jac(2,2) - jac(2,0)*jac(1,2)) +
jac(0,2) * (jac(1,0)*jac(2,1) - jac(2,0)*jac(1,1));
jcL[k] = jacL(0,0) * (jacL(1,1)*jacL(2,2) - jacL(2,1)*jacL(1,2)) -
jacL(0,1) * (jacL(1,0)*jacL(2,2) - jacL(2,0)*jacL(1,2)) +
jacL(0,2) * (jacL(1,0)*jacL(2,1) - jacL(2,0)*jacL(1,1));
}
}
NekDouble scaledJac = Vmath::Vmin(jc.num_elements(),jc,1) /
Vmath::Vmax(jc.num_elements(),jc,1);
Array<OneD, NekDouble> j(pts);
Vmath::Vdiv(jc.num_elements(),jc,1,jcL,1,j,1);
NekDouble scaledJac = Vmath::Vmin(j.num_elements(),j,1) /
Vmath::Vmax(j.num_elements(),j,1);
//cout << scaledJac << endl;
if(scaledJac < thr)
{
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment