Commit 2c38d5c2 authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'feature/qmetric_scaled' into 'master'

scaled jacobian option to qualitymetric module

as titled.
Picked from meshgeneration/feature/varopti_DNC

See merge request !695
parents f238a4b9 49a95466
......@@ -41,6 +41,7 @@ v4.4.0
**FieldConvert:**
- Allow equi-spaced output for 1D and 2DH1D fields (!613)
- Update quality metric to include scaled Jacobian output (!695)
**NekMesh:**
- Modify curve module to allow for spline input (!628)
......
......@@ -427,3 +427,17 @@ year={2011}
volume = {33},
year = {2000}
}
@article{GaRoPeSa15,
title = {Distortion and quality measures for validating and generating
high-order tetrahedral meshes},
author = {Gargallo-Peir{\'o}, Abel and Roca, Xevi and Peraire, Jaime and
Sarrate, Josep},
journal = {Engineering with Computers},
volume = 31,
number = 3,
pages = {423--437},
year = 2015,
publisher = {Springer London}
}
......@@ -122,6 +122,7 @@ possibly also Reynolds stresses) into single file;
\item \inltt{interppoints}: Interpolates a set of points to another, requires fromfld and fromxml to be defined, a line or plane of points can be defined;
\item \inltt{isocontour}: Extract an isocontour of ``fieldid'' variable and at value ``fieldvalue''. Optionally ``fieldstr'' can be specified for a string defiition or ``smooth'' for smoothing;
\item \inltt{jacobianenergy}: Shows high frequency energy of Jacobian;
\item \inltt{qualitymetric}: Evaluate a quality metric of the underlying mesh to show mesh quality;
\item \inltt{meanmode}: Extract mean mode (plane zero) of 3DH1D expansions;
\item \inltt{pointdatatofld}: Given discrete data at quadrature points
project them onto an expansion basis and output fld file;
......@@ -629,6 +630,38 @@ The output file \inltt{jacenergy.fld} can be processed in a similar
way as described in section \ref{s:utilities:fieldconvert:sub:convert}
to visualise the result either in Tecplot, Paraview or VisIt.
\subsection{Calculate mesh quality: \textit{qualitymetric} module}
The \inltt{qualitymetric} module assesses the quality of the mesh by calculating
a per-element quality metric and adding an additional field to any resulting
output. This does not require any field input, therefore an example usage looks
like
\begin{lstlisting}[style=BashInputStyle]
FieldConvert -m qualitymetric mesh.xml mesh-with-quality.dat
\end{lstlisting}
Two quality metrics are implemented that produce scalar fields $Q$:
\begin{itemize}
\item By default a metric outlined in~\cite{GaRoPeSa15} is produced, where all
straight sided elements have quality $Q = 1$ and $Q < 1$ shows the deformation
between the curved element and the straight-sided element. If $Q = 0$ then the
element is invalid. Note that $Q$ varies over the volume of the element but is
not guaranteed to be continuous between elements.
\item Alternatively, if the \inlsh{scaled} option is passed through to the
module, then the scaled Jacobian
\[
J_s =
\frac{\min_{\xi\in\Omega_{\text{st}}}J(\xi)}{\max_{\xi\in\Omega_{\text{st}}}J(\xi)}
\]
(i.e. the ratio of the minimum to maximum Jacobian of each element) is
calculated. Again $Q = 1$ denotes an ideal element, but now invalid elements
are shown by $Q < 0$. Any elements with $Q$ near zero are determined to be low
quality.
\end{itemize}
%
%
%
......
......@@ -60,6 +60,8 @@ ModuleKey ProcessQualityMetric::className =
ProcessQualityMetric::ProcessQualityMetric(FieldSharedPtr f) : ProcessModule(f)
{
m_config["scaled"] =
ConfigOption(true, "", "use scaled jacobian instead");
}
ProcessQualityMetric::~ProcessQualityMetric()
......@@ -85,7 +87,7 @@ void ProcessQualityMetric::Process(po::variables_map &vm)
// copy Jacobian into field
LocalRegions::ExpansionSharedPtr Elmt = m_f->m_exp[0]->GetExp(i);
int offset = m_f->m_exp[0]->GetPhys_Offset(i);
Array<OneD, NekDouble> q = GetQ(Elmt);
Array<OneD, NekDouble> q = GetQ(Elmt,m_config["scaled"].m_beenSet);
Array<OneD, NekDouble> out = phys + offset;
ASSERTL0(q.num_elements() == Elmt->GetTotPoints(),
......@@ -301,8 +303,7 @@ inline vector<DNekMat> MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom,
return ret;
}
Array<OneD, NekDouble> ProcessQualityMetric::GetQ(
LocalRegions::ExpansionSharedPtr e)
Array<OneD, NekDouble> ProcessQualityMetric::GetQ(LocalRegions::ExpansionSharedPtr e, bool s)
{
SpatialDomains::GeometrySharedPtr geom = e->GetGeom();
StdRegions::StdExpansionSharedPtr chi = e->GetGeom()->GetXmap();
......@@ -407,18 +408,40 @@ Array<OneD, NekDouble> ProcessQualityMetric::GetQ(
ASSERTL0(false, "silly exp dim");
}
NekDouble frob = 0.0;
for (int i = 0; i < expDim; ++i)
if(s)
{
for (int j = 0; j < expDim; ++j)
eta[k] = jacDet;
}
else
{
NekDouble frob = 0.0;
for (int i = 0; i < expDim; ++i)
{
frob += jacIdeal(i, j) * jacIdeal(i, j);
for (int j = 0; j < expDim; ++j)
{
frob += jacIdeal(i,j) * jacIdeal(i,j);
}
}
NekDouble sigma = 0.5*(jacDet + sqrt(jacDet*jacDet));
eta[k] = expDim * pow(sigma, 2.0/expDim) / frob;
}
}
NekDouble sigma = 0.5 * (jacDet + sqrt(jacDet * jacDet));
eta[k] = expDim * pow(sigma, 2.0 / expDim) / frob;
if(s)
{
NekDouble mx = -1.0 * numeric_limits<double>::max();
NekDouble mn = numeric_limits<double>::max();
for(int k = 0; k < pts; k++)
{
mx = max(mx,eta[k]);
mn = min(mn,eta[k]);
}
for(int k = 0; k < pts; k++)
{
eta[k] = mn/mx;
}
}
// Project onto output stuff
......
......@@ -66,7 +66,7 @@ public:
}
private:
Array<OneD, NekDouble> GetQ(LocalRegions::ExpansionSharedPtr e);
Array<OneD, NekDouble> GetQ(LocalRegions::ExpansionSharedPtr e, bool s);
};
}
}
......
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