diff --git a/Mesh/meshGFaceRecombine.cpp b/Mesh/meshGFaceRecombine.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..53abddba3d31f5e740b3fce7cc8a61733c14633b
--- /dev/null
+++ b/Mesh/meshGFaceRecombine.cpp
@@ -0,0 +1,863 @@
+// 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
diff --git a/Mesh/meshGFaceRecombine.h b/Mesh/meshGFaceRecombine.h
new file mode 100644
index 0000000000000000000000000000000000000000..619221cec55a15ba016511a4f3a348252219cca0
--- /dev/null
+++ b/Mesh/meshGFaceRecombine.h
@@ -0,0 +1,211 @@
+// 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 carr�ment : arrete d'un triangle = 0, arrete d'un quad = 1/2