From 79b025728c613c5d1e2885a7a2d3cf99a3317236 Mon Sep 17 00:00:00 2001
From: Amaury Johnen <amaury.johnen@uclouvain.be>
Date: Fri, 11 Aug 2017 09:21:28 +0200
Subject: [PATCH] allows to visualize pointwise Jacobian for a certain element
---
Mesh/qualityMeasuresJacobian.cpp | 38 ++++++++----
Mesh/qualityMeasuresJacobian.h | 4 ++
Plugin/AnalyseCurvedMesh.cpp | 99 +++++++++++++++++++++++++++++++-
Plugin/AnalyseCurvedMesh.h | 7 +++
4 files changed, 136 insertions(+), 12 deletions(-)
diff --git a/Mesh/qualityMeasuresJacobian.cpp b/Mesh/qualityMeasuresJacobian.cpp
index 97eb1da16e..b882e54b5c 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;
}
@@ -443,17 +469,7 @@ void sampleIGEMeasure(MElement *el, int deg, double &min, double &max)
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 6c3b62992a..1bb4fc2a35 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 c798ab731f..8343055c84 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -19,6 +19,12 @@
#include <sstream>
#include <fstream>
#include "qualityMeasuresJacobian.h"
+#include "MLine.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "MHexahedron.h"
+#include "MTetrahedron.h"
+#include "MPrism.h"
class bezierBasis;
@@ -29,7 +35,9 @@ StringXNumber CurvedMeshOptions_Number[] = {
{GMSH_FULLRC, "Hiding threshold", NULL, 9},
{GMSH_FULLRC, "Draw PView", NULL, 0},
{GMSH_FULLRC, "Recompute", NULL, 0},
- {GMSH_FULLRC, "Dimension of elements", NULL, -1}
+ {GMSH_FULLRC, "Dimension of elements", NULL, -1},
+ //{GMSH_FULLRC, "Element to print quality", NULL, 2485}
+ {GMSH_FULLRC, "Element to print quality", NULL, 2555}
};
extern "C"
@@ -104,6 +112,11 @@ PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
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);
+ _numElementToScan = static_cast<int>(CurvedMeshOptions_Number[7].def);
+
+ _viewOrder = 10;
+ _elementToScan = NULL;
+ _hoElement = NULL;
if (askedDim < 0 || askedDim > 4) askedDim = _m->getDim();
@@ -163,6 +176,14 @@ PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
if (printStatS) _printStatIGE();
if (printStatI) _printStatICN();
+ if (_hoElement) {
+ std::map<int, std::vector<double>> dataPVelement;
+ dataPVelement[_hoElement->getNum()] = _jacElementToScan;
+ std::stringstream name;
+ name << "Jacobian elem " << _numElementToScan;
+ new PView(name.str().c_str(), "ElementNodeData", _m, dataPVelement, 0, 1);
+ }
+
// Create PView
if (drawPView)
for (int dim = 1; dim <= 3; ++dim) {
@@ -342,6 +363,33 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim)
_data.push_back(data_elementMinMax(el, min, max));
if (min < 0 && max < 0) ++cntInverted;
progress.next();
+
+ if (el->getNum() == _numElementToScan) {
+ _elementToScan = el;
+ fullMatrix<double> points =
+ _elementToScan->getFunctionSpace(_viewOrder)->points;
+ int tag = ElementType::getTag(_elementToScan->getType(), _viewOrder);
+ std::vector<MVertex *> v;
+ for (int k = 0; k < points.size1(); k++) {
+ double t1 = points(k, 0);
+ double t2 = points(k, 1);
+ double t3 = points(k, 2);
+ SPoint3 pos;
+ _elementToScan->pnt(t1, t2, t3, pos);
+ v.push_back(new MVertex(pos.x(), pos.y(), pos.z()));
+ }
+ MElementFactory factory;
+ _hoElement = factory.create(tag, v);
+
+ _addElementInEntity(_hoElement, entity);
+
+ fullVector<double> jac;
+ jacobianBasedQuality::sampleJacobian(_elementToScan, _viewOrder,
+ jac, normals);
+ for (int j = 0; j < jac.size(); ++j) {
+ _jacElementToScan.push_back(jac(j));
+ }
+ }
}
delete normals;
}
@@ -515,4 +563,53 @@ void GMSH_AnalyseCurvedMeshPlugin::_printStatICN()
infminI, avgminI, supminI);
}
+void GMSH_AnalyseCurvedMeshPlugin::_addElementInEntity(MElement *element,
+ GEntity *entity)
+{
+ int dim = entity->dim();
+ switch (dim) {
+ case 1:
+ {
+ GEdge *gedge = dynamic_cast<GEdge*>(entity);
+ int type = element->getType();
+ switch (type) {
+ case TYPE_LIN:
+ gedge->lines.push_back(dynamic_cast<MLine*>(element));
+ }
+ }
+ break;
+ case 2:
+ {
+ GFace *gface = dynamic_cast<GFace*>(entity);
+ int type = element->getType();
+ switch (type) {
+ case TYPE_TRI:
+ gface->triangles.push_back(dynamic_cast<MTriangle*>(element));
+ break;
+ case TYPE_QUA:
+ gface->quadrangles.push_back(dynamic_cast<MQuadrangle*>(element));
+ break;
+ }
+ }
+ break;
+ case 3:
+ {
+ GRegion *gregion = dynamic_cast<GRegion*>(entity);
+ int type = element->getType();
+ switch (type) {
+ case TYPE_TET:
+ gregion->tetrahedra.push_back(dynamic_cast<MTetrahedron*>(element));
+ break;
+ case TYPE_HEX:
+ gregion->hexahedra.push_back(dynamic_cast<MHexahedron*>(element));
+ break;
+ case TYPE_PRI:
+ gregion->prisms.push_back(dynamic_cast<MPrism*>(element));
+ break;
+ }
+ }
+ break;
+ }
+}
+
#endif
diff --git a/Plugin/AnalyseCurvedMesh.h b/Plugin/AnalyseCurvedMesh.h
index d1aac032a5..7722695b87 100644
--- a/Plugin/AnalyseCurvedMesh.h
+++ b/Plugin/AnalyseCurvedMesh.h
@@ -42,6 +42,12 @@ private :
GModel *_m;
double _threshold;
+ // Element to scan
+ int _numElementToScan;
+ MElement *_elementToScan, *_hoElement;
+ int _viewOrder;
+ std::vector<double> _jacElementToScan;
+
// for 1d, 2d, 3d
bool _computedJac[3], _computedIGE[3], _computedICN[3];
bool _PViewJac[3], _PViewIGE[3], _PViewICN[3];
@@ -81,6 +87,7 @@ private :
void _printStatJacobian();
void _printStatIGE();
void _printStatICN();
+ void _addElementInEntity(MElement*, GEntity*);
};
#endif
--
GitLab