give interpolation matrices

Plugin/AnalyseCurvedMesh.cpp

@@ -25,6 +25,7 @@
 #include "MHexahedron.h"
 #include "MTetrahedron.h"
 #include "MPrism.h"
+#include "BasisFactory.h"
 
 class bezierBasis;
 
@@ -34,10 +35,10 @@ StringXNumber CurvedMeshOptions_Number[] = {
   {GMSH_FULLRC, "ICN measure", NULL, 1},
   {GMSH_FULLRC, "Hiding threshold", NULL, 9},
   {GMSH_FULLRC, "Draw PView", NULL, 0},
-  {GMSH_FULLRC, "Recompute", NULL, 0},
+  {GMSH_FULLRC, "Recompute", 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}
+  {GMSH_FULLRC, "Element to print quality", NULL, 2901}
 };
 
 extern "C"
@@ -181,7 +182,15 @@ PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
     dataPVelement[_hoElement->getNum()] = _jacElementToScan;
     std::stringstream name;
     name << "Jacobian elem " << _numElementToScan;
-    new PView(name.str().c_str(), "ElementNodeData", _m, dataPVelement, 0, 1);
+    PView *view = new PView(name.str().c_str(), "ElementNodeData",
+                            _m, dataPVelement, 0, 1);
+    const nodalBasis *fs = BasisFactory::getNodalBasis(_hoElement->getTypeForMSH());
+    const polynomialBasis *pfs = dynamic_cast<const polynomialBasis*>(fs);
+    PViewData *viewData = view->getData();
+    viewData->deleteInterpolationMatrices(_hoElement->getType());
+    viewData->setInterpolationMatrices(_hoElement->getType(),
+                                       pfs->coefficients, pfs->monomials,
+                                       pfs->coefficients, pfs->monomials);
   }
 
   // Create PView
@@ -364,32 +373,7 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim)
       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));
-        }
-      }
+      _computeJacobianToScan(el, normals);
     }
     delete normals;
   }
@@ -563,6 +547,39 @@ void GMSH_AnalyseCurvedMeshPlugin::_printStatICN()
             infminI, avgminI, supminI);
 }
 
+void GMSH_AnalyseCurvedMeshPlugin::_computeJacobianToScan(MElement *el,
+                                                          GEntity *entity,
+                                                          const fullMatrix<double> *normals)
+{
+  if (el->getNum() != _numElementToScan) return;
+
+  _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);
+  _jacElementToScan.clear();
+  for (int j = 0; j < jac.size(); ++j) {
+    _jacElementToScan.push_back(jac(j));
+  }
+}
+
 void GMSH_AnalyseCurvedMeshPlugin::_addElementInEntity(MElement *element,
                                                        GEntity *entity)
 {

Plugin/AnalyseCurvedMesh.h

@@ -87,6 +87,9 @@ private :
   void _printStatJacobian();
   void _printStatIGE();
   void _printStatICN();
+
+  void _computeJacobianToScan(MElement*, GEntity*,
+                              const fullMatrix<double> *normals);
   void _addElementInEntity(MElement*, GEntity*);
 };
 

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