diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 49cdd33cc775db732e6aa5efe2798a9c83e9e91b..108be656e24437a6fade04c53840970a65c5cb49 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -468,11 +468,11 @@ void MElement::getSignedJacobian(fullVector<double> &jacobian, int o)
getJacobianFuncSpace(o)->getSignedJacobian(nodesXYZ,jacobian);
}
-void MElement::getNodesCoord(fullMatrix<double> &nodesXYZ)
+void MElement::getNodesCoord(fullMatrix<double> &nodesXYZ) const
{
const int numNodes = getNumShapeFunctions();
for (int i = 0; i < numNodes; i++) {
- MVertex *v = getShapeFunctionNode(i);
+ const MVertex *v = getShapeFunctionNode(i);
nodesXYZ(i,0) = v->x();
nodesXYZ(i,1) = v->y();
nodesXYZ(i,2) = v->z();
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 1e7bbe3972967716fea3c959d8919dec7c048ea3..db684dbd13b84267e0771d43e201ff9183e12303 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -272,7 +272,7 @@ class MElement
double jac[3][3]; return getJacobian(u, v, w, jac);
}
void getSignedJacobian(fullVector<double> &jacobian, int o = -1);
- void getNodesCoord(fullMatrix<double> &nodesXYZ);
+ void getNodesCoord(fullMatrix<double> &nodesXYZ) const;
virtual int getNumShapeFunctions() const{ return getNumVertices(); }
virtual int getNumPrimaryShapeFunctions() { return getNumPrimaryVertices(); }
virtual const MVertex *getShapeFunctionNode(int i) const{ return getVertex(i); }
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index 6f80a63fadccbe9b5c0c228a054a22d6277ba77f..46b0efee7cc7538710ef7eda5552839ec1b9fc80 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -3,18 +3,20 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
-#include "Gmsh.h"
-#include "GModel.h"
-#include "GmshMessage.h"
-
-#include "polynomialBasis.h"
#include "AnalyseCurvedMesh.h"
+#include "GModel.h"
+#include "OS.h"
#include "Context.h"
-#include <queue>
+
#include <cmath>
-#include<fstream>
-#include "OS.h"
-#include "PView.h"
+#include <queue>
+#include <sstream>
+
+//#include "Gmsh.h"
+//#include "GmshMessage.h"
+//#include "polynomialBasis.h"
+//#include <fstream>
+//#include "PView.h"
#if defined(HAVE_OPENGL)
#include "drawContext.h"
@@ -26,10 +28,10 @@
namespace {
class BezierJacobian;
-struct lessMinB {
+struct lessMinVal {
bool operator()(BezierJacobian*, BezierJacobian*) const;
};
-struct lessMaxB {
+struct lessMax {
bool operator()(BezierJacobian*, BezierJacobian*) const;
};
@@ -37,7 +39,7 @@ class BezierJacobian
{
private:
fullVector<double> _jacBez;
- double _minJ, _maxJ, _minB, _maxB; //Extremum of Jac at corners and of bezier values
+ double _minL, _maxL, _minB, _maxB; //Extremum of Jac at corners and of bezier values
int _depthSub;
const JacobianBasis *_jfs;
@@ -46,28 +48,29 @@ public:
void subDivisions(fullVector<double> &vect) const
{_jfs->subdivideBezierCoeff(_jacBez, vect);}
+ bool boundsOk(double tol, double minL, double maxL) {
+ return (minL <= 0 || _minB > 0) &&
+ minL - _minB <= tol &&
+ _maxB - maxL <= tol;
+ }
+
inline int depth() const {return _depthSub;}
- inline double minJ() const {return _minJ;}
- inline double maxJ() const {return _maxJ;}
+ inline double minL() const {return _minL;}
+ inline double maxL() const {return _maxL;}
inline double minB() const {return _minB;}
inline double maxB() const {return _maxB;}
};
} // namespace
-//#define UNDEF_JAC_TAG -999
-//#define _ANALYSECURVEDMESH_BLAS_
-
-StringXNumber JacobianOptions_Number[] = {
- {GMSH_FULLRC, "Dim", NULL, -1},
- {GMSH_FULLRC, "Analysis", NULL, 2},
- {GMSH_FULLRC, "Effect (1)", NULL, 6},
- {GMSH_FULLRC, "JacBreak (1)", NULL, .0},
- {GMSH_FULLRC, "BezBreak (1)", NULL, .0},
- {GMSH_FULLRC, "MaxDepth (1,2)", NULL, 20},
- {GMSH_FULLRC, "Tolerance (2)", NULL, 1e-3}
+StringXNumber CurvedMeshOptions_Number[] = {
+ {GMSH_FULLRC, "Show: (0) Jacobian, (1) Metric", NULL, 1},
+ {GMSH_FULLRC, "Number of PView", NULL, 2},
+ {GMSH_FULLRC, "Hidding threshold", NULL, .1},
+ {GMSH_FULLRC, "Dimension of elements", NULL, -1},
+ {GMSH_FULLRC, "Recompute bounds", NULL, 0},
+ {GMSH_FULLRC, "Tolerance", NULL, 1e-3}
};
-
extern "C"
{
GMSH_Plugin *GMSH_RegisterAnalyseCurvedMeshPlugin()
@@ -77,318 +80,181 @@ extern "C"
}
int GMSH_AnalyseCurvedMeshPlugin::getNbOptions() const
{
- return sizeof(JacobianOptions_Number) / sizeof(StringXNumber);
+ return sizeof(CurvedMeshOptions_Number) / sizeof(StringXNumber);
}
-
StringXNumber *GMSH_AnalyseCurvedMeshPlugin::getOption(int iopt)
{
- return &JacobianOptions_Number[iopt];
+ return &CurvedMeshOptions_Number[iopt];
}
-
std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const
{
- return "Plugin(AnalyseCurvedMesh) check the jacobian of all elements of dimension 'Dim' or "
- "the greater model dimension if 'Dim' is either <0 or >3.\n\n "
- "Analysis : 0 do nothing\n"
- " +1 find invalid elements (*)\n"
- " +2 compute J_min and J_max of all elements and print some statistics\n\n"
- "Effect (for *) : 0 do nothing\n"
- " +1 print a list of invalid elements\n"
- " +2 print some statistics\n"
- " +4 hide valid elements (for GUI)\n\n"
- "MaxDepth = 0,1,...\n"
- " 0 : only sample the jacobian\n"
- " 1 : compute Bezier coefficients\n"
- " 2+ : execute a maximum of 1+ subdivision(s)\n\n"
- "JacBreak = [0,1[ : if a value of the jacobian <= 'JacBreak' is found, "
- "the element is said to be invalid\n\n"
- "BezBreak = [0,JacBreak[ : if all Bezier coefficients are > 'BezBreak', "
- "the element is said to be valid\n\n"
- "Tolerance = R+ , << 1 : tolerance (relatively to J_min and J_max) used during the computation of J_min and J_max";
+ return "Plugin(AnalyseCurvedMesh) analyse all elements of a given dimension. "
+ "It computes, min(J) where J is the scaled Jacobian determinant. Eventually, "
+ "it computes min(R) where R is the ratio between the smaller and the greater "
+ "of the eigenvalues of the metric. It creates one or more PView and hides "
+ "elements for which min({J, R}) < 'Hidding threshold'.\n"
+ "\n"
+ "Parameters:\n"
+ "\n"
+ "- Show [...] = {0, 1}: If 0, computes Jacobian and shows min(J). If 1, computes "
+ "Jacobian and metric and shows min(R).\n"
+ "\n"
+ "- Number of PView = {0, 1, 2}: If 1, create one PView with all elements. If 2, create "
+ "two PView, one with straight-sided elements and one with curved elements.\n"
+ "\n"
+ "- Hidding threshold = [0,1]: Hides all element for which min(R) or min(J) "
+ "is strictly greater than the threshold. If = 1, no effect, if = 0 hide "
+ "all elements except invalid.\n"
+ "\n"
+ "- Dimension = {-1, 1, 2, 3}: If = -1, analyse element of the greater dimension.\n"
+ "\n"
+ "- Recompute = {0,1}: If the mesh has changed, set to 1 to recompute the bounds.\n"
+ "\n"
+ "- Tolerance = ]0, 1[: Tolerance on the computation of min(R) or min(J)";
}
// Execution
PView *GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
{
_m = GModel::current();
- _dim = (int) JacobianOptions_Number[0].def;
- if (_dim < 0 || _dim > 3)
- _dim = _m->getDim();
- int analysis = (int) JacobianOptions_Number[1].def % 4;
- int toDo = (int) JacobianOptions_Number[2].def % 8;
- _maxDepth = (int) JacobianOptions_Number[5].def;
- _jacBreak = (double) JacobianOptions_Number[3].def;
- _bezBreak = (double) JacobianOptions_Number[4].def;
- _tol = (double) JacobianOptions_Number[6].def;
-
- if (analysis % 2) {
- double t = Cpu();
- Msg::Info("Starting validity check...");
- checkValidity(toDo);
- Msg::Info("Done validity check (%fs)", Cpu()-t);
+ _computeMetric = static_cast<int>(CurvedMeshOptions_Number[0].def);
+ _numPView = static_cast<int>(CurvedMeshOptions_Number[1].def);
+ _threshold = static_cast<double>(CurvedMeshOptions_Number[2].def);
+ _dim = static_cast<int>(CurvedMeshOptions_Number[3].def);
+ _recompute = static_cast<bool>(CurvedMeshOptions_Number[4].def);
+ _tol = static_cast<double>(CurvedMeshOptions_Number[5].def);
+
+ if (_dim < 0 || _dim > 3) _dim = _m->getDim();
+
+ if (_recompute) {
+ _data.clear();
+ _computedJ3D = false;
+ _computedJ2D = false;
+ _computedJ1D = false;
+ _computedR3D = false;
+ _computedR2D = false;
+ _msgHide = true;
+ _1PViewJ = false;
+ _2PViewJ = false;
+ _1PViewR = false;
+ _2PViewR = false;
}
- if (analysis / 2) {
- double t = Cpu();
- Msg::Info("Starting computation J_min, J_max...");
- std::map<int, std::vector<double> > data;
- computeMinMax(&data);
- new PView("Jmin", "ElementData", _m, data);
- Msg::Info("Done computation J_min, J_max (%fs)", Cpu()-t);
- }
- return 0;
-}
-
-void GMSH_AnalyseCurvedMeshPlugin::checkValidity(int toDo)
-{
- std::vector<MElement*> invalids;
- _numAnalysedEl = 0;
- _numInvalid = 0;
- _numValid = 0;
- _numUncertain = 0;
-
- switch (_dim) {
- case 3 :
- for(GModel::riter it = _m->firstRegion(); it != _m->lastRegion(); it++) {
- GRegion *r = *it;
- unsigned int numType[5] = {0, 0, 0, 0, 0};
- r->getNumMeshElements(numType);
-
- for(int type = 0; type < 5; type++) {
- MElement *const *el = r->getStartElementType(type);
- checkValidity(el, numType[type], invalids);
- _numAnalysedEl += numType[type];
- }
- }
- 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);
- checkValidity(el, numType[type], invalids);
- _numAnalysedEl += numType[type];
- }
- }
- 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);
- checkValidity(el, numElement, invalids);
- _numAnalysedEl += numElement;
- }
- break;
-
- default :
- Msg::Error("I can't analyse any element.");
+ if ((_dim == 1 && !_computedJ1D) ||
+ (_dim == 2 && !_computedJ2D) ||
+ (_dim == 3 && !_computedJ3D) ) {
+ double time = Cpu();
+ Msg::StatusBar(true, "Computing Jacobian for %dD elements...", _dim);
+ _computeMinMaxJandValidity();
+ Msg::StatusBar(true, "... Done computing Jacobian (%g seconds)", Cpu()-time);
+ _printStatJacobian();
}
- if (toDo % 2) {
- Msg::Info("Invalids elements :");
- Msg::Info("-------------------");
- for (unsigned int i = 0; i < invalids.size(); ++i) {
- Msg::Info(" %d", invalids[i]->getNum());
- }
- }
- if (toDo / 2 % 2) {
- Msg::Info("Found %d invalid elements and %d valid", _numInvalid, _numValid);
- if (_numUncertain) {
- Msg::Info("%d uncertain elements.", _numUncertain);
- if (_jacBreak < _bezBreak) {
- Msg::Info("'JacBreak' is smaller than 'BezBreak'. Change values in order to decrease the number of uncertain elements.");
- }
- else {
- Msg::Info("Increase MaxDepth in order to decrease the number of uncertain elements.");
- }
- }
- Msg::Info("%d elements analysed", _numAnalysedEl);
+ if (_computeMetric &&
+ ((_dim == 2 && !_computedR2D) ||
+ (_dim == 3 && !_computedR3D) )) {
+ double time = Cpu();
+ Msg::StatusBar(true, "Computing metric for %dD elements...", _dim);
+ _computeMinR();
+ Msg::StatusBar(true, "... Done computing metric (%g seconds)", Cpu()-time);
+ _printStatMetric();
}
- if (toDo / 4 % 2) {
- Msg::Info("Note: Valid elements are hidden. Change 'Effect' to disable this functionality.");
- Msg::Info("(To revert visibility : Tools > Visibility > Interactive > Show All)");
- hideValid_ShowInvalid(invalids);
- CTX::instance()->mesh.changed = (ENT_ALL);
-#if defined(HAVE_FLTK)
- FlGui::instance()->check();
-#endif
-#if defined(HAVE_OPENGL)
- drawContext::global()->draw();
-#endif
- }
-}
-void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
- int numEl,
- std::vector<MElement*> &invalids)
-{
- if (numEl < 1)
- return;
-
- const JacobianBasis *jfs = el[0]->getJacobianFuncSpace();
- if (!jfs) {
- Msg::Error("Jacobian function space not implemented for type of element %d", el[0]->getNum());
- return;
+ if (_computeMetric &&
+ (_dim == 1 ||
+ (_dim == 2 && !_computedR2D) ||
+ (_dim == 3 && !_computedR3D))) {
+ return 0;
}
- const int numSamplingPt = jfs->getNumJacNodes();
- const int numMapNodes = jfs->getNumMapNodes();
-
-#ifdef _ANALYSECURVEDMESH_BLAS_
- fullMatrix<double> jacobianB(numSamplingPt, numEl);
- fullMatrix<double> jacBezB(numSamplingPt, numEl);
- fullVector<double> jacBez, jacobian;
-
- fullMatrix<double> nodesX(numMapNodes, numEl), nodesY(numMapNodes, numEl), nodesZ(numMapNodes, numEl);
- for (int k = 0; k < numEl; ++k)
- for (int i = 0; i < numMapNodes; ++i)
- {
- MVertex *v = (el[k])->getShapeFunctionNode(i);
- nodesX(i,k) = v->x();
- nodesY(i,k) = v->y();
- nodesZ(i,k) = v->z();
- }
-
- jfs->getScaledJacobian(nodesX, nodesY, nodesZ, jacobianB);
- jfs->lag2Bez(jacobianB, jacBezB);
-#else
- fullVector<double> jacobian(numSamplingPt);
- fullVector<double> jacBez(numSamplingPt);
-#endif
-
- for (int k = 0; k < numEl; ++k) {
-
-#ifdef _ANALYSECURVEDMESH_BLAS_
- jacBez.setAsProxy(jacBezB, k);
- jacobian.setAsProxy(jacobianB, k);
-#else
- fullMatrix<double> nodesXYZ(numMapNodes,3);
- el[k]->getNodesCoord(nodesXYZ);
- jfs->getScaledJacobian(nodesXYZ,jacobian);
-#endif
- int i;
- for (i = 0; i < numSamplingPt && jacobian(i) > _jacBreak; ++i);
- if (i < numSamplingPt) {
- invalids.push_back(el[k]);
- ++_numInvalid;
- continue;
+ // Create PView
+
+ std::map<int, std::vector<double> > *dataPV1 = NULL;
+ std::map<int, std::vector<double> > *dataPV2 = NULL;
+ std::stringstream name1, name2;
+ if (_numPView == 1 &&
+ ((_computeMetric && !_1PViewR) || (!_computeMetric && !_1PViewJ))) {
+ dataPV1 = new std::map<int, std::vector<double> >();
+ if (_computeMetric) {
+ name1 << "min R";
+ _1PViewR = true;
}
-
- if (_maxDepth < 1) {
- invalids.push_back(el[k]);
- ++_numUncertain;
- continue;
+ else {
+ name1 << "min J";
+ _1PViewJ = true;
}
-
-#ifndef _ANALYSECURVEDMESH_BLAS_
- jfs->lag2Bez(jacobian, jacBez);
-#endif
-
- for (i = 0; i < jacBez.size() && jacBez(i) > _bezBreak; ++i);
- if (i >= jacBez.size()) {
- ++_numValid;
- continue;
+ for (unsigned int i = 0; i < _data.size(); ++i) {
+ MElement *const el = _data[i].element();
+ if (el->getDim() == _dim) {
+ double val = _computeMetric ? _data[i].minR() : _data[i].minJ();
+ (*dataPV1)[el->getNum()].push_back(val);
+ }
}
-
- if (_maxDepth < 2) {
- invalids.push_back(el[k]);
- ++_numUncertain;
- continue;
+ }
+ else if (_numPView == 2 &&
+ ((_computeMetric && !_2PViewR) || (!_computeMetric && !_2PViewJ))) {
+ dataPV1 = new std::map<int, std::vector<double> >();
+ dataPV2 = new std::map<int, std::vector<double> >();
+ if (_computeMetric) {
+ name1 << "min R straight";
+ name2 << "min R curved";
+ _2PViewR = true;
}
else {
- int result = subDivision(jfs, jacBez, _maxDepth-1);
- if (result < 0) {
- invalids.push_back(el[k]);
- ++_numInvalid;
- continue;
- }
- if (result > 0) {
- ++_numValid;
- continue;
+ name1 << "min J straight";
+ name2 << "min J curved";
+ _2PViewJ = true;
+ }
+ for (unsigned int i = 0; i < _data.size(); ++i) {
+ MElement *const el = _data[i].element();
+ if (el->getDim() == _dim) {
+ double val = _computeMetric ? _data[i].minR() : _data[i].minJ();
+ if (_data[i].maxJ() - _data[i].minJ() < _tol * 1e-2)
+ (*dataPV1)[el->getNum()].push_back(val);
+ else
+ (*dataPV2)[el->getNum()].push_back(val);
}
- invalids.push_back(el[k]);
- ++_numUncertain;
- continue;
}
}
-}
-
-int GMSH_AnalyseCurvedMeshPlugin::subDivision(const JacobianBasis *jb,
- const fullVector<double> &jacobian,
- int depth)
-{
- fullVector<double> newJacobian(jb->getNumSubNodes());
- jb->subdivideBezierCoeff(jacobian, newJacobian);
- for (int i = 0; i < jb->getNumDivisions(); i++)
- for (int j = 0; j < jb->getNumLagCoeff(); j++)
- if (newJacobian(i * jb->getNumJacNodes() + j) <= _jacBreak)
- return -1;
+ if (dataPV1) new PView(name1.str().c_str(), "ElementData", _m, *dataPV1);
+ if (dataPV2) new PView(name2.str().c_str(), "ElementData", _m, *dataPV2);
- int i = 0;
- while (i < newJacobian.size() && newJacobian(i) > _bezBreak)
- ++i;
- if (i >= newJacobian.size()) return 1;
+ //
- if (depth <= 1) {
- return 0;
+#if defined(HAVE_OPENGL)
+ if (_hideWithThreshold() && _msgHide) {
+ _msgHide = false;
+ Msg::Info("Note: Some elements have been hidden (param 'Hidding Threshold').");
+ Msg::Info(" (To revert visibility : Tools > Visibility > Interactive > Show All)");
}
- else{
- fullVector<double> subJacobian;
- std::vector<int> negTag, posTag;
- bool zeroTag = false;
-
- for (int i = 0; i < jb->getNumDivisions(); i++) {
- subJacobian.setAsProxy(newJacobian, i * jacobian.size(), jacobian.size());
- int tag = subDivision(jb, subJacobian, depth-1);
-
- if (tag < 0)
- negTag.push_back(tag);
- else if (tag > 0)
- posTag.push_back(tag);
- else
- zeroTag = true;
- }
-
- if (negTag.size() > 0)
- return *std::min_element(negTag.begin(), negTag.end()) - 1;
-
- if (zeroTag) return 0;
+ CTX::instance()->mesh.changed = (ENT_ALL);
+ drawContext::global()->draw();
+#endif
- return *std::max_element(posTag.begin(), posTag.end()) + 1;
- }
+ return 0;
}
-
-void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(std::map<int, std::vector<double> > *data)
+void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity()
{
- _numAnalysedEl = 0;
- _numInvalid = 0;
- _numValid = 0;
- _numUncertain = 0;
- _min_pJmin = .0, _avg_pJmin = .0;
- _min_ratioJ = .0, _avg_ratioJ = .0;
-
switch (_dim) {
case 3 :
- for(GModel::riter it = _m->firstRegion(); it != _m->lastRegion(); it++) {
+ if (_computedJ3D) break;
+ for (GModel::riter it = _m->firstRegion(); it != _m->lastRegion(); it++) {
GRegion *r = *it;
unsigned int numType[5] = {0, 0, 0, 0, 0};
r->getNumMeshElements(numType);
for(int type = 0; type < 5; type++) {
MElement *const *el = r->getStartElementType(type);
- computeMinMax(el, numType[type], data);
- _numAnalysedEl += numType[type];
+ _computeMinMaxJandValidity(el, numType[type]);
}
}
+ _computedJ3D = true;
break;
case 2 :
+ if (_computedJ2D) break;
for (GModel::fiter it = _m->firstFace(); it != _m->lastFace(); it++) {
GFace *f = *it;
unsigned int numType[3] = {0, 0, 0};
@@ -396,33 +262,30 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(std::map<int, std::vector<doubl
for (int type = 0; type < 3; type++) {
MElement *const *el = f->getStartElementType(type);
- computeMinMax(el, numType[type], data);
- _numAnalysedEl += numType[type];
+ _computeMinMaxJandValidity(el, numType[type]);
}
}
+ _computedJ2D = true;
break;
case 1 :
+ if (_computedJ1D) break;
for (GModel::eiter it = _m->firstEdge(); it != _m->lastEdge(); it++) {
GEdge *e = *it;
unsigned int numElement = e->getNumMeshElements();
MElement *const *el = e->getStartElementType(0);
- computeMinMax(el, numElement, data);
- _numAnalysedEl += numElement;
+ _computeMinMaxJandValidity(el, numElement);
}
+ _computedJ1D = true;
break;
default :
- Msg::Error("I can't analyse any element.");
+ Msg::Error("I can not analyse any element.");
return;
}
- Msg::Info("Minimum of min(~distortion) : %g", _min_pJmin);
- Msg::Info("Average of min(~distortion) : %g", _avg_pJmin / _numAnalysedEl);
- Msg::Info("Minimum of min(J) / max(J) : %g", _min_ratioJ);
- Msg::Info("Average of min(J) / max(J) : %g", _avg_ratioJ / _numAnalysedEl);
}
-void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl, std::map<int, std::vector<double> > *data)
+void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(MElement *const*el, int numEl)
{
if (numEl < 1)
return;
@@ -435,273 +298,249 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
const int numSamplingPt = jfs->getNumJacNodes();
const int numMapNodes = jfs->getNumMapNodes();
-
-#ifdef _ANALYSECURVEDMESH_BLAS_
- fullMatrix<double> jacobianB(numSamplingPt, numEl);
- fullMatrix<double> jacBezB(numSamplingPt, numEl);
- fullVector<double> jacBez, jacobian, jac1;
-
- fullMatrix<double> nodesX(numMapNodes, numEl), nodesY(numMapNodes, numEl), nodesZ(numMapNodes, numEl);
- for (int k = 0; k < numEl; ++k)
- for (int i = 0; i < numMapNodes; ++i)
- {
- MVertex *v = (el[k])->getShapeFunctionNode(i);
- nodesX(i,k) = v->x();
- nodesY(i,k) = v->y();
- nodesZ(i,k) = v->z();
- }
-
- jfs->getScaledJacobian(nodesX, nodesY, nodesZ, jacobianB);
- jfs->lag2Bez(jacobianB, jacBezB);
-#else
fullVector<double> jacobian(numSamplingPt);
fullVector<double> jacBez(numSamplingPt);
-#endif
fullVector<double> subJacBez(jfs->getNumSubNodes());
- _min_pJmin = 1.7e308;
- _min_ratioJ = 1.7e308;
-
for (int k = 0; k < numEl; ++k) {
- if (el[k]->getType() == TYPE_PYR && el[k]->getPolynomialOrder() > 1) {
- static int i = 0;
- if (++i == 1) {
- Msg::Error("High-order pyramids skipped (subdivision not implemented).");
- Msg::Error("Subdivision will come soon :) (before may 2014)");
- }
- continue;
- }
-
-#ifdef _ANALYSECURVEDMESH_BLAS_
- jacBez.setAsProxy(jacBezB, k);
- jacobian.setAsProxy(jacobianB, k);
-#else
fullMatrix<double> nodesXYZ(numMapNodes,3);
el[k]->getNodesCoord(nodesXYZ);
jfs->getScaledJacobian(nodesXYZ,jacobian);
jfs->lag2Bez(jacobian, jacBez);
-#endif
- double minJ, maxJ = minJ = jacobian(0);
+ BezierJacobian *bj = new BezierJacobian(jacBez, jfs, 0);
+ std::vector<BezierJacobian*> heap;
+ heap.push_back(bj);
- for (int i = 1; i < numSamplingPt; ++i) {
- if (jacobian(i) < minJ) minJ = jacobian(i);
- if (jacobian(i) > maxJ) maxJ = jacobian(i);
- }
+ double minL = bj->minL(), maxL = bj->maxL();
+ int currentDepth = 0;
- double minB, maxB = minB = jacBez(0);
+ while (!heap[0]->boundsOk(_tol, minL, maxL)) {
+ bj = heap[0];
+ pop_heap(heap.begin(), heap.end(), lessMinVal());
+ heap.pop_back();
- for (int i = 1; i < numSamplingPt; ++i) {
- if (jacBez(i) < minB) minB = jacBez(i);
- if (jacBez(i) > maxB) maxB = jacBez(i);
- }
+ bj->subDivisions(subJacBez);
+ currentDepth = bj->depth() + 1;
+ if (currentDepth > 20) Msg::Error("depth is too damn high !");
- if (_maxDepth > 1 &&
- (minJ - minB > _tol * (std::abs(minJ) + std::abs(minB)) / 2 ||
- maxB - maxJ > _tol * (std::abs(maxJ) + std::abs(maxB)) / 2 )) {
-
- BezierJacobian *bj = new BezierJacobian(jacBez, jfs, 0);
- std::set<BezierJacobian*> setBJ;
- std::priority_queue<BezierJacobian*, std::vector<BezierJacobian*>, lessMinB> pqMin;
- std::priority_queue<BezierJacobian*, std::vector<BezierJacobian*>, lessMaxB> pqMax;
- setBJ.insert(bj);
- pqMin.push(bj);
-
- int currentDepth = 0;
- while(minJ - minB > _tol * (std::abs(minJ) + std::abs(minB)) / 2 &&
- pqMin.top()->depth() < _maxDepth-1) {
- bj = pqMin.top();
- bj->subDivisions(subJacBez);
- currentDepth = bj->depth() + 1;
- setBJ.erase(bj);
- pqMin.pop();
- delete bj;
-
- for (int i = 0; i < jfs->getNumDivisions(); i++) {
- jacBez.setAsProxy(subJacBez, i * numSamplingPt, numSamplingPt);
- bj = new BezierJacobian(jacBez, jfs, currentDepth);
- pqMin.push(bj);
- setBJ.insert(bj);
- minJ = std::min(minJ, bj->minJ());
- maxJ = std::max(maxJ, bj->maxJ());
- }
+ for (int i = 0; i < jfs->getNumDivisions(); i++) {
+ jacBez.setAsProxy(subJacBez, i * numSamplingPt, numSamplingPt);
+ bj = new BezierJacobian(jacBez, jfs, currentDepth);
+ minL = std::min(minL, bj->minL());
+ maxL = std::max(maxL, bj->maxL());
- minB = minJ;
- maxB = maxJ;
- std::set<BezierJacobian*>::iterator it;
- for (it = setBJ.begin(); it != setBJ.end(); ++it) {
- minB = std::min(minB, (*it)->minB());
- maxB = std::max(maxB, (*it)->maxB());
- }
+ heap.push_back(bj);
+ push_heap(heap.begin(), heap.end(), lessMinVal());
}
+ }
- while (pqMin.size() > 0) {
- bj = pqMin.top();
- pqMin.pop();
- pqMax.push(bj);
- }
+ make_heap(heap.begin(), heap.end(), lessMax());
- while(maxB - maxJ > _tol * (std::abs(maxJ) + std::abs(maxB)) / 2 &&
- pqMax.top()->depth() < _maxDepth-1) {
- bj = pqMax.top();
- bj->subDivisions(subJacBez);
- currentDepth = bj->depth() + 1;
- setBJ.erase(bj);
- pqMax.pop();
- delete bj;
-
- for (int i = 0; i < jfs->getNumDivisions(); i++) {
- jacBez.setAsProxy(subJacBez, i * numSamplingPt, numSamplingPt);
- bj = new BezierJacobian(jacBez, jfs, currentDepth);
- pqMax.push(bj);
- setBJ.insert(bj);
- minJ = std::min(minJ, bj->minJ());
- maxJ = std::max(maxJ, bj->maxJ());
- }
+ while (!heap[0]->boundsOk(_tol, minL, maxL)) {
+ bj = heap[0];
+ pop_heap(heap.begin(), heap.end(), lessMax());
+ heap.pop_back();
- minB = minJ;
- maxB = maxJ;
- std::set<BezierJacobian*>::iterator it;
- for (it = setBJ.begin(); it != setBJ.end(); ++it) {
- minB = std::min(minB, (*it)->minB());
- maxB = std::max(maxB, (*it)->maxB());
- }
- }
+ bj->subDivisions(subJacBez);
+ currentDepth = bj->depth() + 1;
+ if (currentDepth > 20) Msg::Error("depth is too damn high !");
+
+ for (int i = 0; i < jfs->getNumDivisions(); i++) {
+ jacBez.setAsProxy(subJacBez, i * numSamplingPt, numSamplingPt);
+ bj = new BezierJacobian(jacBez, jfs, currentDepth);
+ minL = std::min(minL, bj->minL());
+ maxL = std::max(maxL, bj->maxL());
- while (pqMax.size() > 0) {
- bj = pqMax.top();
- pqMax.pop();
- delete bj;
+ heap.push_back(bj);
+ push_heap(heap.begin(), heap.end(), lessMax());
}
}
- if (data){
- if (1-minB <= _tol * minJ && maxB-1 <= _tol * maxB)
- (*data)[el[k]->getNum()].push_back(1.);
- else if (1-minB <= 1e-8)
- (*data)[el[k]->getNum()].push_back(1.);
- else
- (*data)[el[k]->getNum()].push_back(minB);
+ double minB = minL;
+ double maxB = maxL;
+ for (unsigned int i = 0; i < heap.size(); ++i) {
+ minB = std::min(minB, heap[i]->minB());
+ maxB = std::max(maxB, heap[i]->maxB());
}
- _min_pJmin = std::min(_min_pJmin, minB);
- _avg_pJmin += minB;
- _min_ratioJ = std::min(_min_ratioJ, minB/maxB);
- _avg_ratioJ += minB/maxB;
+
+ // TODO REMOVE #include "BasisFactory.h"
+
+ /*if (el[k]->getNum() == 26776) {
+ Msg::Info("%g %g %g %g", minB, minL, maxL, maxB);
+
+ /*const nodalBasis* nb = BasisFactory::getNodalBasis(el[k]->getTypeForMSH());
+ int tag1 = ElementType::getTag(el[k]->getType(), 1);
+ const nodalBasis* nb1 = BasisFactory::getNodalBasis(tag1);
+
+ //double *f;
+ //f = new double[nb1->points.size1()];
+ fullMatrix<double> f;
+ nb1->f(nb->points, f);
+ Msg::Info("SIZES [%d %d] (%d, %d)", f.size1(), f.size2(), nb->points.size1(), nb1->points.size1());
+
+ for (int i = 0; i < nb->points.size1(); ++i) {
+ double X[3] = {0, 0, 0};
+ for (int j = 0; j < f.size2(); ++j) {
+ X[0] += el[k]->getVertex(j)->x() * f(i, j);
+ X[1] += el[k]->getVertex(j)->y() * f(i, j);
+ X[2] += el[k]->getVertex(j)->z() * f(i, j);
+ }
+ if (std::abs(X[0] - el[k]->getVertex(i)->x()) < 1e-12 &&
+ std::abs(X[1] - el[k]->getVertex(i)->y()) < 1e-12 &&
+ std::abs(X[2] - el[k]->getVertex(i)->z()) < 1e-12);
+ else {
+ Msg::Info("aaarg %d/%d : [%g %g %g] ([%g %g %g])", i+1, nb->points.size1(),
+ el[k]->getVertex(i)->x(), el[k]->getVertex(i)->y(), el[k]->getVertex(i)->z(),
+ X[0], X[1], X[2]);
+ }
+ }* /
+ }*/
+
+ _data.push_back(CurvedMeshPluginData(el[k], minB, maxB));
}
}
-void GMSH_AnalyseCurvedMeshPlugin::hideValid_ShowInvalid(std::vector<MElement*> &invalids)
+void GMSH_AnalyseCurvedMeshPlugin::_computeMinR()
{
- unsigned int current = 0;
- invalids.push_back(NULL);
+ if (_dim == 1 ||
+ (_dim == 2 && _computedR2D) ||
+ (_dim == 3 && _computedR3D)) {
+ return;
+ }
- switch (_dim) {
- case 3 :
- for(GModel::riter it = _m->firstRegion(); it != _m->lastRegion(); it++) {
- GRegion *r = *it;
- unsigned int numType[5] = {0, 0, 0, 0, 0};
- r->getNumMeshElements(numType);
-
- for(int type = 0; type < 5; type++) {
- MElement *const *el = r->getStartElementType(type);
- for (unsigned int i = 0; i < numType[type]; ++i) {
- if (el[i] == invalids[current]) {
- ++current;
- el[i]->setVisibility(1);
- }
- else
- el[i]->setVisibility(0);
- }
- }
+ MetricBasis::setTol(_tol);
+
+ double initial, time = initial = Cpu();
+ int percentage = 0, nextCheck = 0;
+
+ for (unsigned int i = 0; i < _data.size(); ++i) {
+ MElement *const el = _data[i].element();
+ if (el->getDim() == _dim) {
+ if (_data[i].minJ() <= 0)
+ _data[i].setMinR(0);
+ else if (_data[i].maxJ() - _data[i].minJ() < _tol * 1e-2)
+ _data[i].setMinR(MetricBasis::sampleR(el, 0));
+ else
+ _data[i].setMinR(MetricBasis::boundMinR(el));
}
- 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);
- for (unsigned int i = 0; i < numType[type]; ++i) {
- if (el[i] == invalids[current]) {
- ++current;
- el[i]->setVisibility(1);
- }
- else
- el[i]->setVisibility(0);
- }
+ if (i >= nextCheck) {
+ nextCheck += _data.size() / 100;
+ if (Cpu() - time > 10. && i*100/_data.size() > percentage + 4) {
+ percentage = i*100/_data.size();
+ time = Cpu();
+ Msg::StatusBar(true, "%d%% (remaining time ~%g secondes)",
+ percentage, (Cpu()-initial) / (i+1) * (_data.size() - i-1));
}
}
- 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);
- for (unsigned int i = 0; i < numElement; ++i) {
- if (el[i] == invalids[current]) {
- ++current;
- el[i]->setVisibility(1);
- }
- else
- el[i]->setVisibility(0);
+ }
+
+ if (_dim == 2) _computedR2D = true;
+ if (_dim == 3) _computedR3D = true;
+}
+
+bool GMSH_AnalyseCurvedMeshPlugin::_hideWithThreshold()
+{
+ if (_threshold >= 1.) return false;
+
+ bool ans = false;
+ for (unsigned int i = 0; i < _data.size(); ++i) {
+ MElement *const el = _data[i].element();
+ if (el->getDim() == _dim) {
+ if (( _computeMetric && _data[i].minR() > _threshold) ||
+ (!_computeMetric && _data[i].minJ() > _threshold) ) {
+ el->setVisibility(0);
+ ans = true;
+ }
+ else {
+ el->setVisibility(1);
}
}
- break;
-
- default :
- break;
}
+ return ans;
+}
- invalids.pop_back();
-
- switch (_dim) {
- case 3 :
- for (GModel::fiter it = _m->firstFace(); it != _m->lastFace(); it++)
- (*it)->setVisibility(0);
+void GMSH_AnalyseCurvedMeshPlugin::_printStatMetric()
+{
+ if (_data.empty()) {
+ Msg::Info("No stat to print...");
+ return;
+ }
+ double infminR, supminR, avgminR;
+ infminR = supminR = avgminR = _data[0].minR();
- case 2 :
- for (GModel::eiter it = _m->firstEdge(); it != _m->lastEdge(); it++)
- (*it)->setVisibility(0);
+ 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();
- case 1 :
- for (GModel::viter it = _m->firstVertex(); it != _m->lastVertex(); it++)
- (*it)->setVisibility(0);
+ Msg::Info("Minimum of R: in [%g, %g], avg=%g (R ~= measure of anisotropy)",
+ infminR, supminR, avgminR);
+}
- default :
- break;
+void GMSH_AnalyseCurvedMeshPlugin::_printStatJacobian()
+{
+ if (_data.empty()) {
+ Msg::Info("No stat to print...");
+ return;
}
+ double infminJ, supminJ, avgminJ;
+ double infratJ, supratJ, avgratJ;
+ int count = 0;
+
+ infminJ = infratJ = 1e10;
+ supminJ = supratJ = -1e10;
+ avgminJ = avgratJ = 0;
+
+ for (unsigned int i = 0; i < _data.size(); ++i) {
+ if (_data[i].maxJ() - _data[i].minJ() > _tol * 1e-2) {
+ infminJ = std::min(infminJ, _data[i].minJ());
+ supminJ = std::max(supminJ, _data[i].minJ());
+ avgminJ += _data[i].minJ();
+ infratJ = std::min(infratJ, _data[i].minJ()/_data[i].maxJ());
+ supratJ = std::max(supratJ, _data[i].minJ()/_data[i].maxJ());
+ avgratJ += _data[i].minJ()/_data[i].maxJ();
+ ++count;
+ }
+ }
+ avgminJ /= count;
+ avgratJ /= count;
+
+ Msg::Info("Minimum of Jacobian: in [%g, %g], avg=%g (only the %d curved elem.)",
+ infminJ, supminJ, avgminJ, count);
+ Msg::Info("Ratio minJ/maxJ : in [%g, %g], avg=%g (only the %d curved elem.)",
+ infratJ, supratJ, avgratJ, count);
}
-BezierJacobian::BezierJacobian(fullVector<double> &v, const JacobianBasis *jfs, int depth)
+BezierJacobian::BezierJacobian(fullVector<double> &v,
+ const JacobianBasis *jfs, int depth)
{
_jacBez = v;
_depthSub = depth;
_jfs = jfs;
- _minJ = _maxJ = v(0);
+ _minL = _maxL = v(0);
int i = 1;
for (; i < jfs->getNumLagCoeff(); i++) {
- if (_minJ > v(i)) _minJ = v(i);
- if (_maxJ < v(i)) _maxJ = v(i);
+ if (_minL > v(i)) _minL = v(i);
+ if (_maxL < v(i)) _maxL = v(i);
}
- _minB = _minJ;
- _maxB = _maxJ;
+ _minB = _minL;
+ _maxB = _maxL;
for (; i < v.size(); i++) {
if (_minB > v(i)) _minB = v(i);
if (_maxB < v(i)) _maxB = v(i);
}
}
-bool lessMinB::operator()(BezierJacobian *bj1, BezierJacobian *bj2) const
+bool lessMinVal::operator()(BezierJacobian *bj1, BezierJacobian *bj2) const
{
return bj1->minB() > bj2->minB();
}
-bool lessMaxB::operator()(BezierJacobian *bj1, BezierJacobian *bj2) const
+bool lessMax::operator()(BezierJacobian *bj1, BezierJacobian *bj2) const
{
return bj1->maxB() < bj2->maxB();
}
diff --git a/Plugin/AnalyseCurvedMesh.h b/Plugin/AnalyseCurvedMesh.h
index b9153a68bde5d4e82cb3cfda94ebf3ce42da6e84..0fe626a4a055e333a1196e74c648d5a6b0ea880b 100644
--- a/Plugin/AnalyseCurvedMesh.h
+++ b/Plugin/AnalyseCurvedMesh.h
@@ -7,7 +7,10 @@
#define _ANALYSECURVEDMESH_H_
#include "Plugin.h"
+#include "JacobianBasis.h"
+#include "MetricBasis.h"
#include "MElement.h"
+
#include <vector>
extern "C"
@@ -15,41 +18,71 @@ extern "C"
GMSH_Plugin *GMSH_RegisterAnalyseCurvedMeshPlugin();
}
+class CurvedMeshPluginData
+{
+private:
+ MElement *_el;
+ double _minJ, _maxJ, _minR;
+
+public:
+ CurvedMeshPluginData(MElement *e,
+ double minJ = 2,
+ double maxJ = 0,
+ double minR = -1)
+ : _el(e), _minJ(minJ), _maxJ(maxJ), _minR(minR) {}
+
+ void setMinR(double r) {_minR = r;}
+ MElement* element() {return _el;}
+ double minJ() {return _minJ;}
+ double maxJ() {return _maxJ;}
+ double minR() {return _minR;}
+};
+
class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
{
- private :
- int _dim;
- GModel *_m;
- int _maxDepth;
- double _jacBreak, _bezBreak, _tol;
-
- int _numAnalysedEl;
- int _numInvalid, _numValid, _numUncertain;
- double _min_pJmin, _avg_pJmin;
- double _min_ratioJ, _avg_ratioJ;
-
- public :
- GMSH_AnalyseCurvedMeshPlugin(){}
- std::string getName() const { return "AnalyseCurvedMesh"; }
- std::string getShortHelp() const {
- return "Check validity of elements and/or compute min&max of the jacobian";
- }
- std::string getHelp() const;
- std::string getAuthor() const { return "Amaury Johnen"; }
- int getNbOptions() const;
- StringXNumber *getOption(int);
- PView *execute(PView *);
- void checkValidity(MElement *const *, int numEl, std::vector<MElement*> &invalids);
-
- //void storeJMin(MElement *const *, int numEl, std::map<int, std::vector<double> > &data);
-
-
- private :
- void checkValidity(int toDo);
- void computeMinMax(std::map<int, std::vector<double> > *data = 0);
- void computeMinMax(MElement *const *, int numEl, std::map<int, std::vector<double> > *data = 0);
- int subDivision(const JacobianBasis *, const fullVector<double>&, int depth);
- void hideValid_ShowInvalid(std::vector<MElement*> &invalids);
+private :
+ int _dim;
+ GModel *_m;
+ double _threshold, _tol;
+ int _numPView, _computeMetric;
+ bool _recompute;
+ bool _computedR3D, _computedR2D;
+ bool _computedJ3D, _computedJ2D, _computedJ1D;
+ bool _1PViewJ, _2PViewJ, _1PViewR, _2PViewR;
+ bool _msgHide;
+
+ std::vector<CurvedMeshPluginData> _data;
+
+public :
+ GMSH_AnalyseCurvedMeshPlugin() {
+ _computedR3D = false;
+ _computedR2D = false;
+ _computedJ3D = false;
+ _computedJ2D = false;
+ _computedJ1D = false;
+ _msgHide = true;
+ _1PViewJ = false;
+ _2PViewJ = false;
+ _1PViewR = false;
+ _2PViewR = false;
+ }
+ std::string getName() const { return "AnalyseCurvedMesh"; }
+ std::string getShortHelp() const {
+ return "Compute bounds on Jacobian and metric quality.";
+ }
+ std::string getHelp() const;
+ std::string getAuthor() const { return "Amaury Johnen"; }
+ int getNbOptions() const;
+ StringXNumber *getOption(int);
+ PView *execute(PView *);
+
+private :
+ void _computeMinMaxJandValidity();
+ void _computeMinMaxJandValidity(MElement *const *, int numEl);
+ void _computeMinR();
+ bool _hideWithThreshold();
+ void _printStatMetric();
+ void _printStatJacobian();
};
#endif