Merge branch 'master' of https://gitlab.onelab.info/gmsh/gmsh

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?

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);

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

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

@@ -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

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

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);

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);