From fe63d33594a551e8e6690528cd1de800a46dc80b Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Wed, 18 Jan 2012 16:53:40 +0000
Subject: [PATCH]

---
 Mesh/meshGFaceRecombine.cpp | 1590 +++++++++++++++++++----------------
 Mesh/meshGFaceRecombine.h   |  215 +++--
 2 files changed, 1038 insertions(+), 767 deletions(-)

diff --git a/Mesh/meshGFaceRecombine.cpp b/Mesh/meshGFaceRecombine.cpp
index af69ab0b1a..b1b38c7347 100644
--- a/Mesh/meshGFaceRecombine.cpp
+++ b/Mesh/meshGFaceRecombine.cpp
@@ -7,6 +7,12 @@
 //   Amaury Johnen (a.johnen@ulg.ac.be)
 //
 
+#define REC2D_EDGE_BASE 1
+#define REC2D_EDGE_QUAD 1
+#define REC2D_ALIGNMENT .5
+#define REC2D_WAIT_TIME .01
+#define REC2D_NUM_ACTIO 550
+
 #include "meshGFaceRecombine.h"
 #include "MElement.h"
 #include "MTriangle.h"
@@ -18,10 +24,6 @@
 #include "Context.h"
 #include "OS.h"
 
-#define REC2D_EDGE_BASE 1
-#define REC2D_EDGE_QUAD 1
-#define REC2D_WAIT_TIME 1
-
 // #define REC2D_SMOOTH
  #define REC2D_DRAW
 
@@ -30,20 +32,28 @@ bool Recombine2D::bo = 0;
 double **Rec2DVertex::_qualVSnum = NULL;
 double **Rec2DVertex::_gains = NULL;
 
+int otherParity(int a) {
+  if (a % 2)
+    return a - 1;
+  return a + 1;
+}
+
+/*
+  TODO :
+    - map -> myset<pair> (mapCornerVert, mapVert, mapEdge)
+    - set -> myset
+    (_vertices, _edges, _elements) -> set
+    (_vertices, _edges, _elements) -> set
+*/
+
 void Recombine2D::associateParity(int pOld, int pNew)
 {
   std::set<Rec2DVertex*>::iterator it = _current->_parities[pOld].begin();
   for (; it != _current->_parities[pOld].end(); ++it)
     (*it)->setParity(pNew);
   _current->_parities.erase(pOld);
-  if (pOld % 2)
-    pOld -= 1;
-  else
-    pOld +=1;
-  if (pNew % 2)
-    pNew -= 1;
-  else
-    pNew +=1;
+  pOld = otherParity(pOld);
+  pNew = otherParity(pNew);
   it = _current->_parities[pOld].begin();
   for (; it != _current->_parities[pOld].end(); ++it)
     (*it)->setParity(pNew);
@@ -55,21 +65,22 @@ void Recombine2D::associateParity(int pOld, int pNew)
 Recombine2D::Recombine2D(GFace *gf) : _gf(gf)
 {
   if (Recombine2D::_current != NULL) {
-    Msg::Info("An instance of recombination is already in execution");
+    Msg::Warning("An instance of recombination is already in execution");
     return;
   }
   Recombine2D::_current = this;
-  _numChange = 0;
+  
   backgroundMesh::set(gf);
   _bgm = backgroundMesh::current();
-  Rec2DVertex::initStaticTable();
+  _numChange = 0;
   _numEdge = _numVert = 0;
   _valEdge = _valVert = .0;
+  Rec2DVertex::initStaticTable();
   
     _tri = _gf->triangles;
   
   // Be able to compute geometrical angle at corners
-  std::map<MVertex*, std::set<GEdge*> > mapCornerVert;
+  std::map<MVertex*, std::set<GEdge*> > mapCornerVert; // TOBE myset
   {
     std::list<GEdge*> listge = gf->edges();
     std::list<GEdge*>::iterator itge = listge.begin();
@@ -81,11 +92,13 @@ Recombine2D::Recombine2D(GFace *gf) : _gf(gf)
     }
   }
   
-  // Find all Vertices and edges
-  std::map<MVertex*, std::list<MTriangle*> > mapVert;
-  std::map<MEdge, std::list<MTriangle*>, Less_Edge> mapEdge;
+  // Find all Vertices and edges and create the 'Rec2DElement'
+  std::map<MVertex*, std::vector<MTriangle*> > mapVert; // TOBE myset
+  std::map<MEdge, std::vector<MTriangle*>, Less_Edge> mapEdge; // TOBE myset
   for (unsigned int i = 0; i < gf->triangles.size(); ++i) {
     MTriangle *t = gf->triangles[i];
+    Rec2DElement *el = new Rec2DElement(t);
+    _elements[t] = el; // ISITTOKEEP keep only list of relem, map local
     for (int j = 0; j < 3; ++j) {
       mapVert[t->getVertex(j)].push_back(t);
       mapEdge[t->getEdge(j)].push_back(t);
@@ -93,34 +106,31 @@ Recombine2D::Recombine2D(GFace *gf) : _gf(gf)
   }
   
   // Create the 'Rec2DVertex' and store iterator to vertices which have degree 4
-  std::list<std::map<MVertex*, std::list<MTriangle*> >::iterator> fourTri;
+  std::list<std::map<MVertex*, std::vector<MTriangle*> >::iterator> fourTri; // TOBE vector
   {
     mapofVertices::iterator itV2N = _vertices.begin();
-    std::map<MVertex*, std::list<MTriangle*> >::iterator itvert = mapVert.begin();
+    std::map<MVertex*, std::vector<MTriangle*> >::iterator itvert = mapVert.begin();
     for (; itvert != mapVert.end(); ++itvert) {
       MVertex *v = itvert->first;
       int onWhat = 1;
       if (mapCornerVert.find(v) != mapCornerVert.end())
         onWhat = -1;
       Rec2DVertex *rv = new Rec2DVertex(v, itvert->second, onWhat, mapCornerVert);
-      itV2N = _vertices.insert(itV2N, std::make_pair(v,rv));
-      ++_numVert;
+      itV2N = _vertices.insert(itV2N, std::make_pair(v,rv)); // ISITTOKEEP keep only list of rvert, map local
       if (itvert->second.size() == 4)
         fourTri.push_back(itvert);
     }
   }
   
   // Create the 'Rec2DEdge' and store boundary edges
-  std::map<MVertex*, std::list<MEdge> > boundV2E;
+  std::map<MVertex*, std::list<MEdge> > boundV2E; // TOBE myset
   {
     mapofEdges::iterator itE2E = _edges.begin();
-    std::map<MEdge, std::list<MTriangle*> >::iterator itedge = mapEdge.begin();
+    std::map<MEdge, std::vector<MTriangle*> >::iterator itedge = mapEdge.begin();
     for (; itedge != mapEdge.end(); ++itedge) {
       MEdge e = itedge->first;
       Rec2DEdge *re = new Rec2DEdge(e, _vertices, itedge->second);
       itE2E = _edges.insert(itE2E, std::make_pair(e,re));
-      _vertices[e.getVertex(0)]->add(re);
-      _vertices[e.getVertex(1)]->add(re);
       _numEdge += REC2D_EDGE_BASE;
       _valEdge += REC2D_EDGE_BASE * re->getQual();
       if (itedge->second.size() == 1) {
@@ -137,8 +147,10 @@ Recombine2D::Recombine2D(GFace *gf) : _gf(gf)
   // We know now on what are vertices, compute reward
   {
     mapofVertices::iterator itV2N = _vertices.begin();
-    for (; itV2N != _vertices.end(); ++itV2N)
+    for (; itV2N != _vertices.end(); ++itV2N) {
+      ++_numVert;
       _valVert += itV2N->second->getQual();
+    }
   }
   
   // Be dealing with "parity" on boundaries
@@ -186,14 +198,14 @@ Recombine2D::Recombine2D(GFace *gf) : _gf(gf)
   // Create the 'Rec2DTwoTri2Quad' and 'Rec2DCollapseTri'
   {
     mapofEdges::iterator itE2E = _edges.begin();
-    std::map<MEdge, std::list<MTriangle*> >::iterator ite = mapEdge.begin();
+    std::map<MEdge, std::vector<MTriangle*> >::iterator ite = mapEdge.begin();
     for (; ite != mapEdge.end(); ++ite)
     if (ite->second.size() == 2) {
       Rec2DTwoTri2Quad *q = new Rec2DTwoTri2Quad(_edges[ite->first], ite->second);
       //Rec2DCollapseTri *c = new Rec2DCollapseTri(_edges[ite->first], ite->second);
       _actions.push_back(q);
       //_actions.insert(c);
-      std::list<MTriangle*>::iterator it = ite->second.begin();
+      std::vector<MTriangle*>::iterator it = ite->second.begin();
       _mea[*it].insert(q);
       //_mea[*it].insert(c);
       ++it;
@@ -202,58 +214,24 @@ Recombine2D::Recombine2D(GFace *gf) : _gf(gf)
     }
   }
   
-  // Rec2DTwoTri2Quad
-  /*{
-    std::list<std::map<MVertex*, std::list<MTriangle*> >::iterator>::iterator it4;
-    for (it4 = fourTri.begin(); it4 != fourTri.end(); it4++) {
-      Rec2DVertex *rv = _vertices[(*it4)->first];
-      std::set<MEdge, Less_Edge> setEdge;
-      {
-        std::list<MTriangle*>::iterator it = (*it4)->second.begin();
-        for (; it != (*it4)->second.end(); ++it) {
-          setEdge.insert((*it)->getEdge(0));
-          setEdge.insert((*it)->getEdge(1));
-          setEdge.insert((*it)->getEdge(2));
-        }
-      }
-      
-      std::set<MEdge>::iterator it = setEdge.begin();
-      for (; it != setEdge.end(); ++it) {
-        std::map<MEdge, std::list<MTriangle*> >::iterator it2;
-        it2 = mapEdge.find(*it);
-        if (it2 != mapEdge.end())
-        if (it2->second.size() == 2) {
-          Rec2DTwoTri2Quad *q = new Rec2DTwoTri2Quad(_edges[it2->first], it2->second);
-          _actions.insert(q);
-          std::list<MTriangle*>::iterator it = it2->second.begin();
-          _mea[*it].insert(q);
-          ++it;
-          _mea[*it].insert(q);
-        }
-        else
-          Msg::Error("fjdaks;lfjkdsa;");
-      }
-    }
-  }*/
-  
   // Create the 'Rec2DFourTri2Quad' and 'Rec2DCollapseTri'
-  {
-    std::list<std::map<MVertex*, std::list<MTriangle*> >::iterator>::iterator it4;
-    for (it4 = fourTri.begin(); it4 != fourTri.end(); it4++) {
-      Rec2DVertex *rv = _vertices[(*it4)->first];
-      if (!rv->getIsOnBoundary()) {
-        Rec2DFourTri2Quad *q = new Rec2DFourTri2Quad(rv, (*it4)->second);
-        //Rec2DCollapseTri *c = new Rec2DCollapseTri(rv, (*it4)->second);
-        _actions.push_back(q);
-        //_actions.insert(c);
-        std::list<MTriangle*>::iterator it = (*it4)->second.begin();
-        for (; it != (*it4)->second.end(); ++it) {
-          _mea[*it].insert(q);
-          //_mea[*it].insert(c);
-        }
-      }
-    }
-  }
+  //{
+  //  std::list<std::map<MVertex*, std::vector<MTriangle*> >::iterator>::iterator it4;
+  //  for (it4 = fourTri.begin(); it4 != fourTri.end(); it4++) {
+  //    Rec2DVertex *rv = _vertices[(*it4)->first];
+  //    if (!rv->getIsOnBoundary()) {
+  //      Rec2DFourTri2Quad *q = new Rec2DFourTri2Quad(rv, (*it4)->second);
+  //      //Rec2DCollapseTri *c = new Rec2DCollapseTri(rv, (*it4)->second);
+  //      _actions.push_back(q);
+  //      //_actions.insert(c);
+  //      std::vector<MTriangle*>::iterator it = (*it4)->second.begin();
+  //      for (; it != (*it4)->second.end(); ++it) {
+  //        _mea[*it].insert(q);
+  //        //_mea[*it].insert(c);
+  //      }
+  //    }
+  //  }
+  //}
   
   Msg::Info("State");
   Msg::Info("-----");
@@ -278,34 +256,43 @@ Recombine2D::Recombine2D(GFace *gf) : _gf(gf)
 Recombine2D::~Recombine2D()
 {
   Recombine2D::_current = NULL;
-  //delete vertices, edges, actions;
+  //delete vertices, edges, elements, actions;
 }
 
 bool Recombine2D::recombine()
 {
   Rec2DAction *nextAction;
-  while (_actions.size() > 0) {
+  int i = -1;
+  while (_actions.size() > 0 && ++i < REC2D_NUM_ACTIO) {
     double t = Cpu();
     nextAction = *std::max_element(_actions.begin(), _actions.end(),
                                    lessRec2DAction()                );
-    if (nextAction->getReward() > 0) {
+    nextAction->color(0, 0, 200);
+    
+    CTX::instance()->mesh.changed = (ENT_ALL);
+    FlGui::instance()->check();
+    drawContext::global()->draw();
+    //if (nextAction->getReward() > 0) {
       if (_isAllQuad(nextAction))
         _choose(nextAction);
       else {
         nextAction->color(190, 0, 0);
         _removeReference(nextAction);
       }
-    }
+    /*}
     else {
       nextAction->color(0, 100, 0);
       _removeReference(nextAction);
-    }
+    }*/
     
     _gf->triangles = _tri;
     _gf->quadrangles = _quad;
 #ifdef REC2D_DRAW
+    //Msg::Info("%d", REC2D_NUM_ACTIO-i);
+    //Msg::GetAnswer("Continue ?", 1, "no", "yes");
     CTX::instance()->mesh.changed = (ENT_ALL);
     FlGui::instance()->check();
+    if (i%5 == 0)
     drawContext::global()->draw();
     while (Cpu()-t < REC2D_WAIT_TIME)
       FlGui::instance()->check();
@@ -314,47 +301,64 @@ bool Recombine2D::recombine()
     laplaceSmoothing(_gf,1);
 #endif
 
-  std::map<int, std::set<Rec2DVertex*> >::iterator it = _parities.begin();
-  for (; it != _parities.end(); ++it) {
-    Msg::Info("par %d, #%d", it->first, it->second.size());
-  }
-  }
-  
-  Msg::Info("State");
-  Msg::Info("-----");
-  Msg::Info("numEdge %d (%d), valEdge %g", _numEdge, _edges.size(), _valEdge);
-  Msg::Info("numVert %d (%d), valVert %g", _numVert, _vertices.size(), _valVert);
-  Msg::Info("global Value %g", Recombine2D::getGlobalValue(0,0,0,0));
-  Msg::Info("num action %d", _actions.size());
-  std::map<int, std::set<Rec2DVertex*> >::iterator it = _parities.begin();
-  for (; it != _parities.end(); ++it) {
-    Msg::Info("par %d, #%d", it->first, it->second.size());
-  }
-  mapofElementActions::iterator it2 = _mea.begin();
-  int a = 0;
-  for (; it2 != _mea.end(); ++it2) {
-    a += it2->second.size();
+    std::map<int, std::set<Rec2DVertex*> >::iterator it = _parities.begin();
+    for (; it != _parities.end(); ++it) {
+      //Msg::Info("par %d, #%d", it->first, it->second.size());
+    }
   }
-  if (a > 0)
-    Msg::Info("%d elements in mapof elem to action, with average action %g", _mea.size(), (double)a/_mea.size());
-  Msg::Info(" ");
-  
+//  
+//  Msg::Info("State");
+//  Msg::Info("-----");
+//  Msg::Info("numEdge %d (%d), valEdge %g", _numEdge, _edges.size(), _valEdge);
+//  Msg::Info("numVert %d (%d), valVert %g", _numVert, _vertices.size(), _valVert);
+//  Msg::Info("global Value %g", Recombine2D::getGlobalValue(0,0,0,0));
+//  Msg::Info("num action %d", _actions.size());
+//  std::map<int, std::set<Rec2DVertex*> >::iterator it = _parities.begin();
+//  for (; it != _parities.end(); ++it) {
+//    Msg::Info("par %d, #%d", it->first, it->second.size());
+//  }
+//  mapofElementActions::iterator it2 = _mea.begin();
+//  int a = 0;
+//  for (; it2 != _mea.end(); ++it2) {
+//    a += it2->second.size();
+//  }
+//  if (a > 0)
+//    Msg::Info("%d elements in mapof elem to action, with average action %g", _mea.size(), (double)a/_mea.size());
+//  Msg::Info(" ");
+//  
   return 1;
 }
 
 bool Recombine2D::_isAllQuad(Rec2DAction *action)
 {
+  //Msg::Info("    ---new isAllQuad---");
+  
   int p[4];
   action->getParities(p);
   
-  if (p[0] < 0 && p[1] < 0 && p[2] < 0 && p[3] < 0)
+  if (p[0] < 0 && p[1] < 0 && p[2] < 0 && p[3] < 0) {
+    Msg::Info("is isolated");
     return true;
+  }
   
   if (action->isObsolete()) {
-    Msg::Error("[Recombine2D] obsolete action");
+    Msg::Error("[Recombine2D] obsolete action (allQuad)");
     return false;
   }
   
+  std::set<Rec2DVertex*> setRv;
+  std::set<Rec2DElement*> neighbours, checked;
+  for (int i = 0; i < Recombine2D::getNumVertAllQuad(); ++i) {
+    if (Recombine2D::getVertAllQuad(i)->hasTriangle())
+      Recombine2D::getVertAllQuad(i)->getElements(neighbours);
+    else
+      setRv.insert(Recombine2D::getVertAllQuad(i));
+  }
+  std::set<Rec2DVertex*>::iterator itit = setRv.begin();
+  for (; itit != setRv.end(); ++itit)
+    Recombine2D::removeVertAllQuad(*itit);
+  setRv.clear();
+  
   int min = 100, index;
   for (int i = 0; i < 4; ++i) {
     if (p[i] > -1 && min > p[i]) {
@@ -362,254 +366,244 @@ bool Recombine2D::_isAllQuad(Rec2DAction *action)
       index = i;
     }
   }
-  
-  //Quid si p[0..4] = -1 ?
+  //Msg::Info("Passsed through there, [%d %d %d %d] -> min %d", p[0], p[1], p[2], p[3], min);
   
   std::map<Rec2DVertex*, int> newParVert;
+  std::map<Rec2DVertex*, int>::iterator it2;
   std::map<int, int> parAssoc; // first parity value become second parity value
+  std::map<int, int>::iterator itassoc;
   
   for (int i = 0; i < 4; i += 2) {
     int par;
     if ((index/2) * 2 == i)
       par = min;
-    else if (min % 2)
-      par = min - 1;
     else
-      par = min + 1;
+      par = otherParity(min);
     for (int j = 0; j < 2; ++j) {
-      if (p[i+j] == -1)
+      if (p[i+j] == -1) {
         newParVert.insert(std::make_pair(action->getRVertex(i+j), par));
-      else if (p[i+j] != par)
-        parAssoc.insert(std::make_pair(p[i+j], par));
+        action->getRVertex(i+j)->getElements(neighbours);
+      }
+      else if (p[i+j] != par) {
+        if (parAssoc.find(p[i+j]) == parAssoc.end()) {
+          parAssoc[p[i+j]] = par;
+          parAssoc[otherParity(p[i+j])] = otherParity(par);
+          _getElemToLook(p[i+j], neighbours);
+          _getElemToLook(otherParity(p[i+j]), neighbours);
+        }
+      }
     }
   }
   
-  std::set<Rec2DAction*> neighbours;
-  _getNeighbours(action, neighbours);
-  
   while (neighbours.size() > 0) {
-    std::set<Rec2DAction*>::iterator it = neighbours.begin();
-    if ((*it)->getNum() != 1) {
-      neighbours.erase(it);
+    //Msg::Info("num vert %d, num neigh %d", newParVert.size(), neighbours.size());
+    
+    std::set<Rec2DElement*>::iterator itel = neighbours.begin();
+    Rec2DElement *element = *itel;
+    if (element->isQuad()) {
+      neighbours.erase(itel);
       continue;
     }
-    Rec2DVertex *v0 = (*it)->getRVertex(0);
-    Rec2DVertex *v1 = (*it)->getRVertex(1);
-    std::map<Rec2DVertex*, int>::iterator it2;
-    int p0 = v0->getParity(), p1 = v1->getParity();
-    if (p0 < 0 && (it2 = newParVert.find(v0)) == newParVert.end()) {
-      neighbours.erase(it);
-      continue;
+    
+    Rec2DVertex *v[3];
+    element->getVertices(v);
+    int p[3];
+    for (int i = 0; i < 3; ++i) {
+      if ((it2 = newParVert.find(v[i])) != newParVert.end()) {
+        p[i] = it2->second;
+        while ((itassoc = parAssoc.find(p[i])) != parAssoc.end()) {
+          p[i] = itassoc->second;
+        }
+        it2->second = p[i];
+      }
+      else {
+        p[i] = v[i]->getParity();
+        if ((itassoc = parAssoc.find(p[i])) != parAssoc.end()) {
+          do
+            p[i] = itassoc->second;
+          while ((itassoc = parAssoc.find(p[i])) != parAssoc.end());
+          newParVert[v[i]] = p[i];
+        }
+      }
     }
-    if (p0 < 0)
-      p0 = it2->second;
+  //Msg::Info("tri %d [%d %d %d]", (*itel)->getNum(), p[0], p[1], p[2]);
     
-    if (p1 < 0 && (it2 = newParVert.find(v1)) == newParVert.end()) {
-      neighbours.erase(it);
+    bool hasIdentical = false;
+    for (int i = 0; i < 3; ++i) {
+      if (p[i] > -1 && p[i] == p[(i+1)%3]) hasIdentical = true;
+    }
+    if (!hasIdentical) {
+      neighbours.erase(itel);
       continue;
     }
-    if (p1 < 0)
-      p1 = it2->second;
+    if (p[0] == p[1] && p[0] == p[2]) {
+      Msg::Info("3 identical par");
+      return false;
+    }
+  //Msg::Info("has identical");
     
-    std::map<int, int>::iterator it3;
-    if ((it3 = parAssoc.find(p0)) != parAssoc.end())
-      p0 = it3->second;
-    if ((it3 = parAssoc.find(p1)) != parAssoc.end())
-      p1 = it3->second;
+    bool hasAction = false;
+    std::map<Rec2DVertex*, std::vector<int> > suggestions;
+    for (int i = 0; i < element->getNumActions(); ++i) {
+      if (element->getAction(i)->whatImply(v, p, suggestions))
+        hasAction = true;
+    }
+  //Msg::Info("suggest %d", suggestions.size());
     
-    if (p0 != p1) {
-      neighbours.erase(it);
-      continue;
+    if (!hasAction) {
+      Msg::Info("No action %d", (*itel)->getNum());
+      return false;
     }
     
-    int par;
-    bool hasChange = false;
-    if (p0 % 2)
-      par = p0 - 1;
-    else
-      par = p0 + 1;
-    (*it)->getParities(p);
-    for (int j = 2; j < 4; ++j) {
-      if (p[j] == -1) {
-        if ((it2 = newParVert.find((*it)->getRVertex(j))) == newParVert.end()) {
-          newParVert.insert(std::make_pair(action->getRVertex(j), par));
-          hasChange = true;
+    std::map<Rec2DVertex*, std::vector<int> >::iterator itsug;
+    itsug = suggestions.begin();
+    for (; itsug != suggestions.end(); ++itsug) {
+      int par = itsug->second.front();
+      bool onlyOnePar = true;
+      for (unsigned int i = 1; i < itsug->second.size(); ++i) {
+        if (itsug->second[i] != par) {
+          onlyOnePar = false;
+          break;
         }
-        else if (it2->second != par)
-          return false;
       }
-      else if (p[j] != par) {
-        if ((it3 = parAssoc.find(p[j])) == parAssoc.end()) {
-          parAssoc.insert(std::make_pair(p[j], par));
-          hasChange = true;
+      if (onlyOnePar) {
+        if ((itassoc = parAssoc.find(par)) != parAssoc.end()) {
+          do
+            par = itassoc->second;
+          while ((itassoc = parAssoc.find(par)) != parAssoc.end());
+        }
+        
+        int oldPar;
+        Rec2DVertex *v = itsug->first;
+        
+        if ((it2 = newParVert.find(v)) != newParVert.end()) {
+          oldPar = it2->second;
+          while ((itassoc = parAssoc.find(oldPar)) != parAssoc.end()) {
+            oldPar = itassoc->second;
+          }
+          it2->second = oldPar;
         }
         else {
-          int a = it3->second;
-          while (a != par && (it3 = parAssoc.find(a)) != parAssoc.end()) {
-            a = it3->second;
+          oldPar = v->getParity();
+          if ((itassoc = parAssoc.find(oldPar)) != parAssoc.end()) {
+            do
+              oldPar = itassoc->second;
+            while ((itassoc = parAssoc.find(oldPar)) != parAssoc.end());
+            newParVert[v] = oldPar;
+          }
+        }
+        
+        //Msg::Info("a %d, %d", par, oldPar);
+        if (oldPar == -1) {
+        //Msg::Info("b");
+          newParVert[v] = par;
+          v->getElements(neighbours);
+        }
+        else if ((par/2)*2 != (oldPar/2)*2) {
+        //Msg::Info("c");
+          if (oldPar < par) {
+            int a = oldPar;
+            oldPar = par;
+            par = a;
           }
-          if (a != par)
-            return false;
+          parAssoc[oldPar] = par;
+          parAssoc[otherParity(oldPar)] = otherParity(par);
+          _getElemToLook(oldPar, neighbours);
+          _getElemToLook(otherParity(oldPar), neighbours);
+        }
+        else if (par%2 != oldPar%2) {
+        //Msg::Info("d");
+          Msg::Info("not all quad");
+          return false;
         }
       }
     }
-    Rec2DAction *ra = *it;
-    neighbours.erase(it);
-    if (hasChange)
-      _getNeighbours(ra, neighbours);
+    neighbours.erase(itel);
   }
+  Msg::Info("all quad");
   return true;
 }
 
 void Recombine2D::_choose(Rec2DAction *action)
 {
-  action->apply();
-  _removeIncompatible(action);
-  ++_numChange;
+  action->apply(); //to check
+  _removeIncompatible(action); //to check
+  ++_numChange; //to check
 }
 
-void Recombine2D::_removeIncompatible(Rec2DAction *action)
+void Recombine2D::_removeIncompatible(Rec2DAction *action) //to check
 {
-  std::set<MTriangle*> touchedTri;
-  action->getTriangles(touchedTri);
-  std::set<MTriangle*>::iterator itTouched;
+  std::set<Rec2DElement*> touchedTri;
+  action->getRTriangles(touchedTri);
+  std::set<Rec2DElement*>::iterator itTri;
   
   std::set<Rec2DAction*> incomp;
   {
-    itTouched = touchedTri.begin();
-    for (; itTouched != touchedTri.end(); ++itTouched) {
-      mapofElementActions::iterator ite2a = _mea.find(*itTouched);
-      if (ite2a != _mea.end()) {
-        std::set<Rec2DAction*>::iterator ita = ite2a->second.begin();
-        for (; ita != ite2a->second.end(); ++ita)
-          incomp.insert(*ita);
-      }
-      else
-        Msg::Error("[Recombine2D] didn't found triangle %d (rmIncomp1)", (*itTouched)->getNum());
-    }
+    itTri = touchedTri.begin();
+    for (; itTri != touchedTri.end(); ++itTri)
+      (*itTri)->getActions(incomp);
   }
   
-  std::set<Rec2DAction*>::iterator itIncomp = incomp.begin();
-  for (; itIncomp != incomp.end(); ++itIncomp) {
-    std::set<MTriangle*> touched;
-    (*itIncomp)->getTriangles(touched);
-    for (itTouched = touched.begin(); itTouched != touched.end(); ++itTouched) {
-      mapofElementActions::iterator ite2a = _mea.find(*itTouched);
-      if (ite2a != _mea.end()) {
-        bool b = ite2a->second.erase(*itIncomp);
-        if (!b)
-          Msg::Error("[Recombine2D] Tri->Action not found (rmIncomp)");
-      }
-      else
-        Msg::Error("[Recombine2D] didn't found triangle %d (rmIncomp)", (*itTouched)->getNum());
-    }
-    //bool b = _actions.erase(*itIncomp);
-    setofRec2DAction::iterator it2 = _actions.begin();
-    bool b = false;
-    while (it2 != _actions.end()) {
-      if (*it2 == *itIncomp) {
-        _actions.erase(it2);
-        b = true;
-        break;
-      }
-      ++it2;
-    }
-    if (!b)
-      Msg::Error("[Recombine2D] action not found (rmIncomp2)");
-    delete *itIncomp;
-  }
+  //Msg::Warning("size %d", incomp.size());
   
-  itTouched = touchedTri.begin();
-  for (; itTouched != touchedTri.end(); ++itTouched) {
-    bool b = _mea.erase(*itTouched);
-    if (!b)
-      Msg::Error("[Recombine2D] triangle not found %d", (*itTouched)->getNum());
+  std::set<Rec2DAction*>::iterator itAction = incomp.begin();
+  for (; itAction != incomp.end(); ++itAction) {
+    delete *itAction;
   }
 }
 
 void Recombine2D::_removeReference(Rec2DAction *action)
 {
-  std::set<MTriangle*> touchedTri;
-  action->getTriangles(touchedTri);
-  std::set<MTriangle*>::iterator it = touchedTri.begin();
-  for (; it != touchedTri.end(); ++it) {
-    mapofElementActions::iterator ite2a = _mea.find(*it);
-    if (ite2a != _mea.end()) {
-      bool b = ite2a->second.erase(action);
-      if (!b)
-        Msg::Error("[Recombine2D] Tri->Action not found (rmRef)");
-    }
-    else
-      Msg::Error("[Recombine2D] didn't found triangle %d (rmRef)", (*it)->getNum());
-  }
-  //bool b = _actions.erase(action);
-  setofRec2DAction::iterator it2 = _actions.begin();
-  bool b = false;
-  while (it2 != _actions.end()) {
-    if (*it2 == action) {
-      _actions.erase(it2);
-      b = true;
-      break;
-    }
-    ++it2;
-  }
-  if (!b)
-    Msg::Error("[Recombine2D] action not found (rmRef)");
   delete action;
 }
 
-void Recombine2D::_getNeighbours(Rec2DAction *ra, std::set<Rec2DAction*> &actions)
+//void Recombine2D::_getNeighbours(Rec2DAction *ra, std::set<Rec2DAction*> &actions)
+//{
+//  std::set<Rec2DElement*> touchedTri, neighbourTri;
+//  ra->getRTriangles(touchedTri);
+//  std::set<Rec2DElement*>::iterator itTri;
+//  
+//  {
+//    itTri = touchedTri.begin();
+//    for (; itTri != touchedTri.end(); ++itTri)
+//      (*itTri)->getNeighbours(neighbourTri);
+//    //itTri = touchedTri.begin();
+//    //for (; itTri != touchedTri.end(); ++itTri)
+//    //  neighbourTri.erase(*itTri);
+//  }
+//  
+//  std::set<Rec2DAction*> incomp, neighbours;
+//  {
+//    itTri = touchedTri.begin();
+//    for (; itTri != touchedTri.end(); ++itTri)
+//      (*itTri)->getActions(neighbours);
+//    itTri = neighbourTri.begin();
+//    for (; itTri != neighbourTri.end(); ++itTri)
+//      (*itTri)->getActions(neighbours);
+//  }
+//  
+//  std::set<Rec2DAction*>::iterator itAction = incomp.begin();
+//  for (; itAction != incomp.end(); ++itAction) {
+//    bool b = neighbours.erase(*itAction);
+//    if (!b)
+//      Msg::Error("[Recombine2D] didn't found incomp action");
+//  }
+//  
+//  for (itAction = neighbours.begin(); itAction != neighbours.end(); ++itAction) {
+//    actions.insert(*itAction);
+//  }
+//}
+
+void Recombine2D::_getElemToLook(int par, std::set<Rec2DElement*> &elem)
 {
-  std::set<MTriangle*> touchedTri;
-  ra->getTriangles(touchedTri);
-  std::set<MTriangle*>::iterator itTouched;
-  
-  std::set<Rec2DAction*> incomp;
-  {
-    itTouched = touchedTri.begin();
-    for (; itTouched != touchedTri.end(); ++itTouched) {
-      mapofElementActions::iterator ite2a = _mea.find(*itTouched);
-      if (ite2a != _mea.end()) {
-        std::set<Rec2DAction*>::iterator ita = ite2a->second.begin();
-        for (; ita != ite2a->second.end(); ++ita)
-          incomp.insert(*ita);
-      }
-      else
-        Msg::Error("[Recombine2D.] didn't found triangle %d (rmIncomp1)", (*itTouched)->getNum());
-    }
-  }
-  
-  std::set<Rec2DAction*>::iterator it = incomp.begin();
-  for (; it != incomp.end(); ++it) {
-    (*it)->getTriangles(touchedTri);
-  }
-  
-  std::set<Rec2DAction*> neighbours;
-  {
-    itTouched = touchedTri.begin();
-    for (; itTouched != touchedTri.end(); ++itTouched) {
-      mapofElementActions::iterator ite2a = _mea.find(*itTouched);
-      if (ite2a != _mea.end()) {
-        std::set<Rec2DAction*>::iterator ita = ite2a->second.begin();
-        for (; ita != ite2a->second.end(); ++ita)
-          neighbours.insert(*ita);
-      }
-      else
-        Msg::Error("[Recombine2D.] didn't found triangle %d (rmIncomp1)", (*itTouched)->getNum());
-    }
-    std::set<Rec2DAction*>::iterator itIncomp = incomp.begin();
-    for (; itIncomp != incomp.end(); ++itIncomp) {
-      bool b = neighbours.erase(*itIncomp);
-      if (!b)
-        Msg::Error("[Recombine2D] didn't found incomp action");
-    }
-  }
-  
-  for (it = neighbours.begin(); it != neighbours.end(); ++it) {
-    actions.insert(*it);
+  std::set<Rec2DVertex*>::iterator it = _parities[par].begin();
+  std::set<Rec2DVertex*>::iterator itend = _parities[par].end();
+  for (; it != itend; ++it) {
+    (*it)->getElements(elem);
   }
 }
 
+
 double Recombine2D::getGlobalValue(int numEdge, double valEdge,
                                    int numVert, double valVert )
 {
@@ -617,6 +611,15 @@ double Recombine2D::getGlobalValue(int numEdge, double valEdge,
   return (_current->_valEdge + valEdge) / (_current->_numEdge + numEdge) * a * a;
 }
 
+//Rec2DVertex* Recombine2D::getRVertex(MVertex *v)
+//{
+//  mapofVertices::iterator it = _current->_vertices.find(v);
+//  if (it != _current->_vertices.end())
+//    return it->second;
+//  Msg::Error("[Recombine2D::getRVertex] vertex not found");
+//  return NULL;
+//}
+
 Rec2DEdge* Recombine2D::getREdge(MEdge e)
 {
   mapofEdges::iterator it = _current->_edges.find(e);
@@ -626,12 +629,12 @@ Rec2DEdge* Recombine2D::getREdge(MEdge e)
   return NULL;
 }
 
-Rec2DVertex* Recombine2D::getRVertex(MVertex *v)
+Rec2DElement* Recombine2D::getRElement(MElement *el)
 {
-  mapofVertices::iterator it = _current->_vertices.find(v);
-  if (it != _current->_vertices.end())
+  mapofElements::iterator it = _current->_elements.find(el);
+  if (it != _current->_elements.end())
     return it->second;
-  Msg::Error("[Recombine2D::getRVertex] vertex not found");
+  Msg::Error("[Recombine2D::getRElement] element not found");
   return NULL;
 }
 
@@ -653,16 +656,38 @@ void Recombine2D::remove(MVertex *v)
     Msg::Error("[Recombine2D::getRVertex] vertex not found");
 }
 
-void Recombine2D::remove(MTriangle *tri)
+void Recombine2D::remove(Rec2DElement *tri)
 {
+  Rec2DVertex *v[3];
+  tri->getVertices(v);
+  for (int i = 0; i < 3; ++i) {
+    v[i]->remove(tri);
+  }
+  MTriangle *t = tri->getMTriangle();
   std::vector<MTriangle*>::iterator it = _current->_tri.begin();
   for (; it != _current->_tri.end(); ++it) {
-    if (*it == tri) {
+    if (*it == t) {
       _current->_tri.erase(it);
       return;
     }
   }
-  Msg::Error("[Recombine2D] triangle %d was not there", tri->getNum());
+  Msg::Error("[Recombine2D] triangle was not there");
+}
+
+void Recombine2D::remove(Rec2DAction *ra)
+{
+  setofRec2DAction::iterator itAction = _current->_actions.begin();
+  bool b = false;
+  while (itAction != _current->_actions.end()) {
+    if (*itAction == ra) {
+      _current->_actions.erase(itAction);
+      b = true;
+      break;
+    }
+    ++itAction;
+  }
+  if (!b)
+    Msg::Error("[Recombine2D] action to delete not found");
 }
 
 void Recombine2D::add(MQuadrangle *quad)
@@ -673,6 +698,12 @@ void Recombine2D::add(MQuadrangle *quad)
 
 /**  Rec2DAction  **/
 /*******************/
+bool lessRec2DAction::operator()(Rec2DAction *ra1, Rec2DAction *ra2) const
+{
+  return *ra1 < *ra2;
+}
+
+
 bool Rec2DAction::operator<(Rec2DAction &other)
 {
   return getReward() < other.getReward(); 
@@ -680,6 +711,7 @@ bool Rec2DAction::operator<(Rec2DAction &other)
 
 MQuadrangle* Rec2DAction::_createQuad(std::vector<Rec2DEdge*> &boundEdge) const
 {
+  //Msg::Info("Imhere");
   if (boundEdge.size() != 4) {
     Msg::Error("[Rec2DAction] Not 4 edges for quad creation");
     return new MQuadrangle(NULL, NULL, NULL, NULL);
@@ -717,16 +749,17 @@ MQuadrangle* Rec2DAction::_createQuad(std::vector<Rec2DEdge*> &boundEdge) const
 
 /**  Rec2DTwoTri2Quad  **/
 /************************/
-Rec2DTwoTri2Quad::Rec2DTwoTri2Quad(Rec2DEdge *re, std::list<MTriangle*> &tri)
+Rec2DTwoTri2Quad::Rec2DTwoTri2Quad(Rec2DEdge *re, std::vector<MTriangle*> &tri)
 {
   int i;
   std::set<MEdge, Less_Edge> extEdges;
   {
-    std::list<MTriangle*>::iterator it = tri.begin();
+    std::vector<MTriangle*>::iterator it = tri.begin();
     for (i = 0; it != tri.end() && i < 2; ++it, ++i) {
+      _rtriangles[i] = Recombine2D::getRElement(*it);
+      _rtriangles[i]->add(this);
       for (int j = 0; j < 3; ++j)
         extEdges.insert((*it)->getEdge(j));
-      _triangles[i] = *it;
     }
     if (it != tri.end() || i < 2)
       Msg::Error("[Rec2DTwoTri2Quad] Wrong number of triangles");
@@ -742,34 +775,39 @@ Rec2DTwoTri2Quad::Rec2DTwoTri2Quad(Rec2DEdge *re, std::list<MTriangle*> &tri)
   
   _vertices[0] = re->getRVertex(0);
   _vertices[1] = re->getRVertex(1);
-  MVertex *v2, *v3;
-  v2 = _triangles[0]->getOtherVertex(_vertices[0]->getMVertex(),
-                                    _vertices[1]->getMVertex());
-  v3 = _triangles[1]->getOtherVertex(_vertices[0]->getMVertex(),
-                                    _vertices[1]->getMVertex());
-  _vertices[2] = Recombine2D::getRVertex(v2);
-  _vertices[3] = Recombine2D::getRVertex(v3);
+  _vertices[2] = _rtriangles[0]->getOtherVertex(_vertices[0], _vertices[1]);
+  _vertices[3] = _rtriangles[1]->getOtherVertex(_vertices[0], _vertices[1]);
   
   _computeGlobVal();
 }
 
+Rec2DTwoTri2Quad::~Rec2DTwoTri2Quad()
+{
+  _rtriangles[0]->remove(this);
+  _rtriangles[1]->remove(this);
+}
+
 bool Rec2DTwoTri2Quad::isObsolete()
 {
-  int p0 = _vertices[0]->getParity();
-  int p1 = _vertices[1]->getParity();
-  int p2 = _vertices[2]->getParity();
-  int p3 = _vertices[3]->getParity();
-  if (p0>-1 && p1>-1 && p0/2 == p1/2 && p0 % 2 != p1 % 2 ||
-      p2>-1 && p3>-1 && p2/2 == p3/2 && p2 % 2 != p3 % 2 ||
-      p0>-1 && p0 == p2                                  ||
-      p0>-1 && p0 == p3                                  ||
-      p1>-1 && p1 == p2                                  ||
-      p1>-1 && p1 == p3                                    )
+  int p[4];
+  p[0] = _vertices[0]->getParity();
+  p[1] = _vertices[1]->getParity();
+  p[2] = _vertices[2]->getParity();
+  p[3] = _vertices[3]->getParity();
+  return Rec2DTwoTri2Quad::isObsolete(p);
+}
+
+bool Rec2DTwoTri2Quad::isObsolete(int *p)
+{
+  if (p[0]>-1 && p[1]>-1 && p[0]/2 == p[1]/2 && p[0]%2 != p[1]%2 ||
+      p[2]>-1 && p[3]>-1 && p[2]/2 == p[3]/2 && p[2]%2 != p[3]%2 ||
+      p[0]>-1 && (p[0] == p[2] || p[0] == p[3])                  ||
+      p[1]>-1 && (p[1] == p[2] || p[1] == p[3])                    )
     return true;
   return false;
 }
 
-void Rec2DTwoTri2Quad::apply()
+void Rec2DTwoTri2Quad::apply() //to check
 {
   if (isObsolete()) {
     Msg::Error("[Rec2DTwoTri2Quad] Applying obsolet action");
@@ -778,8 +816,8 @@ void Rec2DTwoTri2Quad::apply()
   
   int min = 100, index = -1;
   for (int i = 0; i < 4; ++i) {
-    if (_vertices[0]->getParity() > -1 && min > _vertices[0]->getParity()) {
-      min = _vertices[0]->getParity();
+    if (_vertices[i]->getParity() > -1 && min > _vertices[i]->getParity()) {
+      min = _vertices[i]->getParity();
       index = i;
     }
   }
@@ -789,21 +827,23 @@ void Rec2DTwoTri2Quad::apply()
     _vertices[1]->setParity(par);
     _vertices[2]->setParity(par+1);
     _vertices[3]->setParity(par+1);
+    Recombine2D::addVertAllQuad(_vertices[0]);
+    Recombine2D::addVertAllQuad(_vertices[1]);
+    Recombine2D::addVertAllQuad(_vertices[2]);
+    Recombine2D::addVertAllQuad(_vertices[3]);
   }
-  else { 
+  else {
     for (int i = 0; i < 4; i += 2) {
       int par;
       if ((index/2) * 2 == i)
         par = min;
-      else if (min % 2)
-        par = min - 1;
       else
-        par = min + 1;
-      Msg::Info("%d", par);
+        par = otherParity(min);
       for (int j = 0; j < 2; ++j) {
-        Msg::Info(" %d %d", _vertices[i+j]->getParity(), par);
-        if (_vertices[i+j]->getParity() == -1)
+        if (_vertices[i+j]->getParity() == -1) {
           _vertices[i+j]->setParity(par);
+          Recombine2D::addVertAllQuad(_vertices[i+j]);
+        }
         else if (_vertices[i+j]->getParity() != par)
           Recombine2D::associateParity(_vertices[0]->getParity(), par);
       }
@@ -818,8 +858,8 @@ void Rec2DTwoTri2Quad::apply()
   Recombine2D::addValEdge(valEdge);
   //Recombine2D::addValVert(_vertices[0]->getGain(-1) + _vertices[1]->getGain(-1));
   
-  Recombine2D::remove(_triangles[0]); // FIXME not performant at all !!!
-  Recombine2D::remove(_triangles[1]);
+  Recombine2D::remove(_rtriangles[0]); // FIXME not performant at all !!!
+  Recombine2D::remove(_rtriangles[1]);
   std::vector<Rec2DEdge*> edges;
   for (int i = 1; i < 5; ++i)
     edges.push_back(_edges[i]);
@@ -855,414 +895,479 @@ void Rec2DTwoTri2Quad::_computeGlobVal()
 void Rec2DTwoTri2Quad::color(int a, int b, int c)
 {
   unsigned int col = CTX::instance()->packColor(a, b, c, 255);
-  _triangles[0]->setCol(col);
-  _triangles[1]->setCol(col);
+  _rtriangles[0]->getMElement()->setCol(col);
+  _rtriangles[1]->getMElement()->setCol(col);
 }
 
 
-/**  Rec2DFourTri2Quad  **/
-/*************************/
-Rec2DFourTri2Quad::Rec2DFourTri2Quad(Rec2DVertex *rv, std::list<MTriangle*> &tri)
+bool Rec2DTwoTri2Quad::whatImply(
+                        Rec2DVertex **rv, int *par,
+                        std::map<Rec2DVertex*, std::vector<int> > &suggestions)
 {
-  int i, j;
-  std::set<MEdge, Less_Edge> edges;
-  {
-    std::list<MTriangle*>::iterator it = tri.begin();
-    for (i = 0; it != tri.end() && i < 4; ++it, ++i) {
-      for (j = 0; j < 3; ++j)
-        edges.insert((*it)->getEdge(j));
-      _triangles[i] = *it;
+  int num = 0, p[4];
+  for (int i = 0; i < 4; ++i) {
+    bool b = false;
+    for (int j = 0; j < 3; ++j) {
+      if (_vertices[i] == rv[j]) {
+        ++num;
+        p[i] = par[j];
+        b = true;
+      }
     }
-    if (it != tri.end() || i < 4)
-      Msg::Error("[Rec2DFourTri2Quad] Wrong number of triangles");
+    if (!b)
+      p[i] = _vertices[i]->getParity();
   }
+  if (num != 3)
+    Msg::Error("[Rec2DTwoTri2Quad] %d corresponding vertices instead of 3", num);
   
-  std::set<MEdge>::iterator ite = edges.begin();
-  for (i = 0, j = 4; ite != edges.end() && (i < 4 || j < 8); ++ite) {
-    if ((*ite).getVertex(0) == rv->getMVertex() ||
-        (*ite).getVertex(1) == rv->getMVertex()   )
-      _edges[i++] = Recombine2D::getREdge(*ite);
-    else if (j < 8)
-      _edges[j++] = Recombine2D::getREdge(*ite);
-    else
-      Msg::Error("[Rec2DFourTri2Quad] Too much exterior edges");
+  if (Rec2DTwoTri2Quad::isObsolete(p)) {
+    //Msg::Info("  obsolete action");
+    return false;
   }
-  if (edges.size() > 8 || ite != edges.end() || i < 4 || i > 4 || j < 8)
-    Msg::Error("[Rec2DFourTri2Quad] Wrong number of edges");
   
-  _vertices[4] = rv;
-  // the 4 other must be in order : 2 non adjacent + last 2
-  _vertices[1] = _edges[4]->getRVertex(0);
-  _vertices[3] = _edges[4]->getRVertex(1);
-  for (int i = 5; i < 8; ++i) {
-    if (_edges[i]->getRVertex(0) == _vertices[1])
-      _vertices[0] = _edges[i]->getRVertex(1);
-    if (_edges[i]->getRVertex(1) == _vertices[1])
-      _vertices[0] = _edges[i]->getRVertex(0);
-    if (_edges[i]->getRVertex(0) == _vertices[3])
-      _vertices[2] = _edges[i]->getRVertex(1);
-    if (_edges[i]->getRVertex(1) == _vertices[3])
-      _vertices[2] = _edges[i]->getRVertex(0);
-  }
-  
-  _computeGlobVal();
-}
-
-bool Rec2DFourTri2Quad::isObsolete()
-{
-  int p0 = _vertices[0]->getParity();
-  int p1 = _vertices[1]->getParity();
-  int p2 = _vertices[2]->getParity();
-  int p3 = _vertices[3]->getParity();
-  if (p0>-1 && p1>-1 && p0/2 == p1/2 && p0 % 2 != p1 % 2 ||
-      p2>-1 && p3>-1 && p2/2 == p3/2 && p2 % 2 != p3 % 2 ||
-      p0>-1 && p0 == p2                                  ||
-      p0>-1 && p0 == p3                                  ||
-      p1>-1 && p1 == p2                                  ||
-      p1>-1 && p1 == p3                                    )
-    return true;
-  return false;
-}
-
-void Rec2DFourTri2Quad::apply()
-{
-  double valEdge = 0;
-  for (int i = 0; i < 4; ++i)
-    valEdge -= REC2D_EDGE_BASE * _edges[i]->getQual();
-  for (int i = 4; i < 8; ++i)
-    valEdge += REC2D_EDGE_QUAD * _edges[i]->getQual();
-    
-  double valVert = - _vertices[4]->getQual();
-  for (int i = 0; i < 4; ++i)
-    valEdge += _vertices[i]->getGain(-1);
-  
-  Recombine2D::addNumEdge(4*REC2D_EDGE_QUAD - 4*REC2D_EDGE_BASE);
-  Recombine2D::addNumVert(-1);
-  Recombine2D::addValEdge(valEdge);
-  Recombine2D::addValVert(valVert);
-  
-  Recombine2D::remove(_triangles[0]); // FIXME not performant at all !!!
-  Recombine2D::remove(_triangles[1]);
-  Recombine2D::remove(_triangles[2]);
-  Recombine2D::remove(_triangles[3]);
-  std::vector<Rec2DEdge*> edges;
-  for (int i = 4; i < 8; ++i)
-    edges.push_back(_edges[i]);
-  Recombine2D::add(_createQuad(edges));
+  //Msg::Info("  [%d %d %d %d]", p[0], p[1], p[2], p[3]);
   
+  int min = 100, index = -1;
   for (int i = 0; i < 4; ++i) {
-    _edges[i]->getRVertex(0)->remove(_edges[i]);
-    _edges[i]->getRVertex(1)->remove(_edges[i]);
-    Recombine2D::remove(_edges[i]->getMEdge());
-    delete _edges[i];
-  }
-  Recombine2D::remove(_vertices[4]->getMVertex());
-  delete _vertices[4];
-}
-
-void Rec2DFourTri2Quad::getParities(int *par)
-{
-  par[0] = _vertices[0]->getParity();
-  par[1] = _vertices[1]->getParity();
-  par[2] = _vertices[2]->getParity();
-  par[3] = _vertices[3]->getParity();
-}
-
-void Rec2DFourTri2Quad::_computeGlobVal()
-{
-  double valEdge = 0;
-  for (int i = 0; i < 4; ++i)
-    valEdge -= REC2D_EDGE_BASE * _edges[i]->getQual();
-  for (int i = 4; i < 8; ++i)
-    valEdge += REC2D_EDGE_QUAD * _edges[i]->getQual();
-    
-  double valVert = - _vertices[4]->getQual();
-  for (int i = 0; i < 4; ++i)
-    valEdge += _vertices[i]->getGain(-1);
-  
-  _globValIfExecuted = Recombine2D::getGlobalValue(4*REC2D_EDGE_QUAD - 4*REC2D_EDGE_BASE,
-                                                   valEdge, -1, valVert      );
-  _lastUpdate = Recombine2D::getNumChange();
-}
-
-void Rec2DFourTri2Quad::color(int a, int b, int c)
-{
-  unsigned int col = CTX::instance()->packColor(a, b, c, 255);
-  _triangles[0]->setCol(col);
-  _triangles[1]->setCol(col);
-  _triangles[2]->setCol(col);
-  _triangles[3]->setCol(col); 
-}
-
-
-/**  Rec2DCollapseTri  **/
-/*********************/
-Rec2DCollapseTri::Rec2DCollapseTri(Rec2DVertex *rv, std::list<MTriangle*> &tri)
-{
-  int i;
-  std::set<MEdge, Less_Edge> extEdges;
-  {
-    std::list<MTriangle*>::iterator it = tri.begin();
-    for (i = 0; it != tri.end() && i < 4; ++it, ++i) {
-      for (int j = 0; j < 3; ++j)
-      if ((*it)->getEdge(j).getVertex(0) != rv->getMVertex() &&
-          (*it)->getEdge(j).getVertex(1) != rv->getMVertex()   )
-        extEdges.insert((*it)->getEdge(j));
+    if (p[i] > -1 && min > p[i]) {
+      min = p[i];
+      index = i;
     }
-    if (it != tri.end() || i < 4)
-      Msg::Error("[Rec2DFourTri2Quad] Wrong number of triangles");
   }
-  
-  _vertices[4] = rv;
-  // the 4 other must be in order : 2 non adjacent + last 2
-  std::set<MEdge>::iterator ite = extEdges.begin();
-  Rec2DEdge *re = Recombine2D::getREdge(*ite);
-  _vertices[1] = re->getRVertex(0);
-  _vertices[3] = re->getRVertex(1);
-  ++ite;
-  for (; ite != extEdges.end(); ++ite) {
-    re = Recombine2D::getREdge(*ite);
-    if (re->getRVertex(0) == _vertices[1])
-      _vertices[0] = re->getRVertex(1);
-    if (re->getRVertex(1) == _vertices[1])
-      _vertices[0] = re->getRVertex(0);
-    if (re->getRVertex(0) == _vertices[3])
-      _vertices[2] = re->getRVertex(1);
-    if (re->getRVertex(1) == _vertices[3])
-      _vertices[2] = re->getRVertex(0);
-  }
-  
-  _computeGlobVal();
-}
-
-Rec2DCollapseTri::Rec2DCollapseTri(Rec2DEdge *re, std::list<MTriangle*> &tri)
-: _edge(re)
-{
-  int i;
-  std::set<MEdge, Less_Edge> extEdges;
-  {
-    std::list<MTriangle*>::iterator it = tri.begin();
-    for (i = 0; it != tri.end() && i < 2; ++it, ++i) {
-      for (int j = 0; j < 3; ++j)
-        extEdges.insert((*it)->getEdge(j));
-      _triangles[i] = *it;
-    }
-    if (it != tri.end() || i < 2)
-      Msg::Error("[Rec2DTwoTri2Quad] Wrong number of triangles");
-    extEdges.erase(re->getMEdge());
-  }
-  _triangles[2] = _triangles[3] = NULL;
-  
-  _vertices[0] = re->getRVertex(0);
-  _vertices[1] = re->getRVertex(1);
-  MVertex *v2, *v3;
-  v2 = _triangles[0]->getOtherVertex(_vertices[0]->getMVertex(),
-                                    _vertices[1]->getMVertex());
-  v3 = _triangles[1]->getOtherVertex(_vertices[0]->getMVertex(),
-                                    _vertices[1]->getMVertex());
-  _vertices[2] = Recombine2D::getRVertex(v2);
-  _vertices[3] = Recombine2D::getRVertex(v3);
-  _vertices[4] = NULL;
-  
-  _computeGlobVal();
-}
-
-bool Rec2DCollapseTri::isObsolete()
-{
-  int p0 = _vertices[0]->getParity();
-  int p1 = _vertices[1]->getParity();
-  int p2 = _vertices[2]->getParity();
-  int p3 = _vertices[3]->getParity();
-  int b0 = _vertices[0]->getIsOnBoundary();
-  int b1 = _vertices[1]->getIsOnBoundary();
-  int b2 = _vertices[2]->getIsOnBoundary();
-  int b3 = _vertices[3]->getIsOnBoundary();
-  if ((b0 && b1) || (p0>-1 && p1>-1 && p0/2 == p1/2 && p0 % 2 != p1 % 2) &&
-      (b2 && b3) || (p2>-1 && p3>-1 && p2/2 == p3/2 && p2 % 2 != p3 % 2)   )
-    return true;
-  return false;
-}
-
-void Rec2DCollapseTri::apply()
-{
-  return;
-}
-
-void Rec2DCollapseTri::getParities(int *par)
-{
-  par[0] = _vertices[0]->getParity();
-  par[1] = _vertices[1]->getParity();
-  par[2] = _vertices[2]->getParity();
-  par[3] = _vertices[3]->getParity();
-}
-
-void Rec2DCollapseTri::_computeGlobVal()
-{
-  int b0 = _vertices[0]->getIsOnBoundary();
-  int b1 = _vertices[1]->getIsOnBoundary();
-  int b2 = _vertices[2]->getIsOnBoundary();
-  int b3 = _vertices[3]->getIsOnBoundary();
-  
-  std::map<Rec2DVertex*, int> neighbourVert;
-  std::set<Rec2DEdge*> neighbourEdge; 
-  if (!b0 || !b1) {
-    _vertices[0]->getNeighbours(neighbourVert, neighbourEdge);
-    _vertices[1]->getNeighbours(neighbourVert, neighbourEdge);
-    if (_vertices[4]) {
-      neighbourVert.erase(_vertices[4]);
-      _vertices[4]->getNeighbourEdges(neighbourEdge);
-    }
-    //else
-    //  neighbourEdge.insert(_edge);
-    
-    int numEdge = 0, numVert = 0;
-    double valEdge = .0, valVert = .0;
-    {
-      std::set<Rec2DEdge*>::iterator it = neighbourEdge.begin();
-      for (; it != neighbourEdge.end(); ++it) {
-        valEdge -= (*it)->getWeightedQual();
-        numEdge -= (*it)->getWeight();
-      }
-      valVert -= _vertices[0]->getQual();
-      valVert -= _vertices[1]->getQual();
-      if (_vertices[4]) {
-        valVert -= _vertices[4]->getQual();
-        numVert -= 1;
-        valVert += _vertices[2]->getGain(-2);
-        valVert += _vertices[3]->getGain(-2);
-      }
-      else {
-        valVert += _vertices[2]->getGain(-1);
-        valVert += _vertices[3]->getGain(-1);
-      }
-    }
-    if (b0)
-      _qualCavity(valVert, valEdge, numEdge, neighbourVert, _vertices[0]);
-    else if (b1)
-      _qualCavity(valVert, valEdge, numEdge, neighbourVert, _vertices[1]);
-    else
-      _qualCavity(valVert, valEdge, numEdge, neighbourVert);
-    numVert -= 1;
-    
-    _globValIfExecuted =
-        Recombine2D::getGlobalValue(numEdge, valEdge, numVert, valVert);
+  if (index == -1) {
+    Msg::Error("[Rec2DTwoTri2Quad] no parities");
+    return false;
   }
-  else
-    _globValIfExecuted = -1.;
-  
   
-  neighbourVert.clear();
-  neighbourEdge.clear();
-  if (!b2 || !b3) {
-    _vertices[2]->getNeighbours(neighbourVert, neighbourEdge);
-    _vertices[3]->getNeighbours(neighbourVert, neighbourEdge);
-    if (_vertices[4]) {
-      neighbourVert.erase(_vertices[4]);
-      _vertices[4]->getNeighbourEdges(neighbourEdge);
-    }
-    else
-      neighbourEdge.insert(_edge);
-    
-    int numEdge = 0, numVert = 0;
-    double valEdge = .0, valVert = .0;
-    {
-      std::set<Rec2DEdge*>::iterator it = neighbourEdge.begin();
-      for (; it != neighbourEdge.end(); ++it) {
-        valEdge -= (*it)->getWeightedQual();
-        numEdge -= (*it)->getWeight();
-      }
-      valVert -= _vertices[2]->getQual();
-      valVert -= _vertices[3]->getQual();
-      if (_vertices[4]) {
-        valVert -= _vertices[4]->getQual();
-        numVert -= 1;
-        valVert += _vertices[0]->getGain(-2);
-        valVert += _vertices[1]->getGain(-2);
-      }
-      else {
-        valVert += _vertices[0]->getGain(-2);
-        valVert += _vertices[1]->getGain(-2);
-      }
-    }
-    if (b2)
-      _qualCavity(valVert, valEdge, numEdge, neighbourVert, _vertices[2]);
-    else if (b3)
-      _qualCavity(valVert, valEdge, numEdge, neighbourVert, _vertices[3]);
+  for (int i = 0; i < 4; i += 2) {
+    int par;
+    if ((index/2)*2 == i)
+      par = min;
     else
-      _qualCavity(valVert, valEdge, numEdge, neighbourVert);
-    numVert -= 1;
-    
-    _globValIfExecuted2 =
-        Recombine2D::getGlobalValue(numEdge, valEdge, numVert, valVert);
-  }
-  else
-    _globValIfExecuted2 = -1.;
-    
-  _lastUpdate = Recombine2D::getNumChange();
-}
-
-void Rec2DCollapseTri::_qualCavity(double &valVert, double &valEdge,
-                                   int &numEdge,
-                                   std::map<Rec2DVertex*, int> &vert,
-                                   Rec2DVertex *imposedVert          )
-{
-  Recombine2D::bo = true;
-  std::map<Rec2DVertex*, int>::iterator it;
-  Rec2DVertex *rv = imposedVert;
-  if (rv == NULL) {
-    double u, v = u = .0;
-    it = vert.begin();
-    for (; it != vert.end(); ++it) {
-      u += it->first->u();
-      v += it->first->v();
+      par = otherParity(min);
+    for (int j = 0; j < 2; ++j) {
+      if (p[i+j] != par)
+        suggestions[_vertices[i+j]].push_back(par);
     }
-    u /= vert.size();
-    v /= vert.size();
-    
-    rv = new Rec2DVertex(u, v);
-  }
-  valVert = rv->getQual(vert.size());
-  it = vert.begin();
-  for (; it != vert.end(); ++it) {
-    Rec2DEdge re(rv, it->first);
-    valEdge += it->second * re.getQual();
-    numEdge += it->second;
-  }
-  if (imposedVert == NULL) {
-    rv->deleteMVertex();
-    delete rv;
-  }
-}
-
-void Rec2DCollapseTri::color(int a, int b, int c)
-{
-  unsigned int col = CTX::instance()->packColor(a, b, c, 255);
-  _triangles[0]->setCol(col);
-  _triangles[1]->setCol(col);
-  if (_triangles[2]) {
-    _triangles[2]->setCol(col);
-    _triangles[3]->setCol(col);
   }
+  //Msg::Info("  ok");
+  return true;
 }
 
 
+///**  Rec2DFourTri2Quad  **/
+///*************************/
+//Rec2DFourTri2Quad::Rec2DFourTri2Quad(Rec2DVertex *rv, std::vector<MTriangle*> &tri)
+//{
+//  int i, j;
+//  std::set<MEdge, Less_Edge> edges;
+//  {
+//    std::vector<MTriangle*>::iterator it = tri.begin();
+//    for (i = 0; it != tri.end() && i < 4; ++it, ++i) {
+//      for (j = 0; j < 3; ++j)
+//        edges.insert((*it)->getEdge(j));
+//      _rtriangles[i] = Recombine2D::getRElement(*it);
+//      _rtriangles[i]->add(this);
+//    }
+//    if (it != tri.end() || i < 4)
+//      Msg::Error("[Rec2DFourTri2Quad] Wrong number of triangles");
+//  }
+//  
+//  std::set<MEdge>::iterator ite = edges.begin();
+//  for (i = 0, j = 4; ite != edges.end() && (i < 4 || j < 8); ++ite) {
+//    if ((*ite).getVertex(0) == rv->getMVertex() ||
+//        (*ite).getVertex(1) == rv->getMVertex()   )
+//      _edges[i++] = Recombine2D::getREdge(*ite);
+//    else if (j < 8)
+//      _edges[j++] = Recombine2D::getREdge(*ite);
+//    else
+//      Msg::Error("[Rec2DFourTri2Quad] Too much exterior edges");
+//  }
+//  if (edges.size() > 8 || ite != edges.end() || i < 4 || i > 4 || j < 8)
+//    Msg::Error("[Rec2DFourTri2Quad] Wrong number of edges");
+//  
+//  _vertices[4] = rv;
+//  // the 4 other must be in order : 2 non adjacent + last 2
+//  _vertices[1] = _edges[4]->getRVertex(0);
+//  _vertices[3] = _edges[4]->getRVertex(1);
+//  for (int i = 5; i < 8; ++i) {
+//    if (_edges[i]->getRVertex(0) == _vertices[1])
+//      _vertices[0] = _edges[i]->getRVertex(1);
+//    if (_edges[i]->getRVertex(1) == _vertices[1])
+//      _vertices[0] = _edges[i]->getRVertex(0);
+//    if (_edges[i]->getRVertex(0) == _vertices[3])
+//      _vertices[2] = _edges[i]->getRVertex(1);
+//    if (_edges[i]->getRVertex(1) == _vertices[3])
+//      _vertices[2] = _edges[i]->getRVertex(0);
+//  }
+//  
+//  _computeGlobVal();
+//}
+//
+//Rec2DFourTri2Quad::~Rec2DFourTri2Quad()
+//{
+//  _rtriangles[0]->remove(this);
+//  _rtriangles[1]->remove(this);
+//  _rtriangles[2]->remove(this);
+//  _rtriangles[3]->remove(this);
+//}
+//
+//bool Rec2DFourTri2Quad::isObsolete()
+//{
+//  int p[4];
+//  p[0] = _vertices[0]->getParity();
+//  p[1] = _vertices[1]->getParity();
+//  p[2] = _vertices[2]->getParity();
+//  p[3] = _vertices[3]->getParity();
+//  return Rec2DTwoTri2Quad::isObsolete(p);
+//}
+//
+//void Rec2DFourTri2Quad::apply()
+//{
+//  double valEdge = 0;
+//  for (int i = 0; i < 4; ++i)
+//    valEdge -= REC2D_EDGE_BASE * _edges[i]->getQual();
+//  for (int i = 4; i < 8; ++i)
+//    valEdge += REC2D_EDGE_QUAD * _edges[i]->getQual();
+//    
+//  double valVert = - _vertices[4]->getQual();
+//  for (int i = 0; i < 4; ++i)
+//    valEdge += _vertices[i]->getGain(-1);
+//  
+//  Recombine2D::addNumEdge(4*REC2D_EDGE_QUAD - 4*REC2D_EDGE_BASE);
+//  Recombine2D::addNumVert(-1);
+//  Recombine2D::addValEdge(valEdge);
+//  Recombine2D::addValVert(valVert);
+//  
+//  Recombine2D::remove(_rtriangles[0]); // FIXME not performant at all !!!
+//  Recombine2D::remove(_rtriangles[1]);
+//  Recombine2D::remove(_rtriangles[2]);
+//  Recombine2D::remove(_rtriangles[3]);
+//  std::vector<Rec2DEdge*> edges;
+//  for (int i = 4; i < 8; ++i)
+//    edges.push_back(_edges[i]);
+//  Recombine2D::add(_createQuad(edges));
+//  
+//  for (int i = 0; i < 4; ++i) {
+//    _edges[i]->getRVertex(0)->remove(_edges[i]);
+//    _edges[i]->getRVertex(1)->remove(_edges[i]);
+//    Recombine2D::remove(_edges[i]->getMEdge());
+//    delete _edges[i];
+//  }
+//  Recombine2D::remove(_vertices[4]->getMVertex());
+//  delete _vertices[4];
+//}
+//
+//void Rec2DFourTri2Quad::getParities(int *par)
+//{
+//  par[0] = _vertices[0]->getParity();
+//  par[1] = _vertices[1]->getParity();
+//  par[2] = _vertices[2]->getParity();
+//  par[3] = _vertices[3]->getParity();
+//}
+//
+//void Rec2DFourTri2Quad::_computeGlobVal()
+//{
+//  double valEdge = 0;
+//  for (int i = 0; i < 4; ++i)
+//    valEdge -= REC2D_EDGE_BASE * _edges[i]->getQual();
+//  for (int i = 4; i < 8; ++i)
+//    valEdge += REC2D_EDGE_QUAD * _edges[i]->getQual();
+//    
+//  double valVert = - _vertices[4]->getQual();
+//  for (int i = 0; i < 4; ++i)
+//    valEdge += _vertices[i]->getGain(-1);
+//  
+//  _globValIfExecuted = Recombine2D::getGlobalValue(4*REC2D_EDGE_QUAD - 4*REC2D_EDGE_BASE,
+//                                                   valEdge, -1, valVert      );
+//  _lastUpdate = Recombine2D::getNumChange();
+//}
+//
+//void Rec2DFourTri2Quad::color(int a, int b, int c)
+//{
+//  unsigned int col = CTX::instance()->packColor(a, b, c, 255);
+//  _rtriangles[0]->getMElement()->setCol(col);
+//  _rtriangles[1]->getMElement()->setCol(col);
+//  _rtriangles[2]->getMElement()->setCol(col);
+//  _rtriangles[3]->getMElement()->setCol(col); 
+//}
+//
+//
+///**  Rec2DCollapseTri  **/
+///*********************/
+//Rec2DCollapseTri::Rec2DCollapseTri(Rec2DVertex *rv, std::vector<MTriangle*> &tri)
+//{
+//  int i;
+//  std::set<MEdge, Less_Edge> extEdges;
+//  {
+//    std::vector<MTriangle*>::iterator it = tri.begin();
+//    for (i = 0; it != tri.end() && i < 4; ++it, ++i) {
+//      for (int j = 0; j < 3; ++j) {
+//        if ((*it)->getEdge(j).getVertex(0) != rv->getMVertex() &&
+//            (*it)->getEdge(j).getVertex(1) != rv->getMVertex()   )
+//          extEdges.insert((*it)->getEdge(j));
+//      }
+//      _rtriangles[i] = Recombine2D::getRElement(*it);
+//      _rtriangles[i]->add(this);
+//    }
+//    if (it != tri.end() || i < 4)
+//      Msg::Error("[Rec2DFourTri2Quad] Wrong number of triangles");
+//  }
+//  
+//  _vertices[4] = rv;
+//  // the 4 other must be in order : 2 non adjacent + last 2
+//  std::set<MEdge>::iterator ite = extEdges.begin();
+//  Rec2DEdge *re = Recombine2D::getREdge(*ite);
+//  _vertices[1] = re->getRVertex(0);
+//  _vertices[3] = re->getRVertex(1);
+//  ++ite;
+//  for (; ite != extEdges.end(); ++ite) {
+//    re = Recombine2D::getREdge(*ite);
+//    if (re->getRVertex(0) == _vertices[1])
+//      _vertices[0] = re->getRVertex(1);
+//    if (re->getRVertex(1) == _vertices[1])
+//      _vertices[0] = re->getRVertex(0);
+//    if (re->getRVertex(0) == _vertices[3])
+//      _vertices[2] = re->getRVertex(1);
+//    if (re->getRVertex(1) == _vertices[3])
+//      _vertices[2] = re->getRVertex(0);
+//  }
+//  
+//  _computeGlobVal();
+//}
+//
+//Rec2DCollapseTri::Rec2DCollapseTri(Rec2DEdge *re, std::vector<MTriangle*> &tri)
+//: _edge(re)
+//{
+//  int i;
+//  std::set<MEdge, Less_Edge> extEdges;
+//  {
+//    std::vector<MTriangle*>::iterator it = tri.begin();
+//    for (i = 0; it != tri.end() && i < 2; ++it, ++i) {
+//      for (int j = 0; j < 3; ++j)
+//        extEdges.insert((*it)->getEdge(j));
+//      _rtriangles[i] = Recombine2D::getRElement(*it);
+//      _rtriangles[i]->add(this);
+//    }
+//    if (it != tri.end() || i < 2)
+//      Msg::Error("[Rec2DTwoTri2Quad] Wrong number of triangles");
+//    extEdges.erase(re->getMEdge());
+//  }
+//  _rtriangles[2] = _rtriangles[3] = NULL;
+//  
+//  _vertices[0] = re->getRVertex(0);
+//  _vertices[1] = re->getRVertex(1);
+//  _vertices[2] = _rtriangles[0]->getOtherVertex(_vertices[0], _vertices[1]);
+//  _vertices[3] = _rtriangles[1]->getOtherVertex(_vertices[0], _vertices[1]);
+//  _vertices[4] = NULL;
+//  
+//  _computeGlobVal();
+//}
+//
+//Rec2DCollapseTri::~Rec2DCollapseTri()
+//{
+//  //_rtriangles[0]->remove(this);
+//  //_rtriangles[1]->remove(this);
+//}
+//
+//bool Rec2DCollapseTri::isObsolete()
+//{
+//  int p0 = _vertices[0]->getParity();
+//  int p1 = _vertices[1]->getParity();
+//  int p2 = _vertices[2]->getParity();
+//  int p3 = _vertices[3]->getParity();
+//  int b0 = _vertices[0]->getIsOnBoundary();
+//  int b1 = _vertices[1]->getIsOnBoundary();
+//  int b2 = _vertices[2]->getIsOnBoundary();
+//  int b3 = _vertices[3]->getIsOnBoundary();
+//  if ((b0 && b1) || (p0>-1 && p1>-1 && p0/2 == p1/2 && p0 % 2 != p1 % 2) &&
+//      (b2 && b3) || (p2>-1 && p3>-1 && p2/2 == p3/2 && p2 % 2 != p3 % 2)   )
+//    return true;
+//  return false;
+//}
+//
+//void Rec2DCollapseTri::apply()
+//{
+//  return;
+//}
+//
+//void Rec2DCollapseTri::getParities(int *par)
+//{
+//  par[0] = _vertices[0]->getParity();
+//  par[1] = _vertices[1]->getParity();
+//  par[2] = _vertices[2]->getParity();
+//  par[3] = _vertices[3]->getParity();
+//}
+//
+//void Rec2DCollapseTri::_computeGlobVal()
+//{
+//  int b0 = _vertices[0]->getIsOnBoundary();
+//  int b1 = _vertices[1]->getIsOnBoundary();
+//  int b2 = _vertices[2]->getIsOnBoundary();
+//  int b3 = _vertices[3]->getIsOnBoundary();
+//  
+//  std::map<Rec2DVertex*, int> neighbourVert;
+//  std::set<Rec2DEdge*> neighbourEdge; 
+//  if (!b0 || !b1) {
+//    _vertices[0]->getNeighbours(neighbourVert, neighbourEdge);
+//    _vertices[1]->getNeighbours(neighbourVert, neighbourEdge);
+//    if (_vertices[4]) {
+//      neighbourVert.erase(_vertices[4]);
+//      _vertices[4]->getNeighbourEdges(neighbourEdge);
+//    }
+//    //else
+//    //  neighbourEdge.insert(_edge);
+//    
+//    int numEdge = 0, numVert = 0;
+//    double valEdge = .0, valVert = .0;
+//    {
+//      std::set<Rec2DEdge*>::iterator it = neighbourEdge.begin();
+//      for (; it != neighbourEdge.end(); ++it) {
+//        valEdge -= (*it)->getWeightedQual();
+//        numEdge -= (*it)->getWeight();
+//      }
+//      valVert -= _vertices[0]->getQual();
+//      valVert -= _vertices[1]->getQual();
+//      if (_vertices[4]) {
+//        valVert -= _vertices[4]->getQual();
+//        numVert -= 1;
+//        valVert += _vertices[2]->getGain(-2);
+//        valVert += _vertices[3]->getGain(-2);
+//      }
+//      else {
+//        valVert += _vertices[2]->getGain(-1);
+//        valVert += _vertices[3]->getGain(-1);
+//      }
+//    }
+//    if (b0)
+//      _qualCavity(valVert, valEdge, numEdge, neighbourVert, _vertices[0]);
+//    else if (b1)
+//      _qualCavity(valVert, valEdge, numEdge, neighbourVert, _vertices[1]);
+//    else
+//      _qualCavity(valVert, valEdge, numEdge, neighbourVert);
+//    numVert -= 1;
+//    
+//    _globValIfExecuted =
+//        Recombine2D::getGlobalValue(numEdge, valEdge, numVert, valVert);
+//  }
+//  else
+//    _globValIfExecuted = -1.;
+//  
+//  
+//  neighbourVert.clear();
+//  neighbourEdge.clear();
+//  if (!b2 || !b3) {
+//    _vertices[2]->getNeighbours(neighbourVert, neighbourEdge);
+//    _vertices[3]->getNeighbours(neighbourVert, neighbourEdge);
+//    if (_vertices[4]) {
+//      neighbourVert.erase(_vertices[4]);
+//      _vertices[4]->getNeighbourEdges(neighbourEdge);
+//    }
+//    else
+//      neighbourEdge.insert(_edge);
+//    
+//    int numEdge = 0, numVert = 0;
+//    double valEdge = .0, valVert = .0;
+//    {
+//      std::set<Rec2DEdge*>::iterator it = neighbourEdge.begin();
+//      for (; it != neighbourEdge.end(); ++it) {
+//        valEdge -= (*it)->getWeightedQual();
+//        numEdge -= (*it)->getWeight();
+//      }
+//      valVert -= _vertices[2]->getQual();
+//      valVert -= _vertices[3]->getQual();
+//      if (_vertices[4]) {
+//        valVert -= _vertices[4]->getQual();
+//        numVert -= 1;
+//        valVert += _vertices[0]->getGain(-2);
+//        valVert += _vertices[1]->getGain(-2);
+//      }
+//      else {
+//        valVert += _vertices[0]->getGain(-2);
+//        valVert += _vertices[1]->getGain(-2);
+//      }
+//    }
+//    if (b2)
+//      _qualCavity(valVert, valEdge, numEdge, neighbourVert, _vertices[2]);
+//    else if (b3)
+//      _qualCavity(valVert, valEdge, numEdge, neighbourVert, _vertices[3]);
+//    else
+//      _qualCavity(valVert, valEdge, numEdge, neighbourVert);
+//    numVert -= 1;
+//    
+//    _globValIfExecuted2 =
+//        Recombine2D::getGlobalValue(numEdge, valEdge, numVert, valVert);
+//  }
+//  else
+//    _globValIfExecuted2 = -1.;
+//    
+//  _lastUpdate = Recombine2D::getNumChange();
+//}
+//
+//void Rec2DCollapseTri::_qualCavity(double &valVert, double &valEdge,
+//                                   int &numEdge,
+//                                   std::map<Rec2DVertex*, int> &vert,
+//                                   Rec2DVertex *imposedVert          )
+//{
+//  Recombine2D::bo = true;
+//  std::map<Rec2DVertex*, int>::iterator it;
+//  Rec2DVertex *rv = imposedVert;
+//  if (rv == NULL) {
+//    double u, v = u = .0;
+//    it = vert.begin();
+//    for (; it != vert.end(); ++it) {
+//      u += it->first->u();
+//      v += it->first->v();
+//    }
+//    u /= vert.size();
+//    v /= vert.size();
+//    
+//    rv = new Rec2DVertex(u, v);
+//  }
+//  valVert = rv->getQual(vert.size());
+//  it = vert.begin();
+//  for (; it != vert.end(); ++it) {
+//    Rec2DEdge re(rv, it->first);
+//    valEdge += it->second * re.getQual();
+//    numEdge += it->second;
+//  }
+//  if (imposedVert == NULL) {
+//    rv->deleteMVertex();
+//    delete rv;
+//  }
+//}
+//
+//void Rec2DCollapseTri::color(int a, int b, int c)
+//{
+//  unsigned int col = CTX::instance()->packColor(a, b, c, 255);
+//  _rtriangles[0]->getMElement()->setCol(col);
+//  _rtriangles[1]->getMElement()->setCol(col);
+//  if (_rtriangles[2]) {
+//    _rtriangles[2]->getMElement()->setCol(col);
+//    _rtriangles[3]->getMElement()->setCol(col);
+//  }
+//}
+//
+//
 /**  Rec2DVertex  **/
 /*******************/
-Rec2DVertex::Rec2DVertex(MVertex *v, std::list<MTriangle*> tri, int onWhat,
+Rec2DVertex::Rec2DVertex(MVertex *v, std::vector<MTriangle*> tri, int onWhat,
                          std::map<MVertex*, std::set<GEdge*> > gEdges)
 : _v(v), _onWhat(onWhat), _parity(-1), _lastChange(-1), _numEl(tri.size())
 {
+  for (unsigned int i = 0; i < tri.size(); ++i)
+    _elements.push_back(Recombine2D::getRElement(tri[i]));
+  
   if (_onWhat == -1) {
     std::map<MVertex*, std::set<GEdge*> >::iterator it = gEdges.find(_v);
-    if (it != gEdges.end()) {
+    if (it != gEdges.end())
       _angle = _computeAngle(_v, tri, it->second);
-    }
     else {
       Msg::Warning("[meshRecombine2D] Can't compute corner angle, setting angle to pi/2");
       _angle = M_PI / 2.;
     }
   }
-  reparamMeshVertexOnFace(_v, Recombine2D::getGFace(), _param); 
+  reparamMeshVertexOnFace(_v, Recombine2D::getGFace(), _param);
   _computeQual();
 }
 
@@ -1280,7 +1385,7 @@ Rec2DVertex::Rec2DVertex(double u, double v)
 }
 
 double Rec2DVertex::_computeAngle(MVertex *v,
-                                  std::list<MTriangle*> &listTri,
+                                  std::vector<MTriangle*> &tri,
                                   std::set<GEdge*> &setEdge)
 {
   if (setEdge.size() != 2) {
@@ -1319,8 +1424,8 @@ double Rec2DVertex::_computeAngle(MVertex *v,
   double angle2 = 2. * M_PI - angle1;
   
   double angleMesh = .0;
-  std::list<MTriangle*>::iterator it2 = listTri.begin();
-  for (; it2 != listTri.end(); ++it2) {
+  std::vector<MTriangle*>::iterator it2 = tri.begin();
+  for (; it2 != tri.end(); ++it2) {
     MTriangle *t = *it2;
     int k = 0;
     while (t->getVertex(k) != v && k < 3)
@@ -1357,6 +1462,14 @@ double Rec2DVertex::getQual(int numEdge)
   return _qual;
 }
 
+bool Rec2DVertex::hasTriangle()
+{
+  bool b = false;
+  for (int i = 0; i < _elements.size(); ++i)
+    if (!_elements[i]->isQuad()) b = true;
+  return b;
+}
+
 void Rec2DVertex::_computeQual()
 {
   if (_onWhat > -1)
@@ -1450,15 +1563,31 @@ void Rec2DVertex::remove(Rec2DEdge *re)
     Msg::Error("[Rec2DVertex] Edge not found");
 }
 
+void Rec2DVertex::getElements(std::set<Rec2DElement*> &elem)
+{
+  for (unsigned int i = 0; i < _elements.size(); ++i)
+    elem.insert(_elements[i]);
+}
+
 
 /**  Rec2DEdge  **/
 /*****************/
-Rec2DEdge::Rec2DEdge(MEdge e, mapofVertices &vert, std::list<MTriangle*> &tri)
+Rec2DEdge::Rec2DEdge(MEdge e, mapofVertices &vert, std::vector<MTriangle*> &tri)
 : _lastChange(-1), _qual(.0), _weight(REC2D_EDGE_BASE),
   _rv1(vert[e.getVertex(0)]), _rv2(vert[e.getVertex(1)])
 {
   _rv1->add(this);
   _rv2->add(this);
+  if (tri.size() == 2) {
+    Rec2DElement *t0 = Recombine2D::getRElement(tri[0]);
+    Rec2DElement *t1 = Recombine2D::getRElement(tri[1]);
+    t0->add(this, t1);
+    t1->add(this, t0);
+  }
+  else if (tri.size() == 1)
+    Recombine2D::getRElement(tri[0])->add(this);
+  else
+    Msg::Error("[Rec2DEdge] Wrong number of triangles");
 }
 
 Rec2DEdge::Rec2DEdge(Rec2DVertex *rv1, Rec2DVertex *rv2)
@@ -1475,6 +1604,8 @@ double Rec2DEdge::getQual()
 
 void Rec2DEdge::_computeQual() //*
 {
+  //if (!_rv1->update() && !_rv2->update())
+  //  return;
   static int i = 0; if (++i < 2) Msg::Warning("FIXME compute qual edge and update of vertex");
   //_rv1->update();
   //_rv2->update();
@@ -1482,7 +1613,7 @@ void Rec2DEdge::_computeQual() //*
   double alignment = _straightAlignment();
   if (adimLength > 1)
     adimLength = 1/adimLength;
-  _qual = Recombine2D::edgeReward(adimLength, alignment);
+  _qual = adimLength * ((1-REC2D_ALIGNMENT) + REC2D_ALIGNMENT * alignment);
   _lastChange = Recombine2D::getNumChange();
 }
 
@@ -1527,7 +1658,26 @@ double Rec2DEdge::_straightAlignment()
   return (alpha1/lc1 + alpha2/lc2) / (1/lc1 + 1/lc2);
 }
 
-bool lessRec2DAction::operator()(Rec2DAction *ra1, Rec2DAction *ra2) const
+
+/**  Rec2DElement  **/
+/********************/
+void Rec2DElement::add(Rec2DEdge *re)
 {
-  return *ra1 < *ra2;
+  int i = -1;
+  while (_redges[++i] && i < 4);
+  if (_redges[i])
+    Msg::Error("[Rec2DElement] too much edges");
+  _redges[i] = re;
 }
+
+void Rec2DElement::add(Rec2DEdge *re, Rec2DElement *el)
+{
+  int i = -1;
+  while (_redges[++i] && i < 4);
+  if (_redges[i])
+    Msg::Error("[Rec2DElement] too much edges");
+  _redges[i] = re;
+  _relements[i] = el;
+}
+
+
diff --git a/Mesh/meshGFaceRecombine.h b/Mesh/meshGFaceRecombine.h
index 868fff33f4..ebd9f872ad 100644
--- a/Mesh/meshGFaceRecombine.h
+++ b/Mesh/meshGFaceRecombine.h
@@ -15,13 +15,15 @@
 #include "MEdge.h"
 #include "BackgroundMesh.h"
 
-class Rec2DEdge;
 class Rec2DVertex;
+class Rec2DEdge;
+class Rec2DElement;
 class Rec2DAction;
 
 typedef std::list<Rec2DAction*> setofRec2DAction;
 typedef std::map<MVertex*, Rec2DVertex*> mapofVertices;
 typedef std::map<MEdge, Rec2DEdge*, Less_Edge> mapofEdges;
+typedef std::map<MElement*, Rec2DElement*> mapofElements;
 typedef std::map<MElement*, std::set<Rec2DAction*> > mapofElementActions;
 typedef std::map<MQuadrangle*, std::set<MElement*> > mapofAdjacencies;
 
@@ -35,8 +37,9 @@ class Recombine2D {
     double _valEdge, _valVert;
     GFace *_gf;
     backgroundMesh *_bgm;
-    mapofEdges _edges;
     mapofVertices _vertices;
+    mapofEdges _edges;
+    mapofElements _elements;
     setofRec2DAction _actions;
     mapofElementActions _mea;
     mapofAdjacencies _ma;
@@ -45,7 +48,8 @@ class Recombine2D {
     static Recombine2D *_current;
     
       std::vector<MTriangle*> _tri;
-      std::vector<MQuadrangle*> _quad; 
+      std::vector<MQuadrangle*> _quad;
+      std::vector<Rec2DVertex*> _checkAllQuad;
     
   public :
     static bool bo;
@@ -57,14 +61,30 @@ class Recombine2D {
     bool recombine();
     void apply(bool){return;}
     
-    static inline double edgeReward(double adimLength, double alignment)
-      {return adimLength * (.5 + .5 * alignment);}
+    static int getNumVertAllQuad() {return _current->_checkAllQuad.size();}
+    static Rec2DVertex* getVertAllQuad(int i) {return _current->_checkAllQuad[i];}
+    static void removeVertAllQuad(Rec2DVertex *rv) {
+      std::vector<Rec2DVertex*>::iterator it = _current->_checkAllQuad.begin();
+      while (it != _current->_checkAllQuad.end()) {
+        if (*it == rv)
+          it = _current->_checkAllQuad.erase(it);
+        else
+          ++it;
+      }
+    }
+    static void addVertAllQuad(Rec2DVertex *rv) {
+      _current->_checkAllQuad.push_back(rv);
+    }
+    
+    //static inline double edgeReward(double adimLength, double alignment)
+    //  {return adimLength * ((1-REC2D_ALIGNMENT) + REC2D_ALIGNMENT * alignment);}
     static inline int getNumChange() {return _current->_numChange;}
     static inline backgroundMesh* bgm() {return _current->_bgm;}
-    static inline double getGlobalValue(int numEdge, double valEdge,
-                                        int numVert, double valVert );
+    static double getGlobalValue(int numEdge, double valEdge,
+                                 int numVert, double valVert );
+    //static Rec2DVertex* getRVertex(MVertex*);
     static Rec2DEdge* getREdge(MEdge);
-    static Rec2DVertex* getRVertex(MVertex*);
+    static Rec2DElement* getRElement(MElement*);
     static void remove(MEdge);
     static void remove(MVertex*);
     static inline GFace* getGFace() {return _current->_gf;}
@@ -76,7 +96,8 @@ class Recombine2D {
     static inline void addValEdge(double a) {_current->_valEdge += a;}
     static inline void addValVert(double a) {_current->_valVert += a;}
     
-    static void remove(MTriangle*);
+    static void remove(Rec2DElement*);
+    static void remove(Rec2DAction*);
     static void add(MQuadrangle*);
     
     static inline void addVert2Par(Rec2DVertex *rv, int a)
@@ -86,7 +107,7 @@ class Recombine2D {
       std::map<int, std::set<Rec2DVertex*> >::iterator it;
       it = _current->_parities.end();
       --it;
-      return it->first;
+      return it->first + 1;
     }
     
   private :
@@ -95,6 +116,7 @@ class Recombine2D {
     void _removeIncompatible(Rec2DAction*);
     void _removeReference(Rec2DAction*);
     void _getNeighbours(Rec2DAction*, std::set<Rec2DAction*>&);
+    void _getElemToLook(int par, std::set<Rec2DElement*>&);
 };
 
 class Rec2DAction {
@@ -102,10 +124,11 @@ class Rec2DAction {
     int _position;
     
   protected :
-    int _lastUpdate;
+    int _lastUpdate; // + _lastCheck
     double _globValIfExecuted;
     
   public :
+    virtual ~Rec2DAction() {Recombine2D::remove(this);}
     bool operator<(Rec2DAction&);
     inline int getVectPosition() const {return _position;}
     inline void setVectPosition(int p) {_position = p;}
@@ -117,12 +140,14 @@ class Rec2DAction {
       return _globValIfExecuted - Recombine2D::getGlobalValue(0, .0, 0, .0);
     }
     virtual bool isObsolete() = 0;
-    virtual void getTriangles(std::set<MTriangle*>&) = 0;
-    virtual MTriangle* getTriangle(int) = 0;
+    virtual void getRTriangles(std::set<Rec2DElement*>&) = 0;
+    virtual Rec2DElement* getRTriangle(int) = 0;
     virtual void apply() = 0;
     virtual void getParities(int*) = 0;
     virtual Rec2DVertex* getRVertex(int) = 0;
     virtual void color(int, int, int) = 0;
+    virtual bool whatImply(Rec2DVertex**, int*,
+                           std::map<Rec2DVertex*, std::vector<int> >&) = 0;
   
   protected :
     MQuadrangle* _createQuad(std::vector<Rec2DEdge*>&) const;
@@ -133,52 +158,60 @@ class Rec2DAction {
 
 class Rec2DTwoTri2Quad : public Rec2DAction {
   private :
-    MTriangle *_triangles[2];
+    Rec2DElement *_rtriangles[2];
     Rec2DEdge *_edges[5]; // 1 embedded, 4 boundary
     Rec2DVertex *_vertices[4]; // 4 boundary (2 on embedded edge + 2)
     
   public :
-    Rec2DTwoTri2Quad(Rec2DEdge*, std::list<MTriangle*>&);
+    Rec2DTwoTri2Quad(Rec2DEdge*, std::vector<MTriangle*>&);
+    ~Rec2DTwoTri2Quad();
     bool isObsolete();
-    void getTriangles(std::set<MTriangle*> &tri) {
-      tri.insert(_triangles[0]);
-      tri.insert(_triangles[1]);
+    void getRTriangles(std::set<Rec2DElement*> &tri) {
+      tri.insert(_rtriangles[0]);
+      tri.insert(_rtriangles[1]);
     }
-    MTriangle* getTriangle(int a) {
-      return _triangles[a];
+    Rec2DElement* getRTriangle(int a) {
+      return _rtriangles[a];
     }
     void apply();
     void getParities(int*);
     Rec2DVertex* getRVertex(int i) {return _vertices[i];}
     void color(int, int, int);
+    virtual bool whatImply(Rec2DVertex**, int*,
+                           std::map<Rec2DVertex*, std::vector<int> >&);
+    static bool isObsolete(int*);
     
     int getNum() {return 1;}
   private :
     void _computeGlobVal();
 };
 
-class Rec2DFourTri2Quad : public Rec2DAction {
+/*class Rec2DFourTri2Quad : public Rec2DAction {
   private :
-    MTriangle *_triangles[4];
+    Rec2DElement *_rtriangles[4];
     Rec2DEdge *_edges[8]; // 4 embedded, 4 boundary
     Rec2DVertex *_vertices[5]; // 4 boundary (2 non adjacent + 2), 1 embedded
     
   public :
-    Rec2DFourTri2Quad(Rec2DVertex*, std::list<MTriangle*>&);
+    Rec2DFourTri2Quad(Rec2DVertex*, std::vector<MTriangle*>&);
+    ~Rec2DFourTri2Quad();
     bool isObsolete();
-    void getTriangles(std::set<MTriangle*> &tri) {
-      tri.insert(_triangles[0]);
-      tri.insert(_triangles[1]);
-      tri.insert(_triangles[2]);
-      tri.insert(_triangles[3]);
+    void getRTriangles(std::set<Rec2DElement*> &tri) {
+      tri.insert(_rtriangles[0]);
+      tri.insert(_rtriangles[1]);
+      tri.insert(_rtriangles[2]);
+      tri.insert(_rtriangles[3]);
     }
-    MTriangle* getTriangle(int a) {
-      return _triangles[a];
+    Rec2DElement* getRTriangle(int a) {
+      return _rtriangles[a];
     }
     void apply();
     void getParities(int*);
     Rec2DVertex* getRVertex(int i) {return _vertices[i];}
     void color(int, int, int);
+    virtual void whatImply(Rec2DVertex**, int*,
+                           std::map<Rec2DVertex*, std::vector<int> >&)
+    {return;}
     
     int getNum() {return 2;}
   private :
@@ -187,34 +220,38 @@ class Rec2DFourTri2Quad : public Rec2DAction {
 
 class Rec2DCollapseTri : public Rec2DAction {
   private :
-    MTriangle *_triangles[4]; // either 2 or 4 triangles
+    Rec2DElement *_rtriangles[4]; // either 2 or 4 triangles
     Rec2DEdge *_edge; // 1 embedded or NULL
     Rec2DVertex *_vertices[5]; // 4 boundary (2 on embedded edge + 2), 1 embedded or NULL
     double _globValIfExecuted2;
     
   public :
-    Rec2DCollapseTri(Rec2DVertex*, std::list<MTriangle*>&);
-    Rec2DCollapseTri(Rec2DEdge*, std::list<MTriangle*>&);
+    Rec2DCollapseTri(Rec2DVertex*, std::vector<MTriangle*>&);
+    Rec2DCollapseTri(Rec2DEdge*, std::vector<MTriangle*>&);
+    ~Rec2DCollapseTri();
     double getReward() {
       return std::max(_globValIfExecuted, _globValIfExecuted2)
                     - Recombine2D::getGlobalValue(0, .0, 0, .0);
     }
     bool isObsolete();
-    void getTriangles(std::set<MTriangle*> &tri) {
-      tri.insert(_triangles[0]);
-      tri.insert(_triangles[1]);
-      if (_triangles[2]) {
-        tri.insert(_triangles[2]);
-        tri.insert(_triangles[3]);
+    void getRTriangles(std::set<Rec2DElement*> &tri) {
+      tri.insert(_rtriangles[0]);
+      tri.insert(_rtriangles[1]);
+      if (_rtriangles[2]) {
+        tri.insert(_rtriangles[2]);
+        tri.insert(_rtriangles[3]);
       }
     }
-    MTriangle* getTriangle(int a) {
-      return _triangles[a];
+    Rec2DElement* getRTriangle(int a) {
+      return _rtriangles[a];
     }
     void apply();
     void getParities(int*);
     Rec2DVertex* getRVertex(int i) {return _vertices[i];}
     void color(int, int, int);
+    virtual void whatImply(Rec2DVertex**, int*,
+                           std::map<Rec2DVertex*, std::vector<int> >&)
+    {return;}
     
     int getNum() {return 3;}
   private :
@@ -222,7 +259,8 @@ class Rec2DCollapseTri : public Rec2DAction {
     void _qualCavity(double &valVert, double &valEdge, int &numEdge,
                      std::map<Rec2DVertex*, int> &vertices,
                      Rec2DVertex *imposedVert = NULL                );
-};
+};*/
+//
 
 class Rec2DListAction {
   private :
@@ -236,7 +274,8 @@ class Rec2DVertex {
   private :
     MVertex *_v;
     std::vector<Rec2DEdge*> _edges;
-    int _onWhat, _numEl, _parity, _lastChange;
+    std::vector<Rec2DElement*> _elements;
+    int _onWhat, _numEl, _parity, _lastChange; //lastMove
     // _onWhat={-1:corner,0:edge,1:face,(2:volume)}
     double _angle, _qual;
     SPoint2 _param;
@@ -245,10 +284,11 @@ class Rec2DVertex {
     static double **_gains;
     
   public :
-    Rec2DVertex(MVertex*, std::list<MTriangle*>, int onWhat,
+    Rec2DVertex(MVertex*, std::vector<MTriangle*>, int onWhat,
                 std::map<MVertex*, std::set<GEdge*> >);
     Rec2DVertex(double u, double v);
     
+    bool hasTriangle();
     inline void add(Rec2DEdge *re) {_edges.push_back(re);}
     void remove(Rec2DEdge*);
     inline void setOnBoundary() {if (_onWhat > 0) _onWhat = 0;
@@ -270,22 +310,33 @@ class Rec2DVertex {
     double getGain(int);
     void getNeighbours(std::map<Rec2DVertex*, int>&, std::set<Rec2DEdge*>&);
     void getNeighbourEdges(std::set<Rec2DEdge*>&);
+    void getElements(std::set<Rec2DElement*>&);
+    
+    void remove(Rec2DElement *rel) {
+      std::vector<Rec2DElement*>::iterator it = _elements.begin();
+      while (it != _elements.end()) {
+        if (*it == rel)
+          it = _elements.erase(it);
+        else
+          ++it;
+      }
+    }
     
     static void initStaticTable();
     
   private :
     void _computeQual();
-    double _computeAngle(MVertex*, std::list<MTriangle*>&, std::set<GEdge*>&);
+    double _computeAngle(MVertex*, std::vector<MTriangle*>&, std::set<GEdge*>&);
 };
 
 class Rec2DEdge {
   private :
     Rec2DVertex *_rv1, *_rv2;
     double _qual;
-    int _lastChange, _weight;
+    int _lastChange, _weight; //lastUpdate
     
   public :
-    Rec2DEdge(MEdge, mapofVertices&, std::list<MTriangle*>&);
+    Rec2DEdge(MEdge, mapofVertices&, std::vector<MTriangle*>&);
     Rec2DEdge(Rec2DVertex*, Rec2DVertex*);
     
     double getQual();
@@ -300,4 +351,74 @@ class Rec2DEdge {
     double _straightAlignment();
 };
 
+class Rec2DElement {
+  private :
+    Rec2DEdge *_redges[4]; // if triangle, _redges[3] == NULL
+    Rec2DElement *_relements[4]; // NULL if no neighbour
+    MElement *_mEl; // can be NULL
+    std::vector<Rec2DAction*> _actions;
+    
+  public :
+    Rec2DElement(MTriangle *t) : _mEl((MElement *)t) {
+      for (int i = 0; i < 4; ++i) {
+        _redges[i] = NULL;
+      }
+    }
+    
+    void add(Rec2DEdge*);
+    void add(Rec2DEdge*, Rec2DElement*);
+    inline MElement* getMElement() {return _mEl;}
+    inline MTriangle* getMTriangle()
+      {if (_redges[3])return NULL; return (MTriangle*) _mEl;}
+    inline bool isQuad() {return _redges[3];}
+    inline void add(Rec2DAction *ra) {_actions.push_back(ra);}
+    void remove(Rec2DAction *ra) {
+      std::vector<Rec2DAction*>::iterator it = _actions.begin();
+      while (it != _actions.end()) {
+        if (*it == ra) it = _actions.erase(it);
+        else ++it;
+      }
+    }
+    inline int getNum() {return _mEl->getNum();}
+    inline int getNumActions() {return _actions.size();}
+    inline Rec2DAction* getAction(int i) {return _actions[i];}
+    void getVertices(Rec2DVertex **v) {
+      v[0] = _redges[0]->getRVertex(0);
+      v[1] = _redges[0]->getRVertex(1);
+      if (_redges[1]->getRVertex(0) == v[0] ||
+          _redges[1]->getRVertex(0) == v[1]   )
+        v[2] = _redges[1]->getRVertex(1);
+      else
+        v[2] = _redges[1]->getRVertex(0);
+    }
+    void getActions(std::set<Rec2DAction*> &actions) {
+      for (unsigned int i = 0; i < _actions.size(); ++i)
+        actions.insert(_actions[i]);
+    }
+    void getNeighbours(std::set<Rec2DElement*> &elem) {
+      for (int i = 0; i < 4; ++i) {
+        if (_relements[i])
+          elem.insert(_relements[i]);
+      }
+    }
+    Rec2DVertex* getOtherVertex(Rec2DVertex *rv1, Rec2DVertex *rv2) {
+      if (_redges[3]) {
+        Msg::Error("[Rec2DElement] I'm not a triangle");
+        return NULL;
+      }
+      std::set<Rec2DVertex*> vert;
+      for (int i = 0; i < 3; ++i) {
+        vert.insert(_redges[i]->getRVertex(0));
+        vert.insert(_redges[i]->getRVertex(1));
+      }
+      std::set<Rec2DVertex*>::iterator it = vert.begin();
+      for (; it != vert.end(); ++it) {
+        if (*it != rv1 && *it != rv2)
+          return *it;
+      }
+      Msg::Error("[Rec2DElement] I should not be here... Why this happen to me ?");
+      return NULL;
+    }
+};
+
 #endif
\ No newline at end of file
-- 
GitLab