No commit message

No commit message

Mesh/meshRecombine2D.cpp

@@ -28,8 +28,8 @@ Recombine2D::Recombine2D(GFace *gf, int horizon)
   
   //return;
   
-  laplaceSmoothing(gf,100);
-  gf->model()->writeMSH("befSquare.msh");
+  //laplaceSmoothing(gf,100);
+  //gf->model()->writeMSH("befSquare.msh");
   
   Msg::Info("Branching with horizon %d", _horizon = HORIZ);
   
@@ -120,7 +120,7 @@ void Recombine2D::_buildEdgeToElement(std::vector<E*> &elements, e2t_cont &adj)
 
 void Recombine2D::_recombine()
 {
-  SetBoundingBox();
+  //SetBoundingBox();
   
   /*Map_Tri_Recomb::iterator itt;
   

Mesh/meshRecombine2D.h

@@ -20,12 +20,14 @@
 
 class RecombTriangle;
 class Recomb2D_Node;
-class Rec2d_node;
-class Rec2d_edge;
-class Rec2d_cycle;
+class TrianglesUnion;
+class Rec2d_vertex;
 struct lessRecombTri {
   bool operator()(RecombTriangle *rt1, RecombTriangle *rt2) const;
 };
+struct lessTriUnion {
+  bool operator()(TrianglesUnion *, TrianglesUnion *) const;
+};
 
 typedef std::set<RecombTriangle*,lessRecombTri> Set_Recomb;
 typedef std::map<MElement*,std::set<RecombTriangle*> > Map_Tri_Recomb;
@@ -55,13 +57,15 @@ class Recombine2D {
     ~Recombine2D();
     
     int apply();
-    double getBenef() const {return _benef;}
-    int numTriangle() const {return _isolated.size();}
+    int apply(bool);
+    inline double getBenef() const {return _benef;}
+    inline int numTriangle() const {return _isolated.size();}
     
   private :
-    std::set<Rec2d_node*> _nodes;
-    std::set<Rec2d_edge*> _edges;
-    std::set<Rec2d_cycle*> _cycles;
+    //std::set<TrianglesUnion*, lessTriUnion> _setQuads;
+    std::list<TrianglesUnion*> _setQuads;
+    void _removeImpossible(std::list<TrianglesUnion*>::iterator);
+    void _recombine(bool);
 };
 
 class RecombTriangle {
@@ -82,15 +86,15 @@ class RecombTriangle {
     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();}
+    inline void triangles(std::vector<MElement *> &v) const {v = _t;}
+    inline 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;}
+    inline bool isQuad() const {return _formingQuad;}
+    inline double getBenef() const {return _benefit;}
     
     double compute_alignment(const MEdge&, MElement*, MElement*);
 };
@@ -114,51 +118,71 @@ class Recomb2D_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;}
-    int getnSkip() {return _nskip;}
-    double getBound() {return _totBenef + _benef * _depth;}
+    inline bool isBetter() {return _blocking.size() < 2;}
+    inline Set_Recomb::iterator getItRecomb() {return _recomb;}
+    inline double getTotBenef() {return _totBenef;}
+    inline int getnSkip() {return _nskip;}
+    inline double getBound() {return _totBenef + _benef * _depth;}
 };
 
 
-class Rec2d_node {
+class Rec2d_vertex {
   private :
-    std::set<Rec2d_edge*> _freeEdges;
-    std::set<Rec2d_cycle*> _cycle3;
-    std::set<Rec2d_cycle*> _cycle4;
+    // _onWhat={-1:corner,0:edge,1:face,(2:volume)}
+    int _onWhat, _numEl;
+    double _angle;
+    // GEdge : virtual SVector3 firstDer(double par) const = 0;
+    // GEntity : virtual Range<double> parBounds(int i)
+    static double **_Vvalues;
+    static double **_Vbenefs;
     
-  public :
-    Rec2d_node(){}
+    void _initStaticTable();
+    double _computeAngle(MVertex *, std::list<MTriangle*> &, std::set<GEdge*> &);
   
-    void addCycle(Rec2d_cycle*, int n = 0);
-    void addEdge(Rec2d_edge*);
+  public :
+    Rec2d_vertex(MVertex *, std::list<MTriangle*> &,
+                 std::map<MVertex*, std::set<GEdge*> > &);
     
-    void print();
+    inline void changeNumEl(int c) {_numEl += c;}
+    double getReward();
+    double getReward(int);
 };
-class Rec2d_edge {
-  private :
-    Rec2d_node *_nodes[2];
-    std::set<Rec2d_cycle*> _cycles;
 
+class TrianglesUnion {
+  private :
+    int _numEmbEdge, _numEmbVert, _numBoundVert, _numTri, _computation;
+    double _embEdgeValue, _embVertValue, _globValIfExecuted;
+    Rec2d_vertex **_vertices;
+    MTriangle **_triangles;
+    static int _RECOMPUTE; //incremented when a recombination is executed
     public:
-    Rec2d_edge(Rec2d_node*, Rec2d_node*);
-    
-    void addCycle(Rec2d_cycle*);
-    
-    void print();
-};
-class Rec2d_cycle {
+    static int _NumEdge, _NumVert;
+    static double _ValEdge, _ValVert;
     private:
-    std::set<Rec2d_edge*> _edges;
     
-  public :
-    Rec2d_cycle(){}
-    void addEdge(Rec2d_edge*);
-    int size() {return _edges.size();}
+    double _computeEdgeLength(GFace *, MVertex **, double *u, double *v,
+                              int numIntermedPoint= 1);
+    double _computeAlignment(MEdge &, std::list<MTriangle*> &);
+    void _update();
     
-    void print();
+  public :
+    TrianglesUnion(GFace *, std::list<MTriangle*> &, std::list<MEdge> &,
+                   std::list<Rec2d_vertex*> &, std::map<MVertex*,double*> &);
+    
+    bool operator<(TrianglesUnion &other);
+    inline double getEdgeValue() {return _embEdgeValue;}
+    inline double getNumTriangles() {return _numTri;}
+    inline MTriangle *getTriangle(int num) {return _triangles[num];}
+    void select() {
+      _ValEdge -= _embEdgeValue;
+      _NumEdge -= _numEmbEdge;
+      _ValVert -= _embVertValue;
+      _NumVert -= _numEmbVert;
+      for (int i = 0; i < _numBoundVert; i++) {
+        _ValVert += _vertices[i]->getReward(-1);
+      }
+      _RECOMPUTE++;
+    }
+    MQuadrangle *createQuad() const;
 };
-
-
 #endif

Mesh/meshRecombine2D_2.cpp

@@ -8,6 +8,7 @@
 //
 
 #include "meshRecombine2D.h"
+#include "BackgroundMesh.h"
 #include "GFace.h"
 #include <cmath>
 #include <FL/Fl.H>
@@ -18,102 +19,541 @@
 
 static int HORIZ = 20;
 
+double **Rec2d_vertex::_Vvalues = NULL;
+double **Rec2d_vertex::_Vbenefs = NULL;
+int TrianglesUnion::_RECOMPUTE = 0;
+int TrianglesUnion::_NumEdge = 0;
+int TrianglesUnion::_NumVert = 0;
+double TrianglesUnion::_ValEdge = .0;
+double TrianglesUnion::_ValVert = .0;
+
 Recombine2D::Recombine2D(GFace *gf)
 : _horizon(0), _gf(gf), _benef(.0), _applied(true)
 {
-  std::map<MEdge, Rec2d_edge*, Less_Edge> mapEdge;
-  std::map<MVertex*, Rec2d_node*> mapNode;
+  Msg::Warning("[meshRecombine2D] Alignement computation ok uniquely for xy or yz plane mesh !");
+  
+  // 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++) {
+      if ((*itge)->getBeginVertex()->getNumMeshElements() > 1 ||
+          (*itge)->getEndVertex()->getNumMeshElements() > 1     )
+        Msg::Warning("[meshRecombine2D] Why more than one MPoint, what should I do ?");
+      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++) {
-    Rec2d_cycle *rc = new Rec2d_cycle();
-    _cycles.insert(rc);
     MTriangle *t = gf->triangles[i];
-    
     for (int j = 0; j < 3; j++) {
-      MVertex *v = t->getVertex(j);
-      std::map<MVertex*, Rec2d_node*>::iterator itn = mapNode.find(v);
-      if (itn == mapNode.end()) {
-        Rec2d_node *rn = new Rec2d_node();
-        mapNode[v] = rn;
-        _nodes.insert(rn);
-        rn->addCycle(rc, 3);
+      mapVert[t->getVertex(j)].push_back(t);
+      mapEdge[t->getEdge(j)].push_back(t);
+    }
+  }
+  
+  // Create the 'Rec2d_vertex' and store iterator to vertices which have degree 4
+  std::map<MVertex*, std::list<MTriangle*> >::iterator itvert = mapVert.begin();
+  std::map<MVertex*, Rec2d_vertex*> mapV2N;
+  std::map<MVertex*, Rec2d_vertex*>::iterator itV2N = mapV2N.begin();
+  std::list<std::map<MVertex*, std::list<MTriangle*> >::iterator> fourTri;
+  for (; itvert != mapVert.end(); itvert++) {
+    if ((*itvert).second.size() == 4)
+      fourTri.push_back(itvert);
+    Rec2d_vertex *n = new Rec2d_vertex((*itvert).first, (*itvert).second, mapCornerVert);
+    itV2N = mapV2N.insert(itV2N, std::make_pair((*itvert).first,n));
+    TrianglesUnion::_NumVert++;
+    TrianglesUnion::_ValVert += n->getReward();
+  }
+  
+  // store mesh size and parametric coordinates for better performance
+  std::map<MVertex*,double*> mapV2LcUV;
+  
+  // Create 'TrianglesUnion' for pattern of 4 triangles
+/*
+  +-----+
+  |\   /|
+  | \ / |
+  |  +  |
+  | / \ |
+  |/   \|
+  +-----+
+*/
+  std::list<std::map<MVertex*, std::list<MTriangle*> >::iterator>::iterator it4;
+  int j = 0;
+  for (it4 = fourTri.begin(); it4 != fourTri.end(); it4++) {
+    std::list<MEdge> listEmbEdge;
+    std::list<Rec2d_vertex*> listVert;
+    {
+      std::set<MVertex*> setVert;
+      std::multiset<MEdge, Less_Edge> msetEdge;
+      std::list<MTriangle*>::iterator itTri = (*(*it4)).second.begin();
+      for (; itTri != (*(*it4)).second.end(); itTri++) {
+        MVertex *vt;
+        for (int i = 0; i < 3; i++) {
+          msetEdge.insert((*itTri)->getEdge(i));
+          vt = (*itTri)->getVertex(i);
+          if (vt != (*(*it4)).first)
+            setVert.insert(vt);
+        }
+      }
+      std::multiset<MEdge>::iterator itEdge = msetEdge.begin();
+      MEdge me = *(itEdge++);
+      for (; itEdge != msetEdge.end(); itEdge++) {
+        if (*itEdge == me)
+          listEmbEdge.push_back(*itEdge);
+        me = *itEdge;
+      }
+      std::set<MVertex*>::iterator itsVert = setVert.begin();
+      for (; itsVert != setVert.end(); itsVert++)
+        listVert.push_back(mapV2N[*itsVert]);
+      listVert.push_back(mapV2N[(*(*it4)).first]);
+    }
+    Msg::Info("%d",++j);
+    _setQuads.push_back(new TrianglesUnion(gf, (*(*it4)).second, listEmbEdge, listVert, mapV2LcUV));
+  }
+  
+  // Create 'TrianglesUnion' for pattern of 2 triangles
+/*
+  +---+
+  |\  |
+  | \ |
+  |  \|
+  +---+
+*/
+  /*std::map<MEdge, std::list<MTriangle*> >::iterator itedge = mapEdge.begin();
+  for (; itedge != mapEdge.end(); itedge++) {
+    if ((*itedge).second.size() == 2) {
+      std::list<MEdge> listEmbEdge;
+      listEmbEdge.push_back((*itedge).first);
+      std::list<Rec2d_vertex*> listVert;
+      listVert.push_back(mapV2N[(*itedge).first.getVertex(0)]);
+      listVert.push_back(mapV2N[(*itedge).first.getVertex(1)]);
+      TrianglesUnion *tu = new TrianglesUnion(gf, (*itedge).second, listEmbEdge, listVert, mapV2LcUV);
+      _setQuads.push_back(tu);
+      TrianglesUnion::_NumEdge++;
+      TrianglesUnion::_ValEdge += tu->getEdgeValue();
     }
     else
-        (*itn).second->addCycle(rc, 3);
+      Msg::Info("[bord] %d", (*itedge).second.size());
   }
+  _setQuads.sort(lessTriUnion());
   
-    for (int j = 0; j < 3; j++) {
-      MEdge e = t->getEdge(j);
-      std::map<MEdge, Rec2d_edge*>::iterator ite = mapEdge.find(e);
-      if (ite == mapEdge.end()) {
-        Rec2d_node *n0 = mapNode[e.getVertex(0)];
-        Rec2d_node *n1 = mapNode[e.getVertex(1)];
-        Rec2d_edge *re = new Rec2d_edge(n0, n1);
-        n0->addEdge(re);
-        n1->addEdge(re);
-        mapEdge[e] = re;
-        _edges.insert(re);
-        rc->addEdge(re);
-        re->addCycle(rc);
+  _recombine(true);
+  _applied = true;*/
+}
+
+
+void Recombine2D::_recombine(bool a)
+{ 
+  int i = 0;
+  while (!_setQuads.empty()) {
+    Msg::Info("%d",i);
+    std::list<TrianglesUnion*>::iterator it = _setQuads.begin();
+    
+    (*it)->select();
+    
+    _quads.push_back((*it)->createQuad());
+    
+    _removeImpossible(it);
+    
+    _setQuads.sort(lessTriUnion());
+    i++;
+  }
+  Msg::Info("Done %d (%d)", i, _pairs.size());
+}
+
+MQuadrangle *TrianglesUnion::createQuad() const
+{
+  std::list<MEdge> listBoundEdge;
+  {
+    std::multiset<MEdge, Less_Edge> msetEdge;
+    for (int i = 0; i < _numTri; i++) {
+      MTriangle *t = _triangles[i];
+      for (int i = 0; i < 3; i++) {
+        msetEdge.insert(t->getEdge(i));
+      }
+    }
+    std::multiset<MEdge>::iterator itEdge = msetEdge.begin();
+    MEdge me = *(itEdge++);
+    bool precWasSame = false;
+    for (; itEdge != msetEdge.end(); itEdge++) {
+      if (*itEdge == me)
+        precWasSame = true;
+      else if (precWasSame)
+        precWasSame = false;
+      else
+        listBoundEdge.push_back(me);
+      me = *itEdge;
+    }
+    if (!precWasSame)
+      listBoundEdge.push_back(me);
+  }
+  if (listBoundEdge.size() != 4)
+    Msg::Warning("[meshRecombine2D] Not 4 edges !");
+  
+  MVertex *v1, *v2, *v3, *v4;
+  v1 = listBoundEdge.begin()->getMinVertex();
+  v2 = listBoundEdge.begin()->getMinVertex();
+  std::list<MEdge>::iterator it = listBoundEdge.begin();
+  for (it++; it != listBoundEdge.end(); it++) {
+    if (v2 == (*it).getMinVertex()) {
+      v3 = (*it).getMaxVertex();
+      listBoundEdge.erase(it);
+      break;
+    }
+    if (v2 == (*it).getMaxVertex()) {
+      v3 = (*it).getMinVertex();
+      listBoundEdge.erase(it);
+      break;
+    }
+  }
+  for (it = listBoundEdge.begin(); it != listBoundEdge.end(); it++) {
+    if (v3 == (*it).getMinVertex()) {
+      v4 = (*it).getMaxVertex();
+      break;
+    }
+    if (v3 == (*it).getMaxVertex()) {
+      v4 = (*it).getMinVertex();
+      break;
+    }
+  }
+  
+  return new MQuadrangle(v1, v2, v3, v4);
+}
+
+
+int Recombine2D::apply(bool a)
+{
+  if (!_haveParam || _applied) return 0;
+  
+  std::vector<MTriangle*> triangles2;
+  for (int i = 0; i < _gf->triangles.size(); i++) {
+    delete _gf->triangles[i];
+  }
+  _gf->triangles = triangles2;
+  _gf->quadrangles = _quads;
+  
+  _applied = true;
+  return 1;
+}
+
+void Recombine2D::_removeImpossible(std::list<TrianglesUnion*>::iterator it)
+{
+  std::set<MTriangle *> touched;
+  
+  for (int i = 0; i < (*it)->getNumTriangles(); i++) {
+    touched.insert((*it)->getTriangle(i));
+  }
+  
+  for (it = _setQuads.begin(); it != _setQuads.end();) {
+    bool b = false;
+    for (int i = 0; i < (*it)->getNumTriangles(); i++) {
+      if (touched.find((*it)->getTriangle(i)) != touched.end())
+        b = true;
+    }
+    if (b)
+      _setQuads.erase(it++);
+    else
+      it++;
+  }
+}
+
+
+
+Rec2d_vertex::Rec2d_vertex(MVertex *v,
+                           std::list<MTriangle*> &setTri,
+                           std::map<MVertex*, std::set<GEdge*> > &mapVert)
+: _numEl(setTri.size()), _angle(.0)
+{
+  _onWhat = v->onWhat()->dim() - 1;
+  
+  if (_onWhat < 0) {
+    std::map<MVertex*, std::set<GEdge*> >::iterator it = mapVert.find(v);
+    if (it != mapVert.end()) {
+      _angle = _computeAngle(v, setTri, (*it).second);
     }
     else {
-        rc->addEdge((*ite).second);
-        (*ite).second->addCycle(rc);
+      Msg::Warning("[meshRecombine2D] Can't compute corner angle, setting angle to pi/2");
+      _angle = M_PI / 2.;
     }
   }
+  
+  if (_Vvalues == NULL)
+    _initStaticTable();
+}
+
+void Rec2d_vertex::_initStaticTable()
+{
+  if (_Vvalues == NULL) {
+    // _Vvalues[onWhat][nEl]; onWhat={0:edge,1:face} / nEl=_numEl-1
+    // _Vbenefs[onWhat][nEl]; onWhat={0:edge,1:face} / nEl=_numEl-1 (benef of adding an element)
+    int nE = 5, nF = 10;
+    
+    _Vvalues = new double*[2];
+    _Vvalues[0] = new double[nE]; //AVERIFIER SI CA MARCHE
+    for (int i = 0; i < nE; i++)
+      _Vvalues[0][i] = 1. - fabs(2. / (i+1) - 1.);
+    _Vvalues[1] = new double[nF];
+    for (int i = 0; i < nF; i++)
+      _Vvalues[1][i] = 1. - fabs(4. / (i+1) - 1.);
+      
+    _Vbenefs = new double*[2];
+    _Vbenefs[0] = new double[nE-1];
+    for (int i = 0; i < nE-1; i++)
+      _Vbenefs[0][i] = _Vvalues[0][i+1] - _Vvalues[0][i];
+    _Vbenefs[1] = new double[nF-1];
+    for (int i = 0; i < nF-1; i++)
+      _Vbenefs[1][i] = _Vvalues[1][i+1] - _Vvalues[1][i];
   }
-  }ZzZ
 }
 
-void Rec2d_node::addCycle(Rec2d_cycle *c, int n)
+
+double Rec2d_vertex::_computeAngle(MVertex *v,
+                                   std::list<MTriangle*> &setTri,
+                                   std::set<GEdge*> &setEdge)
 {
-  if (n == 3) {
-    _cycle3.insert(c);
-    return;
+  if (setEdge.size() != 2) {
+    Msg::Warning("[meshRecombine2D] Wrong number of edge : %d, returning pi/2", setEdge.size());
+    return M_PI / 2.;
   }
-  if (n == 4) {
-    _cycle4.insert(c);
-    return;
+  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 = vectub;
+    else if (normub > prec * normlb)
+      firstDer1 = vectlb;
+    else {
+      Msg::Warning("[meshRecombine2D] Bad precision, returning pi/2");
+      return M_PI / 2.;
     }
-  if (c->size() == 3) {
-    _cycle3.insert(c);
-    return;
+    if (k == 0)
+      firstDer0 = firstDer1;
   }
-  if (c->size() == 4) {
-    _cycle4.insert(c);
-    return;
+  
+  Msg::Info("-- %lf, %lf, %lf", norm(firstDer0), norm(firstDer1), dot(firstDer0, firstDer1));
+  double n = 1 / norm(firstDer0);
+  firstDer0 = firstDer0 * n;
+  n = 1 / norm(firstDer1);
+  firstDer1 = firstDer1 * n;
+  Msg::Info("--Angles : %lf, %lf, %lf", norm(firstDer0), norm(firstDer1), dot(firstDer0, firstDer1));
+  
+  double angle1 = acos(dot(firstDer0, firstDer1));
+  double angle2 = 2. * M_PI - angle1;
+  
+  double angleMesh = .0;
+  std::list<MTriangle*>::iterator it2 = setTri.begin();
+  for (; it2 != setTri.end(); it2++) {
+    MTriangle *t = *it2;
+    int k = 0;
+    while (t->getVertex(k) != v && k < 3)
+      k++;
+    if (k == 3) {
+      Msg::Warning("[meshRecombine2D] Wrong triangle, returning pi/2");
+      return M_PI / 2.;
+    }
+    angleMesh += angle3Vertices(t->getVertex((k+1)%3), v, t->getVertex((k+2)%3));
   }
-  Msg::Error("[Recombine2D] Can't tell if it's a cycle3 or a cycle4");
+  
+  Msg::Info("Angles : %lf, %lf, %lf", angle1, angleMesh, angle2);
+  
+  return M_PI / 2.;
+  
+  if (angleMesh < M_PI)
+    return angle1;
+  return angle2;
 }
-void Rec2d_node::addEdge(Rec2d_edge *e)
+
+double Rec2d_vertex::getReward()
 {
-  _freeEdges.insert(e); //que faire ?
+  if (_onWhat > -1)
+    return _Vvalues[_onWhat][_numEl-1];
+  else
+    return 1. - fabs(2. / M_PI * _angle / _numEl - 1.);
 }
-void Rec2d_node::print()
+
+double Rec2d_vertex::getReward(int n)
 {
-  Msg::Info("Node : %d cycle3, %d cycle4, %d edges", _cycle3.size(), _cycle4.size(), _freeEdges.size());
+  if (n == 0) return .0;
+  if (_onWhat > -1) {
+    if (n > 0)
+      return _Vbenefs[_onWhat][_numEl-1];
+    else
+      return - _Vbenefs[_onWhat][_numEl-2];
   }
+  else
+    return fabs(2./M_PI*_angle/_numEl - 1.) - fabs(2./M_PI*_angle/(_numEl + n) - 1.);
+}
+
 
-Rec2d_edge::Rec2d_edge(Rec2d_node *n0, Rec2d_node *n1)
+TrianglesUnion::TrianglesUnion(GFace *gf,
+                               std::list<MTriangle*> &t,
+                               std::list<MEdge> &e,
+                               std::list<Rec2d_vertex*> &v,
+                               std::map<MVertex*,double*> &v2lcUV)
 {
-  _nodes[0] = n0;
-  _nodes[1] = n1;
+  _numTri = t.size();
+  _numEmbEdge = e.size();
+  _numBoundVert = v.size() == 2 ? 2 : 4;
+  _numEmbVert = v.size() - _numBoundVert;
+  _globValIfExecuted = _embEdgeValue = _embVertValue = .0;
+  
+  _triangles = new MTriangle*[_numTri];
+  std::list<MTriangle*>::iterator itt = t.begin();
+  for (int k = 0; itt != t.end(); itt++, k++) {
+    _triangles[k] = *itt;
+  }
+  
+  std::list<MEdge>::iterator ite = e.begin();
+  for (; ite != e.end(); ite++) {
+    double sumlc = .0, u[2], v[2];
+    MVertex *vert[2];
+    for (int i = 0; i < 2; i++) {
+      vert[i] = (*ite).getVertex(i);
+      std::map<MVertex*,double*>::iterator itlc;
+      if ((itlc = v2lcUV.find(vert[i])) == v2lcUV.end()) {
+        double *a = new double[3];
+        gf->XYZtoUV(vert[i]->x(), vert[i]->y(), vert[i]->z(), a[1], a[2], 1.);
+        sumlc += a[0] = BGM_MeshSize(gf, a[1], a[2], vert[i]->x(), vert[i]->y(), vert[i]->z());
+        u[i] = a[1];
+        v[i] = a[2];
+        v2lcUV[vert[i]] = a;
+      }
+      else {
+        sumlc += (*itlc).second[0];
+        u[i] = (*itlc).second[1];
+        v[i] = (*itlc).second[2];
+      }
     }
-void Rec2d_edge::addCycle(Rec2d_cycle *c)
+
+    double length = _computeEdgeLength(gf, vert, u, v, 0);
+    _embEdgeValue += length / sumlc * (1 + _computeAlignment(*ite, t));
+    Msg::Info("Edge a : %lf/%lf = %lf <-> %lf", length, sumlc / 2, 2 * length / sumlc, _computeAlignment(*ite, t));
+  }
+  
+  _vertices = new Rec2d_vertex*[_numBoundVert];
+  std::list<Rec2d_vertex*>::iterator itv = v.begin();
+  for (int k = 0; itv != v.end() && k < _numBoundVert; itv++, k++) {
+    _globValIfExecuted += (*itv)->getReward(-1);
+    _vertices[k] = *itv;
+  }
+  for (_numEmbVert = 0; itv != v.end(); itv++, _numEmbVert++)
+    _embVertValue += (*itv)->getReward();
+  
+  _computation = _RECOMPUTE - 1;
+  _update();
+}
+
+
+double TrianglesUnion::_computeEdgeLength(GFace *gf, MVertex **vert,
+                                          double *u, double *v, int n)
 {
-  _cycles.insert(c);
+  double dx, dy, dz, l = .0;
+  if (n == 0) {
+    dx = vert[1]->x() - vert[0]->x();
+    dy = vert[1]->y() - vert[0]->y();
+    dz = vert[1]->z() - vert[0]->z();
+    l = sqrt(dx * dx + dy * dy + dz * dz);
   }
-void Rec2d_edge::print()
+  else if (n == 1) {
+    double edgeCenter[2] ={(u[0] + u[1]) * .5, (v[0] + v[1]) * .5};
+    GPoint GP = gf->point (edgeCenter[0], edgeCenter[1]);
+    dx = vert[1]->x() - GP.x();
+    dy = vert[1]->y() - GP.y();
+    dz = vert[1]->z() - GP.z();
+    l += sqrt(dx * dx + dy * dy + dz * dz);
+    dx = vert[0]->x() - GP.x();
+    dy = vert[0]->y() - GP.y();
+    dz = vert[0]->z() - GP.z();
+    l += sqrt(dx * dx + dy * dy + dz * dz);
+  }
+  else {
+    Msg::Warning("[meshRecombine2D] edge length computation not implemented for n=%d, returning 0", n);
+    return .0;
+  }
+  return l;
+}
+double TrianglesUnion::_computeAlignment(MEdge &e, std::list<MTriangle*> &lt)
 {
-  Msg::Info("Edge : %d cycles", _cycles.size());
+  std::list<MTriangle*>::iterator itlt = lt.begin();
+  if (lt.size() == 2)
+    return Temporary::compute_alignment(e, *itlt, *(++itlt));
+  
+  MVertex *v0 = e.getVertex(0);
+  MVertex *v1 = e.getVertex(1);
+  MTriangle *tris[2];
+  int k = 0;
+  for (; itlt != lt.end(); itlt++) {
+    MTriangle *t = *itlt;
+    if ((t->getVertex(0) == v0 || t->getVertex(1) == v0 || t->getVertex(2) == v0) &&
+        (t->getVertex(0) == v1 || t->getVertex(1) == v1 || t->getVertex(2) == v1)   )
+      tris[k++] = t;
+  }
+  if (k < 2 || k > 2) {
+    Msg::Warning("[meshRecombine2D] error in alignment computation, returning 0");
+    return .0;
   }
+  return Temporary::compute_alignment(e, tris[0], tris[1]);
+}
+void TrianglesUnion::_update()
+{
+  if (_computation >= _RECOMPUTE)
+    return;
+  double alpha = _ValVert - _embVertValue;
+  for (int i = 0; i < _numBoundVert; i++) {
+    alpha += _vertices[i]->getReward(-1);
+  }
+  alpha /= _NumVert - _numEmbVert;
+  double beta = (_ValEdge - _embEdgeValue) / (_NumEdge - _numEmbEdge);
+  _globValIfExecuted = alpha * beta;
   
-void Rec2d_cycle::addEdge(Rec2d_edge *e)
+  _computation = _RECOMPUTE;
+}
+bool TrianglesUnion::operator<(TrianglesUnion &other)
 {
-  _edges.insert(e);
+  _update();
+  other._update();
+  return _globValIfExecuted < other._globValIfExecuted; 
 }
-void Rec2d_cycle::print()
+
+bool lessTriUnion::operator()(TrianglesUnion *tu1, TrianglesUnion *tu2) const
 {
-  Msg::Info("Cycle : %d edges", _edges.size());
+  return *tu1 < *tu2;
 }
 
+/*map tri to recomb
+//map edge to double value (constructor)
+set Recomb sorted
+function best(copy vertex*[], copy set Recomb sorted);
+
+construction :
+ - list of vertex
+ - set of recomb
+ - map tri to recomb
+ 
+destruction :
+ - vertices, recomb
+ 
+execution :
+ - take best recomb
+ - determine recomb to execute
+ - maj E N e n
+ - reduce maptrirecomb, setrecomb
+ - sort setrecomb*/