rend les elements positif

Plugin/AnalyseCurvedMesh.cpp

@@ -26,12 +26,12 @@
 
 StringXNumber JacobianOptions_Number[] = {
   {GMSH_FULLRC, "Dim", NULL, -1},
-  /*{GMSH_FULLRC, "Analysis", NULL, 1},*/
+  {GMSH_FULLRC, "Analysis", NULL, 1},
   {GMSH_FULLRC, "Effect", NULL, 6},
   {GMSH_FULLRC, "MaxDepth", NULL, 10},
   {GMSH_FULLRC, "JacBreak", NULL, .0},
   {GMSH_FULLRC, "BezBreak", NULL, .0},
-  /*{GMSH_FULLRC, "Tolerance", NULL, 1e-5}*/
+  {GMSH_FULLRC, "Tolerance", NULL, 1e-5}
 };
 
 extern "C"
@@ -55,15 +55,10 @@ std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const
 {
   return "Plugin(AnalyseCurvedMesh) check the jacobian of all elements of dimension 'Dim' or "
     "the greater model dimension if 'Dim' is either <0 or >3.\n\n "
-    //"Elements for which we are absolutely certain that they are good (positive jacobian) are ignored. "
-    //"Others are wrong or suppose to be wrong.\n\n"
-    //"Plugin(AnalyseCurvedMesh) write in the console which elements are wrong. "
-    //"(if labels of analysed type of elements are set visible, only wrong elements will be visible)\n\n"
-    /*"Analysis : 0 do nothing\n"
+    "Analysis : 0 do nothing\n"
     "          +1 find invalid elements (*)\n"
-    "          +2 compute min and max of all elements and print some statistics\n\n"*/
-    /*"Effect (for *) : 0 do nothing\n"*/
-    "Effect : 0 do nothing\n"
+    "          +2 compute J_min and J_max of all elements and print some statistics\n\n"
+    "Effect (for *) : 0 do nothing\n"
     "                +1 print a list of invalid elements\n"
     "                +2 print some statistics\n"
     "                +4 hide valid elements (for GUI)\n\n"
@@ -75,7 +70,7 @@ std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const
     "the element is said to be invalid\n\n"
     "BezBreak = [0,JacBreak[ : if all Bezier coefficients are > 'BezBreak', "
     "the element is said to be valid\n\n"
-    /*"Tolerance = R+ , << 1 : tolerance (relatively to J_avg.) used during the computation of min and max"*/;
+    "Tolerance = R+ , << 1 : tolerance (relatively to J_min and J_max) used during the computation of J_min and J_max";
 }
 
 // Miscellaneous method
@@ -221,7 +216,7 @@ static void setJacobian(MElement *el, const JacobianBasis *jfs, fullVector<doubl
 
 
 
-void setJacobian(MElement *const *el, const JacobianBasis *jfs, fullMatrix<double> &jacobian)
+static void setJacobian(MElement *const *el, const JacobianBasis *jfs, fullMatrix<double> &jacobian)
 {
   int numEl = jacobian.size2();
   int numVertices = el[0]->getNumVertices();
@@ -334,13 +329,12 @@ PView *GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
   _dim = (int) JacobianOptions_Number[0].def;
   if (_dim < 0 || _dim > 3)
     _dim = _m->getDim();
-  /*int analysis = (int) JacobianOptions_Number[1].def % 2;*/
-  int analysis = 1;
+  int analysis = (int) JacobianOptions_Number[1].def % 2;
   int toDo = (int) JacobianOptions_Number[1].def % 8;
   _maxDepth = (int) JacobianOptions_Number[2].def;
   _jacBreak = (double) JacobianOptions_Number[3].def;
   _bezBreak = (double) JacobianOptions_Number[4].def;
-  _tol = 1;/*(double) JacobianOptions_Number[5].def;*/
+  _tol = (double) JacobianOptions_Number[5].def;
   
   if (analysis % 2) {
     double t = Cpu();
@@ -349,7 +343,10 @@ PView *GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
     Msg::Info("Done validity check (%fs)", Cpu()-t);
   }
   if (analysis / 2) {
+    double t = Cpu();
+    Msg::Info("Strating computation J_min / J_max...");
     computeMinMax();
+    Msg::Info("Done computation J_min / J_max (%fs)", Cpu()-t);
   }
 }
   
@@ -518,6 +515,19 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(int toDo)
     
     case 2 :
       Msg::Warning("2D elements must be in a z=cst plane ! If they aren't, results won't be correct.");
+      
+				for (GModel::fiter it = _m->firstFace(); it != _m->lastFace(); it++) {
+					GFace *f = *it;
+					
+					unsigned int numType[3] = {0, 0, 0};
+					f->getNumMeshElements(numType);
+					
+					for (int type = 0; type < 3; type++) {
+						MElement *const *el = f->getStartElementType(type);
+						for (int jo = 0; jo < numType[type]; jo++)
+							el[jo]->setVolumePositive();
+					}
+				}		
       for (GModel::fiter it = _m->firstFace(); it != _m->lastFace(); it++) {
         GFace *f = *it;
         unsigned int numType[3] = {0, 0, 0};
@@ -648,73 +658,74 @@ void GMSH_AnalyseCurvedMeshPlugin::hideValid_ShowInvalid(std::vector<MElement*>
   invalids.pop_back();
 }
 
-void GMSH_AnalyseCurvedMeshPlugin::computeMinMax()
+void GMSH_AnalyseCurvedMeshPlugin::
+  checkValidity(MElement *const*el, int numEl, std::vector<MElement*> &invalids)
 {
-  _numAnalysedEl = 0;
-  _numInvalid = 0;
-  _numValid = 0;
-  _numUncertain = 0;
-  _min_Javg = .0, _max_Javg = .0;
-  _min_pJmin = .0, _avg_pJmin = .0;
-  _min_ratioJ = .0, _avg_ratioJ = .0;
+  if (numEl < 1)
+    return;
   
-  switch (_dim) {
-    case 3 :
-      for(GModel::riter it = _m->firstRegion(); it != _m->lastRegion(); it++) {
-        GRegion *r = *it;
-        unsigned int numType[5] = {0, 0, 0, 0, 0};
-        r->getNumMeshElements(numType);
-        
-        for(int type = 0; type < 5; type++) {
-          MElement *const *el = r->getStartElementType(type);
-          computeMinMax(el, numType[type]);
-          _numAnalysedEl += numType[type];
-        }
-      }
-      break;
+  const JacobianBasis *jfs = el[0]->getJacobianFuncSpace(-1);
+  if (!jfs) {
+    Msg::Error("Jacobian function space not implemented for type of element %d", el[0]->getNum());
+    return;
+  }
+  const bezierBasis *bb = jfs->bezier;
+  
+  int numSamplingPt = bb->points.size1();
+  int dim = bb->points.size2();
+  fullVector<double> jacobian(numSamplingPt);
+  fullVector<double> jacBez(numSamplingPt);
+  
+  for (int k = 0; k < numEl; ++k) {
+    setJacobian(el[k], jfs, jacobian);
     
-    case 2 :
-      Msg::Warning("2D elements must be in a z=cst plane ! If they aren't, results won't be correct.");
-      for (GModel::fiter it = _m->firstFace(); it != _m->lastFace(); it++) {
-        GFace *f = *it;
-        unsigned int numType[3] = {0, 0, 0};
-        f->getNumMeshElements(numType);
-        
-        for (int type = 0; type < 3; type++) {
-          MElement *const *el = f->getStartElementType(type);
-          computeMinMax(el, numType[type]);
-          _numAnalysedEl += numType[type];
-        }
-      }
-      break;
+    int i;
+    for (i = 0; i < numSamplingPt && jacobian(i) > _jacBreak; ++i);
+    if (i < numSamplingPt) {
+      invalids.push_back(el[k]);
+      ++_numInvalid;
+      continue;
+    }
     
-    case 1 :
-      Msg::Warning("1D elements must be on a y=cst & z=cst line ! If they aren't, results won't be correct.");
-      for (GModel::eiter it = _m->firstEdge(); it != _m->lastEdge(); it++) {
-        GEdge *e = *it;
-        unsigned int numElement = e->getNumMeshElements();
-        MElement *const *el = e->getStartElementType(0);
-        computeMinMax(el, numElement);
-        _numAnalysedEl += numElement;
-      }
-      break;
+    if (_maxDepth < 1) {
+      invalids.push_back(el[k]);
+      ++_numUncertain;
+      continue;
+    }
     
-    default :
-      Msg::Error("I can't analyse any element.");
-      return;
+    bb->matrixLag2Bez.mult(jacobian, jacBez);
+    
+    for (i = 0; i < jacBez.size() && jacBez(i) > _bezBreak; ++i);
+    if (i >= jacBez.size()) {
+      ++_numValid;
+      continue;
+    }
+    
+    if (_maxDepth < 2) {
+      invalids.push_back(el[k]);
+      ++_numUncertain;
+      continue;
+    }
+    else {
+      int result = subDivision(bb, jacBez, _maxDepth-1);
+      if (result < 0) {
+        invalids.push_back(el[k]);
+        ++_numInvalid;
+        continue;
+      }
+      if (result > 0) {
+        ++_numValid;
+        continue;
+      }
+      invalids.push_back(el[k]);
+      ++_numUncertain;
+      continue;
+    }
   }
-  Msg::Info("Extrema of J_avg : %g, %g", _min_Javg, _max_Javg);
-  Msg::Info("Minimum of min(distortion) : %g", _min_pJmin);
-  Msg::Info("Average of min(distortion) : %g", _avg_pJmin / _numAnalysedEl);
-  Msg::Info("Minimum of min(J) / max(J) : %g", _min_ratioJ);
-  Msg::Info("Average of min(J) / max(J) : %g", _avg_ratioJ / _numAnalysedEl);
 }
 
-
-
-
 void GMSH_AnalyseCurvedMeshPlugin::
-  checkValidity(MElement *const*el, int numEl, std::vector<MElement*> &invalids)
+  checkValidity_BLAS(MElement *const*el, int numEl, std::vector<MElement*> &invalids)
 {
   if (numEl < 1)
     return;
@@ -728,15 +739,14 @@ void GMSH_AnalyseCurvedMeshPlugin::
   
   int numSamplingPt = bb->points.size1();
   int dim = bb->points.size2();
-  fullVector<double> jacobian(numSamplingPt);
-  fullVector<double> jacBez(numSamplingPt);
+  fullMatrix<double> jacobian(numSamplingPt, numEl);
+  fullVector<double> jacBez(numSamplingPt), proxJac(numSamplingPt);
+  setJacobian(el, jfs, jacobian);
   
   for (int k = 0; k < numEl; ++k) {
-    setJacobian(el[k], jfs, jacobian);
-    
     int i;
-    for (i = 0; i < jacobian.size() && jacobian(i) > _jacBreak; ++i);
-    if (i < jacobian.size()) {
+    for (i = 0; i < numSamplingPt && jacobian(i, k) > _jacBreak; ++i);
+    if (i < numSamplingPt) {
       invalids.push_back(el[k]);
       ++_numInvalid;
       continue;
@@ -748,7 +758,8 @@ void GMSH_AnalyseCurvedMeshPlugin::
       continue;
     }
     
-    bb->matrixLag2Bez.mult(jacobian, jacBez);
+    proxJac.setAsProxy(jacobian, k);
+    bb->matrixLag2Bez.mult(proxJac, jacBez);
     
     for (i = 0; i < jacBez.size() && jacBez(i) > _bezBreak; ++i);
     if (i >= jacBez.size()) {
@@ -825,9 +836,71 @@ int GMSH_AnalyseCurvedMeshPlugin::
 }
 
 
+void GMSH_AnalyseCurvedMeshPlugin::computeMinMax()
+{/*
+  _numAnalysedEl = 0;
+  _numInvalid = 0;
+  _numValid = 0;
+  _numUncertain = 0;
+  _min_Javg = .0, _max_Javg = .0;
+  _min_pJmin = .0, _avg_pJmin = .0;
+  _min_ratioJ = .0, _avg_ratioJ = .0;
+  
+  switch (_dim) {
+    case 3 :
+      for(GModel::riter it = _m->firstRegion(); it != _m->lastRegion(); it++) {
+        GRegion *r = *it;
+        unsigned int numType[5] = {0, 0, 0, 0, 0};
+        r->getNumMeshElements(numType);
+        
+        for(int type = 0; type < 5; type++) {
+          MElement *const *el = r->getStartElementType(type);
+          computeMinMax(el, numType[type]);
+          _numAnalysedEl += numType[type];
+        }
+      }
+      break;
+    
+    case 2 :
+      Msg::Warning("2D elements must be in a z=cst plane ! If they aren't, results won't be correct.");
+      for (GModel::fiter it = _m->firstFace(); it != _m->lastFace(); it++) {
+        GFace *f = *it;
+        unsigned int numType[3] = {0, 0, 0};
+        f->getNumMeshElements(numType);
+        
+        for (int type = 0; type < 3; type++) {
+          MElement *const *el = f->getStartElementType(type);
+          computeMinMax(el, numType[type]);
+          _numAnalysedEl += numType[type];
+        }
+      }
+      break;
+    
+    case 1 :
+      Msg::Warning("1D elements must be on a y=cst & z=cst line ! If they aren't, results won't be correct.");
+      for (GModel::eiter it = _m->firstEdge(); it != _m->lastEdge(); it++) {
+        GEdge *e = *it;
+        unsigned int numElement = e->getNumMeshElements();
+        MElement *const *el = e->getStartElementType(0);
+        computeMinMax(el, numElement);
+        _numAnalysedEl += numElement;
+      }
+      break;
+    
+    default :
+      Msg::Error("I can't analyse any element.");
+      return;
+  }
+  Msg::Info("Extrema of J_avg : %g, %g", _min_Javg, _max_Javg);
+  Msg::Info("Minimum of min(~distortion) : %g", _min_pJmin);
+  Msg::Info("Average of min(~distortion) : %g", _avg_pJmin / _numAnalysedEl);
+  Msg::Info("Minimum of min(J) / max(J) : %g", _min_ratioJ);
+  Msg::Info("Average of min(J) / max(J) : %g", _avg_ratioJ / _numAnalysedEl);
+ */}
+
 void GMSH_AnalyseCurvedMeshPlugin::
   computeMinMax(MElement *const*el, int numEl)
-{
+{/*
   if (numEl < 1)
     return;
   
@@ -843,73 +916,124 @@ void GMSH_AnalyseCurvedMeshPlugin::
   
   fullMatrix<double> jacobian(numSamplingPt, numEl);
   fullMatrix<double> jacBez(numSamplingPt, numEl);
+  fullVector<double> proxBez(numSamplingPt);
   setJacobian(el, jfs, jacobian);
   bb->matrixLag2Bez.mult(jacobian, jacBez);
   
-  fullVector<double> prox(numSamplingPt);
-  
   for (int k = 0; k < numEl; ++k) {
-    prox.setAsProxy(jacBez, k);
+    proxBez.setAsProxy(jacBez, k);
+    double minJ, maxJ = minJ = jacobian(0,k);
     
-    double minJ, maxJ = minJ = jacobian(0,k);//, avgJac = 0;
     for (int i = 1; i < jacobian.size1(); ++i) {
       if (jacobian(i,k) < minJ) minJ = jacobian(i,k);
       if (jacobian(i,k) > maxJ) maxJ = jacobian(i,k);
-      //avgJac += prox(i);
     }
-    //avgJac /= jacobian.size1();
     
+    double minB, maxB = minB = jacBez(0,k), avgJ = .0;
     
+    for (int i = 1; i < proxBez.size1(); ++i) {
+      if (proxBez(i,k) < minB) minB = proxBez(i,k);
+      if (proxBez(i,k) > maxB) maxB = proxBez(i,k);
+      avgJ += jacBez(i,k);
+    }
+    avgJ /= proxBez.size1();
+    _min_Javg = std::min(_min_Javg, avgJ);
+    _max_Javg = std::max(_max_Javg, avgJ);
     
-    //int minJ, maxJ = minJ = prox(0);
-    //
-    //for (int i = 1; i < prox.size(); ++i) {
-    //  if (prox(i) > maxJ) maxJ = prox(i);
-    //  if (prox(i) < minJ) minJ = prox(i);
-    //}
-    
-    //for (int i = 0; i < jacobian.size(); i++)
-    //if (jacobian(i) <= _jacBreak) {
-    //  invalids.push_back(el[k]);
-    //  ++_numInvalid;
-    //  continue;
-    //}
-    //
-    //if (_maxDepth < 1) {
-    //  ++_numUncertain;
-    //  continue;
-    //}
-    //
-    //bb->matrixLag2Bez.mult(jacobian, jacBez);
-    //
-    //int i = 0;
-    //while (i < jacBez.size() && jacBez(i) > _bezBreak)
-    //  ++i;
-    //if (i >= jacBez.size()) {
-    //  ++_numValid;
-    //  continue;
-    //}
-    //
-    //if (_maxDepth < 2) {
-    //  ++_numUncertain;
-    //  continue;
-    //}
-    //else {
-    //  int result = subDivision(bb, jacBez, _maxDepth-1);
-    //  if (result < 0) {
-    //    invalids.push_back(el[k]);
-    //    ++_numInvalid;
-    //    continue;
-    //  }
-    //  if (result > 0) {
-    //    ++_numValid;
-    //    continue;
-    //  }
-    //  ++_numUncertain;
-    //  continue;
-    //}
+    if (_maxDepth > 1 &&
+        (minJ - minB > _tol * (std::abs(minJ) + std::abs(minB)) / 2 ||
+         maxB - maxJ > _tol * (std::abs(maxJ) + std::abs(maxB)) / 2   )) {
+      
+      
+      
+      
+      
+      
+      
+      // Set and map creation :
+      // The set is sorted in such a way that the first key should be the 
+      // better to partition.
+      // The map keep in memory the Bezier minimum and maximum for all part of the element.
+      BezierJacobian *bj = new BezierJacobian(proxBez, jfs, 0);
+      std::list<BezierJacobian> myset;
+      myset.insert(*bj);
+      std::pair<double,double> minmaxB(bj->minB(),bj->maxB());
+      std::map< int, std::pair<double,double> > map_minmaxB;
+      map_minmaxB.insert(std::pair< int, std::pair<double,double> >(bj->num(),minmaxB));
+      
+      bj->setMinMaxBounds(minLB, minUB, maxLB, maxUB);
+      if (min < minUB) minUB = min;
+      if (max > maxLB) maxLB = max;
+      
+      fullVector<double> newJacBez(jfs->divisor.size1());
+      fullVector<double> prox(numSamplingPt);
+      
+      int currentDepth, maxDepth = 0, k = 0;
+      
+      while ((minUB-minLB > tol || maxUB-maxLB > tol)) {
+        // Note that it is not optimized.
+        // When one interval reach tolerance, myset is certainly
+        // no more the more efficiently sorted 
+        // (the fault of BezierJacobian.operator<(...))
+        std::set<BezierJacobian>::iterator it = myset.begin();
+        (*it).partition(newJacBez, jfs);
+        currentDepth = (*it).depth() + 1;
+        if (maxDepth < currentDepth) maxDepth = currentDepth;
+        map_minmaxB.erase((*it).num());
+        myset.erase(it);
+        k++;
+        
+        for (int i = 0; i < jfs->numDivisions; i++) {
+          prox.setAsProxy(newJacBez, i * jacBez.size(), jacBez.size());
+          bj = new BezierJacobian(prox, jfs, currentDepth);
+          if (maxLB < bj->maxJ()) maxLB = bj->maxJ();
+          if (minUB > bj->minJ()) minUB = bj->minJ();
+          myset.insert(*bj);
+          minmaxB = std::make_pair(bj->minB(),bj->maxB());
+          map_minmaxB.insert(std::pair< int, std::pair<double,double> >(bj->num(),minmaxB));
+        }
+        
+        // Deletion of obsolet BezierJacobian
+        it = myset.begin();
+        if ((*it).isObsolet(maxLB, minUB, tol))
+          myset.erase(it++);
+        else ++it;
+        while (it != myset.end()) {
+          if ((*it).isObsolet(maxLB, minUB, tol))
+            myset.erase(it++);
+          else
+            ++it;
+        }
+        
+        // Update of Bezier bounds
+        minLB = minUB;
+        maxUB = maxLB;
+        std::map< int, std::pair<double,double> >::iterator mit;
+        for (mit = map_minmaxB.begin(); mit != map_minmaxB.end(); mit++) {
+          if (minLB > mit->second.first) minLB = mit->second.first;
+          if (maxUB < mit->second.second) maxUB = mit->second.second;
+        }
+      }
+      
+      double *a = new double[6];
+      a[0] = minLB, a[1] = minUB, a[2] = maxLB, a[3] = maxUB;
+      a[4] = (double)k; a[5] = (double)maxDepth;
+      return a;
+      
+      
+      
+      
+      
+      
+      
+      
+    }
+    _min_pJmin = std::min(_min_pJmin, (minJ+minB)/2);
+    _avg_pJmin += (minJ+minB)/2;
+    _min_ratioJ = std::min(_min_ratioJ, (minJ+minB)/(maxB+maxJ));
+    _avg_ratioJ += (minJ+minB)/(maxB+maxJ);
   }
-}
+ */}
 
 //int *GMSH_AnalyseCurvedMeshPlugin::checkJacobian
 //(MElement *const *el, int numEl, int maxDepth, int method)

Plugin/AnalyseCurvedMesh.h

@@ -28,6 +28,7 @@ class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
     double _min_Javg, _max_Javg;
     double _min_pJmin, _avg_pJmin;
     double _min_ratioJ, _avg_ratioJ;
+    //
     
   public :
     GMSH_AnalyseCurvedMeshPlugin(){}
@@ -42,6 +43,7 @@ class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
     StringXNumber *getOption(int);  
     PView *execute(PView *);
     void checkValidity(MElement *const *, int numEl, std::vector<MElement*> &invalids);
+    void checkValidity_BLAS(MElement *const *, int numEl, std::vector<MElement*> &invalids);
     
   private :
     void checkValidity(int toDo);
@@ -62,4 +64,37 @@ class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
     */
 };
 
+class BezierJacobian
+{
+private:
+  fullVector<double> _jacBez;
+  double _minJ, _maxJ, _minB, _maxB; //Extremum of Jac at corners and of bezier values
+  int _depthSub;
+  int _num; // Used for map of minmaxB
+  static int _globalNum;
+  
+public:
+  BezierJacobian(fullVector<double> &, const JacobianBasis *, int depth);
+  inline bool operator<(const BezierJacobian &other) const
+    {return other._maxB - _maxB - other._minB + _minB < 0;}
+  void partition(fullVector<double> &, const JacobianBasis *) const;
+  inline bool isObsolet(double maxLB, double minUB, double tol) const
+    {return _maxB - maxLB < tol &&  minUB - _minB < tol;}
+  
+  int num() const {return _num;}
+  int depth() const {return _depthSub;}
+  double minJ() const {return _minJ;}
+  double minB() const {return _minB;}
+  double maxJ() const {return _maxJ;}
+  double maxB() const {return _maxB;}
+  double setMinMaxBounds(double &minLB, double &minUB, double &maxLB, double &maxUB) const
+  {
+    /*Msg::Info("setminmax : %g = %g - %g",_minJ-_minB,_minJ,_minB);
+    Msg::Info("setminmax : %g = %g - %g",_maxB-_maxJ,_maxB,_maxJ);Msg::Info(" ");
+    
+    int g;
+    std::cin >> g;*/
+      minUB = _minJ; minLB = _minB; maxUB = _maxB; maxLB = _maxJ;}
+};
+
 #endif