Skip to content
Snippets Groups Projects
Commit 4f297927 authored by Amaury Johnen's avatar Amaury Johnen
Browse files

update of the plugin to take into account the 2 quality measures (in addition to jacobian)

parent 5f51edba
No related branches found
No related tags found
No related merge requests found
...@@ -24,13 +24,13 @@ ...@@ -24,13 +24,13 @@
class bezierBasis; class bezierBasis;
StringXNumber CurvedMeshOptions_Number[] = { StringXNumber CurvedMeshOptions_Number[] = {
{GMSH_FULLRC, "Show: 0:J, 1:R, 2:J&&R", NULL, 1}, {GMSH_FULLRC, "Jacobian", NULL, 1},
{GMSH_FULLRC, "Draw PView", NULL, 1}, {GMSH_FULLRC, "Scaled Jacobian", NULL, 0},
{GMSH_FULLRC, "Hidding threshold", NULL, 10}, {GMSH_FULLRC, "Isotropy", NULL, 1},
{GMSH_FULLRC, "Dimension of elements", NULL, -1}, {GMSH_FULLRC, "Hidding threshold", NULL, 1},
{GMSH_FULLRC, "Draw PView", NULL, 0},
{GMSH_FULLRC, "Recompute bounds", NULL, 0}, {GMSH_FULLRC, "Recompute bounds", NULL, 0},
{GMSH_FULLRC, "Tolerance", NULL, 1e-3} {GMSH_FULLRC, "Dimension of elements", NULL, -1}
// tolerance: To be removed when MetricBasis => qualityMeasuresJacobian
}; };
extern "C" extern "C"
...@@ -86,60 +86,77 @@ std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const ...@@ -86,60 +86,77 @@ std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const
"the validity of the mesh."; "the validity of the mesh.";
} }
PView *GMSH_AnalyseCurvedMeshPlugin::execute(PView *v) PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
{ {
_m = GModel::current(); _m = GModel::current();
_computeMetric = static_cast<int>(CurvedMeshOptions_Number[0].def); bool computeJac = static_cast<int>(CurvedMeshOptions_Number[0].def);
bool drawPView = static_cast<int>(CurvedMeshOptions_Number[1].def); bool computeScaledJac = static_cast<int>(CurvedMeshOptions_Number[1].def);
_threshold = static_cast<double>(CurvedMeshOptions_Number[2].def); bool computeIsotropy = static_cast<int>(CurvedMeshOptions_Number[2].def);
int askedDim = static_cast<int>(CurvedMeshOptions_Number[3].def); _threshold = static_cast<double>(CurvedMeshOptions_Number[3].def);
bool recompute = static_cast<bool>(CurvedMeshOptions_Number[4].def); bool drawPView = static_cast<int>(CurvedMeshOptions_Number[4].def);
_tol = static_cast<double>(CurvedMeshOptions_Number[5].def); bool recompute = static_cast<bool>(CurvedMeshOptions_Number[5].def);
int askedDim = static_cast<int>(CurvedMeshOptions_Number[6].def);
if (askedDim < 0 || askedDim > 4) askedDim = _m->getDim(); if (askedDim < 0 || askedDim > 4) askedDim = _m->getDim();
if (recompute) { if (recompute) {
_data.clear(); _data.clear();
if (askedDim < 4) { if (askedDim < 4) {
_computedR[askedDim-1] = false;
_computedJ[askedDim-1] = false; _computedJ[askedDim-1] = false;
_computedS[askedDim-1] = false;
_computedI[askedDim-1] = false;
_PViewJ[askedDim-1] = false; _PViewJ[askedDim-1] = false;
_PViewR[askedDim-1] = false; _PViewS[askedDim-1] = false;
_PViewI[askedDim-1] = false;
} }
else { else {
_computedR[1] = _computedR[2] = false;
_computedJ[1] = _computedJ[2] = false; _computedJ[1] = _computedJ[2] = false;
_computedS[1] = _computedS[2] = false;
_computedI[1] = _computedI[2] = false;
_PViewJ[1] = _PViewJ[2] = false; _PViewJ[1] = _PViewJ[2] = false;
_PViewR[1] = _PViewR[2] = false; _PViewS[1] = _PViewS[2] = false;
_PViewI[1] = _PViewI[2] = false;
} }
_msgHide = true;
} }
// Compute J and/or R // Compute what have to
bool printStatJ = false;
bool printStatS = false;
bool printStatI = false;
for (int dim = 1; dim <= 3 && dim <= _m->getDim(); ++dim) { for (int dim = 1; dim <= 3 && dim <= _m->getDim(); ++dim) {
if ((askedDim == 4 && dim > 1) || dim == askedDim) { if ((askedDim == 4 && dim > 1) || dim == askedDim) {
if (!_computedJ[dim-1]) { if (!_computedJ[dim-1]) {
double time = Cpu(); double time = Cpu();
Msg::StatusBar(true, "Computing Jacobian for %dD elements...", dim); Msg::StatusBar(true, "Computing Jacobian for %dD elements...", dim);
_computeMinMaxJandValidity(dim); _computeMinMaxJandValidity(dim);
Msg::StatusBar(true, "... Done computing Jacobian (%g seconds)", Cpu()-time); Msg::StatusBar(true, "Done in %g seconds", Cpu()-time);
_printStatJacobian(); printStatJ = true;
} }
if (_computeMetric && !_computedR[dim-1]) { if (computeScaledJac && !_computedS[dim-1]) {
double time = Cpu(); double time = Cpu();
Msg::StatusBar(true, "Computing metric for %dD elements...", dim); Msg::StatusBar(true, "Computing scaled Jacobian for %dD elements...", dim);
_computeMinR(dim); _computeMinScaledJac(dim);
Msg::StatusBar(true, "... Done computing metric (%g seconds)", Cpu()-time); Msg::StatusBar(true, "Done in %g seconds", Cpu()-time);
_printStatMetric(); 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);
printStatI = true;
} }
} }
} }
if (printStatJ) _printStatJacobian();
if (printStatS) _printStatScaledJac();
if (printStatI) _printStatIsotropy();
// Create PView // Create PView
if (drawPView) if (drawPView)
for (int dim = 1; dim <= 3; ++dim) { for (int dim = 1; dim <= 3; ++dim) {
if ((askedDim == 4 && dim > 1) || dim == askedDim) { if ((askedDim == 4 && dim > 1) || dim == askedDim) {
if (!_PViewJ[dim-1] && _computeMetric != 1) { if (!_PViewJ[dim-1] && computeJac) {
_PViewJ[dim-1] = true; _PViewJ[dim-1] = true;
std::map<int, std::vector<double> > dataPV; std::map<int, std::vector<double> > dataPV;
for (unsigned int i = 0; i < _data.size(); ++i) { for (unsigned int i = 0; i < _data.size(); ++i) {
...@@ -156,48 +173,48 @@ PView *GMSH_AnalyseCurvedMeshPlugin::execute(PView *v) ...@@ -156,48 +173,48 @@ PView *GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
new PView(name.str().c_str(), "ElementData", _m, dataPV); new PView(name.str().c_str(), "ElementData", _m, dataPV);
} }
} }
if (!_PViewR[dim-1] && _computeMetric) { if (!_PViewS[dim-1] && computeScaledJac) {
_PViewR[dim-1] = true; _PViewS[dim-1] = true;
std::map<int, std::vector<double> > dataPV; std::map<int, std::vector<double> > dataPV;
for (unsigned int i = 0; i < _data.size(); ++i) { for (unsigned int i = 0; i < _data.size(); ++i) {
MElement *const el = _data[i].element(); MElement *const el = _data[i].element();
if (el->getDim() == dim) if (el->getDim() == dim)
dataPV[el->getNum()].push_back(_data[i].minR()); dataPV[el->getNum()].push_back(_data[i].minS());
} }
if (dataPV.size()) { if (dataPV.size()) {
std::stringstream name; std::stringstream name;
name << "min R " << dim << "D"; name << "Scaled Jacobian " << dim << "D";
new PView(name.str().c_str(), "ElementData", _m, dataPV); new PView(name.str().c_str(), "ElementData", _m, dataPV);
} }
}
/*std::map<int, std::vector<double> > dataPV2; if (!_PViewI[dim-1] && computeIsotropy) {
_PViewI[dim-1] = true;
std::map<int, std::vector<double> > dataPV;
for (unsigned int i = 0; i < _data.size(); ++i) { for (unsigned int i = 0; i < _data.size(); ++i) {
MElement *const el = _data[i].element(); MElement *const el = _data[i].element();
if (el->getDim() == dim) if (el->getDim() == dim)
dataPV2[el->getNum()].push_back(_data[i].minR()/_data[i].minRR()); dataPV[el->getNum()].push_back(_data[i].minI());
} }
if (dataPV.size()) { if (dataPV.size()) {
std::stringstream name; std::stringstream name;
name << "min RR " << dim << "D"; name << "Isotropy " << dim << "D";
new PView(name.str().c_str(), "ElementData", _m, dataPV2); new PView(name.str().c_str(), "ElementData", _m, dataPV);
}*/ }
} }
} }
} }
// Hide elements // Hide elements
#if defined(HAVE_OPENGL) #if defined(HAVE_OPENGL)
if (_hideWithThreshold(askedDim) && _msgHide) { _hideWithThreshold(askedDim,
_msgHide = false; computeIsotropy ? 2 :
Msg::Info("Note: Some elements have been hidden (param 'Hidding Threshold')."); computeScaledJac ? 1 :
Msg::Info(" (To revert: set 'Hidding Threshold' to 1 and rerun"); computeJac ? 0 : -1);
Msg::Info(" or, Tools > Visibility > Interactive > Show All)");
}
CTX::instance()->mesh.changed = (ENT_ALL); CTX::instance()->mesh.changed = (ENT_ALL);
drawContext::global()->draw(); drawContext::global()->draw();
#endif #endif
return 0; return NULL;
} }
void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim) void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim)
...@@ -208,7 +225,7 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim) ...@@ -208,7 +225,7 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim)
case 3 : case 3 :
for (GModel::riter it = _m->firstRegion(); it != _m->lastRegion(); it++) { for (GModel::riter it = _m->firstRegion(); it != _m->lastRegion(); it++) {
GRegion *r = *it; GRegion *r = *it;
unsigned int numType[5] = {0, 0, 0, 0, 0}; unsigned int numType[6] = {0, 0, 0, 0, 0, 0};
r->getNumMeshElements(numType); r->getNumMeshElements(numType);
for(int type = 0; type < 5; type++) { for(int type = 0; type < 5; type++) {
...@@ -254,34 +271,27 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(MElement *const*el ...@@ -254,34 +271,27 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(MElement *const*el
double min, max; double min, max;
jacobianBasedQuality::minMaxJacobianDeterminant(el[k], min, max); jacobianBasedQuality::minMaxJacobianDeterminant(el[k], min, max);
_data.push_back(data_elementMinMax(el[k], min, max)); _data.push_back(data_elementMinMax(el[k], min, max));
if (min < 0 && max < 0) {
Msg::Warning("Element %d is completely inverted", el[k]->getNum());
}
} }
} }
void GMSH_AnalyseCurvedMeshPlugin::_computeMinR(int dim) void GMSH_AnalyseCurvedMeshPlugin::_computeMinScaledJac(int dim)
{ {
if (_computedR[dim-1]) return; if (_computedS[dim-1]) return;
MetricBasis::setTol(_tol);
double initial, time = initial = Cpu(); double initial, time = initial = Cpu();
unsigned int percentage = 0, nextCheck = 0; unsigned int percentage = 0, nextCheck = 0;
for (unsigned int i = 0; i < _data.size(); ++i) { for (unsigned int i = 0; i < _data.size(); ++i) {
MElement *const el = _data[i].element(); MElement *const el = _data[i].element();
//if (el->getNum() != 10535) continue; if (el->getDim() != dim) continue;
/*if (el->getNum() != 2712 && if (_data[i].minJ() <= 0) {
el->getNum() != 3378 && _data[i].setMinS(0);
el->getNum() != 3997 && }
el->getNum() != 4314) continue;*/ else {
if (el->getDim() == dim) { _data[i].setMinS(jacobianBasedQuality::minScaledJacobian(el, true));
if (_data[i].minJ() <= 0)
_data[i].setMinR(0);
else {
_data[i].setMinR(MetricBasis::getMinR(el));
//_data[i].setMinR(jacobianBasedQuality::minScaledJacobian(el));
//_data[i].setMinR(MetricBasis::getMinSampledR(el, 10));
}
//_data[i].setMinRR(MetricBasis::getMaxSampledR(el, 10));
} }
if (i >= nextCheck) { if (i >= nextCheck) {
nextCheck += _data.size() / 100; nextCheck += _data.size() / 100;
...@@ -292,7 +302,7 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinR(int dim) ...@@ -292,7 +302,7 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinR(int dim)
percentage = curPercentage; percentage = curPercentage;
time = curTime; time = curTime;
const double remaining = (time-initial) / (i+1) * (_data.size() - i-1); const double remaining = (time-initial) / (i+1) * (_data.size() - i-1);
if (remaining < 300) if (remaining < 60*2)
Msg::StatusBar(true, "%d%% (remaining time ~%g seconds)", Msg::StatusBar(true, "%d%% (remaining time ~%g seconds)",
percentage, remaining); percentage, remaining);
else if (remaining < 60*60*2) else if (remaining < 60*60*2)
...@@ -305,68 +315,80 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinR(int dim) ...@@ -305,68 +315,80 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinR(int dim)
} }
} }
/*std::fstream file; _computedS[dim-1] = true;
file.open("comparisonMyMetric_and_Jacobian_vs_minr_maxr.txt", std::fstream::out); }
file << _data.size() << std::endl;
void GMSH_AnalyseCurvedMeshPlugin::_computeMinIsotropy(int dim)
{
if (_computedI[dim-1]) return;
double initial, time = initial = Cpu();
unsigned int percentage = 0, nextCheck = 0;
int n = 0, next = _data.size()/100;
double tm = Cpu();
for (unsigned int i = 0; i < _data.size(); ++i) { for (unsigned int i = 0; i < _data.size(); ++i) {
if (++n > next) {
Msg::Info("%d%% %g seconds", n*100/_data.size(), Cpu()-tm);
next += _data.size()/100;
}
MElement *const el = _data[i].element(); MElement *const el = _data[i].element();
file << _data[i].minJ() << " "; if (el->getDim() != dim) continue;
file << _data[i].maxJ() << " "; if (_data[i].minJ() <= 0) {
file << _data[i].minR() << " "; _data[i].setMinI(0);
file << MetricBasis::minSampledGlobalRatio(el, 7) << std::endl; }
else {
_data[i].setMinI(jacobianBasedQuality::minIsotropyMeasure(el, true));
}
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);
}
}
} }
file.close();*/
_computedR[dim-1] = true; _computedI[dim-1] = true;
} }
bool GMSH_AnalyseCurvedMeshPlugin::_hideWithThreshold(int askedDim) int GMSH_AnalyseCurvedMeshPlugin::_hideWithThreshold(int askedDim, int whichMeasure)
{ {
if (_threshold >= 1.) return false; int nHidden = 0;
bool ans = false;
for (unsigned int i = 0; i < _data.size(); ++i) { for (unsigned int i = 0; i < _data.size(); ++i) {
MElement *const el = _data[i].element(); MElement *const el = _data[i].element();
const int dim = el->getDim(); const int dim = el->getDim();
if ((askedDim == 4 && dim > 1) || dim == askedDim) { if ((askedDim == 4 && dim > 1) || dim == askedDim) {
if (( _computeMetric && _data[i].minR() > _threshold) || bool toHide = false;
(!_computeMetric && _data[i].minJ() > _threshold) ) { switch (whichMeasure) {
case 2:
toHide = _data[i].minI() > _threshold;
break;
case 1:
toHide = _data[i].minS() > _threshold;
break;
case 0:
toHide = _data[i].minJ() > _threshold;
break;
}
if (toHide) {
el->setVisibility(0); el->setVisibility(0);
ans = true; ++nHidden;
} }
else { else {
el->setVisibility(1); el->setVisibility(1);
} }
} }
} }
return ans; return nHidden;
}
void GMSH_AnalyseCurvedMeshPlugin::_printStatMetric()
{
if (_data.empty()) {
Msg::Info("No stat to print.");
return;
}
double infminR, supminR, avgminR;
infminR = supminR = avgminR = _data[0].minR();
for (unsigned int i = 1; i < _data.size(); ++i) {
infminR = std::min(infminR, _data[i].minR());
supminR = std::max(supminR, _data[i].minR());
avgminR += _data[i].minR();
}
avgminR /= _data.size();
Msg::Info("Minimum of R: in [%g, %g], avg=%g (R ~= measure of anisotropy)",
infminR, supminR, avgminR);
} }
void GMSH_AnalyseCurvedMeshPlugin::_printStatJacobian() void GMSH_AnalyseCurvedMeshPlugin::_printStatJacobian()
...@@ -377,32 +399,79 @@ void GMSH_AnalyseCurvedMeshPlugin::_printStatJacobian() ...@@ -377,32 +399,79 @@ void GMSH_AnalyseCurvedMeshPlugin::_printStatJacobian()
} }
double infminJ, supminJ, avgminJ; double infminJ, supminJ, avgminJ;
double infratJ, supratJ, avgratJ; double infratJ, supratJ, avgratJ;
int count = 0; double avgratJc;
int count = 0, countc = 0;
infminJ = infratJ = 1e10; infminJ = infratJ = 1e10;
supminJ = supratJ = -1e10; supminJ = supratJ = -1e10;
avgminJ = avgratJ = 0; avgminJ = avgratJ = avgratJc = 0;
for (unsigned int i = 0; i < _data.size(); ++i) { for (unsigned int i = 0; i < _data.size(); ++i) {
double q = 0; double q = 0;
if ( _data[i].minJ() > 0) q = _data[i].minJ() / _data[i].maxJ(); if ( _data[i].minJ() > 0) q = _data[i].minJ() / _data[i].maxJ();
if (q < 1-1e-3) { infminJ = std::min(infminJ, _data[i].minJ());
infminJ = std::min(infminJ, _data[i].minJ()); supminJ = std::max(supminJ, _data[i].minJ());
supminJ = std::max(supminJ, _data[i].minJ()); avgminJ += _data[i].minJ();
avgminJ += _data[i].minJ(); infratJ = std::min(infratJ, _data[i].minJ()/_data[i].maxJ());
infratJ = std::min(infratJ, _data[i].minJ()/_data[i].maxJ()); supratJ = std::max(supratJ, _data[i].minJ()/_data[i].maxJ());
supratJ = std::max(supratJ, _data[i].minJ()/_data[i].maxJ()); avgratJ += _data[i].minJ()/_data[i].maxJ();
avgratJ += _data[i].minJ()/_data[i].maxJ(); ++count;
++count; if (q < 1-1e-5) {
avgratJc += _data[i].minJ()/_data[i].maxJ();
++countc;
} }
} }
avgminJ /= count; avgminJ /= count;
avgratJ /= count; avgratJ /= count;
avgratJc /= countc;
Msg::Info("Minimum of Jacobian : in [%g, %g], avg=%g",
infminJ, supminJ, avgminJ);
Msg::Info("Ratio minJ/maxJ : in [%g, %g], avg=%g",
infratJ, supratJ, avgratJ);
Msg::Info(" avg=%g"
" (on the %d non-constant elements)",
avgratJc, count);
}
void GMSH_AnalyseCurvedMeshPlugin::_printStatScaledJac()
{
if (_data.empty()) {
Msg::Info("No stat to print.");
return;
}
double infminS, supminS, avgminS;
infminS = supminS = avgminS = _data[0].minS();
for (unsigned int i = 1; i < _data.size(); ++i) {
infminS = std::min(infminS, _data[i].minS());
supminS = std::max(supminS, _data[i].minS());
avgminS += _data[i].minS();
}
avgminS /= _data.size();
Msg::Info("Minimum of scaled Jacobian : in [%g, %g], avg=%g",
infminS, supminS, avgminS);
}
void GMSH_AnalyseCurvedMeshPlugin::_printStatIsotropy()
{
if (_data.empty()) {
Msg::Info("No stat to print.");
return;
}
double infminI, supminI, avgminI;
infminI = supminI = avgminI = _data[0].minI();
for (unsigned int i = 1; i < _data.size(); ++i) {
infminI = std::min(infminI, _data[i].minI());
supminI = std::max(supminI, _data[i].minI());
avgminI += _data[i].minI();
}
avgminI /= _data.size();
Msg::Info("Minimum of Jacobian: in [%g, %g], avg=%g (only the %d curved elem.)", Msg::Info("Minimum of isotropy : in [%g, %g], avg=%g",
infminJ, supminJ, avgminJ, count); infminI, supminI, avgminI);
Msg::Info("Ratio minJ/maxJ : in [%g, %g], avg=%g (only the %d curved elem.)",
infratJ, supratJ, avgratJ, count);
} }
// For testing // For testing
...@@ -412,14 +481,14 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinR(MElement *const *el, ...@@ -412,14 +481,14 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinR(MElement *const *el,
bool *straight) bool *straight)
{ {
_computedJ[el[0]->getDim()-1] = false; _computedJ[el[0]->getDim()-1] = false;
_computedR[el[0]->getDim()-1] = false; _computedI[el[0]->getDim()-1] = false;
_data.clear(); _data.clear();
_computeMinMaxJandValidity(el, numEl); _computeMinMaxJandValidity(el, numEl);
_computeMinR(el[0]->getDim()); _computeMinIsotropy(el[0]->getDim());
if (minR) { if (minR) {
for (unsigned int i = 0; i < _data.size(); ++i) { for (unsigned int i = 0; i < _data.size(); ++i) {
minR[i] = _data[i].minR(); minR[i] = _data[i].minI();
} }
} }
if (straight) { if (straight) {
......
...@@ -19,33 +19,32 @@ class data_elementMinMax ...@@ -19,33 +19,32 @@ class data_elementMinMax
{ {
private: private:
MElement *_el; MElement *_el;
double _minJ, _maxJ, _minR, _minRR; double _minJ, _maxJ, _minS, _minI;
public: public:
data_elementMinMax(MElement *e, data_elementMinMax(MElement *e,
double minJ = 2, double minJ = 2,
double maxJ = 0, double maxJ = 0,
double minR = -1, double minS = -1,
double minRR = -1) double minI = -1)
: _el(e), _minJ(minJ), _maxJ(maxJ), _minR(minR), _minRR(minRR) {} : _el(e), _minJ(minJ), _maxJ(maxJ), _minS(minS), _minI(minI) {}
void setMinR(double r) { _minR = r; } void setMinS(double r) { _minS = r; }
void setMinRR(double r) { _minRR = r; } void setMinI(double r) { _minI = r; }
MElement* element() { return _el; } MElement* element() { return _el; }
double minJ() { return _minJ; } double minJ() { return _minJ; }
double maxJ() { return _maxJ; } double maxJ() { return _maxJ; }
double minR() { return _minR; } double minS() { return _minS; }
double minRR() { return _minRR; } double minI() { return _minI; }
}; };
class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
{ {
private : private :
GModel *_m; GModel *_m;
double _threshold, _tol; double _threshold;
int _computeMetric;
// for 1d, 2d, 3d // for 1d, 2d, 3d
bool _computedR[3], _computedJ[3], _PViewJ[3], _PViewR[3]; bool _computedJ[3], _computedS[3], _computedI[3];
bool _msgHide; bool _PViewJ[3], _PViewS[3], _PViewI[3];
std::vector<data_elementMinMax> _data; std::vector<data_elementMinMax> _data;
...@@ -53,15 +52,15 @@ public : ...@@ -53,15 +52,15 @@ public :
GMSH_AnalyseCurvedMeshPlugin() GMSH_AnalyseCurvedMeshPlugin()
{ {
_m = NULL; _m = NULL;
_threshold = _tol = -1; _threshold = -1;
_computeMetric = -1;
for (int i = 0; i < 3; ++i) { for (int i = 0; i < 3; ++i) {
_computedR[i] = false;
_computedJ[i] = false; _computedJ[i] = false;
_computedS[i] = false;
_computedI[i] = false;
_PViewJ[i] = false; _PViewJ[i] = false;
_PViewR[i] = false; _PViewS[i] = false;
_PViewI[i] = false;
} }
_msgHide = true;
} }
std::string getName() const { return "AnalyseCurvedMesh"; } std::string getName() const { return "AnalyseCurvedMesh"; }
std::string getShortHelp() const std::string getShortHelp() const
...@@ -71,9 +70,8 @@ public : ...@@ -71,9 +70,8 @@ public :
std::string getHelp() const; std::string getHelp() const;
std::string getAuthor() const { return "Amaury Johnen"; } std::string getAuthor() const { return "Amaury Johnen"; }
int getNbOptions() const; int getNbOptions() const;
StringXNumber *getOption(int); StringXNumber* getOption(int);
PView *execute(PView *); PView* execute(PView *);
void setTol(double tol) { _tol = tol; }
// For testing // For testing
void computeMinJ(MElement *const *el, int numEl, double *minJ, bool *straight) void computeMinJ(MElement *const *el, int numEl, double *minJ, bool *straight)
...@@ -88,7 +86,7 @@ public : ...@@ -88,7 +86,7 @@ public :
} }
if (straight) { if (straight) {
for (unsigned int i = 0; i < _data.size(); ++i) { for (unsigned int i = 0; i < _data.size(); ++i) {
straight[i] = _data[i].maxJ() - _data[i].minJ() < _tol * 1e-1; straight[i] = _data[i].maxJ() - _data[i].minJ() < 1e-5;
} }
} }
_data = save; _data = save;
...@@ -96,24 +94,25 @@ public : ...@@ -96,24 +94,25 @@ public :
void computeMinR(MElement *const *el, int numEl, double *minR, bool *straight); void computeMinR(MElement *const *el, int numEl, double *minR, bool *straight);
void test(MElement *const *el, int numEl, int dim) void test(MElement *const *el, int numEl, int dim)
{ {
_tol = 1e-3;
std::vector<data_elementMinMax> save(_data); std::vector<data_elementMinMax> save(_data);
_data.clear(); _data.clear();
_computeMinMaxJandValidity(el, numEl); _computeMinMaxJandValidity(el, numEl);
Msg::Info("aaa"); Msg::Info("aaa");
Msg::Info("aaa"); Msg::Info("aaa");
_computeMinR(dim); _computeMinIsotropy(dim);
_data = save; _data = save;
} }
private : private :
void _computeMinMaxJandValidity(int dim); void _computeMinMaxJandValidity(int dim);
void _computeMinMaxJandValidity(MElement *const *, int numEl); void _computeMinMaxJandValidity(MElement *const *, int numEl);
void _computeMinR(int dim); void _computeMinScaledJac(int dim);
bool _hideWithThreshold(int askedDim); void _computeMinIsotropy(int dim);
void _printStatMetric(); int _hideWithThreshold(int askedDim, int whichMeasure);
void _printStatJacobian(); void _printStatJacobian();
void _printStatScaledJac();
void _printStatIsotropy();
}; };
#endif #endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment