diff --git a/Common/Context.h b/Common/Context.h
index 8b43003b95de128629e48eda5a5f908615779484..e9121396c3ead20d5a6057fc4162cfd22a45ba16 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -236,8 +236,10 @@ class CTX {
int mouseSelection, mouseHoverMeshes, pickElements;
// disable some warnings for expert users?
int expertMode;
+#if defined(HAVE_VISUDEV)
// Enable heavy visualization capabilities (for development purpose)
int heavyVisu;
+#endif
// dynamic: equal to 1 while gmsh is printing
int printing;
// hide all unselected entities?
diff --git a/Common/Options.h b/Common/Options.h
index 53a35ad8bb72b56ac60a329035d6ee7f62334b34..12b9a3e5e2c366be1bf9cad79154e25d3234dbed 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -285,7 +285,9 @@ double opt_general_zoom_factor(OPT_ARGS_NUM);
double opt_general_expert_mode(OPT_ARGS_NUM);
double opt_general_stereo_mode(OPT_ARGS_NUM);
double opt_general_camera_mode(OPT_ARGS_NUM);
+#if defined(HAVE_VISUDEV)
double opt_general_heavy_visualization(OPT_ARGS_NUM);
+#endif
double opt_general_eye_sep_ratio(OPT_ARGS_NUM);
double opt_general_focallength_ratio(OPT_ARGS_NUM);
double opt_general_camera_aperture(OPT_ARGS_NUM);
diff --git a/Mesh/qualityMeasuresJacobian.cpp b/Mesh/qualityMeasuresJacobian.cpp
index 97eb1da16e508001473a02a6e1d30d1e91510fe8..2742a079780a064ee168a0c08c8d679246bf85de 100644
--- a/Mesh/qualityMeasuresJacobian.cpp
+++ b/Mesh/qualityMeasuresJacobian.cpp
@@ -405,6 +405,32 @@ double minICNMeasure(MElement *el,
}
void sampleIGEMeasure(MElement *el, int deg, double &min, double &max)
+{
+ fullVector<double> ige;
+ sampleIGEMeasure(el, deg, ige);
+
+ min = std::numeric_limits<double>::infinity();
+ max = -min;
+ for (int i = 0; i < ige.size(); ++i) {
+ min = std::min(min, ige(i));
+ max = std::max(max, ige(i));
+ }
+}
+
+void sampleJacobian(MElement *el, int deg, fullVector<double> &jac,
+ const fullMatrix<double> *normals)
+{
+ FuncSpaceData sampleSpace = FuncSpaceData(el, deg);
+ const JacobianBasis *jacBasis = BasisFactory::getJacobianBasis(sampleSpace);
+
+ fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
+ el->getNodesCoord(nodesXYZ);
+
+ jac.resize(jacBasis->getNumJacNodes());
+ jacBasis->getSignedJacobian(nodesXYZ, jac, normals);
+}
+
+void sampleIGEMeasure(MElement *el, int deg, fullVector<double> &ige)
{
fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
el->getNodesCoord(nodesXYZ);
@@ -427,7 +453,7 @@ void sampleIGEMeasure(MElement *el, int deg, double &min, double &max)
jacDetSpace = FuncSpaceData(el, true, deg-1, 1, &serendipFalse);
break;
default:
- Msg::Error("ICN not implemented for type of element %d", el->getType());
+ Msg::Error("IGE not implemented for type of element %d", el->getType());
return;
}
@@ -439,21 +465,12 @@ void sampleIGEMeasure(MElement *el, int deg, double &min, double &max)
fullMatrix<double> coeffMatLag(gradBasis->getNumSamplingPoints(), 9);
gradBasis->getAllGradientsFromNodes(nodesXYZ, coeffMatLag);
+ if (el->getDim() == 2) coeffMatLag.resize(coeffMatLag.size1(), 6, false);
fullMatrix<double> v;
computeCoeffLengthVectors_(coeffMatLag, v, type);
- fullVector<double> ige;
computeIGE_(coeffDeterminant, v, ige, type);
-
- min = std::numeric_limits<double>::infinity();
- max = -min;
- for (int i = 0; i < ige.size(); ++i) {
- min = std::min(min, ige(i));
- max = std::max(max, ige(i));
- }
-
- return;
}
double minSampledICNMeasure(MElement *el, int deg)//fordebug
diff --git a/Mesh/qualityMeasuresJacobian.h b/Mesh/qualityMeasuresJacobian.h
index 6c3b62992a1df10b80eed5c812d02207a1c11103..1bb4fc2a35d7dc8fcf8df15649c6fecf4e46dced 100644
--- a/Mesh/qualityMeasuresJacobian.h
+++ b/Mesh/qualityMeasuresJacobian.h
@@ -24,6 +24,10 @@ double minICNMeasure(MElement *el,
bool knownValid = false,
bool reversedOk = false);
void sampleIGEMeasure(MElement *el, int order, double &min, double &max);
+void sampleJacobian(MElement *el, int order, fullVector<double> &jac,
+ const fullMatrix<double> *normals = NULL);
+void sampleIGEMeasure(MElement *el, int order, fullVector<double> &ige);
+void sampleICNMeasure(MElement *el, int order, fullVector<double> &icn);
double minSampledICNMeasure(MElement *el, int order);//fordebug
double minSampledIGEMeasure(MElement *el, int order);//fordebug
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index c798ab731f8d45ea8a26b9f7c94ef87afa4643af..cc33b71d9966b15b7edbcae1eedab29aba0c26c3 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -7,6 +7,7 @@
#if defined(HAVE_MESH)
+
#include "AnalyseCurvedMesh.h"
#include "OS.h"
#include "Context.h"
@@ -19,6 +20,9 @@
#include <sstream>
#include <fstream>
#include "qualityMeasuresJacobian.h"
+#if defined(HAVE_VISUDEV)
+#include "BasisFactory.h"
+#endif
class bezierBasis;
@@ -30,6 +34,9 @@ StringXNumber CurvedMeshOptions_Number[] = {
{GMSH_FULLRC, "Draw PView", NULL, 0},
{GMSH_FULLRC, "Recompute", NULL, 0},
{GMSH_FULLRC, "Dimension of elements", NULL, -1}
+#if defined(HAVE_VISUDEV)
+ ,{GMSH_FULLRC, "Element to draw quality", NULL, 0}
+#endif
};
extern "C"
@@ -97,13 +104,25 @@ std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const
PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
{
_m = GModel::current();
- bool computeJac = static_cast<int>(CurvedMeshOptions_Number[0].def);
- bool computeIGE = static_cast<int>(CurvedMeshOptions_Number[1].def);
- bool computeICN = static_cast<int>(CurvedMeshOptions_Number[2].def);
- _threshold = static_cast<double>(CurvedMeshOptions_Number[3].def);
- bool drawPView = static_cast<int>(CurvedMeshOptions_Number[4].def);
- bool recompute = static_cast<bool>(CurvedMeshOptions_Number[5].def);
- int askedDim = static_cast<int>(CurvedMeshOptions_Number[6].def);
+ int computeJac = static_cast<int>(CurvedMeshOptions_Number[0].def);
+ int computeIGE = static_cast<int>(CurvedMeshOptions_Number[1].def);
+ int computeICN = static_cast<int>(CurvedMeshOptions_Number[2].def);
+ _threshold = CurvedMeshOptions_Number[3].def;
+ bool drawPView = static_cast<bool>(CurvedMeshOptions_Number[4].def);
+ bool recompute = static_cast<bool>(CurvedMeshOptions_Number[5].def);
+ int askedDim = static_cast<int>(CurvedMeshOptions_Number[6].def);
+
+#if defined(HAVE_VISUDEV)
+ _pwJac = computeJac/2;
+ _pwIGE = computeIGE/2;
+ _pwICN = computeICN/2;
+
+ _numElementToScan = static_cast<int>(CurvedMeshOptions_Number[7].def);
+ _viewOrder = 0;
+ _dataPViewJac.clear();
+ _dataPViewIGE.clear();
+ _dataPViewICN.clear();
+#endif
if (askedDim < 0 || askedDim > 4) askedDim = _m->getDim();
@@ -163,6 +182,10 @@ PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
if (printStatS) _printStatIGE();
if (printStatI) _printStatICN();
+#if defined(HAVE_VISUDEV)
+ _createPViewPointwise();
+#endif
+
// Create PView
if (drawPView)
for (int dim = 1; dim <= 3; ++dim) {
@@ -342,6 +365,10 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim)
_data.push_back(data_elementMinMax(el, min, max));
if (min < 0 && max < 0) ++cntInverted;
progress.next();
+
+#if defined(HAVE_VISUDEV)
+ _computePointwiseQuantities(el, normals);
+#endif
}
delete normals;
}
@@ -515,4 +542,84 @@ void GMSH_AnalyseCurvedMeshPlugin::_printStatICN()
infminI, avgminI, supminI);
}
+#if defined(HAVE_VISUDEV)
+void GMSH_AnalyseCurvedMeshPlugin::_computePointwiseQuantities(MElement *el,
+ const fullMatrix<double> *normals) {
+ if (_numElementToScan != -1 && el->getNum() != _numElementToScan) return;
+
+ if (!_type2tag[el->getType()])
+ _type2tag[el->getType()] = el->getTypeForMSH();
+ else
+ // Skip if element tag is different from previous elements of same type
+ if (_type2tag[el->getType()] != el->getTypeForMSH()) return;
+
+ static fullVector<double> tmpVector;
+ int num = el->getNum();
+
+ if (!_viewOrder) {
+// _viewOrder = std::min(10, 2 * el->getPolynomialOrder());
+ _viewOrder = 9;
+ }
+
+ if (_pwJac) {
+ jacobianBasedQuality::sampleJacobian(el, _viewOrder, tmpVector, normals);
+ std::vector<double> &vec = _dataPViewJac[num];
+ for (int j = 0; j < tmpVector.size(); ++j)
+ vec.push_back(tmpVector(j));
+ }
+
+ if (_pwIGE) {
+ jacobianBasedQuality::sampleIGEMeasure(el, _viewOrder, tmpVector);
+ std::vector<double> &vec = _dataPViewIGE[num];
+ for (int j = 0; j < tmpVector.size(); ++j)
+ vec.push_back(tmpVector(j));
+ }
+
+ if (_pwICN) {
+// jacobianBasedQuality::sampleICNMeasure(el, _viewOrder, tmpVector);
+// std::vector<double> &vec = _dataPViewICN[num];
+// for (int j = 0; j < tmpVector.size(); ++j)
+// vec.push_back(tmpVector(j));
+ }
+}
+
+void GMSH_AnalyseCurvedMeshPlugin::_setInterpolationMatrices(PView *view)
+{
+ PViewData *viewData = view->getData();
+ for (int type = 0; type < 20; ++type) {
+ if (!_type2tag[type]) continue;
+ viewData->deleteInterpolationMatrices(type);
+ const nodalBasis *fsE = BasisFactory::getNodalBasis(_type2tag[type]);
+ const polynomialBasis *pfsE = dynamic_cast<const polynomialBasis *>(fsE);
+ const int hoTag = ElementType::getTag(type, _viewOrder);
+ const nodalBasis *fsV = BasisFactory::getNodalBasis(hoTag);
+ const polynomialBasis *pfsV = dynamic_cast<const polynomialBasis *>(fsV);
+ viewData->setInterpolationMatrices(type,
+ pfsV->coefficients, pfsV->monomials,
+ pfsE->coefficients, pfsE->monomials);
+ }
+}
+
+void GMSH_AnalyseCurvedMeshPlugin::_createPViewPointwise()
+{
+ if (_pwJac) {
+ _setInterpolationMatrices(
+ new PView("Pointwise Jacobian", "ElementNodeData", _m, _dataPViewJac, 0, 1)
+ );
+ }
+
+ if (_pwIGE) {
+ _setInterpolationMatrices(
+ new PView("Pointwise IGE", "ElementNodeData", _m, _dataPViewIGE, 0, 1)
+ );
+ }
+
+ if (_pwICN) {
+ _setInterpolationMatrices(
+ new PView("Pointwise ICN", "ElementNodeData", _m, _dataPViewICN, 0, 1)
+ );
+ }
+}
+#endif
+
#endif
diff --git a/Plugin/AnalyseCurvedMesh.h b/Plugin/AnalyseCurvedMesh.h
index d1aac032a5306c5db3438f49739f82779a8e130d..cddb2796fcbfb1086125fd1c20b0588b447be9f9 100644
--- a/Plugin/AnalyseCurvedMesh.h
+++ b/Plugin/AnalyseCurvedMesh.h
@@ -22,11 +22,9 @@ private:
double _minJac, _maxJac, _minIGE, _minICN;
public:
data_elementMinMax(MElement *e,
- double minJ = 2,
- double maxJ = 0,
- double minIGE = -1,
- double minICN = -1)
- : _el(e), _minJac(minJ), _maxJac(maxJ), _minIGE(minIGE), _minICN(minICN) {}
+ double minJ = 2, double maxJ = 0,
+ double minIGE = -1, double minICN = -1)
+ : _el(e), _minJac(minJ), _maxJac(maxJ), _minIGE(minIGE), _minICN(minICN) {}
void setMinS(double r) { _minIGE = r; }
void setMinI(double r) { _minICN = r; }
MElement* element() { return _el; }
@@ -42,6 +40,17 @@ private :
GModel *_m;
double _threshold;
+#if defined(HAVE_VISUDEV)
+ // Pointwise data
+ int _numElementToScan;
+ bool _pwJac, _pwIGE, _pwICN;
+ std::map<int, std::vector<double>> _dataPViewJac;
+ std::map<int, std::vector<double>> _dataPViewIGE;
+ std::map<int, std::vector<double>> _dataPViewICN;
+ int _type2tag[20] = {0};
+ int _viewOrder = 0;
+#endif
+
// for 1d, 2d, 3d
bool _computedJac[3], _computedIGE[3], _computedICN[3];
bool _PViewJac[3], _PViewIGE[3], _PViewICN[3];
@@ -81,6 +90,12 @@ private :
void _printStatJacobian();
void _printStatIGE();
void _printStatICN();
+
+#if defined(HAVE_VISUDEV)
+ void _computePointwiseQuantities(MElement *, const fullMatrix<double> *normals);
+ void _createPViewPointwise();
+ void _setInterpolationMatrices(PView *);
+#endif
};
#endif
diff --git a/Post/PViewData.cpp b/Post/PViewData.cpp
index 59d56e5155b7a33ff0feac502e4407cc291d5402..0473d24ec2da431bc3d1266b4d80f29a5edf348e 100644
--- a/Post/PViewData.cpp
+++ b/Post/PViewData.cpp
@@ -184,6 +184,11 @@ bool PViewData::haveInterpolationMatrices(int type)
return _interpolation.count(type) ? true : false;
}
+void PViewData::deleteInterpolationMatrices(int type)
+{
+ _interpolation.erase(type);
+}
+
void PViewData::removeInterpolationScheme(const std::string &name)
{
std::map<std::string, interpolationMatrices>::iterator it = _interpolationSchemes.find(name);
diff --git a/Post/PViewData.h b/Post/PViewData.h
index 9b5449ec1227eb610707e2891c246d09dccdd6de..922083f18097416c4b61b40962f604d8b01bfd35 100644
--- a/Post/PViewData.h
+++ b/Post/PViewData.h
@@ -225,6 +225,7 @@ class PViewData {
const fullMatrix<double> &expGeo);
int getInterpolationMatrices(int type, std::vector<fullMatrix<double>*> &p);
bool haveInterpolationMatrices(int type=0);
+ void deleteInterpolationMatrices(int type=0);
// access to global interpolation schemes
static void removeInterpolationScheme(const std::string &name);