Commit d3a92bd9 authored by Michael Turner's avatar Michael Turner

fix

parent 5aa2c58b
......@@ -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");
}
ProcessLinear::~ProcessLinear()
......@@ -234,36 +234,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