allows to visualize pointwise Jacobian for a certain element

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

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
 

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

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