From 8ed03598b4b54d9b6bb20629c1f0c2b3dd4f0dab Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Wed, 22 Jun 2016 09:06:51 +0000
Subject: [PATCH] (1) improved stats for analyseCurvedMesh plugin. (2) add
function that sample quality measures for testing
---
Mesh/qualityMeasuresJacobian.cpp | 123 ++++++++++++++++++++++++++
Mesh/qualityMeasuresJacobian.h | 5 ++
Plugin/AnalyseCurvedMesh.cpp | 143 ++++++++++++++++++++-----------
3 files changed, 223 insertions(+), 48 deletions(-)
diff --git a/Mesh/qualityMeasuresJacobian.cpp b/Mesh/qualityMeasuresJacobian.cpp
index d9c8f500c7..48ed8d53fb 100644
--- a/Mesh/qualityMeasuresJacobian.cpp
+++ b/Mesh/qualityMeasuresJacobian.cpp
@@ -391,6 +391,129 @@ double minIsotropyMeasure(MElement *el,
// return min;
//}
+double minSampledIsotropyMeasure(MElement *el, int deg, bool writeInFile)//fordebug
+{
+ fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
+ el->getNodesCoord(nodesXYZ);
+
+ const JacobianBasis *jacBasis;
+ const GradientBasis *gradBasis;
+
+ const int type = el->getType();
+// const int order = el->getPolynomialOrder();
+// const int jacOrder = order * el->getDim();
+ const bool serendipFalse = false;
+
+ FuncSpaceData jacMatSpace, jacDetSpace;
+
+ switch(type) {
+ case TYPE_TRI:
+ case TYPE_TET:
+ case TYPE_QUA:
+ case TYPE_HEX:
+ case TYPE_PRI:
+ jacMatSpace = FuncSpaceData(el, deg, &serendipFalse);
+ jacDetSpace = FuncSpaceData(el, deg, &serendipFalse);
+ break;
+ case TYPE_PYR:
+// jacMatSpace = FuncSpaceData(el, false, order, order-1, &serendipFalse);
+// jacDetSpace = FuncSpaceData(el, false, jacOrder, jacOrder-3, &serendipFalse);
+ break;
+ default:
+ Msg::Error("Isotropy not implemented for type of element %d", el->getType());
+ return -1;
+ }
+ gradBasis = BasisFactory::getGradientBasis(jacMatSpace);
+ jacBasis = BasisFactory::getJacobianBasis(jacDetSpace);
+
+ fullVector<double> coeffDetLag(jacBasis->getNumJacNodes());
+ jacBasis->getSignedIdealJacobian(nodesXYZ, coeffDetLag);
+
+ fullMatrix<double> coeffMatLag(gradBasis->getNumSamplingPoints(), 9);
+ gradBasis->getAllIdealGradientsFromNodes(nodesXYZ, coeffMatLag);
+
+ double min = std::numeric_limits<double>::infinity();
+ for (int i = 0; i < coeffDetLag.size(); ++i) {
+ double frobNorm = 0;
+ if (el->getDim() == 2) {
+ for (int k = 0; k < 6; ++k)
+ frobNorm += coeffMatLag(i, k) * coeffMatLag(i, k);
+ min = std::min(min, 2*coeffDetLag(i)/frobNorm);
+ }
+ else if (el->getDim() == 3) {
+ for (int k = 0; k < 9; ++k)
+ frobNorm += coeffMatLag(i, k) * coeffMatLag(i, k);
+ min = std::min(min, 3*std::pow(coeffDetLag(i), 2/3.)/frobNorm);
+ }
+ }
+
+ return min;
+}
+
+double minSampledScaledJacobian(MElement *el, int deg, bool writeInFile)//fordebug
+{
+ fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
+ el->getNodesCoord(nodesXYZ);
+
+ const JacobianBasis *jacBasis;
+ const GradientBasis *gradBasis;
+
+ const int type = el->getType();
+// const int order = el->getPolynomialOrder();
+// const int jacOrder = order * el->getDim();
+ const bool serendipFalse = false;
+
+ FuncSpaceData jacMatSpace, jacDetSpace;
+
+ switch(type) {
+ case TYPE_TRI:
+ case TYPE_TET:
+ case TYPE_QUA:
+ case TYPE_HEX:
+ case TYPE_PRI:
+ jacMatSpace = FuncSpaceData(el, deg, &serendipFalse);
+ jacDetSpace = FuncSpaceData(el, deg, &serendipFalse);
+ break;
+ case TYPE_PYR:
+// jacMatSpace = FuncSpaceData(el, false, order, order-1, &serendipFalse);
+// jacDetSpace = FuncSpaceData(el, false, jacOrder, jacOrder-3, &serendipFalse);
+ break;
+ default:
+ Msg::Error("Isotropy not implemented for type of element %d", el->getType());
+ return -1;
+ }
+ gradBasis = BasisFactory::getGradientBasis(jacMatSpace);
+ jacBasis = BasisFactory::getJacobianBasis(jacDetSpace);
+
+ fullVector<double> coeffDetLag(jacBasis->getNumJacNodes());
+ jacBasis->getSignedIdealJacobian(nodesXYZ, coeffDetLag);
+
+ fullMatrix<double> coeffMatLag(gradBasis->getNumSamplingPoints(), 9);
+ gradBasis->getAllIdealGradientsFromNodes(nodesXYZ, coeffMatLag);
+
+ double min = std::numeric_limits<double>::infinity();
+ for (int i = 0; i < coeffDetLag.size(); ++i) {
+ if (el->getDim() == 2) {
+ double v[2] = {0, 0};
+ for (int k = 0; k < 2; ++k) {
+ for (int l = k*3; l < k*3+3; ++l)
+ v[k] += coeffMatLag(i, l) * coeffMatLag(i, l);
+ }
+ min = std::min(min, coeffDetLag(i)/v[0]/v[1]);
+ }
+ else if (el->getDim() == 3) {
+ double v[3] = {0, 0, 0};
+ for (int k = 0; k < 3; ++k) {
+ for (int l = k*3; l < k*3+3; ++l)
+ v[k] += coeffMatLag(i, l) * coeffMatLag(i, l);
+ }
+ min = std::min(min, coeffDetLag(i)/std::sqrt(v[0]*v[1]*v[2]));
+ }
+ }
+
+ return min;
+}
+
// Virtual class _CoeffData
bool _lessMinB::operator()(_CoeffData *cd1, _CoeffData *cd2) const
{
diff --git a/Mesh/qualityMeasuresJacobian.h b/Mesh/qualityMeasuresJacobian.h
index 5c4921c60a..dbf4eae0bf 100644
--- a/Mesh/qualityMeasuresJacobian.h
+++ b/Mesh/qualityMeasuresJacobian.h
@@ -28,6 +28,10 @@ double minIsotropyMeasure(MElement *el,
bool reversedOk = false);
//double minSampledAnisotropyMeasure(MElement *el, int order,//fordebug
// bool writeInFile = false);
+double minSampledIsotropyMeasure(MElement *el, int order,//fordebug
+ bool writeInFile = false);
+double minSampledScaledJacobian(MElement *el, int order,//fordebug
+ bool writeInFile = false);
class _CoeffData
{
@@ -45,6 +49,7 @@ public:
inline double maxL() const {return _maxL;}
inline double minB() const {return _minB;}
inline double maxB() const {return _maxB;}
+ inline int depth() const {return _depth;}
virtual bool boundsOk(double minL, double maxL) const = 0;
virtual void getSubCoeff(std::vector<_CoeffData*>&) const = 0;
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index f4aa42dc40..49771e7dd4 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -27,7 +27,7 @@ StringXNumber CurvedMeshOptions_Number[] = {
{GMSH_FULLRC, "Jacobian determinant", NULL, 1},
{GMSH_FULLRC, "Scaled Jacobian", NULL, 0},
{GMSH_FULLRC, "Isotropy", NULL, 1},
- {GMSH_FULLRC, "Hidding threshold", NULL, 10},
+ {GMSH_FULLRC, "Hidding threshold", NULL, 9},
{GMSH_FULLRC, "Draw PView", NULL, 0},
{GMSH_FULLRC, "Recompute", NULL, 0},
{GMSH_FULLRC, "Dimension of elements", NULL, -1}
@@ -133,21 +133,24 @@ PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
double time = Cpu();
Msg::StatusBar(true, "Computing Jacobian for %dD elements...", dim);
_computeMinMaxJandValidity(dim);
- Msg::StatusBar(true, "Done in %g seconds", Cpu()-time);
+ Msg::StatusBar(true, "Done computing Jacobian for %dD elements (%g s)",
+ dim, Cpu()-time);
printStatJ = true;
}
if (computeScaledJac && !_computedS[dim-1]) {
double time = Cpu();
Msg::StatusBar(true, "Computing scaled Jacobian for %dD elements...", dim);
_computeMinScaledJac(dim);
- Msg::StatusBar(true, "Done in %g seconds", Cpu()-time);
+ Msg::StatusBar(true, "Done computing scaled Jacobian for %dD elements (%g s)",
+ dim, Cpu()-time);
printStatS = true;
}
if (computeIsotropy && !_computedI[dim-1]) {
double time = Cpu();
Msg::StatusBar(true, "Computing isotropy for %dD elements...", dim);
_computeMinIsotropy(dim);
- Msg::StatusBar(true, "Done in %g seconds", Cpu()-time);
+ Msg::StatusBar(true, "Done computing isotropy for %dD elements (%g s)",
+ dim, Cpu()-time);
printStatI = true;
}
}
@@ -222,47 +225,79 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim)
{
if (_computedJ[dim-1]) return;
+ std::set<GEntity*, GEntityLessThan> entities;
switch (dim) {
- case 3 :
- for (GModel::riter it = _m->firstRegion(); it != _m->lastRegion(); it++) {
- GRegion *r = *it;
- unsigned int numType[6] = {0, 0, 0, 0, 0, 0};
- r->getNumMeshElements(numType);
-
- for(int type = 0; type < 5; type++) {
- MElement *const *el = r->getStartElementType(type);
- _computeMinMaxJandValidity(el, numType[type]);
- }
- }
+ case 3:
+ for (GModel::riter it = _m->firstRegion(); it != _m->lastRegion(); it++)
+ entities.insert(*it);
break;
-
- case 2 :
- for (GModel::fiter it = _m->firstFace(); it != _m->lastFace(); it++) {
- GFace *f = *it;
- unsigned int numType[3] = {0, 0, 0};
- f->getNumMeshElements(numType);
-
- for (int type = 0; type < 3; type++) {
- MElement *const *el = f->getStartElementType(type);
- _computeMinMaxJandValidity(el, numType[type]);
- }
- }
+ case 2:
+ for (GModel::fiter it = _m->firstFace(); it != _m->lastFace(); it++)
+ entities.insert(*it);
break;
-
- case 1 :
- for (GModel::eiter it = _m->firstEdge(); it != _m->lastEdge(); it++) {
- GEdge *e = *it;
- unsigned int numElement = e->getNumMeshElements();
- MElement *const *el = e->getStartElementType(0);
- _computeMinMaxJandValidity(el, numElement);
- }
+ case 1:
+ for (GModel::eiter it = _m->firstEdge(); it != _m->lastEdge(); it++)
+ entities.insert(*it);
break;
-
- default :
+ default:
Msg::Fatal("This should not happen.");
return;
}
+ std::set<GEntity*, GEntityLessThan>::iterator it;
+ for (it = entities.begin(); it != entities.end(); ++it) {
+ GEntity *entity = *it;
+ unsigned num = entity->getNumMeshElements();
+ switch (dim) {
+ case 3:
+ Msg::StatusBar(true, "Volume %d: checking the Jacobian of %d elements",
+ entity->tag(), num);
+ break;
+ case 2:
+ Msg::StatusBar(true, "Surface %d: checking the Jacobian of %d elements",
+ entity->tag(), num);
+ break;
+ case 1:
+ Msg::StatusBar(true, "Line %d: checking the Jacobian of %d elements",
+ entity->tag(), num);
+ break;
+ default: break;
+ }
+
+ double initial, time = initial = Cpu();
+ unsigned int percentage = 0, nextCheck = 0;
+
+ for (unsigned i = 0; i < num; ++i) {
+ MElement *el = entity->getMeshElement(i);
+ double min, max;
+ jacobianBasedQuality::minMaxJacobianDeterminant(el, min, max);
+ _data.push_back(data_elementMinMax(el, min, max));
+ if (min < 0 && max < 0) {
+ Msg::Warning("Element %d is completely inverted", el->getNum());
+ }
+
+ if (i >= nextCheck) {
+ nextCheck += _data.size() / 100;
+ double curTime = Cpu();
+ unsigned int curPercentage = i*100/_data.size();
+ if ((curTime - time > 10. && curPercentage > percentage + 4) ||
+ (curTime - time > 15. && curPercentage < 5)) {
+ percentage = curPercentage;
+ time = curTime;
+ const double remaining = (time-initial) / (i+1) * (_data.size() - i-1);
+ if (remaining < 60*2)
+ Msg::StatusBar(true, "%d%% (remaining time ~%g seconds)",
+ percentage, remaining);
+ else if (remaining < 60*60*2)
+ Msg::StatusBar(true, "%d%% (remaining time ~%g minutes)",
+ percentage, remaining/60);
+ else
+ Msg::StatusBar(true, "%d%% (remaining time ~%g hours)",
+ percentage, remaining/3600);
+ }
+ }
+ }
+ }
_computedJ[dim-1] = true;
}
@@ -294,6 +329,12 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinScaledJac(int dim)
else {
_data[i].setMinS(jacobianBasedQuality::minScaledJacobian(el, true));
}
+// Msg::Info("Scaled Jac");
+// Msg::Info("==========");
+// for (int k = 1; k < 30; ++k) {
+// Msg::Info("%.10g", jacobianBasedQuality::minSampledScaledJacobian(el, k));
+// }
+// Msg::Info(" ");
if (i >= nextCheck) {
nextCheck += _data.size() / 100;
double curTime = Cpu();
@@ -335,6 +376,12 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinIsotropy(int dim)
else {
_data[i].setMinI(jacobianBasedQuality::minIsotropyMeasure(el, true));
}
+// Msg::Info("Isotropy");
+// Msg::Info("========");
+// for (int k = 1; k < 30; ++k) {
+// Msg::Info("%.10g", jacobianBasedQuality::minSampledIsotropyMeasure(el, k));
+// }
+// Msg::Info(" ");
if (i >= nextCheck) {
nextCheck += _data.size() / 100;
double curTime = Cpu();
@@ -362,7 +409,7 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinIsotropy(int dim)
int GMSH_AnalyseCurvedMeshPlugin::_hideWithThreshold(int askedDim, int whichMeasure)
{
- // only hide for wuality measures
+ // only hide for quality measures
if (_threshold > 1 || whichMeasure == 0) return 0;
int nHidden = 0;
@@ -429,12 +476,12 @@ void GMSH_AnalyseCurvedMeshPlugin::_printStatJacobian()
avgratJ /= count;
avgratJc /= countc;
- Msg::Info("Jacobian determinant: in [%.3g, %.3g] (avg=%.3g)",
- infminJ, supminJ, avgminJ);
- Msg::Info("Ratio minJ/maxJ : in [%.3g, %.3g] (avg=%.3g)",
- infratJ, supratJ, avgratJ);
- if (countc < count)
- Msg::Info(" (avg=%.3g"
+ Msg::Info("Min Jac. det. (minJ): %.3g, %.3g, %.3g (= min, avg, max)",
+ infminJ, avgminJ, supminJ);
+ Msg::Info("Ratio minJ/maxJ : %.3f, %.3f, %.3f (= worst, avg, best)",
+ infratJ, avgratJ, supratJ);
+ if (countc && countc < count)
+ Msg::Info(" (avg=%3.3g"
" on the %d non-constant elements)",
avgratJc, countc);
}
@@ -455,8 +502,8 @@ void GMSH_AnalyseCurvedMeshPlugin::_printStatScaledJac()
}
avgminS /= _data.size();
- Msg::Info("Scaled Jacobian : in [%.3g, %.3g] (avg=%.3g)",
- infminS, supminS, avgminS);
+ Msg::Info("Scaled Jacobian : %.3f, %.3f, %.3f (= worst, avg, best)",
+ infminS, avgminS, supminS);
}
void GMSH_AnalyseCurvedMeshPlugin::_printStatIsotropy()
@@ -475,8 +522,8 @@ void GMSH_AnalyseCurvedMeshPlugin::_printStatIsotropy()
}
avgminI /= _data.size();
- Msg::Info("Isotropy : in [%.3g, %.3g] (avg=%.3g)",
- infminI, supminI, avgminI);
+ Msg::Info("Isotropy : %.3f, %.3f, %.3f (= worst, avg, best)",
+ infminI, avgminI, supminI);
}
// For testing
--
GitLab