diff --git a/Mesh/meshRecombine2D.cpp b/Mesh/meshRecombine2D.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fc535af1096ab06552d44253bda56b5e1115f7bc
--- /dev/null
+++ b/Mesh/meshRecombine2D.cpp
@@ -0,0 +1,437 @@
+// 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 (amjohnen@gmail.com) adapted from meshGFaceOptimize
+//
+
+#include "meshRecombine2D.h"
+#include "GFace.h"
+
+static int HORIZ = 25;
+
+static int NUM_NEW = 0;
+static int NEW_REC = 0;
+static int NEW_NODE = 0;
+static int NEW_QUAD = 0;
+
+static int NUM_DEL = 0;
+static int DEL_REC = 0;
+static int DEL_NODE = 0;
+static int DEL_QUAD = 0;
+
+Recombine2D::Recombine2D(GFace *gf, int horizon)
+ : _horizon(horizon), _gf(gf), _benef(.0), _applied(false)
+{
+ if (_horizon < 1) _horizon = 1;
+
+ Msg::Info("Branching with horizon %d", _horizon = HORIZ);
+
+ _haveParam = gf->geomType() != GEntity::DiscreteSurface || gf->getCompound();
+ if(!_haveParam) {
+ Msg::Warning("Surfaces must have prametrization in order to recombine triangles");
+ //return;
+ _haveParam=true;
+ }
+
+ e2t_cont adj;
+ _buildEdgeToElement(gf->triangles, adj);
+
+ std::map<MElement*, int> boundaryTriangle;
+ std::pair<Set_Recomb::iterator, bool> ret;
+
+ for (e2t_cont::iterator it = adj.begin(); it != adj.end(); ++it) {
+ if (it->second.second) {
+ RecombTriangle *rt = new RecombTriangle(it->first, it->second.first, it->second.second);
+NUM_NEW++;
+NEW_REC++;
+ do {
+ rt->touch();
+ ret = _pairs.insert(rt);
+ } while (ret.second == false);
+ _possibleRecomb[it->second.first].insert(rt);
+ _possibleRecomb[it->second.second].insert(rt);
+ }
+ else {
+ if (boundaryTriangle.find(it->second.first) == boundaryTriangle.end())
+ boundaryTriangle[it->second.first] = 1;
+ else
+ boundaryTriangle[it->second.first] = 2;
+ }
+ }
+
+ std::map<MElement*, int>::iterator itt;
+
+ for (itt = boundaryTriangle.begin(); itt != boundaryTriangle.end(); ++itt) {
+ RecombTriangle *rt;
+ if (itt->second == 1)
+ rt = new RecombTriangle(itt->first, -100);
+ else
+ rt = new RecombTriangle(itt->first, -200);
+NUM_NEW++;
+NEW_REC++;
+
+ ret = _pairs.insert(rt);
+ while (!ret.second) {
+ rt->touch();
+ ret = _pairs.insert(rt);
+ }
+ _possibleRecomb[itt->first].insert(rt);
+ }
+
+ Msg::Info("poss %d", _possibleRecomb.size());
+ Msg::Info("pairs %d", _pairs.size());
+ Msg::Info("tri %d", gf->triangles.size());
+
+ _recombine();
+
+ Msg::Info("new %d = %d + %d + %d", NUM_NEW, NEW_REC, NEW_NODE, NEW_QUAD);
+ Msg::Info("del %d = %d + %d + %d", NUM_DEL, DEL_REC, DEL_NODE, DEL_QUAD);
+}
+
+Recombine2D::~Recombine2D()
+{
+ if (!_applied)
+ for (int i = 0; i < _quads.size(); i++) {
+ delete _quads[i];
+NUM_DEL++;
+DEL_QUAD++;
+ }
+}
+
+template <class E>
+void Recombine2D::_buildEdgeToElement(std::vector<E*> &elements, e2t_cont &adj)
+{
+ for (unsigned int i = 0; i < elements.size(); i++) {
+ E *el = elements[i];
+ for (int j = 0; j < el->getNumEdges(); j++) {
+ MEdge e = el->getEdge(j);
+ e2t_cont::iterator it = adj.find(e);
+ if (it == adj.end()) {
+ std::pair<MElement*, MElement*> one = std::make_pair(el, (MElement*)0);
+ adj[e] = one;
+ }
+ else
+ it->second.second = el;
+ }
+ }
+}
+
+void Recombine2D::_recombine()
+{
+ int i = 0;
+ while ((!_pairs.empty() || !_lastRecomb.empty()) && i < 1000) {
+ Msg::Info("%d",i);
+ Set_Recomb::iterator it = _bestNextRecombination();
+ _benef += (*it)->getBenef();
+ if ((*it)->isQuad())
+ _quads.push_back((*it)->createQuad());
+ else
+ for (int i = 0; i < (*it)->numTriangles(); i++)
+ _isolated.insert((*it)->triangle(i));
+ _removeImpossible(it);
+ i++;
+ }
+ Msg::Info("Done %d (%d)", i, _pairs.size());
+}
+
+Set_Recomb::iterator Recombine2D::_bestNextRecombination()
+{
+ if (!_lastRecomb.empty())
+ return _lastRecomb.begin();
+
+ int h = 0, nSkip = 0;
+ double maxBenef = _horizon * -200.;
+
+
+ Map_Tri_Node touched;
+ std::vector<Node*> nodes(_horizon);
+ Set_Recomb::iterator it = _pairs.begin(), itmax = _pairs.begin();
+
+ bool end = false;
+
+ int numb = 0;
+ while (!end) {
+ numb++;
+ while (h < _horizon && it != _pairs.end()) {
+ bool possible = true;
+
+ Map_Tri_Node::iterator itm;
+ for (int i = 0; i < (*it)->numTriangles(); i++) {
+ itm = touched.find((*it)->triangle(i));
+ if (itm != touched.end()) {
+ itm->second->blocking(itm);
+ possible = false;
+ }
+ }
+
+ if (possible) {
+ if (h == 0) nodes[h] = new Node(it, touched, nSkip);
+ else nodes[h] = new Node(it, touched, nSkip, nodes[h-1]->getTotBenef());
+NUM_NEW++;
+NEW_NODE++;
+ h++;
+ nSkip = 0;
+ }
+ else
+ nSkip++;
+
+ it++;
+ }
+
+ if (nodes[--h]->getTotBenef() > maxBenef) {
+ maxBenef = nodes[h]->getTotBenef();
+ itmax = nodes[0]->getItRecomb();
+ }
+
+ /*nSkip = 0;
+ while (h >= 0 && nSkip < 2) {
+ nodes[h]->erase();
+ nSkip += nodes[h--]->getnSkip();
+
+ if (h == _horizon - 2 && nodes[h]->getnSkip() >= 3) nodes[h--]->erase(); //probleme si h=-1
+ while (h >= 0 &&
+ nodes[h]->beSkipped() &&
+ nodes[h]->getnSkip() >= 2 * (_horizon - h))
+ nodes[h--]->erase();*/
+ delete nodes[h];
+NUM_DEL++;
+DEL_NODE++;
+ while (--h >= 0 &&
+ (nodes[h]->isBetter() ||
+ nodes[h]->getnSkip() >= 2 * (_horizon - h))){
+NUM_DEL++;
+DEL_NODE++;
+ delete nodes[h];
+ }
+
+ if (h < 0) end = true;
+ else {
+ it = ++nodes[h]->getItRecomb();
+ nSkip = nodes[h]->getnSkip() + 1;
+ delete nodes[h];
+NUM_DEL++;
+DEL_NODE++;
+ }
+ }
+
+ Msg::Info("num test %d", numb);
+
+ return itmax;
+}
+
+void Recombine2D::_removeImpossible(Set_Recomb::iterator it)
+{
+ std::vector<MElement *> t;
+ (*it)->triangles(t);
+
+ Map_Tri_Recomb::iterator itm;
+ std::set<RecombTriangle*>::iterator itrt;
+
+ for (int i = 0; i < t.size(); i++) {
+
+ itm = _possibleRecomb.find(t[i]);
+ if (itm == _possibleRecomb.end())
+ Msg::Info("not _possibleRecomb");
+ else {
+ itrt = itm->second.begin();
+ while (itrt != itm->second.end()) { // loop over unusable RecombTriangle
+ if (*itrt != *it)
+ _rmRT(*itrt, t[i]);
+ itrt++;
+ }
+ }
+ _possibleRecomb.erase(itm);
+ }
+
+ RecombTriangle *rt = *it;
+ if (!_lastRecomb.erase(*it))
+ _pairs.erase(it);
+ delete rt;
+NUM_DEL++;
+DEL_REC++;
+}
+
+void Recombine2D::_rmRT(RecombTriangle *rt, MElement *el)
+{
+ MElement *tri;
+ Map_Tri_Recomb::iterator itmtri;
+
+ for (int i = 0; i < rt->numTriangles(); i++)
+ if ((tri = rt->triangle(i)) != el) {
+ itmtri = _possibleRecomb.find(tri);
+ if (itmtri == _possibleRecomb.end())
+ Msg::Info("-not _possibleRecomb");
+ else {
+ itmtri->second.erase(rt);
+ switch (itmtri->second.size()) {
+ case 1 :
+ _pairs.erase(*itmtri->second.begin());
+ _lastRecomb.insert(*itmtri->second.begin());
+ break;
+ case 0 :
+ _isolated.insert(tri);
+ _possibleRecomb.erase(itmtri);
+ default :;
+ }
+ }
+ }
+ _pairs.erase(rt);
+ _lastRecomb.erase(rt);
+ delete rt;
+NUM_DEL++;
+DEL_REC++;
+}
+
+int Recombine2D::apply()
+{
+ if (!_haveParam && !_applied) return 0;
+
+ std::vector<MTriangle*> triangles2;
+ for (int i = 0; i < _gf->triangles.size(); i++) {
+ if (_isolated.find(_gf->triangles[i]) == _isolated.end())
+ delete _gf->triangles[i];
+ else
+ triangles2.push_back(_gf->triangles[i]);
+ }
+ _gf->triangles = triangles2;
+ _gf->quadrangles = _quads;
+
+ _applied = true;
+ return 1;
+}
+
+RecombTriangle::RecombTriangle(const MEdge &me, MElement *t1, MElement *t2)
+ : _formingQuad(true)
+{
+ _t.push_back(t1);
+ _t.push_back(t2);
+ _n1 = me.getVertex(0);
+ _n2 = me.getVertex(1);
+
+ if(t1->getVertex(0) != _n1 && t1->getVertex(0) != _n2) _n3 = t1->getVertex(0);
+ else if(t1->getVertex(1) != _n1 && t1->getVertex(1) != _n2) _n3 = t1->getVertex(1);
+ else if(t1->getVertex(2) != _n1 && t1->getVertex(2) != _n2) _n3 = t1->getVertex(2);
+ if(t2->getVertex(0) != _n1 && t2->getVertex(0) != _n2) _n4 = t2->getVertex(0);
+ else if(t2->getVertex(1) != _n1 && t2->getVertex(1) != _n2) _n4 = t2->getVertex(1);
+ else if(t2->getVertex(2) != _n1 && t2->getVertex(2) != _n2) _n4 = t2->getVertex(2);
+
+ double a1 = 180 * angle3Vertices(_n1, _n4, _n2) / M_PI;
+ double a2 = 180 * angle3Vertices(_n4, _n2, _n3) / M_PI;
+ double a3 = 180 * angle3Vertices(_n2, _n3, _n1) / M_PI;
+ double a4 = 180 * angle3Vertices(_n3, _n1, _n4) / M_PI;
+ _angle = fabs(90. - a1);
+ _angle = std::max(fabs(90. - a2), _angle);
+ _angle = std::max(fabs(90. - a3), _angle);
+ _angle = std::max(fabs(90. - a4), _angle);
+
+ _benefit = std::max(1. - _angle/90., .0);
+ //_benefit = compute_alignment(me,t1,t2);
+
+ //cost_alignment = Temporary::compute_alignment(me,_t1,_t2); //addition for class Temporary
+ //total_cost = Temporary::compute_total_cost(cost_quality,cost_alignment); //addition for class Temporary
+ //total_cost = 100.0*cost_quality; //addition for class Temporary
+}
+
+RecombTriangle::RecombTriangle(MElement *t, double a)
+ : _benefit(a), _formingQuad(false) {
+ _t.push_back(t);
+}
+
+double RecombTriangle::touch()
+{
+ double prec = 1e-10;
+
+ double a = (double)(rand() - RAND_MAX / 2) / RAND_MAX * prec;
+ _benefit += a;
+ return a;
+}
+
+MQuadrangle *RecombTriangle::createQuad() const
+{
+NUM_NEW++;
+NEW_QUAD++;
+ switch (_t.size()) {
+ case 1 :
+ {
+ Msg::Error("Can't make a quad from only one triangle");
+ return new MQuadrangle(_t[0]->getVertex(0),
+ _t[0]->getVertex(1),
+ _t[0]->getVertex(2),
+ _t[0]->getVertex(0));
+ }
+ case 2 :
+ {
+ int orientation = 0;
+ for (int i = 0; i < 3; i++) {
+ if (_t[0]->getVertex(i) == _n1) {
+ if (_t[0]->getVertex((i + 1) % 3) == _n2)
+ orientation = 1;
+ else
+ orientation = -1;
+ break;
+ }
+ }
+ if (orientation < 0)
+ return new MQuadrangle(_n1, _n3, _n2, _n4);
+ else
+ return new MQuadrangle(_n1, _n4, _n2, _n3);
+ }
+ default :
+ {
+ Msg::Warning("Quad creation not yet implemented for more than two triangles");
+ return new MQuadrangle(_n1, _n2, _n3, _n4);
+ }
+ }
+}
+
+double RecombTriangle::compute_alignment(const MEdge &edge, MElement *e1, MElement* e2)
+{
+ return 1 - Temporary::compute_alignment(edge, e1, e2);
+}
+
+bool RecombTriangle::operator<(const RecombTriangle &other) const
+{
+ return _benefit < other._benefit;
+}
+
+bool lessRecombTri::operator()(RecombTriangle *rt1, RecombTriangle *rt2) const
+{
+ return *rt2 < *rt1;
+}
+
+Node::Node(Set_Recomb::iterator it,
+ Map_Tri_Node &touched,
+ int nskp, double ben) : _nskip(nskp), _recomb(it)
+{
+ _benef = (*it)->getBenef();
+ _totBenef = ben + _benef;
+ for (int i = 0; i < (*it)->numTriangles(); i++)
+ _vit.push_back(touched.insert( touched.begin(), std::make_pair((*it)->triangle(i), this) ));
+ _touch = &touched;
+}
+
+Node::~Node()
+{
+ for (int i = 0; i < _vit.size(); i++)
+ _touch->erase(_vit[i]);
+}
+
+void Node::erase()
+{
+ for (int i = 0; i < _vit.size(); i++)
+ _touch->erase(_vit[i]);
+}
+
+void Node::blocking(const Map_Tri_Node::iterator &it)
+{
+ for (int i = 0; i < _vit.size(); i++) {
+ if (it == _vit[i]) {
+ _blocking.insert(i);
+ }
+ }
+}
+
+//
diff --git a/Mesh/meshRecombine2D.h b/Mesh/meshRecombine2D.h
new file mode 100644
index 0000000000000000000000000000000000000000..dc216e6f79e2225171c7e0e386ca355a8d5953c5
--- /dev/null
+++ b/Mesh/meshRecombine2D.h
@@ -0,0 +1,111 @@
+// 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 (amjohnen@gmail.com)
+//
+
+#ifndef _MESH_RECOMBINE_2D_H_
+#define _MESH_RECOMBINE_2D_H_
+
+#include <map>
+#include <set>
+#include <vector>
+#include "MVertex.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "meshGFaceOptimize.h"
+
+class RecombTriangle;
+class Node;
+struct lessRecombTri {
+ bool operator()(RecombTriangle *rt1, RecombTriangle *rt2) const;
+};
+
+typedef std::set<RecombTriangle*,lessRecombTri> Set_Recomb;
+typedef std::map<MElement*,std::set<RecombTriangle*> > Map_Tri_Recomb;
+typedef std::map<MElement*,Node*> Map_Tri_Node;
+
+class Recombine2D {
+ private :
+ GFace *_gf;
+ int _horizon;
+ double _benef;
+ bool _haveParam, _applied;
+ Set_Recomb _pairs, _lastRecomb;
+ Map_Tri_Recomb _possibleRecomb;
+ std::set<MElement*> _isolated;
+ std::vector<MQuadrangle*> _quads;
+
+ template <class E>
+ void _buildEdgeToElement(std::vector<E*> &, e2t_cont &);
+ void _recombine();
+ Set_Recomb::iterator _bestNextRecombination();
+ void _removeImpossible(Set_Recomb::iterator);
+ void _rmRT(RecombTriangle *, MElement *);
+
+ public :
+ Recombine2D(GFace*, int horizon);
+ ~Recombine2D();
+
+ int apply();
+ double getBenef() const {return _benef;}
+ int numTriangle() const {return _isolated.size();}
+};
+
+class RecombTriangle {
+ private :
+ std::vector<MElement *> _t;
+ MVertex *_n1, *_n2, *_n3, *_n4;
+ double _angle;
+ double _benefit;
+ bool _formingQuad;
+
+ public :
+ RecombTriangle(const MEdge &, MElement *, MElement *);
+ RecombTriangle(MElement *, double);
+
+ double touch();
+
+ MQuadrangle *createQuad() const;
+ bool operator<(const RecombTriangle &other) const;
+
+ void triangles(std::vector<MElement *> &v) const {v = _t;}
+ int numTriangles() const {return (int) _t.size();}
+ MElement *triangle(int i) const {
+ if (i >= 0 && i < _t.size())
+ return _t[i];
+ return NULL;
+ }
+ bool isQuad() const {return _formingQuad;}
+ double getBenef() const {return _benefit;}
+
+ double compute_alignment(const MEdge&, MElement*, MElement*);
+};
+
+class Node {
+ private :
+ double _benef, _totBenef;
+ int _nskip;
+ Set_Recomb::iterator _recomb;
+ Map_Tri_Node *_touch;
+ std::vector<Map_Tri_Node::iterator> _vit;
+ std::set<int> _blocking;
+
+ public :
+ //Node(){}
+ Node(Set_Recomb::iterator, Map_Tri_Node &, int, double ben = .0);
+ ~Node();
+
+ void erase();
+ void blocking(const Map_Tri_Node::iterator &);
+ bool isBetter() {return _blocking.size() < 2;}
+ Set_Recomb::iterator getItRecomb() {return _recomb;}
+ double getTotBenef() {return _totBenef;}
+ double getnSkip() {return _nskip;}
+};
+
+
+#endif