No commit message

No commit message

Mesh/meshRecombine2D.cpp

@@ -326,9 +326,9 @@ int Recombine2D::apply()
   _gf->quadrangles = _quads;
   
   _applied = true;
-  _gf->model()->writeMSH("recSquare.msh");
-  laplaceSmoothing(_gf,100);
-  _gf->model()->writeMSH("aftSquare.msh");
+  //_gf->model()->writeMSH("recSquare.msh");
+  laplaceSmoothing(_gf,10);
+  //_gf->model()->writeMSH("aftSquare.msh");
   return 1;
 }
 

Mesh/meshRecombine2D.h

@@ -22,16 +22,184 @@ class RecombTriangle;
 class Recomb2D_Node;
 class TrianglesUnion;
 class Rec2d_vertex;
+class list_actions;
+class Recombine2D;
+//class list_actions;
+//class list_actions::iterator;
 struct lessRecombTri {
   bool operator()(RecombTriangle *rt1, RecombTriangle *rt2) const;
 };
 struct lessTriUnion {
   bool operator()(TrianglesUnion *, TrianglesUnion *) const;
 };
-
+struct lessTriUnionInv {
+  bool operator()(TrianglesUnion *, TrianglesUnion *) const;
+};
 typedef std::set<RecombTriangle*,lessRecombTri> Set_Recomb;
 typedef std::map<MElement*,std::set<RecombTriangle*> > Map_Tri_Recomb;
 typedef std::map<MElement*,Recomb2D_Node*> Map_Tri_Node;
+typedef std::map<MTriangle*, std::set<TrianglesUnion*> > mapTriUnion;
+typedef std::set<TrianglesUnion*, lessTriUnionInv> setTriUnion;
+typedef std::vector<setTriUnion> vectSetTriUnion;
+
+class list_actions {
+  private :
+    vectSetTriUnion _cont;
+    
+  public :
+    list_actions(int h) {_cont.reserve(h);}
+    
+    class iterator {
+      private :
+        int pos;
+        list_actions *la;
+        setTriUnion::iterator it;
+      
+      public :
+        void print() {
+          Msg::Info("iterator(pos %d, la %d, tu %d)", pos, la, *it);
+        }
+        iterator(list_actions *l, int i, setTriUnion::iterator t) {
+          pos = i;
+          la = l;
+          it = t;
+        }
+        iterator() {
+          pos = -1;
+        }
+        bool operator==(iterator a) const {
+          if (pos != a.pos)
+            return false;
+          if (pos < 0)
+            return true;
+          return it == a.it;
+        }
+        bool operator!=(iterator a) const {
+          if (pos != a.pos)
+            return true;
+          if (pos < 0)
+            return false;
+          return it != a.it;
+        }
+        iterator& operator=(const iterator &a) {
+          if (this != &a) {
+            pos = a.pos;
+            la = a.la;
+            it = a.it;
+          }
+          return *this;
+        }
+        TrianglesUnion* operator*() const {
+          //Msg::Info("pos %d", pos);
+          //Msg::Info(" ");
+          if (pos < 0)
+            return NULL;
+          return *it;
+        }
+        iterator& operator++() {
+          int p = pos;
+          ++it;
+          while (it == la->_end(p) && ++p != la->_size()) {
+            ++pos;
+            it = la->_begin(pos);
+            if (p != pos)
+              Msg::Error(" ");
+          }
+          return *this;
+        }
+        iterator operator++(int n) {
+          iterator t = *this;
+          int p = pos;
+          ++it;
+          while (it == la->_end(p) && ++p != la->_size()) {  
+            ++pos;
+            it = la->_begin(pos);
+          }
+          return t;
+        }
+        inline int getPos() {return pos;}
+    };
+    
+    int add(setTriUnion &s) {
+      setTriUnion::iterator it = s.begin();
+      while (it != s.end()) {
+        if (find(*it) != end())
+          s.erase(it++);
+        else
+          ++it;
+      }
+      _cont.push_back(s);
+      return _cont.size() - 1;
+    }
+    void remove(int a) {
+      for (; a < _cont.size(); a++)
+        _cont.pop_back();
+    }
+    iterator find(TrianglesUnion* tu) {
+      if (_cont.empty())
+        return end();
+      for (int i = 0; i < _cont.size(); ++i) {
+        setTriUnion::iterator it = _cont[i].find(tu);
+        if (it != _cont[i].end())
+          return iterator(this, i, it);
+      }
+      return end();
+    }
+    iterator begin() {
+      if (_cont.empty())
+        return iterator(this, -1, _begin());
+      return iterator(this, 0, _begin());
+    }
+    iterator end() {
+      return iterator(this, _cont.size() - 1, _end()); 
+    }
+    void pop_back() {
+      _cont.pop_back(); 
+    }
+    
+    void sizes(){
+      switch (_cont.size()) {
+        case 3 :
+          Msg::Info("--Actions : %d + %d + %d", _cont[0].size(), _cont[1].size(), _cont[2].size());
+          break;
+        case 2 :
+          Msg::Info("--Actions : %d + %d", _cont[0].size(), _cont[1].size());
+          break;
+        case 1 :
+          Msg::Info("--Actions : %d", _cont[0].size());
+          break;
+        default :
+          break;
+      }
+      //Msg::Info("%d sets", _cont.size());
+      //for (int i = 0; i < _cont.size(); ++i)
+      //  Msg::Info("[%d] %d", i+1, _cont[i].size());
+    }
+  
+  private :
+    setTriUnion::iterator _end(int pos = -1) {
+      if (_cont.size() == 0) {
+        setTriUnion::iterator it;
+        return it; // Not clean !
+      }
+      if (pos < 0)
+        return _cont.back().end();
+      return _cont[pos].end();
+    }
+    setTriUnion::iterator _begin(int pos = -1) {
+      if (_cont.size() == 0) {
+        setTriUnion::iterator it;
+        return it; // Not clean !
+      }
+      if (pos < 0)
+        return _cont.front().begin();
+      return _cont[pos].begin();
+    }
+    inline int _size() {return _cont.size();}
+    
+};
+
+
 
 class Recombine2D {
   private :
@@ -43,6 +211,7 @@ class Recombine2D {
     Map_Tri_Recomb _possibleRecomb;
     std::set<MElement*> _isolated;
     std::vector<MQuadrangle*> _quads;
+    static int ra;
     
     template <class E>
     void _buildEdgeToElement(std::vector<E*> &, e2t_cont &);
@@ -57,15 +226,27 @@ class Recombine2D {
     ~Recombine2D();
     
     int apply();
-    int apply(bool);
     inline double getBenef() const {return _benef;}
     inline int numTriangle() const {return _isolated.size();}
     
   private :
     //std::set<TrianglesUnion*, lessTriUnion> _setQuads;
     std::list<TrianglesUnion*> _setQuads;
-    void _removeImpossible(std::list<TrianglesUnion*>::iterator);
+    std::list<TrianglesUnion*> _lastQuad;
+    mapTriUnion _possibleActions;
+    mapTriUnion _lessActions;
+    void _removeImpossible(TrianglesUnion*);
     void _recombine(bool);
+    
+    void _lookAhead(std::list<TrianglesUnion*>&, int horiz);
+    void _getIncompatibles(const TrianglesUnion*, setTriUnion&);
+    void _getNeighbors(const TrianglesUnion*, setTriUnion&);
+    void _addNeighbors(int horiz, std::vector<list_actions::iterator>&, int current, list_actions*);
+    void _removeNeighbors(int horiz, int current, list_actions*);
+    int _checkIsolatedTri(std::vector<list_actions::iterator>&, int size, std::set<MTriangle*>&);
+
+  public :
+    int apply(bool);
 };
 
 class RecombTriangle {
@@ -135,9 +316,6 @@ class Rec2d_vertex {
     // GEntity : virtual Range<double> parBounds(int i)
     static double **_Vvalues;
     static double **_Vbenefs;
-    
-    void _initStaticTable();
-    double _computeAngle(MVertex *, std::list<MTriangle*> &, std::set<GEdge*> &);
   
   public :
     Rec2d_vertex(MVertex *, std::list<MTriangle*> &,
@@ -146,6 +324,12 @@ class Rec2d_vertex {
     inline void changeNumEl(int c) {_numEl += c;}
     double getReward();
     double getReward(int);
+    //int onwhat() {return _onWhat;}
+    //int numel() {return _numEl;}
+    
+  private :
+    void _initStaticTable();
+    double _computeAngle(MVertex *, std::list<MTriangle*> &, std::set<GEdge*> &);
 };
 
 class TrianglesUnion {
@@ -155,24 +339,31 @@ class TrianglesUnion {
     Rec2d_vertex **_vertices;
     MTriangle **_triangles;
     static int _RECOMPUTE; //incremented when a recombination is executed
+    
     public:
     static int _NumEdge, _NumVert;
     static double _ValEdge, _ValVert;
     private:
     
-    double _computeEdgeLength(GFace *, MVertex **, double *u, double *v,
-                              int numIntermedPoint= 1);
-    double _computeAlignment(MEdge &, std::list<MTriangle*> &);
-    void _update();
-    
   public :
     TrianglesUnion(GFace *, std::list<MTriangle*> &, std::list<MEdge> &,
                    std::list<Rec2d_vertex*> &, std::map<MVertex*,double> &);
     
     bool operator<(TrianglesUnion &other);
-    inline double getEdgeValue() {return _embEdgeValue;}
-    inline double getNumTriangles() {return _numTri;}
-    inline MTriangle *getTriangle(int num) {return _triangles[num];}
+    void addTri(std::set<MTriangle*>&) const;
+    void removeTri(std::set<MTriangle*>&) const;
+    bool isIn(std::set<MTriangle*>&) const;
+    void addInfluence(int*, double*, std::map<Rec2d_vertex*, int>&) const;
+    void removeInfluence(int*, double*, std::map<Rec2d_vertex*, int>&) const;
+    static double computeReturn(const int*, const double*, const std::map<Rec2d_vertex*, int>&);
+    inline double getEdgeValue() const {return _embEdgeValue;}
+    inline double getVertValue() const {return _embVertValue;}
+    inline int getNumVerts() const {return _numEmbVert;}
+    inline int getNumEdges() const {return _numEmbEdge;}
+    inline int getNumVertices() const {return _numBoundVert;}
+    inline Rec2d_vertex* getVertex(int num) const {return _vertices[num];}
+    inline int getNumTriangles() const {return _numTri;}
+    inline MTriangle* getTriangle(int num) const {return _triangles[num];}
     void select() {
       _ValEdge -= _embEdgeValue;
       _NumEdge -= _numEmbEdge;
@@ -184,7 +375,7 @@ class TrianglesUnion {
       }
       _RECOMPUTE++;
     }
-    MQuadrangle *createQuad() const;
+    MQuadrangle* createQuad() const;
     void print() {
       Msg::Info("Printing TU (%d,%d,%d,%d)", _numEmbVert, _numBoundVert, _numEmbEdge, _numTri);
       for (int i = 0; i < _numTri; i++)
@@ -192,5 +383,11 @@ class TrianglesUnion {
     }
     static inline void addRec() {_RECOMPUTE++;}
     inline double getReturn() {return _globValIfExecuted;}
+  
+  private:  
+    double _computeEdgeLength(GFace *, MVertex **, double *u, double *v,
+                              int numIntermedPoint= 1);
+    double _computeAlignment(MEdge &, std::list<MTriangle*> &);
+    void _update();
 };
 #endif

Mesh/meshRecombine2D_2.cpp

@@ -17,8 +17,8 @@
 #include "Context.h"
 #include "OpenFile.h"//pas propre
 #include "Field.h"
-
-static int HORIZ = 20;
+#include "OS.h"
+#define HORIZ2 7
 
 double **Rec2d_vertex::_Vvalues = NULL;
 double **Rec2d_vertex::_Vbenefs = NULL;
@@ -27,6 +27,16 @@ int TrianglesUnion::_NumEdge = 0;
 int TrianglesUnion::_NumVert = 0;
 double TrianglesUnion::_ValEdge = .0;
 double TrianglesUnion::_ValVert = .0;
+int Recombine2D::ra = 0;
+
+/*
+  a faire :
+  - modifier gain des edges
+  - fonction retournant qualite
+  - verifier la qualite
+  - action supplementaire (faire classe de base)
+  - Faire un super conteneur pour les actions
+*/
 
 Recombine2D::Recombine2D(GFace *gf)
 : _horizon(0), _gf(gf), _benef(.0), _applied(true)
@@ -69,10 +79,10 @@ Recombine2D::Recombine2D(GFace *gf)
   std::map<MVertex*, Rec2d_vertex*>::iterator itV2N = mapV2N.begin();
   std::list<std::map<MVertex*, std::list<MTriangle*> >::iterator> fourTri;
   for (; itvert != mapVert.end(); itvert++) {
-    if ((*itvert).second.size() == 4 && (*itvert).first->onWhat()->dim() == 2)
+    if (itvert->second.size() == 4 && itvert->first->onWhat()->dim() == 2)
       fourTri.push_back(itvert);
-    Rec2d_vertex *n = new Rec2d_vertex((*itvert).first, (*itvert).second, mapCornerVert);
-    itV2N = mapV2N.insert(itV2N, std::make_pair((*itvert).first,n));
+    Rec2d_vertex *n = new Rec2d_vertex(itvert->first, itvert->second, mapCornerVert);
+    itV2N = mapV2N.insert(itV2N, std::make_pair(itvert->first,n));
     TrianglesUnion::_NumVert++;
     TrianglesUnion::_ValVert += n->getReward();
   }
@@ -98,13 +108,13 @@ Recombine2D::Recombine2D(GFace *gf)
     {
       std::set<MVertex*> setVert;
       std::multiset<MEdge, Less_Edge> msetEdge;
-      std::list<MTriangle*>::iterator itTri = (*(*it4)).second.begin();
-      for (; itTri != (*(*it4)).second.end(); itTri++) {
+      std::list<MTriangle*>::iterator itTri = (*it4)->second.begin();
+      for (; itTri != (*it4)->second.end(); itTri++) {
         MVertex *vt;
         for (int i = 0; i < 3; i++) {
           msetEdge.insert((*itTri)->getEdge(i));
           vt = (*itTri)->getVertex(i);
-          if (vt != (*(*it4)).first)
+          if (vt != (*it4)->first)
             setVert.insert(vt);
         }
       }
@@ -118,10 +128,14 @@ Recombine2D::Recombine2D(GFace *gf)
       std::set<MVertex*>::iterator itsVert = setVert.begin();
       for (; itsVert != setVert.end(); itsVert++)
         listVert.push_back(mapV2N[*itsVert]);
-      listVert.push_back(mapV2N[(*(*it4)).first]);
+      listVert.push_back(mapV2N[(*it4)->first]);
     }
     Msg::Info("%d",++j);
-    _setQuads.push_back(new TrianglesUnion(gf, (*(*it4)).second, listEmbEdge, listVert, mapV2LcUV));
+    TrianglesUnion *tu = new TrianglesUnion(gf, (*it4)->second, listEmbEdge, listVert, mapV2LcUV);
+    std::list<MTriangle*>::iterator itTri = (*it4)->second.begin();
+    for (; itTri != (*it4)->second.end(); itTri++)
+      _possibleActions[*itTri].insert(tu);
+    _setQuads.push_back(tu);
   }
   
   // Create 'TrianglesUnion' for pattern of 2 triangles
@@ -134,22 +148,25 @@ Recombine2D::Recombine2D(GFace *gf)
 */
   std::map<MEdge, std::list<MTriangle*> >::iterator itedge = mapEdge.begin();
   for (; itedge != mapEdge.end(); itedge++) {
-    if ((*itedge).second.size() == 2) {
+    if (itedge->second.size() == 2) {
       std::list<MEdge> listEmbEdge;
-      listEmbEdge.push_back((*itedge).first);
+      listEmbEdge.push_back(itedge->first);
       std::list<Rec2d_vertex*> listVert;
-      listVert.push_back(mapV2N[(*itedge).first.getVertex(0)]);
-      listVert.push_back(mapV2N[(*itedge).first.getVertex(1)]);
-      TrianglesUnion *tu = new TrianglesUnion(gf, (*itedge).second, listEmbEdge, listVert, mapV2LcUV);
+      listVert.push_back(mapV2N[itedge->first.getVertex(0)]);
+      listVert.push_back(mapV2N[itedge->first.getVertex(1)]);
+      TrianglesUnion *tu = new TrianglesUnion(gf, itedge->second, listEmbEdge, listVert, mapV2LcUV);
       _setQuads.push_back(tu);
       TrianglesUnion::_NumEdge++;
       TrianglesUnion::_ValEdge += tu->getEdgeValue();
+      std::list<MTriangle*>::iterator itTri = itedge->second.begin();
+      for (; itTri != itedge->second.end(); itTri++)
+        _possibleActions[*itTri].insert(tu);
     }
     else
-      ;//Msg::Info("[bord] %d", (*itedge).second.size());
+      ;//Msg::Info("[bord] %d", itedge->second.size());
   }
   TrianglesUnion::addRec();
-  _setQuads.sort(lessTriUnion());
+  //_setQuads.sort(lessTriUnion());
   
   _recombine(true);
   _applied = false;
@@ -161,31 +178,324 @@ void Recombine2D::_recombine(bool a)
   SetBoundingBox();
   
   int i = 0;
-  while (!_setQuads.empty() && i < 1000) {
-    std::list<TrianglesUnion*>::iterator it = _setQuads.begin();
-        Msg::Info("------------------ %d", _setQuads.size());
-    {
-      Msg::Info("Expected return %lf", (*it)->getReturn());
+  std::list<TrianglesUnion*> vectTu;
+  
+  while (!_setQuads.empty() && i < 2000) {
+    TrianglesUnion *tu;
+    if (_lastQuad.empty()) {
+      if (vectTu.size() < 2)
+        tu = *std::max_element(_setQuads.begin(), _setQuads.end(), lessTriUnion());
+      else {
+        vectTu.pop_front();
+        //tu = *std::max_element(vectTu.begin(),vectTu.end(), lessTriUnion());
+        tu = *vectTu.begin();
+      }
+      vectTu.clear();
+      vectTu.insert(vectTu.begin(), tu);
+      
+            Msg::Info("------------------ %d", _setQuads.size());
+      _lookAhead(vectTu, HORIZ2);
+    }
+    else {
+      tu = *_lastQuad.begin();
+      vectTu.clear();
+      vectTu.insert(vectTu.begin(), tu);
+      std::set<MTriangle*> touched;
+      for (int i = 0; i < tu->getNumTriangles(); i++) {
+        touched.insert(tu->getTriangle(i));
+      }
+      std::list<TrianglesUnion*>::iterator it = _lastQuad.begin();
+      while (it != _lastQuad.end()) {
+        bool toBeErased = false;
+        for (int i = 0; i < (*it)->getNumTriangles(); i++) {
+          if (touched.find((*it)->getTriangle(i)) != touched.end())
+            toBeErased = true;
+        }
+        if (toBeErased)
+          _lastQuad.erase(it++);
+        else
+          ++it;
+      }
+      //_lastQuad.erase(_lastQuad.begin());
+    }
+    tu = *vectTu.begin();
+    tu->select();
+    _quads.push_back(tu->createQuad());
+    _removeImpossible(tu);
+    { // draw Mesh
+      Msg::Info("Expected return %lf", tu->getReturn());
       _applied = false;
-      apply();
+      apply(true);
       _applied = false;
       CTX::instance()->mesh.changed = (ENT_ALL);
       //drawContext::global()->draw();
-      //FlGui::instance()->check();
+      FlGui::instance()->check();
       drawContext::global()->draw();
-      //if (!Msg::GetAnswer("Continue ?", 1, "no", "yes"))
-      //  Msg::Info("I continue anyway :p");
-    }   
-    (*it)->select();
+      double t = Cpu();
+      while (Cpu() - t < .0001);
+      if (i % 1000 == 0)
+      if (!Msg::GetAnswer("Continue ?", 1, "no", "yes"))
+        Msg::Info("I continue anyway :p");
+    }
+    
+    ++i;
+  }
+  Msg::Info("Done %d (%d)", i, _pairs.size());
+}
+
+void printIt(std::vector<list_actions::iterator> &vIt, int size)
+{
+  switch (size) {
+    case 3 :
+      Msg::Info("--Iterators : %d , %d , %d", vIt[0].getPos(), vIt[1].getPos(), vIt[2].getPos());
+      break;
+    case 2 :
+      Msg::Info("--Iterators : %d , %d", vIt[0].getPos(), vIt[1].getPos());
+      break;
+    case 1 :
+      Msg::Info("--Iterators : %d", vIt[0].getPos());
+      break;
+    default :
+      break;
+  }
+}
+
+void Recombine2D::_lookAhead(std::list<TrianglesUnion*> &candidate, int horiz)
+{
+  if (horiz < 1)
+    horiz = 1;
+  _lessActions.clear();
+  
+  double maxReturn = .0;
+  //TrianglesUnion *bestAction;
+  int numbers[2] = {0, 0}; // number of embedded verts & edges
+  double values[2] = {.0, .0}; // embedded vert & edge values
+  std::set<MTriangle*> setTri;
+  std::map<Rec2d_vertex*, int> modifiedVert;
+  std::vector<list_actions::iterator> vecIt(horiz);
+  std::vector<std::pair<TrianglesUnion*, int> > bestSequence;
+  
+  list_actions la(horiz);
+  {
+    setTriUnion candidates;
+    _getIncompatibles(*candidate.begin(), candidates);
+    la.add(candidates);
+  }
+  
+  vecIt[0] = la.begin();
+  (*vecIt[0])->addTri(setTri);
+  (*vecIt[0])->addInfluence(numbers, values, modifiedVert);
+  _addNeighbors(horiz, vecIt, 0, &la);
+  
+  int current = 1;
+  int num=0;
+  
+  while (current > 0) {
+    bool best = false;
+    // set first acceptable action for each step
+    while (current < horiz && vecIt[current-1] != la.end()) {
+      vecIt[current] = vecIt[current-1];
+      while (++vecIt[current] != la.end() && (*vecIt[current])->isIn(setTri));
+      if (vecIt[current] != la.end()) {
+        (*vecIt[current])->addTri(setTri);
+        (*vecIt[current])->addInfluence(numbers, values, modifiedVert);
+        _addNeighbors(horiz, vecIt, current, &la);
+      }
+      ++current;
+    }
+
+    { // color
+      for (int i = 1; i < current &&  vecIt[i] != la.end(); ++i) {
+        TrianglesUnion *tu = *vecIt[i];
+        int pos = vecIt[i].getPos();
+        double frac = (double)pos / (double)(horiz-1);
+        unsigned int col;
+        //la.sizes();
+        //Msg::Info("%d / %d -> %d", pos, horiz, (int)(200.*frac));
+        col = CTX::instance()->packColor(50 - (int)(25.*frac), 200 - (int)(25.*frac), (int)(200.*frac), 255);
+        for (int j = 0; j < tu->getNumTriangles(); ++j) {
+          tu->getTriangle(j)->setCol(col);
+        }
+      }
+    }
+    
+    double expectedReturn = TrianglesUnion::computeReturn(numbers, values, modifiedVert);
+    if (maxReturn < expectedReturn) {
+      int sizeSequence = current;
+      if (vecIt[sizeSequence - 1] == la.end())
+        --sizeSequence;
+      expectedReturn -= .1 * _checkIsolatedTri(vecIt, sizeSequence, setTri);
+      if (maxReturn < expectedReturn) {
+        maxReturn = expectedReturn;
+        //bestAction = *vecIt[0];
+        bestSequence.resize(sizeSequence);
+        for (int i = 0; i < sizeSequence; ++i)
+          bestSequence[i] = std::make_pair(*vecIt[i], vecIt[i].getPos());
+        { //color
+          best = true;
+          TrianglesUnion *tu = *vecIt[0];
+          unsigned int col;
+          col = CTX::instance()->packColor(190, 0, 190, 255);
+          for (int j = 0; j < tu->getNumTriangles(); ++j) {
+            tu->getTriangle(j)->setCol(col);
+          }
+        }
+      }
+    }
+    else { //color
+      TrianglesUnion *tu = *vecIt[0];
+      unsigned int col;
+      col = CTX::instance()->packColor(190, 0, 0, 255);
+      for (int j = 0; j < tu->getNumTriangles(); ++j) {
+        tu->getTriangle(j)->setCol(col);
+      }
+    }
     
-    _quads.push_back((*it)->createQuad());
+    if (true || best) { // draw mesh
+      apply(false);
+      CTX::instance()->mesh.changed = (ENT_ALL);
+      FlGui::instance()->check();
+      drawContext::global()->draw();
+      double t = Cpu();
+      while (Cpu() - t < .0001);
+      //if (i % 1 == 0)
+      if (/*false && best &&*/ !Msg::GetAnswer("Continue ?", 1, "no", "yes"))
+        Msg::Info("I continue anyway :p");
+      best = false;
+    }
     
-    _removeImpossible(it);
+    { // decolor
+      for (int i = 0; i < current &&  vecIt[i] != la.end(); ++i) {
+        TrianglesUnion *tu = *vecIt[i];
+        int pos = vecIt[i].getPos();
+        double frac = (double)pos / (double)horiz;
+        unsigned int col;
+        col = CTX::instance()->packColor(255, 255, 0, 255);
+        for (int j = 0; j < tu->getNumTriangles(); ++j) {
+          tu->getTriangle(j)->setCol(col);
+        }
+      }
+    }
     
-    _setQuads.sort(lessTriUnion());
-    i++;
+    if (vecIt[current-1] != la.end() || --current > 0) {
+      do {
+        (*vecIt[current-1])->removeTri(setTri);
+        (*vecIt[current-1])->removeInfluence(numbers, values, modifiedVert);
+        _removeNeighbors(horiz, current - 1, &la);
+        while (++vecIt[current-1] != la.end() && (*vecIt[current-1])->isIn(setTri));
+      } while (vecIt[current-1] == la.end() && --current > 0);
+      if (current > 0) {
+        (*vecIt[current-1])->addTri(setTri);
+        (*vecIt[current-1])->addInfluence(numbers, values, modifiedVert);
+        _addNeighbors(horiz, vecIt, current - 1, &la);
+      }
+    }
   }
-  Msg::Info("Done %d (%d)", i, _pairs.size());
+  
+  { // color
+    for (int i = 1; i < bestSequence.size(); ++i) {
+      TrianglesUnion *tu = bestSequence[i].first;
+      unsigned int col = CTX::instance()->packColor(50, 200, 0, 255);
+      for (int j = 0; j < tu->getNumTriangles(); ++j) {
+        tu->getTriangle(j)->setCol(col);
+      }
+    }
+  }
+
+  candidate.clear();
+  for (int i = 0; i < bestSequence.size(); ++i) {
+    candidate.insert(candidate.end(), bestSequence[i].first);
+  }
+}
+
+void Recombine2D::_getIncompatibles(const TrianglesUnion *tu, setTriUnion &set)
+{
+  for (int i = 0; i < tu->getNumTriangles(); i++) {
+    mapTriUnion::iterator it = _lessActions.find(tu->getTriangle(i));
+    if (it == _lessActions.end()) {
+      it = _possibleActions.find(tu->getTriangle(i));
+      if (it == _possibleActions.end()) {
+        Msg::Error("[Rec2D] error no triangle !, stopping with partially filled set");
+      }
+      _lessActions.insert(*it);
+    }
+    setTriUnion::iterator it2 = it->second.begin();
+    for (; it2 != it->second.end(); ++it2) {
+      set.insert(*it2);
+    }
+  }
+}
+
+void Recombine2D::_getNeighbors(const TrianglesUnion *tu, setTriUnion &set)
+{
+  setTriUnion incomp;
+  setTriUnion neighbors;
+  _getIncompatibles(tu, incomp);
+  setTriUnion::iterator it = incomp.begin();
+  for (; it != incomp.end(); ++it)
+    _getIncompatibles(*it, neighbors);
+  it = incomp.begin();
+  for (; it != incomp.end(); ++it)
+    neighbors.erase(*it);
+  it = neighbors.begin();
+  for (; it != neighbors.end(); ++it)
+    set.insert(*it);
+}
+
+void Recombine2D::_addNeighbors(int horiz, std::vector<list_actions::iterator> &vecIt, int current, list_actions *la)
+{
+  if (current < horiz - 1 && current > -1) {
+    //Msg::Info("+ra %d", ++ra);
+    setTriUnion neighbors;
+    _getNeighbors(*vecIt[current], neighbors);
+    la->add(neighbors);
+  }
+}
+
+void Recombine2D::_removeNeighbors(int horiz, int current, list_actions *la)
+{
+  if (current < horiz - 1 && current > -1) {
+    //Msg::Info("-ra %d", --ra);
+    la->pop_back();
+  }
+}
+
+
+int Recombine2D::_checkIsolatedTri(std::vector<list_actions::iterator> &vecIt, int size, std::set<MTriangle*> &seqTri)
+{
+  setTriUnion setTu;
+  for (int i = 0; i < size; ++i) {
+    _getIncompatibles(*vecIt[i], setTu);
+  }
+  
+  std::set<MTriangle*> setBoundaryTri;
+  setTriUnion::iterator it = setTu.begin();
+  for (; it != setTu.end(); ++it) {
+    for (int i = 0; i < (*it)->getNumTriangles(); i++)
+    if (seqTri.find((*it)->getTriangle(i)) == seqTri.end())
+      setBoundaryTri.insert((*it)->getTriangle(i));
+  }
+  
+  int num = 0;
+  std::set<MTriangle*>::iterator itTri = setBoundaryTri.begin();
+  for (; itTri != setBoundaryTri.end(); ++itTri) {
+    mapTriUnion::iterator it1 = _lessActions.find(*itTri);
+    if (it1 == _lessActions.end()) {
+      it1 = _possibleActions.find(*itTri);
+      if (it1 == _possibleActions.end()) {
+        Msg::Error("[Rec2D] error no triangle !, stopping with partially filled set");
+      }
+      _lessActions.insert(*it1);
+    }
+    int numPossibleRecombinations = 0;
+    setTriUnion::iterator it2 = it1->second.begin();
+    for (; it2 != it1->second.end(); ++it2)
+    if (setTu.find(*it2) == setTu.end())
+      ++numPossibleRecombinations;
+    if (!numPossibleRecombinations)
+      ++num;
+  }
+  return num;
 }
 
 MQuadrangle *TrianglesUnion::createQuad() const
@@ -222,24 +532,24 @@ MQuadrangle *TrianglesUnion::createQuad() const
   v2 = listBoundEdge.begin()->getMaxVertex();
   std::list<MEdge>::iterator it = listBoundEdge.begin();
   for (it++; it != listBoundEdge.end(); it++) {
-    if (v2 == (*it).getMinVertex()) {
-      v3 = (*it).getMaxVertex();
+    if (v2 == it->getMinVertex()) {
+      v3 = it->getMaxVertex();
       listBoundEdge.erase(it);
       break;
     }
-    if (v2 == (*it).getMaxVertex()) {
-      v3 = (*it).getMinVertex();
+    if (v2 == it->getMaxVertex()) {
+      v3 = it->getMinVertex();
       listBoundEdge.erase(it);
       break;
     }
   }
   for (it = listBoundEdge.begin(); it != listBoundEdge.end(); it++) {
-    if (v3 == (*it).getMinVertex()) {
-      v4 = (*it).getMaxVertex();
+    if (v3 == it->getMinVertex()) {
+      v4 = it->getMaxVertex();
       break;
     }
-    if (v3 == (*it).getMaxVertex()) {
-      v4 = (*it).getMinVertex();
+    if (v3 == it->getMaxVertex()) {
+      v4 = it->getMinVertex();
       break;
     }
   }
@@ -250,50 +560,96 @@ MQuadrangle *TrianglesUnion::createQuad() const
 
 int Recombine2D::apply(bool a)
 {
-  if (!_haveParam || _applied) return 0;
-  
-  std::vector<MTriangle*> triangles2;
-  for (int i = 0; i < _gf->triangles.size(); i++) {
-    if (_isolated.find(_gf->triangles[i]) != _isolated.end())
-      delete _gf->triangles[i];
-    else
-      triangles2.push_back(_gf->triangles[i]);
+  //Msg::Info("(%d, %d, %d)", _quads.size(), _gf->triangles.size(), _isolated.size());
+  if (a) {
+    std::vector<MTriangle*> triangles2;
+    for (int i = 0; i < _gf->triangles.size(); i++) {
+      if (_isolated.find(_gf->triangles[i]) != _isolated.end())
+        delete _gf->triangles[i];
+      else
+        triangles2.push_back(_gf->triangles[i]);
+    }
+    _gf->triangles = triangles2;
+    _gf->quadrangles = _quads;
+    _isolated.clear();
   }
-  _gf->triangles = triangles2;
-  _gf->quadrangles = _quads;
-  _isolated.clear();
   
   _applied = true;
   return 1;
 }
 
-void Recombine2D::_removeImpossible(std::list<TrianglesUnion*>::iterator it)
+void Recombine2D::_removeImpossible(TrianglesUnion *tu)
 {
-  std::set<MTriangle *> touched;
+  for (int i = 0; i < tu->getNumTriangles(); i++) {
+    _isolated.insert(tu->getTriangle(i));
+  }// for test
+  
+  std::set<MTriangle*> touched;
+  for (int i = 0; i < tu->getNumTriangles(); i++) {
+    touched.insert(tu->getTriangle(i));
+  }
+  
+  setTriUnion incomp;
+  std::set<MTriangle*>::iterator itTri = touched.begin();
+  for (; itTri != touched.end(); ++itTri) {
+    mapTriUnion::iterator it = _possibleActions.find(*itTri);
+    if (it != _possibleActions.end()) {
+      setTriUnion::iterator it2 = it->second.begin();
+      for (; it2 != it->second.end(); ++it2)
+        incomp.insert(*it2);
+    }
+  }
+  
+  setTriUnion::iterator itTu = incomp.begin();
+  for (; itTu != incomp.end(); ++itTu) {
+    for (int i = 0; i < (*itTu)->getNumTriangles(); i++) {
+      mapTriUnion::iterator it4 = _possibleActions.find((*itTu)->getTriangle(i));
+      it4->second.erase(*itTu);
+      TrianglesUnion *tu;
+      if (touched.find(it4->first) == touched.end())
+      switch (it4->second.size()) {
+        case 1 :
+          _lastQuad.insert(_lastQuad.begin(), *it4->second.begin());
+          tu = *it4->second.begin();
+          unsigned int col;
+          col = CTX::instance()->packColor(0, 0, 0, 0);
+          for (int j = 0; j < tu->getNumTriangles(); ++j)
+            tu->getTriangle(j)->setCol(col);
+          break;
+        case 0 :
+        default :
+          break;
+      }
+    }
+    //_setQuads.erase(*itTu); //marche pas avec list
+  } 
   
-  for (int i = 0; i < (*it)->getNumTriangles(); i++) {
-    touched.insert((*it)->getTriangle(i));
+  std::list<TrianglesUnion*>::iterator it = _setQuads.begin();
+  while (it != _setQuads.end()) {
+    if (incomp.find(*it) != incomp.end())
+      it = _setQuads.erase(it); // replacement de _setQuads.erase(*itTu);
+    else
+      ++it;
+  }
+  
+  for (int i = 0; i < tu->getNumTriangles(); i++) {
+    _possibleActions.erase(tu->getTriangle(i));
   }
   
-  for (it = _setQuads.begin(); it != _setQuads.end();) {
+  /*std::list<TrianglesUnion*>::iterator it = _setQuads.begin();
+  while (it != _setQuads.end()) {
     bool b = false;
     for (int i = 0; i < (*it)->getNumTriangles(); i++) {
       if (touched.find((*it)->getTriangle(i)) != touched.end())
         b = true;
     }
-    if (b) {
-      for (int i = 0; i < (*it)->getNumTriangles(); i++) {
-        _isolated.insert((*it)->getTriangle(i));
-      }
-      _setQuads.erase(it++);
-    }
+    if (b)
+      it = _setQuads.erase(it);
     else
-      it++;
-  }
+      ++it;
+  }*/
 }
 
-
-
 Rec2d_vertex::Rec2d_vertex(MVertex *v,
                            std::list<MTriangle*> &setTri,
                            std::map<MVertex*, std::set<GEdge*> > &mapVert)
@@ -304,7 +660,7 @@ Rec2d_vertex::Rec2d_vertex(MVertex *v,
   if (_onWhat < 0) {
     std::map<MVertex*, std::set<GEdge*> >::iterator it = mapVert.find(v);
     if (it != mapVert.end()) {
-      _angle = _computeAngle(v, setTri, (*it).second);
+      _angle = _computeAngle(v, setTri, it->second);
     }
     else {
       Msg::Warning("[meshRecombine2D] Can't compute corner angle, setting angle to pi/2");
@@ -324,7 +680,7 @@ void Rec2d_vertex::_initStaticTable()
     int nE = 5, nF = 10;
     
     _Vvalues = new double*[2];
-    _Vvalues[0] = new double[nE]; //AVERIFIER SI CA MARCHE
+    _Vvalues[0] = new double[nE];
     for (int i = 0; i < nE; i++)
       _Vvalues[0][i] = 1. - fabs(2. / (i+1) - 1.);
     _Vvalues[1] = new double[nF];
@@ -345,7 +701,7 @@ void Rec2d_vertex::_initStaticTable()
 double Rec2d_vertex::_computeAngle(MVertex *v,
                                    std::list<MTriangle*> &setTri,
                                    std::set<GEdge*> &setEdge)
-{
+  {
   if (setEdge.size() != 2) {
     Msg::Warning("[meshRecombine2D] Wrong number of edge : %d, returning pi/2", setEdge.size());
     return M_PI / 2.;
@@ -364,9 +720,9 @@ double Rec2d_vertex::_computeAngle(MVertex *v,
     vectub -= vectv;
     double normlb, normub;
     if ((normlb = norm(vectlb)) > prec * (normub = norm(vectub)))
-      firstDer1 = vectub;
+      firstDer1 = ge->firstDer(ge->getUpperBound());
     else if (normub > prec * normlb)
-      firstDer1 = vectlb;
+      firstDer1 = ge->firstDer(ge->getLowerBound());
     else {
       Msg::Warning("[meshRecombine2D] Bad precision, returning pi/2");
       return M_PI / 2.;
@@ -375,12 +731,8 @@ double Rec2d_vertex::_computeAngle(MVertex *v,
       firstDer0 = firstDer1;
   }
   
-  Msg::Info("-- %lf, %lf, %lf", norm(firstDer0), norm(firstDer1), dot(firstDer0, firstDer1));
-  double n = 1 / norm(firstDer0);
-  firstDer0 = firstDer0 * n;
-  n = 1 / norm(firstDer1);
-  firstDer1 = firstDer1 * n;
-  Msg::Info("--Angles : %lf, %lf, %lf", norm(firstDer0), norm(firstDer1), dot(firstDer0, firstDer1));
+  firstDer0 = firstDer0 * (1 / norm(firstDer0));
+  firstDer1 = firstDer1 * (1 / norm(firstDer1));
   
   double angle1 = acos(dot(firstDer0, firstDer1));
   double angle2 = 2. * M_PI - angle1;
@@ -399,10 +751,6 @@ double Rec2d_vertex::_computeAngle(MVertex *v,
     angleMesh += angle3Vertices(t->getVertex((k+1)%3), v, t->getVertex((k+2)%3));
   }
   
-  Msg::Info("Angles : %lf, %lf, %lf", angle1, angleMesh, angle2);
-  
-  return M_PI / 2.;
-  
   if (angleMesh < M_PI)
     return angle1;
   return angle2;
@@ -420,10 +768,14 @@ double Rec2d_vertex::getReward(int n)
 {
   if (n == 0) return .0;
   if (_onWhat > -1) {
-    if (n > 0)
-      return _Vbenefs[_onWhat][_numEl-1];
-    else
-      return - _Vbenefs[_onWhat][_numEl-2];
+    switch (n) {
+      case 1 :
+        return _Vbenefs[_onWhat][_numEl-1];
+      case -1 :
+        return - _Vbenefs[_onWhat][_numEl-2];
+      default :
+        return _Vvalues[_onWhat][_numEl-1+n] - _Vvalues[_onWhat][_numEl-1];
+    }
   }
   else
     return fabs(2./M_PI*_angle/_numEl - 1.) - fabs(2./M_PI*_angle/(_numEl + n) - 1.);
@@ -452,7 +804,7 @@ TrianglesUnion::TrianglesUnion(GFace *gf,
     double sumlc = .0, u[2], v[2];
     MVertex *vert[2];
     for (int i = 0; i < 2; i++) {
-      vert[i] = (*ite).getVertex(i);
+      vert[i] = ite->getVertex(i);
       vert[i]->getParameter(0, u[i]);
       vert[i]->getParameter(1, v[i]); // Warning : should check if vertex on face or on edge
       std::map<MVertex*,double>::iterator itlc;
@@ -460,13 +812,13 @@ TrianglesUnion::TrianglesUnion(GFace *gf,
         sumlc += v2lcUV[vert[i]] = BGM_MeshSize(gf, u[i], v[i], vert[i]->x(), vert[i]->y(), vert[i]->z());
       }
       else
-        sumlc += (*itlc).second;
+        sumlc += itlc->second;
     }
     
     sumlc = .2; // FIXME BGM_MeshSize returns wrong meshsize
 
     double length = _computeEdgeLength(gf, vert, u, v, 0);
-    double a = 0;
+    double a = .0;
     _embEdgeValue += length / sumlc * (a + (2-a) *_computeAlignment(*ite, t));
     //Msg::Info("Edge a : {%lf/.1 = %lf <-> %lf} => %lf", length, 2 * length / sumlc, 1 - _computeAlignment(*ite, t),length / sumlc * (1 + _computeAlignment(*ite, t)));
   }
@@ -536,17 +888,77 @@ void TrianglesUnion::_update()
 {
   if (_computation >= _RECOMPUTE)
     return;
-  double r = _globValIfExecuted;
   double alpha = _ValVert - _embVertValue;
   for (int i = 0; i < _numBoundVert; i++) {
     alpha += _vertices[i]->getReward(-1);
   }
   alpha /= _NumVert - _numEmbVert;
   double beta = (_ValEdge - _embEdgeValue) / (_NumEdge - _numEmbEdge);
-  _globValIfExecuted = alpha * beta;
+  _globValIfExecuted = alpha * alpha * beta;
   
   _computation = _RECOMPUTE;
 }
+
+void TrianglesUnion::addTri(std::set<MTriangle*> &setTri) const
+{
+  for (int i = 0; i < _numTri; ++i)
+    setTri.insert(_triangles[i]);
+}
+
+void TrianglesUnion::removeTri(std::set<MTriangle*> &setTri) const
+{
+  for (int i = 0; i < _numTri; ++i)
+    setTri.erase(_triangles[i]);
+}
+
+bool TrianglesUnion::isIn(std::set<MTriangle*> &setTri) const
+{
+  for (int i = 0; i < _numTri; ++i)
+  if (setTri.find(_triangles[i]) != setTri.end())
+    return true;
+  return false;
+}
+
+void TrianglesUnion::addInfluence(int *num, double *val, std::map<Rec2d_vertex*, int> &mapVert) const
+{
+  num[0] += _numEmbVert;
+  num[1] += _numEmbEdge;
+  val[0] += _embVertValue;
+  val[1] += _embEdgeValue;
+  for (int i = 0; i < _numBoundVert; ++i) {
+    if (mapVert.find(_vertices[i]) == mapVert.end())
+      mapVert[_vertices[i]] = -1;
+    else
+      mapVert[_vertices[i]] -= 1;  
+  }
+}
+
+void TrianglesUnion::removeInfluence(int *num, double *val, std::map<Rec2d_vertex*, int> &mapVert) const
+{
+  num[0] -= _numEmbVert;
+  num[1] -= _numEmbEdge;
+  val[0] -= _embVertValue;
+  val[1] -= _embEdgeValue;
+  for (int i = 0; i < _numBoundVert; ++i) {
+    mapVert[_vertices[i]] += 1;  
+  }
+}
+
+double TrianglesUnion::computeReturn(const int *num,
+                                     const double *val,
+                                     const std::map<Rec2d_vertex*,
+                                     int> &mapVert)
+{
+  double alpha = _ValVert - val[0];
+  std::map<Rec2d_vertex*, int>::const_iterator it = mapVert.begin();
+  for (; it != mapVert.end(); ++it) {
+    alpha += it->first->getReward(it->second);
+  }
+  alpha /= _NumVert - num[0];
+  double beta = (_ValEdge - val[1]) / (_NumEdge - num[1]);
+  return alpha * beta;
+}
+
 bool TrianglesUnion::operator<(TrianglesUnion &other)
 {
   _update();
@@ -555,6 +967,11 @@ bool TrianglesUnion::operator<(TrianglesUnion &other)
 }
 
 bool lessTriUnion::operator()(TrianglesUnion *tu1, TrianglesUnion *tu2) const
+{
+  return *tu1 < *tu2;
+}
+
+bool lessTriUnionInv::operator()(TrianglesUnion *tu1, TrianglesUnion *tu2) const
 {
   return *tu2 < *tu1;
 }