Commit e3022395 authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'feature/qualitythreshold_linearise' into 'master'

feature/qualitythreshold_linearise

add quality threshold to linearise module so that any elements below a certain value will linearised. This works based on the scaled jacobian.

See merge request !690
parents 062312ea 9cddc721
......@@ -52,6 +52,7 @@ v4.4.0
- Add support for Gmsh high-order output (!679)
- Move CAD classes to factory format (!676)
- Add option to `linearise` module to linearise only prisms (!688)
- Add option to `linearise` to use element quality (!690)
**FieldConvert:**
- Move all modules to a new library, FieldUtils, to support post-processing
......
......@@ -52,7 +52,7 @@ ProcessLinear::ProcessLinear(MeshSharedPtr m) : ProcessModule(m)
m_config["all"] =
ConfigOption(true, "0", "remove curve nodes for all elements.");
m_config["invalid"] =
ConfigOption(true, "0", "remove curve nodes if element is invalid.");
ConfigOption(false, "0", "remove curve nodes if element is invalid.");
m_config["prismonly"] =
ConfigOption(false, "", "only acts on prims");
}
......@@ -69,7 +69,8 @@ void ProcessLinear::Process()
}
bool all = m_config["all"].as<bool>();
bool invalid = m_config["invalid"].as<bool>();
bool invalid = m_config["invalid"].beenSet;
NekDouble thr = m_config["invalid"].as<NekDouble>();
ASSERTL0(all || invalid,
"Must specify an option: all (to remove all curvature) or invalid "
......@@ -107,7 +108,8 @@ void ProcessLinear::Process()
map<int,vector<FaceSharedPtr> > eidToFace;
map<int,vector<ElementSharedPtr> > eidToElm;
vector<ElementSharedPtr> el = m_mesh->m_element[m_mesh->m_expDim];
vector<ElementSharedPtr> els = m_mesh->m_element[m_mesh->m_expDim];
vector<ElementSharedPtr> el = els;
for(int i = 0; i < el.size(); i++)
{
......@@ -149,15 +151,7 @@ void ProcessLinear::Process()
}
}
// Create elemental geometry.
SpatialDomains::GeometrySharedPtr geom =
el[i]->GetGeom(m_mesh->m_spaceDim);
// Generate geometric factors.
SpatialDomains::GeomFactorsSharedPtr gfac =
geom->GetGeomFactors();
if (!gfac->IsValid())
if (Invalid(el[i],thr))//(!gfac->IsValid())
{
clearedElmts.insert(el[i]->GetId());;
el[i]->SetVolumeNodes(zeroNodes);
......@@ -201,15 +195,16 @@ void ProcessLinear::Process()
}
}
vector<ElementSharedPtr> tmp = el;
el.clear();
set<int>::iterator it;
for(int i = 0; i < tmp.size(); i++)
set<int>::iterator it1;
boost::unordered_set<int>::iterator it2;
for(int i = 0; i < els.size(); i++)
{
it = neigh.find(tmp[i]->GetId());
if(it != neigh.end())
it1 = neigh.find(els[i]->GetId());
it2 = clearedElmts.find(els[i]->GetId());
if(it1 != neigh.end() && it2 == clearedElmts.end())
{
el.push_back(tmp[i]);
el.push_back(els[i]);
}
}
neigh.clear();
......@@ -223,5 +218,59 @@ void ProcessLinear::Process()
}
}
}
bool ProcessLinear::Invalid(ElementSharedPtr el, NekDouble thr)
{
// Create elemental geometry.
SpatialDomains::GeometrySharedPtr geom =
el->GetGeom(m_mesh->m_spaceDim);
// Generate geometric factors.
SpatialDomains::GeomFactorsSharedPtr gfac =
geom->GetGeomFactors();
if(!gfac->IsValid())
{
return true;
}
LibUtilities::PointsKeyVector p = geom->GetPointsKeys();
SpatialDomains::DerivStorage deriv = gfac->GetDeriv(p);
const int pts = deriv[0][0].num_elements();
Array<OneD,NekDouble> jc(pts);
for (int k = 0; k < pts; ++k)
{
DNekMat jac(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];
}
}
if(m_mesh->m_expDim == 2)
{
jc[k] = jac(0,0) * jac(1,1) - jac(0,1)*jac(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));
}
}
NekDouble scaledJac = Vmath::Vmin(jc.num_elements(),jc,1) /
Vmath::Vmax(jc.num_elements(),jc,1);
if(scaledJac < thr)
{
return true;
}
return false;
}
}
}
......@@ -61,6 +61,8 @@ public:
/// Write mesh to output file.
virtual void Process();
private:
bool Invalid(ElementSharedPtr el, NekDouble 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