Creation of actions

Mesh/meshGFaceRecombine.cpp

+// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+//   Amaury Johnen (a.johnen@ulg.ac.be)
+//
+
+// TODO :
+//   - recombine w. take care of parity, update of quality...
+
+#include "meshGFaceRecombine.h"
+#include "MElement.h"
+#include "MTriangle.h"
+
+#define EDGE_BASE 2 
+#define EDGE_QUAD 1
+
+Recombine2D *Recombine2D::_current = NULL;
+bool Recombine2D::bo = 0;
+double **Rec2DVertex::_qualVSnum = NULL;
+double **Rec2DVertex::_gains = NULL;
+
+/**  Recombine2D  **/
+/*******************/
+Recombine2D::Recombine2D(GFace *gf) : _gf(gf)
+{
+  if (Recombine2D::_current != NULL) {
+    Msg::Info("An instance of recombination is already in execution");
+    return;
+  }
+  Recombine2D::_current = this;
+  _numChange = 0;
+  backgroundMesh::set(gf);
+  _bgm = backgroundMesh::current();
+  Rec2DVertex::initStaticTable();
+  _numEdge = _numVert = 0;
+  _valEdge = _valVert = .0;
+  
+  // Be able to compute geometrical angle at corners
+  std::map<MVertex*, std::set<GEdge*> > mapCornerVert;
+  {
+    std::list<GEdge*> listge = gf->edges();
+    std::list<GEdge*>::iterator itge = listge.begin();
+    for (; itge != listge.end(); ++itge) {
+      mapCornerVert[(*itge)->getBeginVertex()->getMeshElement(0)->getVertex(0)]
+        .insert(*itge);
+      mapCornerVert[(*itge)->getEndVertex()->getMeshElement(0)->getVertex(0)]
+        .insert(*itge);
+    }
+  }
+  
+  // Find all Vertices and edges
+  std::map<MVertex*, std::list<MTriangle*> > mapVert;
+  std::map<MEdge, std::list<MTriangle*>, Less_Edge> mapEdge;
+  for (unsigned int i = 0; i < gf->triangles.size(); ++i) {
+    MTriangle *t = gf->triangles[i];
+    for (int j = 0; j < 3; ++j) {
+      mapVert[t->getVertex(j)].push_back(t);
+      mapEdge[t->getEdge(j)].push_back(t);
+    }
+  }
+  
+  // Create the 'Rec2DVertex' and store iterator to vertices which have degree 4
+  std::list<std::map<MVertex*, std::list<MTriangle*> >::iterator> fourTri;
+  {
+    mapofVertices::iterator itV2N = _vertices.begin();
+    std::map<MVertex*, std::list<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;
+      if (itvert->second.size() == 4)
+        fourTri.push_back(itvert);
+    }
+  }
+  
+  // Create the 'Rec2DEdge' and store boundary edges
+  std::map<MVertex*, std::list<MEdge> > boundV2E;
+  {
+    mapofEdges::iterator itE2E = _edges.begin();
+    std::map<MEdge, std::list<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 += EDGE_BASE;
+      _valEdge += EDGE_BASE * re->getQual();
+      if (itedge->second.size() == 1) {
+        _numEdge += EDGE_QUAD;
+        _valEdge += EDGE_QUAD * re->getQual();
+        _vertices[e.getVertex(0)]->setOnBoundary();
+        _vertices[e.getVertex(1)]->setOnBoundary();
+        boundV2E[e.getVertex(0)].push_back(e);
+        boundV2E[e.getVertex(1)].push_back(e);
+      }
+    }
+  }
+  
+  // We know now on what are vertices, compute reward
+  {
+    mapofVertices::iterator itV2N = _vertices.begin();
+    for (; itV2N != _vertices.end(); ++itV2N)
+      _valVert += itV2N->second->getQual();
+  }
+  
+  // Be dealing with "parity" on boundaries
+  {
+    std::map<MVertex*, std::list<MEdge> >::iterator itBound = boundV2E.begin();
+    int actualParity = 0;
+    bool p = 0;
+    while (itBound != boundV2E.end()) {
+      MVertex *firstV = itBound->first, *actualV = firstV;
+      MEdge actualEdge = *itBound->second.begin();
+      Rec2DVertex *rv = _vertices[firstV];
+      rv->setParity(actualParity+p);
+      _parities[actualParity+p].insert(rv);
+      p = !p;
+      if (actualEdge.getVertex(0) == actualV)
+        actualV = actualEdge.getVertex(1);
+      else
+        actualV = actualEdge.getVertex(0);
+      while (actualV != firstV) {
+        rv = _vertices[actualV];
+        rv->setParity(actualParity+p);
+        _parities[actualParity+p].insert(rv);
+        p = !p;
+        itBound = boundV2E.find(actualV);
+        std::list<MEdge>::iterator it = itBound->second.begin();
+        if (*it == actualEdge)
+          actualEdge = *(++it);
+        else
+          actualEdge = *it;
+        if (actualEdge.getVertex(0) == actualV)
+          actualV = actualEdge.getVertex(1);
+        else
+          actualV = actualEdge.getVertex(0);
+        boundV2E.erase(itBound);
+      }
+      if (_vertices[actualV]->getParity() != actualParity+p)
+        Msg::Error("Not even number of mesh edges on boundary !");
+      
+      boundV2E.erase(boundV2E.begin());
+      itBound = boundV2E.begin();
+      actualParity += 2;
+    }
+  }
+  
+  // Create the 'Rec2DTwoTri2Quad' and 'Rec2DCollapseTri'
+  {
+    mapofEdges::iterator itE2E = _edges.begin();
+    std::map<MEdge, std::list<MTriangle*> >::iterator ite = mapEdge.begin();
+    for (; ite != mapEdge.end(); ++ite)
+    if (ite->second.size() == 2) {
+      _actions.insert(new Rec2DTwoTri2Quad(_edges[ite->first], ite->second));
+      _actions.insert(new Rec2DCollapseTri(_edges[ite->first], ite->second));
+    }
+  }
+  
+  // 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()) {
+        _actions.insert(new Rec2DFourTri2Quad(rv, (*it4)->second));
+        _actions.insert(new Rec2DCollapseTri(rv, (*it4)->second));
+      }
+    }
+  }
+  
+  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(" ");
+}
+
+Recombine2D::~Recombine2D()
+{
+  Recombine2D::_current = NULL;
+  //delete vertices, edges, actions;
+}
+
+double Recombine2D::getGlobalValue(int numEdge, double valEdge,
+                                   int numVert, double valVert )
+{
+  return (_current->_valEdge + valEdge) / (_current->_numEdge + numEdge) *
+         (_current->_valVert + valVert) / (_current->_numVert + numVert);
+}
+
+Rec2DEdge* Recombine2D::getREdge(MEdge e)
+{
+  mapofEdges::iterator it = _current->_edges.find(e);
+  if (it != _current->_edges.end())
+    return it->second;
+  Msg::Error("[Recombine2D::getREdge] edge not found");
+  return NULL;
+}
+
+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;
+}
+
+
+/**  Rec2DTwoTri2Quad  **/
+/************************/
+Rec2DTwoTri2Quad::Rec2DTwoTri2Quad(Rec2DEdge *re, 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 < 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());
+  }
+  
+  _edges[0] = re;
+  std::set<MEdge>::iterator ite = extEdges.begin();
+  for (i = 1; ite != extEdges.end() && i < 5; ++ite, ++i)
+    _edges[i] = Recombine2D::getREdge(*ite);
+  if (ite != extEdges.end() || i < 5)
+    Msg::Error("[Rec2DTwoTri2Quad] Wrong number of edges");
+  
+  _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);
+  
+  _computeGlobVal();
+}
+
+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                                    )
+    return true;
+  return false;
+}
+
+void Rec2DTwoTri2Quad::_computeGlobVal()
+{
+  double valEdge = - EDGE_BASE * _edges[0]->getQual();
+  for (int i = 1; i < 5; ++i)
+    valEdge += EDGE_QUAD * _edges[i]->getQual();
+  
+  double valVert = _vertices[0]->getGain(-1) + _vertices[1]->getGain(-1);
+  
+  _globValIfExecuted = Recombine2D::getGlobalValue(4*EDGE_QUAD - EDGE_BASE,
+                                                   valEdge, 0, valVert     );
+  _lastUpdate = Recombine2D::getNumChange();
+}
+
+/**  Rec2DFourTri2Quad  **/
+/*************************/
+Rec2DFourTri2Quad::Rec2DFourTri2Quad(Rec2DVertex *rv, std::list<MTriangle*> &tri)
+{
+  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;
+    }
+    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();
+}
+
+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::_computeGlobVal()
+{
+  double valEdge = 0;
+  for (int i = 0; i < 4; ++i)
+    valEdge -= EDGE_BASE * _edges[i]->getQual();
+  for (int i = 4; i < 8; ++i)
+    valEdge += EDGE_QUAD * _edges[i]->getQual();
+    
+  double valVert = - _vertices[0]->getQual();
+  for (int i = 1; i < 5; ++i)
+    valEdge += _vertices[i]->getGain(-1);
+  
+  _globValIfExecuted = Recombine2D::getGlobalValue(4*EDGE_QUAD - 4*EDGE_BASE,
+                                                   valEdge, -1, valVert      );
+  _lastUpdate = Recombine2D::getNumChange();
+}
+
+/**  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 (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::_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;
+  }
+}
+
+/**  Rec2DVertex  **/
+/*******************/
+Rec2DVertex::Rec2DVertex(MVertex *v, std::list<MTriangle*> tri, int onWhat,
+                         std::map<MVertex*, std::set<GEdge*> > gEdges)
+: _v(v), _onWhat(onWhat), _parity(-1), _lastChange(-1), _numEl(tri.size())
+{
+  if (_onWhat == -1) {
+    std::map<MVertex*, std::set<GEdge*> >::iterator it = gEdges.find(_v);
+    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); 
+  _computeQual();
+}
+
+Rec2DVertex::Rec2DVertex(double u, double v)
+: _onWhat(1), _parity(-1), _lastChange(-1)
+{
+  _param[0] = u;
+  _param[1] = v;
+  
+  GPoint p = Recombine2D::getGFace()->point(_param);
+  _v = new MVertex(p.x(), p.y(), p.z(), Recombine2D::getGFace());
+  static int i = 0;
+  if (++i < 20)
+    Msg::Info("resulting point = [%g,%g,%g]", p.x(), p.y(), p.z());
+}
+
+double Rec2DVertex::_computeAngle(MVertex *v,
+                                  std::list<MTriangle*> &listTri,
+                                  std::set<GEdge*> &setEdge)
+{
+  if (setEdge.size() != 2) {
+    Msg::Warning("[meshGFaceRecombine] Wrong number of edge : %d, returning pi/2", setEdge.size());
+    return M_PI / 2.;
+  }
+  static const double prec = 100.;
+  
+  SVector3 vectv = SVector3(v->x(), v->y(), v->z());
+  SVector3 firstDer0, firstDer1;
+  
+  std::set<GEdge*>::iterator it = setEdge.begin();
+  for (int k = 0; k < 2; ++it, ++k) {
+    GEdge *ge = *it;
+    SVector3 vectlb = ge->position(ge->getLowerBound());
+    SVector3 vectub = ge->position(ge->getUpperBound());
+    vectlb -= vectv;
+    vectub -= vectv;
+    double normlb, normub;
+    if ((normlb = norm(vectlb)) > prec * (normub = norm(vectub)))
+      firstDer1 = ge->firstDer(ge->getUpperBound());
+    else if (normub > prec * normlb)
+      firstDer1 = ge->firstDer(ge->getLowerBound());
+    else {
+      Msg::Warning("[meshGFaceRecombine] Bad precision, returning pi/2");
+      return M_PI / 2.;
+    }
+    if (k == 0)
+      firstDer0 = firstDer1;
+  }
+  
+  firstDer0 = firstDer0 * (1 / norm(firstDer0));
+  firstDer1 = firstDer1 * (1 / norm(firstDer1));
+  
+  double angle1 = acos(dot(firstDer0, firstDer1));
+  double angle2 = 2. * M_PI - angle1;
+  
+  double angleMesh = .0;
+  std::list<MTriangle*>::iterator it2 = listTri.begin();
+  for (; it2 != listTri.end(); ++it2) {
+    MTriangle *t = *it2;
+    int k = 0;
+    while (t->getVertex(k) != v && k < 3)
+      ++k;
+    if (k == 3) {
+      Msg::Warning("[meshGFaceRecombine] Wrong triangle, returning pi/2");
+      return M_PI / 2.;
+    }
+    angleMesh += angle3Vertices(t->getVertex((k+1)%3), v, t->getVertex((k+2)%3));
+  }
+  
+  if (angleMesh < M_PI)
+    return angle1;
+  return angle2;
+}
+
+double Rec2DVertex::getQual(int numEdge)
+{
+  if (numEdge > 1) {
+    switch (_onWhat) {
+      case 1 :
+        return _qualVSnum[_onWhat][numEdge-1];
+      case 0 :
+        return _qualVSnum[_onWhat][numEdge-2];
+      case -1 :
+        return 1. - fabs(2./M_PI * _angle/(numEdge-1) - 1.);
+      default :
+        Msg::Error("[Rec2DVertex] Don't know on what is the vertex");
+        return -1.;
+    }
+  }
+  if (_lastChange < Recombine2D::getNumChange())
+    _computeQual();
+  return _qual;
+}
+
+void Rec2DVertex::_computeQual()
+{
+  if (_onWhat > -1)
+    _qual = _qualVSnum[_onWhat][_numEl-1];
+  else
+    _qual = 1. - fabs(2./M_PI * _angle/_numEl - 1.);
+  _lastChange = Recombine2D::getNumChange();
+}
+
+double Rec2DVertex::getGain(int n)
+{
+  if (_onWhat > -1) {
+    switch (n) {
+      case 1 :
+        return _gains[_onWhat][_numEl-1];
+      case -1 :
+        return - _gains[_onWhat][_numEl-2];
+      default :
+        return _qualVSnum[_onWhat][_numEl-1+n] - _qualVSnum[_onWhat][_numEl-1];
+    }
+  }
+  else
+    return fabs(2./M_PI * _angle/_numEl - 1.)
+         - fabs(2./M_PI * _angle/(_numEl + n) - 1.);
+}
+
+void Rec2DVertex::getNeighbours(std::map<Rec2DVertex*, int> &vert,
+                                std::set<Rec2DEdge*> &edge        )
+{
+  Rec2DVertex *rv;
+  std::map<Rec2DVertex*, int>::iterator it;
+  for (unsigned int i = 0; i < _edges.size(); ++i){
+    edge.insert(_edges[i]);
+    if (_edges[i]->getRVertex(0) != this)
+      rv = _edges[i]->getRVertex(0);
+    else
+      rv = _edges[i]->getRVertex(1);
+    if ((it = vert.find(rv)) != vert.end())
+      it->second += _edges[i]->getWeight() - EDGE_BASE;
+    else
+      vert[rv] = _edges[i]->getWeight();
+  }
+}
+
+void Rec2DVertex::getNeighbourEdges(std::set<Rec2DEdge*> &edge)
+{
+  for (unsigned int i = 0; i < _edges.size(); ++i)
+    edge.insert(_edges[i]);
+}
+
+void Rec2DVertex::initStaticTable()
+{
+  if (_qualVSnum == NULL) {
+    // _qualVSnum[onWhat][nEl]; onWhat={0:edge,1:face} / nEl=_numEl-1
+    // _gains[onWhat][nEl]; onWhat={0:edge,1:face} / nEl=_numEl-1 (earning of adding an element)
+    int nE = 10, nF = 20;
+    
+    _qualVSnum = new double*[2];
+    _qualVSnum[0] = new double[nE];
+    for (int i = 0; i < nE; ++i)
+      _qualVSnum[0][i] = 1. - fabs(2./(i+1) - 1.);
+    _qualVSnum[1] = new double[nF];
+    for (int i = 0; i < nF; ++i)
+      _qualVSnum[1][i] = 1. - fabs(4./(i+1) - 1.);
+      
+    _gains = new double*[2];
+    _gains[0] = new double[nE-1];
+    for (int i = 0; i < nE-1; ++i)
+      _gains[0][i] = _qualVSnum[0][i+1] - _qualVSnum[0][i];
+    _gains[1] = new double[nF-1];
+    for (int i = 0; i < nF-1; ++i)
+      _gains[1][i] = _qualVSnum[1][i+1] - _qualVSnum[1][i];
+  }
+}
+
+/**  Rec2DEdge  **/
+/*****************/
+Rec2DEdge::Rec2DEdge(MEdge e, mapofVertices &vert, std::list<MTriangle*> &tri)
+: _lastChange(-1), _qual(.0), _weight(EDGE_BASE),
+  _rv1(vert[e.getVertex(0)]), _rv2(vert[e.getVertex(1)])
+{
+  _rv1->add(this);
+  _rv2->add(this);
+}
+
+Rec2DEdge::Rec2DEdge(Rec2DVertex *rv1, Rec2DVertex *rv2)
+: _lastChange(-1), _qual(.0), _weight(EDGE_BASE), _rv1(rv1), _rv2(rv2)
+{
+}
+
+double Rec2DEdge::getQual()
+{
+  if (_lastChange < Recombine2D::getNumChange())
+    _computeQual();
+  return _qual;
+}
+
+void Rec2DEdge::_computeQual()
+{
+  static int i = 0; if (++i < 2) Msg::Warning("FIXME compute qual edge and update of vertex");
+  //_rv1->update();
+  //_rv2->update();
+  double adimLength = _straightAdimLength();
+  double alignment = _straightAlignment();
+  _qual = Recombine2D::edgeReward(adimLength, alignment);
+  _lastChange = Recombine2D::getNumChange();
+}
+
+double Rec2DEdge::_straightAdimLength()
+{
+  double dx, dy, dz, **xyz;
+  xyz = new double*[2];
+  xyz[0] = new double[3];
+  xyz[1] = new double[3];
+  _rv1->getxyz(xyz[0]);
+  _rv2->getxyz(xyz[1]);
+  dx = xyz[0][0] - xyz[1][0];
+  dy = xyz[0][1] - xyz[1][1];
+  dz = xyz[0][2] - xyz[1][2];
+  double length = sqrt(dx * dx + dy * dy + dz * dz);
+  delete[] xyz[0];
+  delete[] xyz[1];
+  delete[] xyz;
+  
+  double lc1 = (*Recombine2D::bgm())(_rv1->u(), _rv1->v(), .0);
+  double lc2 = (*Recombine2D::bgm())(_rv2->u(), _rv2->v(), .0);
+  
+  return length * (1/lc1 + 1/lc2) / 2;
+}
+
+double Rec2DEdge::_straightAlignment()
+{
+  double angle1 = Recombine2D::bgm()->getAngle(_rv1->u(), _rv1->v(), .0);
+  double angle2 = Recombine2D::bgm()->getAngle(_rv1->u(), _rv1->v(), .0);
+  double angleEdge = atan2(_rv1->u()-_rv2->u(), _rv1->v()-_rv2->v());
+  
+  double alpha1 = angleEdge - angle1;
+  double alpha2 = angleEdge - angle2;
+  crossField2d::normalizeAngle(alpha1);
+  crossField2d::normalizeAngle(alpha2);
+  alpha1 = std::min(alpha1, .5 * M_PI - alpha1);
+  alpha2 = std::min(alpha2, .5 * M_PI - alpha2);
+  
+  double lc1 = (*Recombine2D::bgm())(_rv1->u(), _rv1->v(), .0);
+  double lc2 = (*Recombine2D::bgm())(_rv1->u(), _rv1->v(), .0);
+  
+  return (alpha1/lc1 + alpha2/lc2) / (1/lc1 + 1/lc2);
+}
\ No newline at end of file

Mesh/meshGFaceRecombine.h

+// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+//   Amaury Johnen (a.johnen@ulg.ac.be)
+//
+
+#ifndef _MESH_GFACE_RECOMBINE_H_
+#define _MESH_GFACE_RECOMBINE_H_
+
+#include "GFace.h"
+#include "GModel.h"
+#include "MEdge.h"
+#include "BackgroundMesh.h"
+
+class Rec2DEdge;
+class Rec2DVertex;
+class Rec2DAction;
+
+typedef std::set<Rec2DAction*/*, lessRec2DAction*/> setofRec2DAction;
+typedef std::map<MVertex*, Rec2DVertex*> mapofVertices;
+typedef std::map<MEdge, Rec2DEdge*, Less_Edge> mapofEdges;
+typedef std::map<MElement*, std::set<Rec2DAction*> > mapofElementActions;
+typedef std::map<MQuadrangle*, std::set<MElement*> > mapofAdjacencies;
+
+class Recombine2D {
+  private :
+    int _numEdge, _numVert, _numChange;
+    double _valEdge, _valVert;
+    GFace *_gf;
+    backgroundMesh *_bgm;
+    mapofEdges _edges;
+    mapofVertices _vertices;
+    setofRec2DAction _actions;
+    mapofElementActions _mea;
+    mapofAdjacencies _ma;
+    std::map<int, std::set<Rec2DVertex*> > _parities;
+    std::vector<Rec2DAction*> _obsoleteAction;
+    static Recombine2D *_current;
+    
+  public :
+    static bool bo;
+    
+    Recombine2D(GFace*);
+    Recombine2D(GModel*);
+    ~Recombine2D();
+    
+    void apply(bool){return;}
+    
+    static inline double edgeReward(double adimLength, double alignment)
+      {return adimLength * (.5 + .5 * 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 Rec2DEdge* getREdge(MEdge);
+    static Rec2DVertex* getRVertex(MVertex*);
+    static inline GFace* getGFace() {return _current->_gf;}
+    static inline void addObsleteAction(Rec2DAction *a) {
+      _current->_obsoleteAction.push_back(a);
+    }
+    //static inline int* getPtrNumChange() {return &_numChange;}
+};
+
+class Rec2DAction {
+  private :
+    int _position;
+    
+  protected :
+    int _lastUpdate;
+    double _globValIfExecuted;
+    
+  public :
+    inline int getVectPosition() {return _position;}
+    inline void setVectPosition(int p) {_position = p;}
+    //virtual double getReward() = 0;
+    virtual bool isObsolete() = 0;
+  
+  private :
+    virtual void _computeGlobVal() = 0;
+};
+
+class Rec2DTwoTri2Quad : public Rec2DAction {
+  private :
+    MTriangle *_triangles[2];
+    Rec2DEdge *_edges[5]; // 1 embedded, 4 boundary
+    Rec2DVertex *_vertices[4]; // 4 boundary (2 on embedded edge + 2)
+    
+  public :
+    Rec2DTwoTri2Quad(Rec2DEdge*, std::list<MTriangle*>&);
+    //double getReward();
+    bool isObsolete();
+  
+  private :
+    void _computeGlobVal();
+};
+
+class Rec2DFourTri2Quad : public Rec2DAction {
+  private :
+    MTriangle *_triangles[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*>&);
+    //double getReward();
+    bool isObsolete();
+  
+  private :
+    void _computeGlobVal();
+};
+
+class Rec2DCollapseTri : public Rec2DAction {
+  private :
+    MTriangle *_triangles[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*>&);
+    //double getReward();
+    bool isObsolete();
+  
+  private :
+    void _computeGlobVal();
+    void _qualCavity(double &valVert, double &valEdge, int &numEdge,
+                     std::map<Rec2DVertex*, int> &vertices,
+                     Rec2DVertex *imposedVert = NULL                );
+    
+};
+
+class Rec2DListAction {
+  private :
+    
+    
+  public :
+    
+};
+
+class Rec2DVertex {
+  private :
+    MVertex *_v;
+    std::vector<Rec2DEdge*> _edges;
+    int _onWhat, _numEl, _parity, _lastChange;
+    // _onWhat={-1:corner,0:edge,1:face,(2:volume)}
+    double _angle, _qual;
+    SPoint2 _param;
+    
+    static double **_qualVSnum;
+    static double **_gains;
+    
+  public :
+    Rec2DVertex(MVertex*, std::list<MTriangle*>, int onWhat,
+                std::map<MVertex*, std::set<GEdge*> >);
+    Rec2DVertex(double u, double v);
+    
+    void add(Rec2DEdge *re) {_edges.push_back(re);}
+    void setOnBoundary() {if (_onWhat > 0) _onWhat = -1;
+      static bool a=1; if(a){Msg::Error("FIXME boundary");a=0;}}
+    bool getIsOnBoundary() {return _onWhat < 1;}
+    double getQual(int numEdge = -1);
+    inline int getParity() {return _parity;}
+    inline void setParity(int p) {_parity = p;}
+    inline MVertex* getMVertex() {return _v;}
+    void getxyz(double *xyz) {
+      xyz[0] = _v->x(); xyz[1] = _v->y(); xyz[2] = _v->z();}
+    inline double u() {return _param[0];}
+    inline double v() {return _param[1];}
+    inline void deleteMVertex() {delete _v;}
+    
+    double getGain(int);
+    void getNeighbours(std::map<Rec2DVertex*, int>&, std::set<Rec2DEdge*>&);
+    void getNeighbourEdges(std::set<Rec2DEdge*>&);
+    
+    static void initStaticTable();
+    
+  private :
+    void _computeQual();
+    double _computeAngle(MVertex*, std::list<MTriangle*>&, std::set<GEdge*>&);
+};
+
+class Rec2DEdge {
+  private :
+    Rec2DVertex *_rv1, *_rv2;
+    double _qual;
+    int _lastChange, _weight;
+    
+  public :
+    Rec2DEdge(MEdge, mapofVertices&, std::list<MTriangle*>&);
+    Rec2DEdge(Rec2DVertex*, Rec2DVertex*);
+    
+    double getQual();
+    inline MEdge getMEdge() {return MEdge(_rv1->getMVertex(), _rv2->getMVertex());}
+    inline Rec2DVertex* getRVertex(int i) {if (i) return _rv2; return _rv1;}
+    int getWeight() {return _weight;}
+    double getWeightedQual() {return _weight * _qual;}
+    
+  private :
+    void _computeQual();
+    double _straightAdimLength();
+    double _straightAlignment();
+};
+
+#endif
+
+//idee : lors du parcours de larbre avec un horizon : copier les Rec2dVertex et les Rec2DEdge.
+//idee : faire carrment : arrete d'un triangle = 0, arrete d'un quad = 1/2