diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index c15f629ead80b1bf809726cc230510823bc686a0..99b7150cd59253b9efb82c0107f1971c60e65d37 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -489,85 +489,73 @@ static void Mesh3D(GModel *m)
   FindConnectedRegions(delaunay, connected);
 
   // remove quads elements for volumes that are recombined
-  // pragma OMP ICI ??
+  // pragma OMP here ?
   for(unsigned int i = 0; i < connected.size(); i++){
-    for(unsigned j=0;j<connected[i].size();j++){
+    for(unsigned j = 0; j < connected[i].size(); j++){
       GRegion *gr = connected[i][j];
       if(CTX::instance()->mesh.recombine3DAll || gr->meshAttributes.recombine3D){
         std::list<GFace*> f = gr->faces();
-        for (std::list<GFace*>::iterator it = f.begin();
-            it != f.end() ; ++it) quadsToTriangles (*it,1000000);
+        for(std::list<GFace*>::iterator it = f.begin(); it != f.end() ; ++it)
+          quadsToTriangles(*it, 1000000);
       }
     }
   }
 
-  // pragma OMP ICI ??
-  double time_recombination, vol_element_recombination, vol_hexa_recombination;
-  int nb_elements_recombination,nb_hexa_recombination;
-  time_recombination = vol_hexa_recombination = vol_element_recombination = 0.;
-  nb_elements_recombination = nb_hexa_recombination = 0;
+  double time_recombination = 0., vol_element_recombination = 0.;
+  double vol_hexa_recombination = 0.;
+  int nb_elements_recombination = 0, nb_hexa_recombination = 0;
+
+  // pragma OMP here ?
   for(unsigned int i = 0; i < connected.size(); i++){
     MeshDelaunayVolume(connected[i]);
 
-    //Additional code for hex mesh begin
-    for(unsigned j=0;j<connected[i].size();j++){
+    // additional code for experimental hex mesh
+    for(unsigned j = 0; j < connected[i].size(); j++){
       GRegion *gr = connected[i][j];
-      //R-tree
-      bool treat_region_ok=false;
+      bool treat_region_ok = false;
       if(CTX::instance()->mesh.algo3d == ALGO_3D_RTREE){
-        // PEB MODIF
         if (old_algo_hexa()){
           Filler f;
           f.treat_region(gr);
-          treat_region_ok=true;
+          treat_region_ok = true;
         }
         else{
           Filler3D f;
           treat_region_ok = f.treat_region(gr);
         }
-        // END PEB MODIF
       }
-      //Recombine3D into hex
-      if((CTX::instance()->mesh.recombine3DAll || gr->meshAttributes.recombine3D) && treat_region_ok){
-        clock_t a = clock();
-        //        Recombinator rec;
-        //        rec.execute();
+      if(treat_region_ok && (CTX::instance()->mesh.recombine3DAll ||
+                             gr->meshAttributes.recombine3D)){
+        double a = Cpu();
         Recombinator rec;
         rec.execute(gr);
-//                Recombinator_Graph rec(1.e7,"test");
-//                rec.execute(gr);
-        //        Supplementary sup;
-        //        sup.execute();
-        //        PostOp post;
-        //        post.execute(0);
         Supplementary sup;
         sup.execute(gr);
         PostOp post;
         post.execute(gr,0);
-
         nb_elements_recombination += post.get_nb_elements();
         nb_hexa_recombination += post.get_nb_hexahedra();
         vol_element_recombination += post.get_vol_elements();
         vol_hexa_recombination += post.get_vol_hexahedra();
-
-        // -- partial export
+        // partial export
         stringstream ss;
         ss << "yamakawa_part_";
         ss << gr->tag();
         ss << ".msh";
-        export_gregion_mesh(gr,ss.str().c_str());
-        // -- end partial export
-
-        time_recombination += (clock() - a) / (double)CLOCKS_PER_SEC;
+        export_gregion_mesh(gr, ss.str().c_str());
+        time_recombination += (Cpu() - a);
       }
     }
   }
+
   if(CTX::instance()->mesh.recombine3DAll){
-    cout << "RECOMBINATION timing:" << endl;
-    cout << "  ------- CUMULATIVE TIME RECOMBINATION : " << time_recombination << " s." << endl;
-    cout << "RECOMBINATION CUMULATIVE STATISTICS:" << endl;
-    cout << "............................  Percentage of hexahedra   (#) : " << nb_hexa_recombination*100./nb_elements_recombination << endl;
-    cout << "............................  Percentage of hexahedra (Vol) : " << vol_hexa_recombination*100./vol_element_recombination << endl;
+    Msg::Info("RECOMBINATION timing:");
+    Msg::Info(" --- CUMULATIVE TIME RECOMBINATION : %g s.", time_recombination);
+    Msg::Info("RECOMBINATION CUMULATIVE STATISTICS:");
+    Msg::Info(".... Percentage of hexahedra   (#) : %g",
+              nb_hexa_recombination*100./nb_elements_recombination);
+    Msg::Info(".... Percentage of hexahedra (Vol) : %g",
+              vol_hexa_recombination*100./vol_element_recombination);
   }
 
   // Ensure that all volume Jacobians are positive
diff --git a/Mesh/pointInsertionRTreeTools.h b/Mesh/pointInsertionRTreeTools.h
index 32daa4e422b4add31fcfab022c0883bc469ea400..22cf1ffe6aeb252f8290f6b1c3988b26378c49f9 100644
--- a/Mesh/pointInsertionRTreeTools.h
+++ b/Mesh/pointInsertionRTreeTools.h
@@ -192,13 +192,14 @@ class compareSmoothnessVertexPairs
   private:
     const double threshold;
     double roundit(const double &d)const{
-      return (round(d/threshold)*threshold);
+      //return (round(d/threshold)*threshold);
+      return ((int)(d/threshold + 0.5)*threshold);
     }
   public:
     static bool vectorial;
 
     compareSmoothnessVertexPairs():threshold(0.05){};
-    
+
     inline bool operator () (const smoothness_vertex_pair *a, const smoothness_vertex_pair *b)  const
     {
 
diff --git a/Mesh/yamakawa.cpp b/Mesh/yamakawa.cpp
index 5df5c2dd4c7a079effebd5432a8e253531c63c27..24382028e6ec1866640b7b66cc4ead48b25c8d12 100644
--- a/Mesh/yamakawa.cpp
+++ b/Mesh/yamakawa.cpp
@@ -6,13 +6,19 @@
 // Contributor(s):
 //   Tristan Carrier
 
+// FIXME: This code should be cleaned up and made to conform with Gmsh's coding
+// style
+
+#include <numeric>
 #include <iterator>
+#include <fstream>
+#include <sstream>
+#include <math.h>
 #include "yamakawa.h"
 #include "GModel.h"
+#include "OS.h"
 #include "MVertex.h"
 #include "MElement.h"
-#include <fstream>
-#include <sstream>
 #include "MHexahedron.h"
 #include "MQuadrangle.h"
 #include "MTriangle.h"
@@ -20,54 +26,61 @@
 #include "MPrism.h"
 #include "MTetrahedron.h"
 #include "MHexahedron.h"
-#include <numeric>
-#include <math.h>
-
 
-
-
-void export_gregion_mesh(GRegion *gr, string filename){
+void export_gregion_mesh(GRegion *gr, string filename)
+{
+  // FIXME: why not use MElement::writeMSH !?
 
   // create set of all tets
   map<MVertex*,int> vertices;
   int counterv=1;
-  
-  for (vector<MTetrahedron*>::iterator it = gr->tetrahedra.begin();it!=gr->tetrahedra.end();it++){
-      for (int i=0;i<(*it)->getNumVertices();i++){
-        vertices.insert(make_pair((*it)->getVertex(i),counterv));
-        counterv++;
-      }
+
+  for (vector<MTetrahedron*>::iterator it = gr->tetrahedra.begin();
+       it!=gr->tetrahedra.end(); it++){
+    for (int i=0;i<(*it)->getNumVertices();i++){
+      vertices.insert(make_pair((*it)->getVertex(i),counterv));
+      counterv++;
+    }
   }
-  for (vector<MHexahedron*>::iterator it = gr->hexahedra.begin();it!=gr->hexahedra.end();it++){
-      for (int i=0;i<(*it)->getNumVertices();i++){
-        vertices.insert(make_pair((*it)->getVertex(i),counterv));
-        counterv++;
-      }
+  for (vector<MHexahedron*>::iterator it = gr->hexahedra.begin();
+       it!=gr->hexahedra.end(); it++){
+    for (int i=0;i<(*it)->getNumVertices();i++){
+      vertices.insert(make_pair((*it)->getVertex(i),counterv));
+      counterv++;
+    }
   }
-  for (vector<MPrism*>::iterator it = gr->prisms.begin();it!=gr->prisms.end();it++){
-      for (int i=0;i<(*it)->getNumVertices();i++){
-        vertices.insert(make_pair((*it)->getVertex(i),counterv));
-        counterv++;
-      }
+  for (vector<MPrism*>::iterator it = gr->prisms.begin();
+       it!=gr->prisms.end(); it++){
+    for (int i=0;i<(*it)->getNumVertices();i++){
+      vertices.insert(make_pair((*it)->getVertex(i),counterv));
+      counterv++;
+    }
   }
-  for (vector<MPyramid*>::iterator it = gr->pyramids.begin();it!=gr->pyramids.end();it++){
-      for (int i=0;i<(*it)->getNumVertices();i++){
-        vertices.insert(make_pair((*it)->getVertex(i),counterv));
-        counterv++;
-      }
+  for (vector<MPyramid*>::iterator it = gr->pyramids.begin();
+       it!=gr->pyramids.end(); it++){
+    for (int i=0;i<(*it)->getNumVertices();i++){
+      vertices.insert(make_pair((*it)->getVertex(i),counterv));
+      counterv++;
+    }
   }
 
   // export mesh
   ofstream out(filename.c_str());
-  out << "$MeshFormat" << endl << "2.2 0 8" << endl << "$EndMeshFormat" << endl << "$Nodes" << endl << vertices.size() << endl;
+  out << "$MeshFormat" << endl << "2.2 0 8" << endl << "$EndMeshFormat"
+      << endl << "$Nodes" << endl << vertices.size() << endl;
   // write vertices
-  for (map<MVertex*,int>::iterator it = vertices.begin();it!=vertices.end();it++)
-    out << it->second << " " << it->first->x() << " " << it->first->y() << " " << it->first->z() << endl;
-  out << "$EndNodes" << endl << "$Elements" << endl << (gr->tetrahedra.size()+gr->hexahedra.size()+gr->prisms.size()+gr->pyramids.size()) << endl;
-  
+  for (map<MVertex*,int>::iterator it = vertices.begin();
+       it!=vertices.end(); it++)
+    out << it->second << " " << it->first->x() << " " << it->first->y()
+        << " " << it->first->z() << endl;
+  out << "$EndNodes" << endl << "$Elements" << endl
+      << (gr->tetrahedra.size() + gr->hexahedra.size() +
+          gr->prisms.size() + gr->pyramids.size()) << endl;
+
   // write elems
-  int counter=1;
-  for (vector<MTetrahedron*>::iterator it = gr->tetrahedra.begin();it!=gr->tetrahedra.end();it++){
+  int counter = 1;
+  for (vector<MTetrahedron*>::iterator it = gr->tetrahedra.begin();
+       it!=gr->tetrahedra.end();it++){
     out << counter << " 4 2 0 26";
     for (int i=0;i<(*it)->getNumVertices();i++){
       MVertex *v = (*it)->getVertex(i);
@@ -76,7 +89,8 @@ void export_gregion_mesh(GRegion *gr, string filename){
     out << endl;
     counter++;
   }
-  for (vector<MHexahedron*>::iterator it = gr->hexahedra.begin();it!=gr->hexahedra.end();it++){
+  for (vector<MHexahedron*>::iterator it = gr->hexahedra.begin();
+       it!=gr->hexahedra.end();it++){
     out << counter << " 5 2 0 26";
     for (int i=0;i<(*it)->getNumVertices();i++){
       MVertex *v = (*it)->getVertex(i);
@@ -85,7 +99,8 @@ void export_gregion_mesh(GRegion *gr, string filename){
     out << endl;
     counter++;
   }
-  for (vector<MPrism*>::iterator it = gr->prisms.begin();it!=gr->prisms.end();it++){
+  for (vector<MPrism*>::iterator it = gr->prisms.begin();
+       it!=gr->prisms.end();it++){
     out << counter << " 6 2 0 26";
     for (int i=0;i<(*it)->getNumVertices();i++){
       MVertex *v = (*it)->getVertex(i);
@@ -94,7 +109,8 @@ void export_gregion_mesh(GRegion *gr, string filename){
     out << endl;
     counter++;
   }
-  for (vector<MPyramid*>::iterator it = gr->pyramids.begin();it!=gr->pyramids.end();it++){
+  for (vector<MPyramid*>::iterator it = gr->pyramids.begin();
+       it!=gr->pyramids.end();it++){
     out << counter << " 7 2 0 26";
     for (int i=0;i<(*it)->getNumVertices();i++){
       MVertex *v = (*it)->getVertex(i);
@@ -103,20 +119,15 @@ void export_gregion_mesh(GRegion *gr, string filename){
     out << endl;
     counter++;
   }
-  
-  
-  
+
   out << "$EndElements" << endl;
   out.close();
-
-
 }
 
-
-
 template <class T>
-void export_the_clique_graphviz_format(cliques_compatibility_graph<T> &cl, int clique_number,string filename){
-
+void export_the_clique_graphviz_format(cliques_compatibility_graph<T> &cl,
+                                       int clique_number,string filename)
+{
   ofstream out(filename.c_str());
   out << "Graph G {" << endl;
   typename multimap<int,set<T> >::reverse_iterator it_all = cl.allQ.rbegin();
@@ -131,7 +142,6 @@ void export_the_clique_graphviz_format(cliques_compatibility_graph<T> &cl, int c
   //typedef tr1::unordered_multimap<hash_key, T> graph_data;
   //typedef tr1::unordered_multimap<hash_key, pair<T, graph_data > > graph;
 
-
   int counter=1;
   map<T,int> visited_hex;
   multimap<int,int> done;
@@ -155,9 +165,6 @@ void export_the_clique_graphviz_format(cliques_compatibility_graph<T> &cl, int c
     else
       num1 = itfind->second;
 
-
-
-
     itgraphdata = itgraph->second.second.begin();
     for(;itgraphdata!=itgraph->second.second.end();itgraphdata++){
       T secondhex = itgraphdata->second;
@@ -173,7 +180,8 @@ void export_the_clique_graphviz_format(cliques_compatibility_graph<T> &cl, int c
 
       // search if num1 - num2 has already been written...
       bool found=false;
-      pair<multimap<int,int>::iterator, multimap<int,int>::iterator> range = done.equal_range(num1);
+      pair<multimap<int,int>::iterator, multimap<int,int>::iterator> range =
+        done.equal_range(num1);
       for(multimap<int,int>::iterator it = range.first;it!=range.second;it++){
         if (it->second==num2){
           found=true;
@@ -200,33 +208,29 @@ void export_the_clique_graphviz_format(cliques_compatibility_graph<T> &cl, int c
       int num2=0;
       num2=counter;
       visited_hex[itfind->first] = counter++;
-      out << num2 << " [shape=circle, style=filled, fillcolor=red];" << endl; 
+      out << num2 << " [shape=circle, style=filled, fillcolor=red];" << endl;
     }
     else
-      out << itfind->second << " [shape=circle, style=filled, fillcolor=red];" << endl; 
+      out << itfind->second << " [shape=circle, style=filled, fillcolor=red];" << endl;
   }
 
-
-
   out << "}" << endl;
   out.close();
 }
 
-
-
-
-
 template<class T>
-clique_stop_criteria<T>::clique_stop_criteria(map<T, std::set<MElement*> > &_m, int _i):hex_to_tet(_m),total_number_tet(_i){
+clique_stop_criteria<T>::clique_stop_criteria(map<T, std::set<MElement*> > &_m, int _i)
+  : hex_to_tet(_m),total_number_tet(_i)
+{
 };
 
-
 template<class T>
 clique_stop_criteria<T>::~clique_stop_criteria(){};
 
-
 template<class T>
-void clique_stop_criteria<T>::export_corresponding_mesh(const graph_data_no_hash &clique)const{
+void clique_stop_criteria<T>::export_corresponding_mesh
+(const graph_data_no_hash &clique) const
+{
   // NB: fct not often called...
 
   // filename
@@ -271,15 +275,19 @@ void clique_stop_criteria<T>::export_corresponding_mesh(const graph_data_no_hash
     hexs.insert(h);
   }
 
-  // export mesh
+  // export mesh FIXME: why not use MElement::writeMSH
   ofstream out(filename.c_str());
   ofstream outtets(filenametets.c_str());
-  out << "$MeshFormat" << endl << "2.2 0 8" << endl << "$EndMeshFormat" << endl << "$Nodes" << endl << vertices.size() << endl;
-  outtets << "$MeshFormat" << endl << "2.2 0 8" << endl << "$EndMeshFormat" << endl << "$Nodes" << endl << vertices.size() << endl;
+  out << "$MeshFormat" << endl << "2.2 0 8" << endl << "$EndMeshFormat"
+      << endl << "$Nodes" << endl << vertices.size() << endl;
+  outtets << "$MeshFormat" << endl << "2.2 0 8" << endl << "$EndMeshFormat"
+          << endl << "$Nodes" << endl << vertices.size() << endl;
   // write vertices
   for (map<MVertex*,int>::iterator it = vertices.begin();it!=vertices.end();it++){
-    out << it->second << " " << it->first->x() << " " << it->first->y() << " " << it->first->z() << endl;
-    outtets << it->second << " " << it->first->x() << " " << it->first->y() << " " << it->first->z() << endl;
+    out << it->second << " " << it->first->x() << " "
+        << it->first->y() << " " << it->first->z() << endl;
+    outtets << it->second << " " << it->first->x() << " "
+            << it->first->y() << " " << it->first->z() << endl;
   }
   out << "$EndNodes" << endl << "$Elements" << endl << (hexs.size()+tets.size()) << endl;
   outtets << "$EndNodes" << endl << "$Elements" << endl << (hexs.size()+tets.size()) << endl;
@@ -313,13 +321,12 @@ void clique_stop_criteria<T>::export_corresponding_mesh(const graph_data_no_hash
   out.close();
   outtets << "$EndElements" << endl;
   outtets.close();
-
-
 }
 
 
 template<class T>
-bool clique_stop_criteria<T>::stop(const graph_data_no_hash &clique)const{
+bool clique_stop_criteria<T>::stop(const graph_data_no_hash &clique)const
+{
   unsigned int total = 0;
 
   // need to detect slivers !!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -350,7 +357,7 @@ bool clique_stop_criteria<T>::stop(const graph_data_no_hash &clique)const{
 
   }
 
-  // to be sure, adding volume criteria... 
+  // to be sure, adding volume criteria...
   vector<double> volumes;
   for (set<MElement*>::iterator it = thetets.begin();it!=thetets.end();it++){
     volumes.push_back((*it)->getVolume());
@@ -374,21 +381,28 @@ bool clique_stop_criteria<T>::stop(const graph_data_no_hash &clique)const{
 
 
 template<class T>
-cliques_compatibility_graph<T>::cliques_compatibility_graph(graph &_g, const map<T, std::vector<double> > &_hex_ranks, unsigned int _max_nb_cliques, unsigned int _nb_hex_potentiels, clique_stop_criteria<T> *csc, ptrfunction_export fct):debug(false), max_nb_cliques(_max_nb_cliques), nb_hex_potentiels(_nb_hex_potentiels), max_clique_size(0), position(0), total_nodes_number(0), cancel_search(false), hex_ranks(_hex_ranks), G(_g),criteria(csc),export_clique_graph(fct),found_the_ultimate_max_clique(false){
-
+cliques_compatibility_graph<T>::cliques_compatibility_graph
+(graph &_g, const map<T, std::vector<double> > &_hex_ranks,
+ unsigned int _max_nb_cliques, unsigned int _nb_hex_potentiels,
+ clique_stop_criteria<T> *csc, ptrfunction_export fct)
+  : debug(false), max_nb_cliques(_max_nb_cliques),
+    nb_hex_potentiels(_nb_hex_potentiels), max_clique_size(0), position(0),
+    total_nodes_number(0), cancel_search(false), hex_ranks(_hex_ranks),
+    G(_g),criteria(csc),export_clique_graph(fct),
+    found_the_ultimate_max_clique(false)
+{
   total_nb_of_cliques_searched = 0;
   max_nb_of_stored_cliques = 10;
-
 };
 
-
 template<class T>
-cliques_compatibility_graph<T>::~cliques_compatibility_graph(){
+cliques_compatibility_graph<T>::~cliques_compatibility_graph()
+{
 };
 
-
 template<class T>
-void cliques_compatibility_graph<T>::find_cliques(){
+void cliques_compatibility_graph<T>::find_cliques()
+{
   // init
   graph_data s;
   for (typename graph::iterator it = G.begin();it!=G.end();it++){
@@ -402,9 +416,9 @@ void cliques_compatibility_graph<T>::find_cliques(){
 
 };
 
-
 template<class T>
-void cliques_compatibility_graph<T>::export_cliques(){
+void cliques_compatibility_graph<T>::export_cliques()
+{
   typename multimap<int, set<T> >::reverse_iterator itstore = allQ.rbegin();
   for (;itstore!=allQ.rend();itstore++){
     cout << "clique of size " << itstore->first << ": {";
@@ -417,11 +431,10 @@ void cliques_compatibility_graph<T>::export_cliques(){
 
 
 template<class T>
-void cliques_compatibility_graph<T>::store_clique(int n){
-
+void cliques_compatibility_graph<T>::store_clique(int n)
+{
   total_nb_of_cliques_searched++;
 
-
   if(total_nb_of_cliques_searched%10000==0){
     if (max_nb_cliques>0)
       cout << "found " << total_nb_of_cliques_searched << " cliques on " << max_nb_cliques << endl << flush;
@@ -429,13 +442,11 @@ void cliques_compatibility_graph<T>::store_clique(int n){
       cout << "found " << total_nb_of_cliques_searched << " cliques " << endl << flush;
   }
 
-
   if (debug){
     for (int i=0;i<n;i++) cout << " ";
     cout << "entering store_clique." << endl;
   }
 
-
   //  // check if the current clique already exists
   //  typename multimap<int, set<T> >::iterator itall = allQ.begin();
   //  size_t currentsize = Q.size();
@@ -443,7 +454,7 @@ void cliques_compatibility_graph<T>::store_clique(int n){
   //    if (itall->first != currentsize) continue;
   //    typename std::vector<T> common(currentsize);
   //    typename std::vector<T>::iterator itfind = std::set_intersection (itall->second.begin(),itall->second.end(), Q.begin(),Q.end(), common.begin());
-  //    common.resize(itfind-common.begin()); 
+  //    common.resize(itfind-common.begin());
   //    if (common.size()==currentsize){
   //      cout << "deux cliques les mêmes !!!" << endl;
   //      cout << "première: " << endl;
@@ -462,8 +473,6 @@ void cliques_compatibility_graph<T>::store_clique(int n){
   //  }
   //  // END
 
-
-
   bool found_best_clique_so_far=false;
   if (Q.size()>max_clique_size){
     max_clique_size=Q.size();
@@ -484,9 +493,6 @@ void cliques_compatibility_graph<T>::store_clique(int n){
     found_the_ultimate_max_clique=true;
   }
 
-
-
-
   const bool verbose=false;
   if (verbose) cout << "MAX CLIQUE found of size " << Q.size() << " # total cliques searched:" << total_nb_of_cliques_searched << endl;
 
@@ -496,7 +502,7 @@ void cliques_compatibility_graph<T>::store_clique(int n){
   }
 
   // storing the current clique... or not, depending of the number of cliques to store, to reduce memory footprint.
-  bool store_it = true; 
+  bool store_it = true;
   bool delete_worst=false;
   if ((max_nb_of_stored_cliques) && (allQ.size()>=max_nb_of_stored_cliques)){
     // then, compare sizes...
@@ -521,7 +527,7 @@ void cliques_compatibility_graph<T>::store_clique(int n){
     }
   }
 
-  // finally, possibly delete worst clique, to reduce memory footprint 
+  // finally, possibly delete worst clique, to reduce memory footprint
   if (delete_worst) allQ.erase(allQ.begin());
 
   if (debug){
@@ -539,7 +545,8 @@ void cliques_compatibility_graph<T>::store_clique(int n){
 
 
 template<class T>
-void cliques_compatibility_graph<T>::find_cliques(graph_data &s, int n){
+void cliques_compatibility_graph<T>::find_cliques(graph_data &s, int n)
+{
   // possible end
   if (s.size()==0){
     store_clique(n);
@@ -650,8 +657,10 @@ void cliques_compatibility_graph<T>::find_cliques(graph_data &s, int n){
 
 
 template<class T>
-void cliques_compatibility_graph<T>::erase_entry(graph_data &s, T &u, hash_key &key){
-  pair<typename graph_data::iterator, typename graph_data::iterator> range = s.equal_range(key);
+void cliques_compatibility_graph<T>::erase_entry(graph_data &s, T &u, hash_key &key)
+{
+  pair<typename graph_data::iterator, typename graph_data::iterator> range =
+    s.equal_range(key);
 
   typename graph_data::iterator it = range.first;
   for (;it!=range.second;it++){
@@ -664,7 +673,8 @@ void cliques_compatibility_graph<T>::erase_entry(graph_data &s, T &u, hash_key &
 
 
 template<class T>
-void cliques_compatibility_graph<T>::choose_u(const graph_data &s, T &u, hash_key &u_key){
+void cliques_compatibility_graph<T>::choose_u(const graph_data &s, T &u, hash_key &u_key)
+{
   if (s.size()==0){
     u=T();
     u_key=0;
@@ -698,19 +708,17 @@ void cliques_compatibility_graph<T>::choose_u(const graph_data &s, T &u, hash_ke
     }
   }
 
-
-
   //  typename map<T, std::vector<double> >::const_iterator itfind = hex_ranks.find(u);
   //  cout << "u: " << u << " # face tri: " << itfind->second[0] << " quality: " << u->get_quality() << " value max: " << valuemax << endl;
 
-
-
   return;
 }
 
-
 template<class T>
-void cliques_compatibility_graph<T>::split_set_BW(const T &u, const hash_key &u_key, const graph_data &s, graph_data &white, graph_data &black){
+void cliques_compatibility_graph<T>::split_set_BW(const T &u, const hash_key &u_key,
+                                                  const graph_data &s, graph_data &white,
+                                                  graph_data &black)
+{
   // splitting set s into white and black nodes
   white.insert(make_pair(u_key,u));
   typename graph_data::const_iterator it = s.begin();
@@ -726,7 +734,9 @@ void cliques_compatibility_graph<T>::split_set_BW(const T &u, const hash_key &u_
 
 
 template<class T>
-void cliques_compatibility_graph<T>::fill_black_set(const T &u, const hash_key &u_key, const graph_data &s, graph_data &black){
+void cliques_compatibility_graph<T>::fill_black_set(const T &u, const hash_key &u_key,
+                                                    const graph_data &s, graph_data &black)
+{
   // filling black set
   typename graph_data::const_iterator it = s.begin();
   for (;it!=s.end();it++){
@@ -739,7 +749,10 @@ void cliques_compatibility_graph<T>::fill_black_set(const T &u, const hash_key &
 
 // the maximum score (int) will be chosen...
 template<class T>
-double cliques_compatibility_graph<T>::function_to_maximize_for_u(const T &u, const hash_key &u_key, const graph_data &s){
+double cliques_compatibility_graph<T>::function_to_maximize_for_u(const T &u,
+                                                                  const hash_key &u_key,
+                                                                  const graph_data &s)
+{
   typename graph_data::const_iterator it = s.begin();
   int counter=0;
   for (;it!=s.end();it++){
@@ -760,8 +773,9 @@ double cliques_compatibility_graph<T>::function_to_maximize_for_u(const T &u, co
 
 
 template<class T>
-bool cliques_compatibility_graph<T>::compatibility(const T &u, const hash_key &u_key, const T &v, const hash_key &v_key){
-
+bool cliques_compatibility_graph<T>::compatibility(const T &u, const hash_key &u_key,
+                                                   const T &v, const hash_key &v_key)
+{
   // first, find u in graph
   pair<typename graph::const_iterator, typename graph::const_iterator> range_ukey = G.equal_range(u_key);
   typename graph::const_iterator itfind_u = range_ukey.first;
@@ -784,8 +798,9 @@ bool cliques_compatibility_graph<T>::compatibility(const T &u, const hash_key &u
 
 
 template<class T>
-bool cliques_losses_graph<T>::compatibility(const T &u, const hash_key &u_key, const T &v, const hash_key &v_key){
-
+bool cliques_losses_graph<T>::compatibility(const T &u, const hash_key &u_key,
+                                            const T &v, const hash_key &v_key)
+{
   // first, find u in graph
   pair<typename graph::const_iterator, typename graph::const_iterator> range_ukey = G.equal_range(u_key);
   typename graph::const_iterator itfind_u = range_ukey.first;
@@ -808,22 +823,24 @@ bool cliques_losses_graph<T>::compatibility(const T &u, const hash_key &u_key, c
 
 
 template<class T>
-cliques_losses_graph<T>::cliques_losses_graph(graph &_g, const map<T, std::vector<double> > &_hex_ranks, unsigned int _max_nb_cliques, unsigned int _nb_hex_potentiels, clique_stop_criteria<T> *csc, ptrfunction_export fct):cliques_compatibility_graph<T>(_g, _hex_ranks, _max_nb_cliques, _nb_hex_potentiels,csc,fct),G(_g){
+cliques_losses_graph<T>::cliques_losses_graph
+(graph &_g, const map<T, std::vector<double> > &_hex_ranks,
+ unsigned int _max_nb_cliques, unsigned int _nb_hex_potentiels,
+ clique_stop_criteria<T> *csc, ptrfunction_export fct)
+  : cliques_compatibility_graph<T>(_g, _hex_ranks, _max_nb_cliques,
+                                   _nb_hex_potentiels,csc,fct),G(_g)
+{
 };
 
-
 template<class T>
 cliques_losses_graph<T>::~cliques_losses_graph(){
 };
 
 
-
-
-
-
-
-bool compare_hex_ptr_by_quality (Hex *a,Hex *b) { return (a->get_quality()>(b->get_quality())); }
-
+bool compare_hex_ptr_by_quality (Hex *a,Hex *b)
+{
+  return (a->get_quality()>(b->get_quality()));
+}
 
 /*****************************************/
 /****************class Hex****************/
@@ -911,33 +928,33 @@ MVertex* Hex::get_h(){
 MVertex* Hex::getVertex(int n){
   MVertex *v;
   switch (n){
-    case 0:
-      v = get_a();
-      break;
-    case 1:
-      v = get_b();
-      break;
-    case 2:
-      v = get_c();
-      break;
-    case 3:
-      v = get_d();
-      break;
-    case 4:
-      v = get_e();
-      break;
-    case 5:
-      v = get_f();
-      break;
-    case 6:
-      v = get_g();
-      break;
-    case 7:
-      v = get_h();
-      break;
-    default:
-      cout << "Hex: unknown vertex number " << n << endl;
-      throw;
+  case 0:
+    v = get_a();
+    break;
+  case 1:
+    v = get_b();
+    break;
+  case 2:
+    v = get_c();
+    break;
+  case 3:
+    v = get_d();
+    break;
+  case 4:
+    v = get_e();
+    break;
+  case 5:
+    v = get_f();
+    break;
+  case 6:
+    v = get_g();
+    break;
+  case 7:
+    v = get_h();
+    break;
+  default:
+    cout << "Hex: unknown vertex number " << n << endl;
+    throw;
   }
   return v;
 }
@@ -1297,12 +1314,12 @@ void Recombinator::execute(){
   GModel::riter it;
 
   for(it=model->firstRegion();it!=model->lastRegion();it++)
-  {
-    gr = *it;
-    if(gr->getNumMeshElements()>0){
-      execute(gr);
+    {
+      gr = *it;
+      if(gr->getNumMeshElements()>0){
+        execute(gr);
+      }
     }
-  }
 }
 
 void Recombinator::execute(GRegion* gr){
@@ -1880,21 +1897,21 @@ void Recombinator::build_tuples(GRegion* gr){
   faces = gr->faces();
 
   for(it=faces.begin();it!=faces.end();it++)
-  {
-    gf = *it;
+    {
+      gf = *it;
 
-    for(i=0;i<gf->getNumMeshElements();i++){
-      element = gf->getMeshElement(i);
+      for(i=0;i<gf->getNumMeshElements();i++){
+        element = gf->getMeshElement(i);
 
-      if(element->getNumVertices()==3){
-        a = element->getVertex(0);
-        b = element->getVertex(1);
-        c = element->getVertex(2);
+        if(element->getNumVertices()==3){
+          a = element->getVertex(0);
+          b = element->getVertex(1);
+          c = element->getVertex(2);
 
-        tuples.insert(Tuple(a,b,c,element,gf));
+          tuples.insert(Tuple(a,b,c,element,gf));
+        }
       }
     }
-  }
 }
 
 void Recombinator::modify_surfaces(GRegion* gr){
@@ -1933,28 +1950,28 @@ void Recombinator::modify_surfaces(GRegion* gr){
   faces = gr->faces();
 
   for(it=faces.begin();it!=faces.end();it++)
-  {
-    gf = *it;
+    {
+      gf = *it;
 
-    opt.clear();
+      opt.clear();
 
-    for(i=0;i<gf->getNumMeshElements();i++){
-      element = gf->getMeshElement(i);
+      for(i=0;i<gf->getNumMeshElements();i++){
+        element = gf->getMeshElement(i);
 
-      if(element->getNumVertices()==3){
-        it2 = triangles.find(element);
-        if(it2==triangles.end()){
-          opt.push_back(element);
+        if(element->getNumVertices()==3){
+          it2 = triangles.find(element);
+          if(it2==triangles.end()){
+            opt.push_back(element);
+          }
         }
       }
-    }
 
-    gf->triangles.clear();
+      gf->triangles.clear();
 
-    for(i=0;i<opt.size();i++){
-      gf->triangles.push_back((MTriangle*)opt[i]);
+      for(i=0;i<opt.size();i++){
+        gf->triangles.push_back((MTriangle*)opt[i]);
+      }
     }
-  }
 }
 
 void Recombinator::modify_surfaces(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
@@ -2448,7 +2465,7 @@ MVertex* Recombinator::find(MVertex* v1,MVertex* v2,MVertex* v3,MVertex* already
 }
 
 void Recombinator::intersection(const std::set<MVertex*>& bin1,const std::set<MVertex*>& bin2,
-    const std::vector<MVertex*>& already,std::set<MVertex*>& final){
+                                const std::vector<MVertex*>& already,std::set<MVertex*>& final){
   size_t i;
   bool ok;
   std::set<MVertex*> temp;
@@ -2473,7 +2490,7 @@ void Recombinator::intersection(const std::set<MVertex*>& bin1,const std::set<MV
 }
 
 void Recombinator::intersection(const std::set<MVertex*>& bin1,const std::set<MVertex*>& bin2,const std::set<MVertex*>& bin3,
-    const std::vector<MVertex*>& already,std::set<MVertex*>& final){
+                                const std::vector<MVertex*>& already,std::set<MVertex*>& final){
   size_t i;
   bool ok;
   std::set<MVertex*> temp;
@@ -2644,8 +2661,8 @@ bool Recombinator::duplicate(Diagonal diagonal){
 // est "conforme A" une face si ses 4 facets existent dans tableA, ou bien si aucune des ses facets ne se trouve dans table A
 // ça veut dire: un hex est condorme A s'il est tout seul.
 // Sinon, s'il est pas tout seul, s'il a des "voisins de face"... il faut que les faces en contact soient bien en contact... pas juste un triangle en commun, mais il faut un QUAD en commun !
-// ne pouvait-on pas écrire ça avec des quad, du coup ??? 
-// ça veut dire que pour les slivers, qu'ils soient ou pas dans l'hex, dans tous les cas, les hex sont compatibles ! 
+// ne pouvait-on pas écrire ça avec des quad, du coup ???
+// ça veut dire que pour les slivers, qu'ils soient ou pas dans l'hex, dans tous les cas, les hex sont compatibles !
 bool Recombinator::conformityA(Hex &hex){
   bool c1,c2,c3,c4,c5,c6;
   MVertex *a,*b,*c,*d;
@@ -2691,8 +2708,8 @@ bool Recombinator::conformityA(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
 }
 
 // return false si:
-//- une des 12 arrêtes de l'hex se trouve dans tableB !!! (pas C !!!), càd si une arrete a été utilisée comme diagonale d'un autre hex 
-//- (ou bien) si, pour chaque face de l'hex, on a une diagonale dans tableB et pas l'autre 
+//- une des 12 arrêtes de l'hex se trouve dans tableB !!! (pas C !!!), càd si une arrete a été utilisée comme diagonale d'un autre hex
+//- (ou bien) si, pour chaque face de l'hex, on a une diagonale dans tableB et pas l'autre
 // ce test conformityB n'est-il pas redondant avec conformityA ???
 bool Recombinator::conformityB(Hex &hex){
   bool flag1;
@@ -2932,7 +2949,7 @@ void Recombinator::build_vertex_to_vertices(GRegion* gr){
   for(i=0;i<gr->getNumMeshElements();i++){
     if(i%progress==0){
       done += percentage*100;
-      cout << "..." << done << "% " << flush; 
+      cout << "..." << done << "% " << flush;
     }
     element = gr->getMeshElement(i);
     for(j=0;j<element->getNumVertices();j++){
@@ -2967,7 +2984,7 @@ void Recombinator::build_vertex_to_elements(GRegion* gr){
   MVertex* vertex;
   std::set<MElement*> bin;
   std::map<MVertex*,std::set<MElement*> >::iterator it;
-  
+
   int nbElements = gr->getNumMeshElements();
   double percentage=0.05;
   double done=0.;
@@ -2978,7 +2995,7 @@ void Recombinator::build_vertex_to_elements(GRegion* gr){
   for(i=0;i<gr->getNumMeshElements();i++){
     if(i%progress==0){
       done += percentage*100;
-      cout << "..." << done << "% " << flush; 
+      cout << "..." << done << "% " << flush;
     }
     element = gr->getMeshElement(i);
     for(j=0;j<element->getNumVertices();j++){
@@ -2998,7 +3015,7 @@ void Recombinator::build_vertex_to_elements(GRegion* gr){
   cout << endl;
 }
 
-// pour les 6 faces de l'hex, stocke les 4 "facet" (équivalent à PETriangle) possibles dans la table A 
+// pour les 6 faces de l'hex, stocke les 4 "facet" (équivalent à PETriangle) possibles dans la table A
 void Recombinator::build_hash_tableA(Hex hex){
   MVertex *a,*b,*c,*d;
   MVertex *e,*f,*g,*h;
@@ -3054,7 +3071,7 @@ void Recombinator::build_hash_tableA(Facet facet){
   }
 }
 
-// pour les 6 faces de l'hex, stoke les 2 "diagonal" possibles dans la table B 
+// pour les 6 faces de l'hex, stoke les 2 "diagonal" possibles dans la table B
 void Recombinator::build_hash_tableB(Hex hex){
   MVertex *a,*b,*c,*d;
   MVertex *e,*f,*g,*h;
@@ -3215,9 +3232,9 @@ void Recombinator::print_hash_tableA(){
 
 void Recombinator::print_segment(SPoint3 p1,SPoint3 p2,std::ofstream& file){
   file << "SL ("
-    << p1.x() << ", " << p1.y() << ", " << p1.z() << ", "
-    << p2.x() << ", " << p2.y() << ", " << p2.z() << ")"
-    << "{10, 20};\n";
+       << p1.x() << ", " << p1.y() << ", " << p1.z() << ", "
+       << p2.x() << ", " << p2.y() << ", " << p2.z() << ")"
+       << "{10, 20};\n";
 }
 
 double Recombinator::scaled_jacobian(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
@@ -3393,12 +3410,12 @@ void Supplementary::execute(){
   GModel::riter it;
 
   for(it=model->firstRegion();it!=model->lastRegion();it++)
-  {
-    gr = *it;
-    if(gr->getNumMeshElements()>0){
-      execute(gr);
+    {
+      gr = *it;
+      if(gr->getNumMeshElements()>0){
+        execute(gr);
+      }
     }
-  }
 }
 
 void Supplementary::execute(GRegion* gr){
@@ -3701,21 +3718,21 @@ void Supplementary::build_tuples(GRegion* gr){
   faces = gr->faces();
 
   for(it=faces.begin();it!=faces.end();it++)
-  {
-    gf = *it;
+    {
+      gf = *it;
 
-    for(i=0;i<gf->getNumMeshElements();i++){
-      element = gf->getMeshElement(i);
+      for(i=0;i<gf->getNumMeshElements();i++){
+        element = gf->getMeshElement(i);
 
-      if(element->getNumVertices()==3){
-        a = element->getVertex(0);
-        b = element->getVertex(1);
-        c = element->getVertex(2);
+        if(element->getNumVertices()==3){
+          a = element->getVertex(0);
+          b = element->getVertex(1);
+          c = element->getVertex(2);
 
-        tuples.insert(Tuple(a,b,c,element,gf));
+          tuples.insert(Tuple(a,b,c,element,gf));
+        }
       }
     }
-  }
 }
 
 void Supplementary::modify_surfaces(GRegion* gr){
@@ -3749,28 +3766,28 @@ void Supplementary::modify_surfaces(GRegion* gr){
   faces = gr->faces();
 
   for(it=faces.begin();it!=faces.end();it++)
-  {
-    gf = *it;
+    {
+      gf = *it;
 
-    opt.clear();
+      opt.clear();
 
-    for(i=0;i<gf->getNumMeshElements();i++){
-      element = gf->getMeshElement(i);
+      for(i=0;i<gf->getNumMeshElements();i++){
+        element = gf->getMeshElement(i);
 
-      if(element->getNumVertices()==3){
-        it2 = triangles.find(element);
-        if(it2==triangles.end()){
-          opt.push_back(element);
+        if(element->getNumVertices()==3){
+          it2 = triangles.find(element);
+          if(it2==triangles.end()){
+            opt.push_back(element);
+          }
         }
       }
-    }
 
-    gf->triangles.clear();
+      gf->triangles.clear();
 
-    for(i=0;i<opt.size();i++){
-      gf->triangles.push_back((MTriangle*)opt[i]);
+      for(i=0;i<opt.size();i++){
+        gf->triangles.push_back((MTriangle*)opt[i]);
+      }
     }
-  }
 }
 
 void Supplementary::modify_surfaces(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
@@ -4068,7 +4085,7 @@ void Supplementary::find(MVertex* vertex,Prism prism,std::set<MElement*>& final)
 }
 
 void Supplementary::intersection(const std::set<MVertex*>& bin1,const std::set<MVertex*>& bin2,
-    const std::vector<MVertex*>& already,std::set<MVertex*>& final){
+                                 const std::vector<MVertex*>& already,std::set<MVertex*>& final){
   size_t i;
   bool ok;
   std::set<MVertex*> temp;
@@ -4714,12 +4731,12 @@ void PostOp::execute(bool flag){
   GModel::riter it;
 
   for(it=model->firstRegion();it!=model->lastRegion();it++)
-  {
-    gr = *it;
-    if(gr->getNumMeshElements()>0){
-      execute(gr,flag);
+    {
+      gr = *it;
+      if(gr->getNumMeshElements()>0){
+        execute(gr,flag);
+      }
     }
-  }
 }
 
 void PostOp::execute(GRegion* gr,bool flag){
@@ -5232,21 +5249,21 @@ void PostOp::build_tuples(GRegion* gr){
   faces = gr->faces();
 
   for(it=faces.begin();it!=faces.end();it++)
-  {
-    gf = *it;
+    {
+      gf = *it;
 
-    for(i=0;i<gf->getNumMeshElements();i++){
-      element = gf->getMeshElement(i);
+      for(i=0;i<gf->getNumMeshElements();i++){
+        element = gf->getMeshElement(i);
 
-      if(element->getNumVertices()==3){
-        a = element->getVertex(0);
-        b = element->getVertex(1);
-        c = element->getVertex(2);
+        if(element->getNumVertices()==3){
+          a = element->getVertex(0);
+          b = element->getVertex(1);
+          c = element->getVertex(2);
 
-        tuples.insert(Tuple(a,b,c,element,gf));
+          tuples.insert(Tuple(a,b,c,element,gf));
+        }
       }
     }
-  }
 }
 
 void PostOp::modify_surfaces(GRegion* gr){
@@ -5276,28 +5293,28 @@ void PostOp::modify_surfaces(GRegion* gr){
   faces = gr->faces();
 
   for(it=faces.begin();it!=faces.end();it++)
-  {
-    gf = *it;
+    {
+      gf = *it;
 
-    opt.clear();
+      opt.clear();
 
-    for(i=0;i<gf->getNumMeshElements();i++){
-      element = gf->getMeshElement(i);
+      for(i=0;i<gf->getNumMeshElements();i++){
+        element = gf->getMeshElement(i);
 
-      if(element->getNumVertices()==3){
-        it2 = triangles.find(element);
-        if(it2==triangles.end()){
-          opt.push_back(element);
+        if(element->getNumVertices()==3){
+          it2 = triangles.find(element);
+          if(it2==triangles.end()){
+            opt.push_back(element);
+          }
         }
       }
-    }
 
-    gf->triangles.clear();
+      gf->triangles.clear();
 
-    for(i=0;i<opt.size();i++){
-      gf->triangles.push_back((MTriangle*)opt[i]);
+      for(i=0;i<opt.size();i++){
+        gf->triangles.push_back((MTriangle*)opt[i]);
+      }
     }
-  }
 }
 
 void PostOp::modify_surfaces(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
@@ -6049,10 +6066,10 @@ void Recombinator_Graph::fill_tet_to_hex_table(Hex *hex){
 
 
 
-  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  
-  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  
-  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  
-  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  
+  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!
+  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!
+  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!
+  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!
   // le test suivant déconne, cf. cube à 125 hex... sais pas pourquoi, mais ça élimine des hex potentiels tout à fait valides (et nécessaires à la solution optimale)!!!
   // dans le cas cube125: ce test fait passer le #hex potentiels recensés par les patterns de recombinaison d'environ 1000 à 2434 ??? !!!, ça double aussi le nbre d'hex dans le graphe des pertes
   // en fait, si les noeuds par exemple a b c d deviennent b c d a, on a une rotation et un hex tout pourri... mais qui a même hash et mêmes noeuds ?!?!
@@ -6065,10 +6082,10 @@ void Recombinator_Graph::fill_tet_to_hex_table(Hex *hex){
       return;
     }
   }
-  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  
-  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  
-  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  
-  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  
+  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!
+  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!
+  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!
+  // TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!  TODO !!!!
 
 
 
@@ -6126,12 +6143,12 @@ void Recombinator_Graph::buildGraphOnly(unsigned int max_nb_cliques,string filen
   GModel::riter it;
 
   for(it=model->firstRegion();it!=model->lastRegion();it++)
-  {
-    gr = *it;
-    if(gr->getNumMeshElements()>0){
-      buildGraphOnly(gr, max_nb_cliques,filename);
+    {
+      gr = *it;
+      if(gr->getNumMeshElements()>0){
+        buildGraphOnly(gr, max_nb_cliques,filename);
+      }
     }
-  }
 }
 
 
@@ -6172,12 +6189,12 @@ void Recombinator_Graph::execute(){
   model->writeMSH("beforeyamakawa.msh");
 
   for(it=model->firstRegion();it!=model->lastRegion();it++)
-  {
-    gr = *it;
-    if(gr->getNumMeshElements()>0){
-      execute(gr);
+    {
+      gr = *it;
+      if(gr->getNumMeshElements()>0){
+        execute(gr);
+      }
     }
-  }
 }
 
 
@@ -6189,12 +6206,12 @@ void Recombinator_Graph::execute_blossom(unsigned int max_nb_cliques,string file
   model->writeMSH("beforeyamakawa.msh");
 
   for(it=model->firstRegion();it!=model->lastRegion();it++)
-  {
-    gr = *it;
-    if(gr->getNumMeshElements()>0){
-      execute_blossom(gr, max_nb_cliques,filename);
+    {
+      gr = *it;
+      if(gr->getNumMeshElements()>0){
+        execute_blossom(gr, max_nb_cliques,filename);
+      }
     }
-  }
 }
 
 
@@ -6210,7 +6227,7 @@ Recombinator_Graph::~Recombinator_Graph(){
 
 void Recombinator_Graph::createBlossomInfo(){
 
-  throw; 
+  throw;
 
 
   GRegion* gr;
@@ -6218,12 +6235,12 @@ void Recombinator_Graph::createBlossomInfo(){
   GModel::riter it;
 
   for(it=model->firstRegion();it!=model->lastRegion();it++)
-  {
-    gr = *it;
-    //    if(gr->getNumMeshElements()>0){
-    createBlossomInfo(gr);
-    //    }
-  }
+    {
+      gr = *it;
+      //    if(gr->getNumMeshElements()>0){
+      createBlossomInfo(gr);
+      //    }
+    }
 }
 
 
@@ -6350,7 +6367,7 @@ void Recombinator_Graph::execute_blossom(GRegion* gr, unsigned int max_nb_clique
 
   Msg::Info("Building Connectivity...");
 
-  clock_t a=clock();
+  double a = Cpu();
 
   build_vertex_to_vertices(gr);
   build_vertex_to_elements(gr);
@@ -6366,12 +6383,12 @@ void Recombinator_Graph::execute_blossom(GRegion* gr, unsigned int max_nb_clique
   create_losses_graph(gr);
   // add points to potential hexas containing original blossom pairs of triangles
   compute_hex_ranks_blossom();
-  
-  double time_building_graph = (clock() - a) / (double)CLOCKS_PER_SEC;
 
+  double time_building_graph = (Cpu() - a);
 
-  a=clock();
-  // a criteria to stop when the whole domain is exclusively composed of hex 
+
+  a=Cpu();
+  // a criteria to stop when the whole domain is exclusively composed of hex
   clique_stop_criteria<Hex*> criteria(hex_to_tet, gr->tetrahedra.size());
 
   cliques_losses_graph<Hex*> cl(incompatibility_graph, hex_ranks, max_nb_cliques, hex_to_tet.size(),&criteria,export_the_clique_graphviz_format);
@@ -6379,8 +6396,8 @@ void Recombinator_Graph::execute_blossom(GRegion* gr, unsigned int max_nb_clique
   //cl.export_cliques();
 
 
-  double time_cliques = (clock() - a) / (double)CLOCKS_PER_SEC;
-    
+  double time_cliques = (Cpu() - a) ;
+
   cout << "RECOMBINATOR_GRAPH timing:" << endl;
   cout << "  ------- TIME BUILDING GRAPH : " << time_building_graph << " s." << endl;
   cout << "  ------- TIME CLIQUE         : " << time_cliques << " s." << endl;
@@ -6481,7 +6498,7 @@ void Recombinator_Graph::execute(GRegion* gr){
   compute_hex_ranks();
 
 
-  // a criteria to stop when the whole domain is exclusively composed of hex 
+  // a criteria to stop when the whole domain is exclusively composed of hex
   found_the_ultimate_max_clique=false;
   clique_stop_criteria<Hex*> criteria(hex_to_tet, gr->tetrahedra.size());
 
@@ -6892,7 +6909,7 @@ void Recombinator_Graph::export_all_hex(int &file, GRegion *gr){
 }
 
 
-// check if the hex is good enough to be put into the graph. If not in the graph, it cannot be chosen... 
+// check if the hex is good enough to be put into the graph. If not in the graph, it cannot be chosen...
 bool Recombinator_Graph::is_not_good_enough(Hex* hex){
 
   //  // Hack
@@ -6909,701 +6926,703 @@ bool Recombinator_Graph::is_not_good_enough(Hex* hex){
     return true;
   }
   return false;
-  }
+}
 
 
-  // fills incompatibility_graph if two hex share a common (non-sliver!) tet 
-  void Recombinator_Graph::create_indirect_neighbors_graph(){
-    std::map<MElement*, std::set<Hex*> >::iterator it_tet = tet_to_hex.begin();
+// fills incompatibility_graph if two hex share a common (non-sliver!) tet
+void Recombinator_Graph::create_indirect_neighbors_graph()
+{
+  std::map<MElement*, std::set<Hex*> >::iterator it_tet = tet_to_hex.begin();
 
-    for (;it_tet!=tet_to_hex.end();it_tet++){
-      std::set<Hex*>::iterator it_hex1 = it_tet->second.begin();
-      for (;it_hex1!=it_tet->second.end();it_hex1++){
-        Hex *hex = *it_hex1;
+  for (;it_tet!=tet_to_hex.end();it_tet++){
+    std::set<Hex*>::iterator it_hex1 = it_tet->second.begin();
+    for (;it_hex1!=it_tet->second.end();it_hex1++){
+      Hex *hex = *it_hex1;
 
-        //      cout << " create_indirect_neighbors_graph:: treating hex " << hex << " made of ";
-        //      for (int i=0;i<8;i++)
-        //        cout << "  " << hex->getVertex(i)->getNum();
-        //      cout << endl;
+      //      cout << " create_indirect_neighbors_graph:: treating hex " << hex << " made of ";
+      //      for (int i=0;i<8;i++)
+      //        cout << "  " << hex->getVertex(i)->getNum();
+      //      cout << endl;
 
-        if (is_not_good_enough(hex)){
-          continue;
-        }
-        graph::iterator itfind_graph = find_hex_in_graph(hex);
-        if (itfind_graph==incompatibility_graph.end()){
-          incompatibility_graph.insert(make_pair(hex->get_hash(),make_pair(hex,graph_data())));// this creates an entry if none exists. Important ! if hex is good, but surrounded by "not good enough" hex, it has to be in the graph anyway ! A single node in the graph, without any connection... the perfect lonely hex.
-          set_of_all_hex_in_graph.insert(hex);
-        }
-        // check if the hex quality is sufficient to be put into the graph...
-        //cout << "creating graph for hex " << hex << " with quality " << hex->get_quality() << endl;
-        //      if (hex->get_quality()>=0.1)
-        // check if the tet is a sliver or not...
-        // if sliver, the tet should not be taken into account in the graph creation ! 
-        bool is_sliver = sliver(it_tet->first,*hex);
-        if (!is_sliver){
-          //cout << " is not sliver, size it_tet->second:" << it_tet->second.size() << endl;
-          std::set<Hex*>::iterator it_hex2 = it_tet->second.begin();
-          for (;it_hex2!=it_tet->second.end();it_hex2++){
-            if ((hex)!=(*it_hex2)){
-              if (is_not_good_enough(*it_hex2)){
-                continue;
-              }
-              //cout << "   hex put in graph, linked to hex " << *it_hex2 << endl;
-              add_graph_entry(hex,*it_hex2);
-              //cout << "hex's " << hex << " and " << *it_hex2 << " have one tet in common: " << it_tet->first << endl;
-              //            incompatibility_graph[*it_hex2].insert(hex);
+      if (is_not_good_enough(hex)){
+        continue;
+      }
+      graph::iterator itfind_graph = find_hex_in_graph(hex);
+      if (itfind_graph==incompatibility_graph.end()){
+        incompatibility_graph.insert(make_pair(hex->get_hash(),make_pair(hex,graph_data())));// this creates an entry if none exists. Important ! if hex is good, but surrounded by "not good enough" hex, it has to be in the graph anyway ! A single node in the graph, without any connection... the perfect lonely hex.
+        set_of_all_hex_in_graph.insert(hex);
+      }
+      // check if the hex quality is sufficient to be put into the graph...
+      //cout << "creating graph for hex " << hex << " with quality " << hex->get_quality() << endl;
+      //      if (hex->get_quality()>=0.1)
+      // check if the tet is a sliver or not...
+      // if sliver, the tet should not be taken into account in the graph creation !
+      bool is_sliver = sliver(it_tet->first,*hex);
+      if (!is_sliver){
+        //cout << " is not sliver, size it_tet->second:" << it_tet->second.size() << endl;
+        std::set<Hex*>::iterator it_hex2 = it_tet->second.begin();
+        for (;it_hex2!=it_tet->second.end();it_hex2++){
+          if ((hex)!=(*it_hex2)){
+            if (is_not_good_enough(*it_hex2)){
+              continue;
             }
+            //cout << "   hex put in graph, linked to hex " << *it_hex2 << endl;
+            add_graph_entry(hex,*it_hex2);
+            //cout << "hex's " << hex << " and " << *it_hex2 << " have one tet in common: " << it_tet->first << endl;
+            //            incompatibility_graph[*it_hex2].insert(hex);
           }
         }
-        //      else{
-        //        cout << " hex's " << hex << " has a sliver !!!!!!!!!!!!!!! : " << it_tet->first << endl;
-        //      }
       }
+      //      else{
+      //        cout << " hex's " << hex << " has a sliver !!!!!!!!!!!!!!! : " << it_tet->first << endl;
+      //      }
     }
   }
+}
 
 
-  std::multimap<unsigned long long, Hex* >::const_iterator  Recombinator_Graph::find_the_created_potential_hex(Hex *hex, const std::multimap<unsigned long long, Hex*> &list){
-    std::pair<std::multimap<unsigned long long, Hex* >::const_iterator, std::multimap<unsigned long long, Hex* >::const_iterator> range = list.equal_range(hex->get_hash());
-    for (std::multimap<unsigned long long, Hex*>::const_iterator it=range.first;it!=range.second;it++){
-      Hex *candidate = it->second;
-      if (candidate->same_vertices(hex)){
-        return it;
-      }
+std::multimap<unsigned long long, Hex* >::const_iterator  Recombinator_Graph::find_the_created_potential_hex(Hex *hex, const std::multimap<unsigned long long, Hex*> &list){
+  std::pair<std::multimap<unsigned long long, Hex* >::const_iterator, std::multimap<unsigned long long, Hex* >::const_iterator> range = list.equal_range(hex->get_hash());
+  for (std::multimap<unsigned long long, Hex*>::const_iterator it=range.first;it!=range.second;it++){
+    Hex *candidate = it->second;
+    if (candidate->same_vertices(hex)){
+      return it;
     }
-    return list.end();
   }
+  return list.end();
+}
 
 
-  std::multimap<unsigned long long, pair<PETriangle*,int> >::iterator  Recombinator_Graph::find_the_triangle(PETriangle *t, std::multimap<unsigned long long, pair<PETriangle*, int> > &list){
-    std::pair<std::multimap<unsigned long long, pair<PETriangle*,int> >::iterator, std::multimap<unsigned long long, pair<PETriangle*,int> >::iterator> range = list.equal_range(t->get_hash());
-    for (std::multimap<unsigned long long, pair<PETriangle*,int> >::iterator it=range.first;it!=range.second;it++){
-      PETriangle *candidate = it->second.first;
-      if (candidate->same_vertices(t)){
-        it->second.second++;
-        return it;
-      }
+std::multimap<unsigned long long, pair<PETriangle*,int> >::iterator  Recombinator_Graph::find_the_triangle(PETriangle *t, std::multimap<unsigned long long, pair<PETriangle*, int> > &list){
+  std::pair<std::multimap<unsigned long long, pair<PETriangle*,int> >::iterator, std::multimap<unsigned long long, pair<PETriangle*,int> >::iterator> range = list.equal_range(t->get_hash());
+  for (std::multimap<unsigned long long, pair<PETriangle*,int> >::iterator it=range.first;it!=range.second;it++){
+    PETriangle *candidate = it->second.first;
+    if (candidate->same_vertices(t)){
+      it->second.second++;
+      return it;
     }
-    return list.end();
   }
+  return list.end();
+}
 
 
 
-  Recombinator_Graph::citer  Recombinator_Graph::find_the_triangle(PETriangle *t, const trimap &list){
-    std::pair<citer, citer> range = list.equal_range(t->get_hash());
-    for (citer it=range.first;it!=range.second;it++){
-      if (it->second->same_vertices(t)) return it;
-    }
-    return list.end();
+Recombinator_Graph::citer  Recombinator_Graph::find_the_triangle(PETriangle *t, const trimap &list){
+  std::pair<citer, citer> range = list.equal_range(t->get_hash());
+  for (citer it=range.first;it!=range.second;it++){
+    if (it->second->same_vertices(t)) return it;
   }
+  return list.end();
+}
 
 
-  Recombinator_Graph::linemap::const_iterator  Recombinator_Graph::find_the_line(PELine *t, const linemap &list){
-    std::pair<linemap::const_iterator, linemap::const_iterator> range = list.equal_range(t->get_hash());
-    for (linemap::const_iterator it=range.first;it!=range.second;it++){
-      if (it->second->same_vertices(t)) return it;
-    }
-    return list.end();
+Recombinator_Graph::linemap::const_iterator  Recombinator_Graph::find_the_line(PELine *t, const linemap &list){
+  std::pair<linemap::const_iterator, linemap::const_iterator> range = list.equal_range(t->get_hash());
+  for (linemap::const_iterator it=range.first;it!=range.second;it++){
+    if (it->second->same_vertices(t)) return it;
   }
+  return list.end();
+}
 
 
-  void Recombinator_Graph::export_direct_neighbor_table(int max){
-    stringstream ss;
-    ss << "neighbors_table";
-    ofstream out(ss.str().c_str());
+void Recombinator_Graph::export_direct_neighbor_table(int max){
+  stringstream ss;
+  ss << "neighbors_table";
+  ofstream out(ss.str().c_str());
 
-    std::multimap<int,Hex*>::iterator it =  ndegree.begin();
+  std::multimap<int,Hex*>::iterator it =  ndegree.begin();
 
-    int counter=0;
-    out << " n  neighbors_rank hex* quality" << endl;
-    for (;it!=ndegree.end();it++,counter++){
-      if (counter>=max) break;
-      Hex *hex = it->second;
-      out << counter << "  " << it->first << "  " << hex << "  " << hex->get_quality() << endl;
-      stringstream ss2;
-      ss2 << "neighbors_table_hex";
-      char chose[256];
-      sprintf(chose, "_%0*d", 2, counter);
-      ss2 << chose;
-      ss2 << ".pos";
-      ofstream out2(ss2.str().c_str());
-      out2 << "View \"hex\" {" << endl;
-      out2 << "SH(";
-      for (int n=0;n<8;n++){
-        MVertex *v = hex->getVertex(n);
-        out2 << v->x() << "," << v->y() << "," << v->z();
-        if (n!=7) out2 << ",";
-      }
-      out2 << "){";
-      for (int n=0;n<8;n++){
-        out2 << it->first;
-        if (n!=7) out2 << ",";
-      }
-      out2 << "};" << endl;
-      out2 << "};" << endl;
-      out2.close();
+  int counter=0;
+  out << " n  neighbors_rank hex* quality" << endl;
+  for (;it!=ndegree.end();it++,counter++){
+    if (counter>=max) break;
+    Hex *hex = it->second;
+    out << counter << "  " << it->first << "  " << hex << "  " << hex->get_quality() << endl;
+    stringstream ss2;
+    ss2 << "neighbors_table_hex";
+    char chose[256];
+    sprintf(chose, "_%0*d", 2, counter);
+    ss2 << chose;
+    ss2 << ".pos";
+    ofstream out2(ss2.str().c_str());
+    out2 << "View \"hex\" {" << endl;
+    out2 << "SH(";
+    for (int n=0;n<8;n++){
+      MVertex *v = hex->getVertex(n);
+      out2 << v->x() << "," << v->y() << "," << v->z();
+      if (n!=7) out2 << ",";
     }
-
-
-
-
-    out.close();
-  }
-
-
-  // if two hex are not connected in the incompatibility_graph, they are compatible
-  void Recombinator_Graph::create_losses_graph(GRegion *gr){
-    incompatibility_graph.clear();
-    create_indirect_neighbors_graph();
-    create_direct_neighbors_incompatibility_graph();
-
-
-    // stats
-    graph::iterator it =  incompatibility_graph.begin();
-    graph::iterator ite =  incompatibility_graph.end();
-    int total=0;
-    for (;it!=ite;it++){
-      total+=it->second.second.size();
+    out2 << "){";
+    for (int n=0;n<8;n++){
+      out2 << it->first;
+      if (n!=7) out2 << ",";
     }
-    nbhex_in_losses_graph = incompatibility_graph.size();
-    cout << "total=" << total << endl;
-    cout << "nbhex_in_losses_graph=" << nbhex_in_losses_graph << endl;
-    average_connectivity = total/((double)(nbhex_in_losses_graph));
-    cout << "#hex potentiels recensés dans le graphe des pertes: " << nbhex_in_losses_graph << "   #connectivité moyenne: " << average_connectivity<< endl;
-    cout << "#hex potentiels recensés par les patterns de recombinaison: " << hex_to_tet.size() << endl;
-  };
-
+    out2 << "};" << endl;
+    out2 << "};" << endl;
+    out2.close();
+  }
 
-  // TODO: check only the direct neighbors !!! change the algo pour ne as parcourir les hex qui ont un tet commun. cad Filter au niveau des faces... juste les faces extérieures de l'hex !!! et attention aux slivers, du coup... c'est sans doute la source du pb...
 
-  // pour chaque hex, sortir un .pos avec l'hex, ses tets et ses faces extérieures pour checker... 
-  // pourquoi pas sortir aussi tous les slivers à part... ? 
 
 
+  out.close();
+}
 
 
-  // fills incompatibility_graph if two hex are incompatible direct neighbors, 
-  // i.e. (they DO NOT have one tet in common) and (they are neighbors (face(s) in common) but DO NOT pass the compatibility tests)
-  void Recombinator_Graph::create_direct_neighbors_incompatibility_graph(){
+// if two hex are not connected in the incompatibility_graph, they are compatible
+void Recombinator_Graph::create_losses_graph(GRegion *gr)
+{
+  incompatibility_graph.clear();
+  create_indirect_neighbors_graph();
+  create_direct_neighbors_incompatibility_graph();
 
-    vector<Hex*> visited_hex;
-    std::map<Hex*, std::set<MElement*> >::iterator it_hex =  hex_to_tet.begin();
-    for (;it_hex!=hex_to_tet.end();it_hex++){// for all hex
-      Hex *hex = it_hex->first;
-      if (is_not_good_enough(hex))
-        continue;
-      graph::iterator itfind_graph = find_hex_in_graph(hex);
-      if (itfind_graph==incompatibility_graph.end()){
-        incompatibility_graph.insert(make_pair(hex->get_hash(),make_pair(hex,graph_data())));// this creates an entry if none exists. Important ! if hex is good, but surrounded by "not good enough" hex, it has to be in the graph anyway ! A single node in the graph, without any connection... the perfect lonely hex.
-        set_of_all_hex_in_graph.insert(hex);
-      }
 
-      hash_tableA.clear();
-      hash_tableB.clear();
-      hash_tableC.clear();
-      build_hash_tableA(*hex);
-      build_hash_tableB(*hex);
-      build_hash_tableC(*hex);
-
-      visited_hex.clear();
-
-      std::set<PETriangle*>::iterator it_faces = hex_to_faces[hex].begin();
-      std::set<PETriangle*>::iterator it_facesen = hex_to_faces[hex].end();
-      for (;it_faces!=it_facesen;it_faces++){// i.e. for all triangular external face
-        PETriangle *face = *it_faces;
-        std::set<Hex*>::iterator it_neighbors = faces_to_hex[face].begin();
-        std::set<Hex*>::iterator it_neighborsen = faces_to_hex[face].end();
-        for (;it_neighbors!=it_neighborsen;it_neighbors++){// for all its neighbors
-          if ((*it_neighbors)==hex) continue;
-          Hex *other_hex = *it_neighbors;
-
-          vector<Hex*>::iterator itfind_hex = std::find(visited_hex.begin(),visited_hex.end(),other_hex);
-          if (itfind_hex!=visited_hex.end()) continue;// already done
-          else visited_hex.push_back(other_hex);
-
-          evaluate_hex_couple(hex,other_hex);
-        }
-      }
-      // change following...
-      std::set<PELine*>::iterator it_line = hex_to_edges[hex].begin();
-      std::set<PELine*>::iterator it_lineen = hex_to_edges[hex].end();
-      for (;it_line!=it_lineen;it_line++){// i.e. for all edges and diagonals -> check the hex neighbors too
-        PELine *line = *it_line;
-        std::set<Hex*>::iterator it_neighbors = edges_to_hex[line].begin();
-        std::set<Hex*>::iterator it_neighborsen = edges_to_hex[line].end();
-        for (;it_neighbors!=it_neighborsen;it_neighbors++){// for all its neighbors
-          if ((*it_neighbors)==hex) continue;
-          Hex *other_hex = *it_neighbors;
-
-          vector<Hex*>::iterator itfind_hex = std::find(visited_hex.begin(),visited_hex.end(),other_hex);
-          if (itfind_hex!=visited_hex.end()) continue;// already done
-          else visited_hex.push_back(other_hex);
-
-          evaluate_hex_couple(hex,other_hex);
-        }
-      }
-    }
+  // stats
+  graph::iterator it =  incompatibility_graph.begin();
+  graph::iterator ite =  incompatibility_graph.end();
+  int total=0;
+  for (;it!=ite;it++){
+    total+=it->second.second.size();
+  }
+  nbhex_in_losses_graph = incompatibility_graph.size();
+  cout << "total=" << total << endl;
+  cout << "nbhex_in_losses_graph=" << nbhex_in_losses_graph << endl;
+  average_connectivity = total/((double)(nbhex_in_losses_graph));
+  cout << "#hex potentiels recensés dans le graphe des pertes: " << nbhex_in_losses_graph << "   #connectivité moyenne: " << average_connectivity<< endl;
+  cout << "#hex potentiels recensés par les patterns de recombinaison: " << hex_to_tet.size() << endl;
+};
 
 
-  }
+// TODO: check only the direct neighbors !!! change the algo pour ne as parcourir les hex qui ont un tet commun. cad Filter au niveau des faces... juste les faces extérieures de l'hex !!! et attention aux slivers, du coup... c'est sans doute la source du pb...
 
+// pour chaque hex, sortir un .pos avec l'hex, ses tets et ses faces extérieures pour checker...
+// pourquoi pas sortir aussi tous les slivers à part... ?
 
-  void Recombinator_Graph::evaluate_hex_couple(Hex* hex, Hex* other_hex){
 
-    const bool very_verbose=false;
 
 
+// fills incompatibility_graph if two hex are incompatible direct neighbors,
+// i.e. (they DO NOT have one tet in common) and (they are neighbors (face(s) in common) but DO NOT pass the compatibility tests)
+void Recombinator_Graph::create_direct_neighbors_incompatibility_graph()
+{
 
-    if (very_verbose){
-      //    export_single_hex_all(other_hex,"");
-      //    export_single_hex_all(hex,"");
-      cout << " evaluate_hex_couple:: treating hex " << hex << " made of ";
-      for (int i=0;i<8;i++)
-        cout << "  " << hex->getVertex(i)->getNum();
-      cout << endl;
-      cout << " evaluate_hex_couple:: treating other_hex " << other_hex << " made of ";
-      for (int i=0;i<8;i++)
-        cout << "  " << other_hex->getVertex(i)->getNum();
-      cout << endl;
+  vector<Hex*> visited_hex;
+  std::map<Hex*, std::set<MElement*> >::iterator it_hex =  hex_to_tet.begin();
+  for (;it_hex!=hex_to_tet.end();it_hex++){// for all hex
+    Hex *hex = it_hex->first;
+    if (is_not_good_enough(hex))
+      continue;
+    graph::iterator itfind_graph = find_hex_in_graph(hex);
+    if (itfind_graph==incompatibility_graph.end()){
+      incompatibility_graph.insert(make_pair(hex->get_hash(),make_pair(hex,graph_data())));// this creates an entry if none exists. Important ! if hex is good, but surrounded by "not good enough" hex, it has to be in the graph anyway ! A single node in the graph, without any connection... the perfect lonely hex.
+      set_of_all_hex_in_graph.insert(hex);
     }
 
+    hash_tableA.clear();
+    hash_tableB.clear();
+    hash_tableC.clear();
+    build_hash_tableA(*hex);
+    build_hash_tableB(*hex);
+    build_hash_tableC(*hex);
 
+    visited_hex.clear();
 
-    if (is_not_good_enough(other_hex)){
-      if (very_verbose) cout << "hex " << hex << " not good enough" << endl;
-      return;
-    }
+    std::set<PETriangle*>::iterator it_faces = hex_to_faces[hex].begin();
+    std::set<PETriangle*>::iterator it_facesen = hex_to_faces[hex].end();
+    for (;it_faces!=it_facesen;it_faces++){// i.e. for all triangular external face
+      PETriangle *face = *it_faces;
+      std::set<Hex*>::iterator it_neighbors = faces_to_hex[face].begin();
+      std::set<Hex*>::iterator it_neighborsen = faces_to_hex[face].end();
+      for (;it_neighbors!=it_neighborsen;it_neighbors++){// for all its neighbors
+        if ((*it_neighbors)==hex) continue;
+        Hex *other_hex = *it_neighbors;
 
+        vector<Hex*>::iterator itfind_hex = std::find(visited_hex.begin(),visited_hex.end(),other_hex);
+        if (itfind_hex!=visited_hex.end()) continue;// already done
+        else visited_hex.push_back(other_hex);
 
-    if (find_hex_couple_in_graph(hex,other_hex)){
-      if (very_verbose) cout << "already incompatible: hex's " << hex << " and " << other_hex << endl;
-      return;
+        evaluate_hex_couple(hex,other_hex);
+      }
     }
-
-    // count the number of faces in common... 
-    // intuitivement: si !=2, c'est incompatible
-    //  bool two_faces_in_common=true;
-    //  std::vector<PETriangle*> common(12);
-    //  std::vector<PETriangle*>::iterator itfind_face = std::set_intersection (hex_to_faces[other_hex].begin(),hex_to_faces[other_hex].end(), hex_to_faces[hex].begin(), hex_to_faces[hex].end(), common.begin());
-    //  common.resize(itfind_face-common.begin()); 
-    //  if (common.size()!=2){
-    //    two_faces_in_common=false;
-    //  }
-
-    // check compatibility between the two neighbors
-
-
-    // if they do not pass the tests
-    //-> incompatible !
-
-    // the following depends on hex... on "hash_tables's"
-    //        std::set<MElement*> parts = hex_to_tet[other_hex];
-    //  std::set<MElement*> parts;
-    //  find(other_hex->get_a(),*other_hex,parts);
-    //  find(other_hex->get_b(),*other_hex,parts);
-    //  find(other_hex->get_c(),*other_hex,parts);
-    //  find(other_hex->get_d(),*other_hex,parts);
-    //  find(other_hex->get_e(),*other_hex,parts);
-    //  find(other_hex->get_f(),*other_hex,parts);
-    //  find(other_hex->get_g(),*other_hex,parts);
-    //  find(other_hex->get_h(),*other_hex,parts);
-    //  if(!valid(*other_hex,parts)){
-    //    if (very_verbose) cout << "valid incompatible: other_hex " << other_hex << " relative to " << hex << endl;
-    //    add_graph_entry(hex,other_hex);
-    //    return;
-    //  }
-    //  else if (very_verbose) cout << "valid compatible: other_hex " << other_hex << " relative to " << hex << endl;
-
-    if(!conformityA(*other_hex)){
-      if (very_verbose) cout << "conformA incompatible: other_hex " << other_hex << " relative to " << hex << endl;
-      add_graph_entry(hex,other_hex);
-      return;
+    // change following...
+    std::set<PELine*>::iterator it_line = hex_to_edges[hex].begin();
+    std::set<PELine*>::iterator it_lineen = hex_to_edges[hex].end();
+    for (;it_line!=it_lineen;it_line++){// i.e. for all edges and diagonals -> check the hex neighbors too
+      PELine *line = *it_line;
+      std::set<Hex*>::iterator it_neighbors = edges_to_hex[line].begin();
+      std::set<Hex*>::iterator it_neighborsen = edges_to_hex[line].end();
+      for (;it_neighbors!=it_neighborsen;it_neighbors++){// for all its neighbors
+        if ((*it_neighbors)==hex) continue;
+        Hex *other_hex = *it_neighbors;
+
+        vector<Hex*>::iterator itfind_hex = std::find(visited_hex.begin(),visited_hex.end(),other_hex);
+        if (itfind_hex!=visited_hex.end()) continue;// already done
+        else visited_hex.push_back(other_hex);
+
+        evaluate_hex_couple(hex,other_hex);
+      }
     }
-    else if (very_verbose) cout << "conformA compatible: other_hex " << other_hex << " relative to " << hex << endl;
+  }
 
 
-    if(!conformityB(*other_hex)){
-      if (very_verbose) cout << "conformB incompatible: other_hex " << other_hex << " relative to " << hex << endl;
-      add_graph_entry(hex,other_hex);
-      return;
-    }
-    else if (very_verbose) cout << "conformB compatible: other_hex " << other_hex << " relative to " << hex << endl;
+}
 
-    if(!conformityC(*other_hex)){
-      if (very_verbose) cout << "conformC incompatible: other_hex " << other_hex << " relative to " << hex << endl;
-      add_graph_entry(hex,other_hex);
-      return;
-    }
-    else if (very_verbose) cout << "conformC compatible: other_hex " << other_hex << " relative to " << hex << endl;
 
-    if(!faces_statuquo(*other_hex)){
-      if (very_verbose) cout << "faces_statuquo incompatible: other_hex " << other_hex << " relative to " << hex << endl;
-      add_graph_entry(hex,other_hex);
-      return;
-    }
-    else if (very_verbose) cout << "faces_statuquo compatible: other_hex " << other_hex << " relative to " << hex << endl;
+void Recombinator_Graph::evaluate_hex_couple(Hex* hex, Hex* other_hex)
+{
 
-    if (very_verbose) cout << "other_hex " << other_hex << " relative to " << hex << " totaly compatible !!! " << endl;
-    //  if (!two_faces_in_common){
+  const bool very_verbose=false;
+
+  if (very_verbose){
     //    export_single_hex_all(other_hex,"");
     //    export_single_hex_all(hex,"");
-    //
-    //    //cout << "******************************** not two faces in common and compatible !!!********************************************************** " << endl;
-    //  }
-  }
-
-
-  bool Recombinator_Graph::post_check_validation(Hex* current_hex){
-    // post check...
-    //  std::set<MElement*> parts;
-    //  find(current_hex->get_a(),*current_hex,parts);
-    //  find(current_hex->get_b(),*current_hex,parts);
-    //  find(current_hex->get_c(),*current_hex,parts);
-    //  find(current_hex->get_d(),*current_hex,parts);
-    //  find(current_hex->get_e(),*current_hex,parts);
-    //  find(current_hex->get_f(),*current_hex,parts);
-    //  find(current_hex->get_g(),*current_hex,parts);
-    //  find(current_hex->get_h(),*current_hex,parts);
-    //  if(!valid(*current_hex,parts)){
-    //    std::cout << "     not valid ! : hex " << current_hex  << " made of ";
-    //    for (int i=0;i<8;i++)
-    //      cout << "  " << current_hex->getVertex(i)->getNum();
-    //    cout << endl;
-    //    return false;
-    //  }
-
-    if(!conformityA(*current_hex)){
-      std::cout << "     not conform A! : hex " << current_hex << " made of ";
-      for (int i=0;i<8;i++)
-        cout << "  " << current_hex->getVertex(i)->getNum();
-      cout << endl;
-      return false;
-    }
-
-    if(!conformityB(*current_hex)){
-      std::cout << "     not conform B! : hex " << current_hex  << " made of ";
-      for (int i=0;i<8;i++)
-        cout << "  " << current_hex->getVertex(i)->getNum();
-      cout << endl;
-      return false;
-    }
+    cout << " evaluate_hex_couple:: treating hex " << hex << " made of ";
+    for (int i=0;i<8;i++)
+      cout << "  " << hex->getVertex(i)->getNum();
+    cout << endl;
+    cout << " evaluate_hex_couple:: treating other_hex " << other_hex << " made of ";
+    for (int i=0;i<8;i++)
+      cout << "  " << other_hex->getVertex(i)->getNum();
+    cout << endl;
+  }
 
-    if(!conformityC(*current_hex)){
-      std::cout << "     not conform C! : hex " << current_hex  << " made of ";
-      for (int i=0;i<8;i++)
-        cout << "  " << current_hex->getVertex(i)->getNum();
-      cout << endl;
-      return false;
-    }
+  if (is_not_good_enough(other_hex)){
+    if (very_verbose) cout << "hex " << hex << " not good enough" << endl;
+    return;
+  }
 
-    if(!faces_statuquo(*current_hex)){
-      std::cout << "     not ok faces status quo! : hex " << current_hex << " made of ";
-      for (int i=0;i<8;i++)
-        cout << "  " << current_hex->getVertex(i)->getNum();
-      cout << endl;
-      return false;
-    }
+  if (find_hex_couple_in_graph(hex,other_hex)){
+    if (very_verbose) cout << "already incompatible: hex's " << hex << " and " << other_hex << endl;
+    return;
+  }
 
-    // end post check...
+  // count the number of faces in common...
+  // intuitivement: si !=2, c'est incompatible
+  //  bool two_faces_in_common=true;
+  //  std::vector<PETriangle*> common(12);
+  //  std::vector<PETriangle*>::iterator itfind_face = std::set_intersection (hex_to_faces[other_hex].begin(),hex_to_faces[other_hex].end(), hex_to_faces[hex].begin(), hex_to_faces[hex].end(), common.begin());
+  //  common.resize(itfind_face-common.begin());
+  //  if (common.size()!=2){
+  //    two_faces_in_common=false;
+  //  }
 
+  // check compatibility between the two neighbors
+
+
+  // if they do not pass the tests
+  //-> incompatible !
+
+  // the following depends on hex... on "hash_tables's"
+  //        std::set<MElement*> parts = hex_to_tet[other_hex];
+  //  std::set<MElement*> parts;
+  //  find(other_hex->get_a(),*other_hex,parts);
+  //  find(other_hex->get_b(),*other_hex,parts);
+  //  find(other_hex->get_c(),*other_hex,parts);
+  //  find(other_hex->get_d(),*other_hex,parts);
+  //  find(other_hex->get_e(),*other_hex,parts);
+  //  find(other_hex->get_f(),*other_hex,parts);
+  //  find(other_hex->get_g(),*other_hex,parts);
+  //  find(other_hex->get_h(),*other_hex,parts);
+  //  if(!valid(*other_hex,parts)){
+  //    if (very_verbose) cout << "valid incompatible: other_hex " << other_hex << " relative to " << hex << endl;
+  //    add_graph_entry(hex,other_hex);
+  //    return;
+  //  }
+  //  else if (very_verbose) cout << "valid compatible: other_hex " << other_hex << " relative to " << hex << endl;
 
-    return true;
+  if(!conformityA(*other_hex)){
+    if (very_verbose) cout << "conformA incompatible: other_hex " << other_hex << " relative to " << hex << endl;
+    add_graph_entry(hex,other_hex);
+    return;
   }
+  else if (very_verbose) cout << "conformA compatible: other_hex " << other_hex << " relative to " << hex << endl;
 
 
-  void Recombinator_Graph::add_face(const MVertex *a,const MVertex* b,const MVertex *c,std::multimap<unsigned long long, pair<PETriangle*,int> > &f){
-    vector<const MVertex*> v;
-    v.push_back(a);
-    v.push_back(b);
-    v.push_back(c);
-    PETriangle *q = new PETriangle(v);
-    std::multimap<unsigned long long, pair<PETriangle*,int> >::iterator itfind = find_the_triangle(q,f);
-    if (itfind==f.end()){
-      f.insert(make_pair(q->get_hash(),make_pair(q,1)));
-    }
-    else{
-      delete q;
-    }
+  if(!conformityB(*other_hex)){
+    if (very_verbose) cout << "conformB incompatible: other_hex " << other_hex << " relative to " << hex << endl;
+    add_graph_entry(hex,other_hex);
+    return;
   }
+  else if (very_verbose) cout << "conformB compatible: other_hex " << other_hex << " relative to " << hex << endl;
 
-
-  void Recombinator_Graph::add_face(const MVertex *a,const MVertex* b,const MVertex *c,Hex *hex){
-    vector<const MVertex*> v;
-    v.push_back(a);
-    v.push_back(b);
-    v.push_back(c);
-    PETriangle *q = new PETriangle(v);
-    citer itfind = find_the_triangle(q,triangular_faces);
-    if (itfind==triangular_faces.end()){
-      itfind = triangular_faces.insert(make_pair(q->get_hash(),q));
-    }
-    else{
-      delete q;
-      q = itfind->second;
-    }
-
-    hex_to_faces[hex].insert(q);
-    faces_to_hex[q].insert(hex);
+  if(!conformityC(*other_hex)){
+    if (very_verbose) cout << "conformC incompatible: other_hex " << other_hex << " relative to " << hex << endl;
+    add_graph_entry(hex,other_hex);
+    return;
   }
+  else if (very_verbose) cout << "conformC compatible: other_hex " << other_hex << " relative to " << hex << endl;
 
+  if(!faces_statuquo(*other_hex)){
+    if (very_verbose) cout << "faces_statuquo incompatible: other_hex " << other_hex << " relative to " << hex << endl;
+    add_graph_entry(hex,other_hex);
+    return;
+  }
+  else if (very_verbose) cout << "faces_statuquo compatible: other_hex " << other_hex << " relative to " << hex << endl;
 
-  void Recombinator_Graph::add_edges(Hex *hex){
-    MVertex *a,*b,*c,*d;
-    MVertex *e,*f,*g,*h;
+  if (very_verbose) cout << "other_hex " << other_hex << " relative to " << hex << " totaly compatible !!! " << endl;
+  //  if (!two_faces_in_common){
+  //    export_single_hex_all(other_hex,"");
+  //    export_single_hex_all(hex,"");
+  //
+  //    //cout << "******************************** not two faces in common and compatible !!!********************************************************** " << endl;
+  //  }
+}
 
-    a = hex->get_a();
-    b = hex->get_b();
-    c = hex->get_c();
-    d = hex->get_d();
-    e = hex->get_e();
-    f = hex->get_f();
-    g = hex->get_g();
-    h = hex->get_h();
 
-    fill_edges_table(a,b,c,d,hex);
-    fill_edges_table(e,f,g,h,hex);
-    fill_edges_table(a,b,f,e,hex);
-    fill_edges_table(b,c,g,f,hex);
-    fill_edges_table(d,c,g,h,hex);
-    fill_edges_table(d,a,e,h,hex);
-  }
-
-
-  void Recombinator_Graph::fill_edges_table(const MVertex *a, const MVertex *b, const MVertex *c, const MVertex *d, Hex *hex){
-    vector<const MVertex* > v;
-    vector<const MVertex* > u;
-    v.push_back(a);
-    v.push_back(b);
-    v.push_back(c);
-    v.push_back(d);
-
-    vector<const MVertex*>::const_iterator it1 = v.begin();
-    for (;it1!=v.end();it1++){
-      vector<const MVertex*>::const_iterator it2 = it1;
-      for (;it2!=v.end();it2++){
-        if (it1==it2) continue;
-
-        u.clear();
-        u.push_back(*it1);
-        u.push_back(*it2);
-        // see if already exists or not...
-        PELine *l = new PELine(u);
-        linemap::const_iterator itfind = find_the_line(l,edges_and_diagonals);
-        if (itfind==edges_and_diagonals.end()){
-          itfind = edges_and_diagonals.insert(make_pair(l->get_hash(),l));
-        }
-        else{
-          delete l;
-          l = itfind->second;
-        }
+bool Recombinator_Graph::post_check_validation(Hex* current_hex)
+{
+  // post check...
+  //  std::set<MElement*> parts;
+  //  find(current_hex->get_a(),*current_hex,parts);
+  //  find(current_hex->get_b(),*current_hex,parts);
+  //  find(current_hex->get_c(),*current_hex,parts);
+  //  find(current_hex->get_d(),*current_hex,parts);
+  //  find(current_hex->get_e(),*current_hex,parts);
+  //  find(current_hex->get_f(),*current_hex,parts);
+  //  find(current_hex->get_g(),*current_hex,parts);
+  //  find(current_hex->get_h(),*current_hex,parts);
+  //  if(!valid(*current_hex,parts)){
+  //    std::cout << "     not valid ! : hex " << current_hex  << " made of ";
+  //    for (int i=0;i<8;i++)
+  //      cout << "  " << current_hex->getVertex(i)->getNum();
+  //    cout << endl;
+  //    return false;
+  //  }
 
-        hex_to_edges[hex].insert(l);
-        edges_to_hex[l].insert(hex);
-      }
-    }
+  if(!conformityA(*current_hex)){
+    std::cout << "     not conform A! : hex " << current_hex << " made of ";
+    for (int i=0;i<8;i++)
+      cout << "  " << current_hex->getVertex(i)->getNum();
+    cout << endl;
+    return false;
   }
 
-
-  Recombinator_Graph::graph::iterator Recombinator_Graph::find_hex_in_graph(Hex* hex){
-    pair<graph::iterator, graph::iterator> range = incompatibility_graph.equal_range(hex->get_hash());
-    if (range.first==range.second) return incompatibility_graph.end();
-
-    graph::iterator it = range.first;
-    for (;it!=range.second;it++){
-      if (it->second.first==hex){
-        return it;
-      }
-    }
-    return incompatibility_graph.end();
+  if(!conformityB(*current_hex)){
+    std::cout << "     not conform B! : hex " << current_hex  << " made of ";
+    for (int i=0;i<8;i++)
+      cout << "  " << current_hex->getVertex(i)->getNum();
+    cout << endl;
+    return false;
   }
 
+  if(!conformityC(*current_hex)){
+    std::cout << "     not conform C! : hex " << current_hex  << " made of ";
+    for (int i=0;i<8;i++)
+      cout << "  " << current_hex->getVertex(i)->getNum();
+    cout << endl;
+    return false;
+  }
 
-  Recombinator_Graph::graph_data::iterator Recombinator_Graph::find_hex_in_graphrow(Hex* hex, graph_data &row){
-    pair<graph_data::iterator, graph_data::iterator> range = row.equal_range(hex->get_hash());
-    if (range.first==range.second) return row.end();
-
-    graph_data::iterator it = range.first;
-    for (;it!=range.second;it++){
-      if (it->second==hex){
-        return it;
-      }
-    }
-    return row.end();
+  if(!faces_statuquo(*current_hex)){
+    std::cout << "     not ok faces status quo! : hex " << current_hex << " made of ";
+    for (int i=0;i<8;i++)
+      cout << "  " << current_hex->getVertex(i)->getNum();
+    cout << endl;
+    return false;
   }
 
+  // end post check...
 
 
-  bool Recombinator_Graph::find_hex_couple_in_graph(Hex* hex, Hex* other_hex){
-    graph::iterator it = find_hex_in_graph(hex);
-    if (it==incompatibility_graph.end()) return false;
+  return true;
+}
 
-    graph_data::iterator itt = find_hex_in_graphrow(other_hex, it->second.second);
-    if (itt==it->second.second.end()) return false;
-    return true;
 
+void Recombinator_Graph::add_face(const MVertex *a,const MVertex* b,const MVertex *c,std::multimap<unsigned long long, pair<PETriangle*,int> > &f)
+{
+  vector<const MVertex*> v;
+  v.push_back(a);
+  v.push_back(b);
+  v.push_back(c);
+  PETriangle *q = new PETriangle(v);
+  std::multimap<unsigned long long, pair<PETriangle*,int> >::iterator itfind = find_the_triangle(q,f);
+  if (itfind==f.end()){
+    f.insert(make_pair(q->get_hash(),make_pair(q,1)));
   }
+  else{
+    delete q;
+  }
+}
 
 
-  void Recombinator_Graph::add_graph_entry(Hex* hex, Hex* other_hex){
-
-    graph::iterator itfind_graph = find_hex_in_graph(hex);
-
-    if (itfind_graph==incompatibility_graph.end()){
-      itfind_graph = incompatibility_graph.insert(make_pair(hex->get_hash(),make_pair(hex, graph_data())));
-      itfind_graph->second.second.insert(make_pair(other_hex->get_hash(),other_hex));
-      set_of_all_hex_in_graph.insert(hex);
-      return;
-    }
-    else{
-      itfind_graph->second.second.insert(make_pair(other_hex->get_hash(),other_hex));
-    }
-
+void Recombinator_Graph::add_face(const MVertex *a,const MVertex* b,const MVertex *c,Hex *hex)
+{
+  vector<const MVertex*> v;
+  v.push_back(a);
+  v.push_back(b);
+  v.push_back(c);
+  PETriangle *q = new PETriangle(v);
+  citer itfind = find_the_triangle(q,triangular_faces);
+  if (itfind==triangular_faces.end()){
+    itfind = triangular_faces.insert(make_pair(q->get_hash(),q));
+  }
+  else{
+    delete q;
+    q = itfind->second;
   }
 
+  hex_to_faces[hex].insert(q);
+  faces_to_hex[q].insert(hex);
+}
 
-  void Recombinator_Graph::compute_hex_ranks(){
-
-    create_faces_connectivity();
 
-    PETriangle* face;
-    Hex *hex;
-    double boundary_count;
+void Recombinator_Graph::add_edges(Hex *hex)
+{
+  MVertex *a,*b,*c,*d;
+  MVertex *e,*f,*g,*h;
 
-    for (map<Hex*,set<PETriangle*> >::iterator it = hex_to_faces.begin();it!=hex_to_faces.end();it++){
-      hex = it->first;
-      boundary_count=0.;
-      for (set<PETriangle*>::iterator itf = it->second.begin();itf!=it->second.end();itf++){
-        face = *itf;
-        if (faces_connectivity[face]==1) boundary_count+=1.;
+  a = hex->get_a();
+  b = hex->get_b();
+  c = hex->get_c();
+  d = hex->get_d();
+  e = hex->get_e();
+  f = hex->get_f();
+  g = hex->get_g();
+  h = hex->get_h();
+
+  fill_edges_table(a,b,c,d,hex);
+  fill_edges_table(e,f,g,h,hex);
+  fill_edges_table(a,b,f,e,hex);
+  fill_edges_table(b,c,g,f,hex);
+  fill_edges_table(d,c,g,h,hex);
+  fill_edges_table(d,a,e,h,hex);
+}
+
+
+void Recombinator_Graph::fill_edges_table(const MVertex *a, const MVertex *b, const MVertex *c, const MVertex *d, Hex *hex)
+{
+  vector<const MVertex* > v;
+  vector<const MVertex* > u;
+  v.push_back(a);
+  v.push_back(b);
+  v.push_back(c);
+  v.push_back(d);
+
+  vector<const MVertex*>::const_iterator it1 = v.begin();
+  for (;it1!=v.end();it1++){
+    vector<const MVertex*>::const_iterator it2 = it1;
+    for (;it2!=v.end();it2++){
+      if (it1==it2) continue;
+
+      u.clear();
+      u.push_back(*it1);
+      u.push_back(*it2);
+      // see if already exists or not...
+      PELine *l = new PELine(u);
+      linemap::const_iterator itfind = find_the_line(l,edges_and_diagonals);
+      if (itfind==edges_and_diagonals.end()){
+        itfind = edges_and_diagonals.insert(make_pair(l->get_hash(),l));
+      }
+      else{
+        delete l;
+        l = itfind->second;
       }
-      //cout << " --- hex " << hex << " has " << boundary_count << " tri faces on boundaries, " << hex_to_faces[hex].size() << " total tri faces " << endl;
-      map<Hex*,vector<double> >::iterator itfind = hex_ranks.find(hex);
-      if (itfind==hex_ranks.end())
-        hex_ranks.insert(make_pair(hex,vector<double>(1)));
-      hex_ranks[hex][0] = boundary_count;
-      hex_ranks[hex][1] = hex->get_quality();
+
+      hex_to_edges[hex].insert(l);
+      edges_to_hex[l].insert(hex);
     }
   }
+}
+
 
+Recombinator_Graph::graph::iterator Recombinator_Graph::find_hex_in_graph(Hex* hex)
+{
+  pair<graph::iterator, graph::iterator> range = incompatibility_graph.equal_range(hex->get_hash());
+  if (range.first==range.second) return incompatibility_graph.end();
 
-  void Recombinator_Graph::create_faces_connectivity(){
-    for (std::map<MElement*, std::set<Hex*> >::iterator it_tet = tet_to_hex.begin(); it_tet!=tet_to_hex.end();it_tet++){
-      add_face_connectivity(it_tet->first, 0,1,2);
-      add_face_connectivity(it_tet->first, 0,1,3);
-      add_face_connectivity(it_tet->first, 0,2,3);
-      add_face_connectivity(it_tet->first, 1,2,3);
+  graph::iterator it = range.first;
+  for (;it!=range.second;it++){
+    if (it->second.first==hex){
+      return it;
     }
   }
+  return incompatibility_graph.end();
+}
 
+Recombinator_Graph::graph_data::iterator Recombinator_Graph::find_hex_in_graphrow(Hex* hex, graph_data &row)
+{
+  pair<graph_data::iterator, graph_data::iterator> range = row.equal_range(hex->get_hash());
+  if (range.first==range.second) return row.end();
 
-  void Recombinator_Graph::add_face_connectivity(MElement *tet, int i, int j, int k){
-    vector<const MVertex*> v;
-    PETriangle *t;
-    v.push_back(tet->getVertex(i));
-    v.push_back(tet->getVertex(j));
-    v.push_back(tet->getVertex(k));
-    t = new PETriangle(v);
-    citer itfind = find_the_triangle(t,triangular_faces);
-    if (itfind!=triangular_faces.end()){
-      faces_connectivity[itfind->second]++;
+  graph_data::iterator it = range.first;
+  for (;it!=range.second;it++){
+    if (it->second==hex){
+      return it;
     }
-    delete t;
   }
+  return row.end();
+}
+
+bool Recombinator_Graph::find_hex_couple_in_graph(Hex* hex, Hex* other_hex)
+{
+  graph::iterator it = find_hex_in_graph(hex);
+  if (it==incompatibility_graph.end()) return false;
 
+  graph_data::iterator itt = find_hex_in_graphrow(other_hex, it->second.second);
+  if (itt==it->second.second.end()) return false;
+  return true;
 
-  void Recombinator_Graph::compute_hex_ranks_blossom(){
+}
 
-    create_faces_connectivity();
 
-    PETriangle* face;
-    Hex *hex;
-    double boundary_count;
-    MVertex *a,*b,*c,*d;
-    MVertex *e,*f,*g,*h;
+void Recombinator_Graph::add_graph_entry(Hex* hex, Hex* other_hex)
+{
 
-    for (map<Hex*,set<PETriangle*> >::iterator it = hex_to_faces.begin();it!=hex_to_faces.end();it++){
-      hex = it->first;
-      boundary_count=0.;
-      for (set<PETriangle*>::iterator itf = it->second.begin();itf!=it->second.end();itf++){
-        face = *itf;
-        if (faces_connectivity[face]==1) boundary_count+=1.;
-      }
-      //cout << " --- hex " << hex << " has " << boundary_count << " tri faces on boundaries, " << hex_to_faces[hex].size() << " total tri faces " << endl;
-      map<Hex*,vector<double> >::iterator itfind = hex_ranks.find(hex);
-      if (itfind==hex_ranks.end())
-        hex_ranks.insert(make_pair(hex,vector<double>(1)));
-      hex_ranks[hex][0] = boundary_count;
-      hex_ranks[hex][1] = hex->get_quality();
+  graph::iterator itfind_graph = find_hex_in_graph(hex);
+
+  if (itfind_graph==incompatibility_graph.end()){
+    itfind_graph = incompatibility_graph.insert(make_pair(hex->get_hash(),make_pair(hex, graph_data())));
+    itfind_graph->second.second.insert(make_pair(other_hex->get_hash(),other_hex));
+    set_of_all_hex_in_graph.insert(hex);
+    return;
+  }
+  else{
+    itfind_graph->second.second.insert(make_pair(other_hex->get_hash(),other_hex));
+  }
 
+}
 
-      a = hex->get_a();
-      b = hex->get_b();
-      c = hex->get_c();
-      d = hex->get_d();
-      e = hex->get_e();
-      f = hex->get_f();
-      g = hex->get_g();
-      h = hex->get_h();
+void Recombinator_Graph::compute_hex_ranks()
+{
 
-      int count_blossom=0;
-      if (find_face_in_blossom_info(a,b,c,d)) count_blossom++;
-      if (find_face_in_blossom_info(e,f,g,h)) count_blossom++;
-      if (find_face_in_blossom_info(a,b,f,e)) count_blossom++;
-      if (find_face_in_blossom_info(b,c,g,f)) count_blossom++;
-      if (find_face_in_blossom_info(d,c,g,h)) count_blossom++;
-      if (find_face_in_blossom_info(d,a,e,h)) count_blossom++;
+  create_faces_connectivity();
 
-      hex_ranks[hex][2] = count_blossom;
+  PETriangle* face;
+  Hex *hex;
+  double boundary_count;
 
+  for (map<Hex*,set<PETriangle*> >::iterator it = hex_to_faces.begin();it!=hex_to_faces.end();it++){
+    hex = it->first;
+    boundary_count=0.;
+    for (set<PETriangle*>::iterator itf = it->second.begin();itf!=it->second.end();itf++){
+      face = *itf;
+      if (faces_connectivity[face]==1) boundary_count+=1.;
     }
+    //cout << " --- hex " << hex << " has " << boundary_count << " tri faces on boundaries, " << hex_to_faces[hex].size() << " total tri faces " << endl;
+    map<Hex*,vector<double> >::iterator itfind = hex_ranks.find(hex);
+    if (itfind==hex_ranks.end())
+      hex_ranks.insert(make_pair(hex,vector<double>(1)));
+    hex_ranks[hex][0] = boundary_count;
+    hex_ranks[hex][1] = hex->get_quality();
   }
+}
 
 
-  bool Recombinator_Graph::find_face_in_blossom_info(MVertex *a, MVertex *b, MVertex *c, MVertex *d){
-    PETriangle *t1,*t2;
+void Recombinator_Graph::create_faces_connectivity()
+{
+  for (std::map<MElement*, std::set<Hex*> >::iterator it_tet = tet_to_hex.begin(); it_tet!=tet_to_hex.end();it_tet++){
+    add_face_connectivity(it_tet->first, 0,1,2);
+    add_face_connectivity(it_tet->first, 0,1,3);
+    add_face_connectivity(it_tet->first, 0,2,3);
+    add_face_connectivity(it_tet->first, 1,2,3);
+  }
+}
 
-    t1 = get_triangle(a,b,c);
-    t2 = get_triangle(a,c,d);
 
-    if (is_blossom_pair(t1,t2))
-      return true;
+void Recombinator_Graph::add_face_connectivity(MElement *tet, int i, int j, int k)
+{
+  vector<const MVertex*> v;
+  PETriangle *t;
+  v.push_back(tet->getVertex(i));
+  v.push_back(tet->getVertex(j));
+  v.push_back(tet->getVertex(k));
+  t = new PETriangle(v);
+  citer itfind = find_the_triangle(t,triangular_faces);
+  if (itfind!=triangular_faces.end()){
+    faces_connectivity[itfind->second]++;
+  }
+  delete t;
+}
 
+void Recombinator_Graph::compute_hex_ranks_blossom()
+{
 
-    t1 = get_triangle(a,b,d);
-    t2 = get_triangle(b,c,d);
-    if (is_blossom_pair(t1,t2))
-      return true;
+  create_faces_connectivity();
 
-    return false;
-  }
+  PETriangle* face;
+  Hex *hex;
+  double boundary_count;
+  MVertex *a,*b,*c,*d;
+  MVertex *e,*f,*g,*h;
+
+  for (map<Hex*,set<PETriangle*> >::iterator it = hex_to_faces.begin();it!=hex_to_faces.end();it++){
+    hex = it->first;
+    boundary_count=0.;
+    for (set<PETriangle*>::iterator itf = it->second.begin();itf!=it->second.end();itf++){
+      face = *itf;
+      if (faces_connectivity[face]==1) boundary_count+=1.;
+    }
+    //cout << " --- hex " << hex << " has " << boundary_count << " tri faces on boundaries, " << hex_to_faces[hex].size() << " total tri faces " << endl;
+    map<Hex*,vector<double> >::iterator itfind = hex_ranks.find(hex);
+    if (itfind==hex_ranks.end())
+      hex_ranks.insert(make_pair(hex,vector<double>(1)));
+    hex_ranks[hex][0] = boundary_count;
+    hex_ranks[hex][1] = hex->get_quality();
 
 
-  PETriangle* Recombinator_Graph::get_triangle(MVertex*a, MVertex* b, MVertex *c){
-    vector<const MVertex*> v;
-    v.push_back(a);
-    v.push_back(b);
-    v.push_back(c);
-    PETriangle *t = new PETriangle(v);
-    citer it_find_tri = find_the_triangle(t, triangular_faces);
-    delete t;
-    return (it_find_tri->second);
-  }
+    a = hex->get_a();
+    b = hex->get_b();
+    c = hex->get_c();
+    d = hex->get_d();
+    e = hex->get_e();
+    f = hex->get_f();
+    g = hex->get_g();
+    h = hex->get_h();
 
+    int count_blossom=0;
+    if (find_face_in_blossom_info(a,b,c,d)) count_blossom++;
+    if (find_face_in_blossom_info(e,f,g,h)) count_blossom++;
+    if (find_face_in_blossom_info(a,b,f,e)) count_blossom++;
+    if (find_face_in_blossom_info(b,c,g,f)) count_blossom++;
+    if (find_face_in_blossom_info(d,c,g,h)) count_blossom++;
+    if (find_face_in_blossom_info(d,a,e,h)) count_blossom++;
+
+    hex_ranks[hex][2] = count_blossom;
 
-  bool Recombinator_Graph::is_blossom_pair(PETriangle *t1, PETriangle *t2){
-    tripair::iterator itfind = blossom_info.find(t1);
-    if (itfind!=blossom_info.end()){
-      if (t2==itfind->second)
-        return true;
-    }
-    return false;
   }
+}
 
+bool Recombinator_Graph::find_face_in_blossom_info(MVertex *a, MVertex *b,
+                                                   MVertex *c, MVertex *d)
+{
+  PETriangle *t1,*t2;
 
+  t1 = get_triangle(a,b,c);
+  t2 = get_triangle(a,c,d);
 
+  if (is_blossom_pair(t1,t2))
+    return true;
+
+
+  t1 = get_triangle(a,b,d);
+  t2 = get_triangle(b,c,d);
+  if (is_blossom_pair(t1,t2))
+    return true;
 
+  return false;
+}
 
 
+PETriangle* Recombinator_Graph::get_triangle(MVertex*a, MVertex* b, MVertex *c)
+{
+  vector<const MVertex*> v;
+  v.push_back(a);
+  v.push_back(b);
+  v.push_back(c);
+  PETriangle *t = new PETriangle(v);
+  citer it_find_tri = find_the_triangle(t, triangular_faces);
+  delete t;
+  return (it_find_tri->second);
+}
 
 
+bool Recombinator_Graph::is_blossom_pair(PETriangle *t1, PETriangle *t2)
+{
+  tripair::iterator itfind = blossom_info.find(t1);
+  if (itfind!=blossom_info.end()){
+    if (t2==itfind->second)
+      return true;
+  }
+  return false;
+}
 
diff --git a/Mesh/yamakawa.h b/Mesh/yamakawa.h
index be66e4a094aa0bc0d83be21be9fae8bab4d97d94..2fc54050464c7b685083997c6dbd92b1dc13239a 100644
--- a/Mesh/yamakawa.h
+++ b/Mesh/yamakawa.h
@@ -25,234 +25,234 @@ extern void export_gregion_mesh(GRegion *gr, string filename);
 
 
 class PEEntity{
-  protected:
-    vector<const MVertex *> vertices;
-    size_t hash;
-    void compute_hash();
-  public:
-    PEEntity(const vector<const MVertex*> &_v);
-    //PEEntity(size_t l);
-    ~PEEntity();
-    virtual size_t get_max_nb_vertices() const=0;
-    const MVertex* getVertex(size_t n) const;
-    bool hasVertex(const MVertex *v)const;
-    size_t get_hash() const;
-    bool same_vertices(const PEEntity *t)const;
-    bool operator<(const PEEntity&) const;
-    //bool operator==(const PEEntity&) const;
-    //bool operator==(const size_t&) const;
+protected:
+  vector<const MVertex *> vertices;
+  size_t hash;
+  void compute_hash();
+public:
+  PEEntity(const vector<const MVertex*> &_v);
+  //PEEntity(size_t l);
+  ~PEEntity();
+  virtual size_t get_max_nb_vertices() const=0;
+  const MVertex* getVertex(size_t n) const;
+  bool hasVertex(const MVertex *v)const;
+  size_t get_hash() const;
+  bool same_vertices(const PEEntity *t)const;
+  bool operator<(const PEEntity&) const;
+  //bool operator==(const PEEntity&) const;
+  //bool operator==(const size_t&) const;
 };
 
 class PELine : public PEEntity{
-  public:
-    PELine(const vector<const MVertex*> &_v);
-    ~PELine();
-    size_t get_max_nb_vertices() const;
+public:
+  PELine(const vector<const MVertex*> &_v);
+  virtual ~PELine();
+  size_t get_max_nb_vertices() const;
 };
 
 class PETriangle : public PEEntity{
-  public:
-    PETriangle(const vector<const MVertex*> &_v);
-    //PETriangle(size_t l);
-    ~PETriangle();
-    size_t get_max_nb_vertices() const;
+public:
+  PETriangle(const vector<const MVertex*> &_v);
+  //PETriangle(size_t l);
+  virtual ~PETriangle();
+  size_t get_max_nb_vertices() const;
 };
 
 class PEQuadrangle : public PEEntity{
-  public:
-    PEQuadrangle(const vector<const MVertex*> &_v);
-    //PEQuadrangle(size_t l);
-    ~PEQuadrangle();
-    size_t get_max_nb_vertices() const;
+public:
+  PEQuadrangle(const vector<const MVertex*> &_v);
+  //PEQuadrangle(size_t l);
+  virtual ~PEQuadrangle();
+  size_t get_max_nb_vertices() const;
 };
 
 
 template<class T>
 class clique_stop_criteria{
-  public:
-    //typedef tr1::unordered_set<T> graph_data_no_hash;
-    typedef std::set<T> graph_data_no_hash;
-    clique_stop_criteria(map<T, std::set<MElement*> > &_m, int _i);
-    ~clique_stop_criteria();
-    bool stop(const graph_data_no_hash &clique)const;
-    void export_corresponding_mesh(const graph_data_no_hash &clique)const;
-
-  private:
-    const map<T, std::set<MElement*> > &hex_to_tet;
-    const unsigned int total_number_tet;
+public:
+  //typedef tr1::unordered_set<T> graph_data_no_hash;
+  typedef std::set<T> graph_data_no_hash;
+  clique_stop_criteria(map<T, std::set<MElement*> > &_m, int _i);
+  ~clique_stop_criteria();
+  bool stop(const graph_data_no_hash &clique)const;
+  void export_corresponding_mesh(const graph_data_no_hash &clique)const;
+
+private:
+  const map<T, std::set<MElement*> > &hex_to_tet;
+  const unsigned int total_number_tet;
 };
 
 
 template<class T>
 class cliques_compatibility_graph{
-  public:
-    typedef unsigned long long hash_key;
-
-//    typedef  set<T> graph_data;
-//    typedef  map<T, graph_data > graph;
-//    typedef multimap<int,T> ranking_data;
-    //typedef tr1::unordered_set<T> graph_data_no_hash;
-    typedef std::set<T> graph_data_no_hash;
-//    typedef tr1::unordered_multimap<hash_key, T> graph_data;
-//    typedef tr1::unordered_multimap<hash_key, pair<T, graph_data > > graph;
-//    typedef tr1::unordered_map<int,T> ranking_data;
-    typedef std::multimap<hash_key, T> graph_data;
-    typedef std::multimap<hash_key, pair<T, graph_data > > graph;
-    typedef std::map<int,T> ranking_data;
-
-    typedef void (*ptrfunction_export)(cliques_compatibility_graph<T>&, int, string);
-
-    cliques_compatibility_graph(graph &_g, const map<T, std::vector<double> > &_hex_ranks, unsigned int _max_nb_cliques, unsigned int _nb_hex_potentiels, clique_stop_criteria<T> *csc, ptrfunction_export fct);
-    ~cliques_compatibility_graph();
-    void find_cliques();
-    void export_cliques();
-
-    virtual typename graph::const_iterator begin_graph(){return G.begin();};
-    virtual typename graph::const_iterator end_graph(){return G.end();};
-
-    bool found_the_ultimate_max_clique;
-
-    multimap<int, set<T> > allQ;// all cliques
-
-  protected:
-    void erase_entry(graph_data &s, T &u, hash_key &key);
-    void find_cliques(graph_data &s,int n);
-    void split_set_BW(const T &u, const hash_key &u_key,const graph_data &s, graph_data &white, graph_data &black);
-    void fill_black_set(const T &u, const hash_key &u_key, const graph_data &s, graph_data &black);
-    void choose_u(const graph_data &s, T &u, hash_key &u_key);
-    // the maximum score (int) will be chosen...
-    double function_to_maximize_for_u(const T &u, const hash_key &u_key, const graph_data &s);
-    void store_clique(int n);
-    // returns true if two nodes are connected in the compatibility graph
-    virtual bool compatibility(const T &u, const hash_key &u_key, const T &v, const hash_key &v_key);
-
-    ptrfunction_export export_clique_graph;
-
-    const bool debug;
-    unsigned int max_nb_cliques;
-    unsigned int nb_hex_potentiels;
-    unsigned int max_clique_size;
-    unsigned int position;
-    unsigned int total_nodes_number;
-    unsigned int total_nb_of_cliques_searched;
-    unsigned int max_nb_of_stored_cliques;// to reduce memory footprint (set to zero if no limit)
-    
-    clique_stop_criteria<T>* criteria;
-
-    bool cancel_search;
-    const map<T, std::vector<double> > &hex_ranks;
-    graph &G;
-    graph_data_no_hash Q;// the current clique
+public:
+  typedef unsigned long long hash_key;
+
+  //    typedef  set<T> graph_data;
+  //    typedef  map<T, graph_data > graph;
+  //    typedef multimap<int,T> ranking_data;
+  //typedef tr1::unordered_set<T> graph_data_no_hash;
+  typedef std::set<T> graph_data_no_hash;
+  //    typedef tr1::unordered_multimap<hash_key, T> graph_data;
+  //    typedef tr1::unordered_multimap<hash_key, pair<T, graph_data > > graph;
+  //    typedef tr1::unordered_map<int,T> ranking_data;
+  typedef std::multimap<hash_key, T> graph_data;
+  typedef std::multimap<hash_key, pair<T, graph_data > > graph;
+  typedef std::map<int,T> ranking_data;
+
+  typedef void (*ptrfunction_export)(cliques_compatibility_graph<T>&, int, string);
+
+  cliques_compatibility_graph(graph &_g, const map<T, std::vector<double> > &_hex_ranks, unsigned int _max_nb_cliques, unsigned int _nb_hex_potentiels, clique_stop_criteria<T> *csc, ptrfunction_export fct);
+  ~cliques_compatibility_graph();
+  void find_cliques();
+  void export_cliques();
+
+  virtual typename graph::const_iterator begin_graph(){return G.begin();};
+  virtual typename graph::const_iterator end_graph(){return G.end();};
+
+  bool found_the_ultimate_max_clique;
+
+  multimap<int, set<T> > allQ;// all cliques
+
+protected:
+  void erase_entry(graph_data &s, T &u, hash_key &key);
+  void find_cliques(graph_data &s,int n);
+  void split_set_BW(const T &u, const hash_key &u_key,const graph_data &s, graph_data &white, graph_data &black);
+  void fill_black_set(const T &u, const hash_key &u_key, const graph_data &s, graph_data &black);
+  void choose_u(const graph_data &s, T &u, hash_key &u_key);
+  // the maximum score (int) will be chosen...
+  double function_to_maximize_for_u(const T &u, const hash_key &u_key, const graph_data &s);
+  void store_clique(int n);
+  // returns true if two nodes are connected in the compatibility graph
+  virtual bool compatibility(const T &u, const hash_key &u_key, const T &v, const hash_key &v_key);
+
+  ptrfunction_export export_clique_graph;
+
+  const bool debug;
+  unsigned int max_nb_cliques;
+  unsigned int nb_hex_potentiels;
+  unsigned int max_clique_size;
+  unsigned int position;
+  unsigned int total_nodes_number;
+  unsigned int total_nb_of_cliques_searched;
+  unsigned int max_nb_of_stored_cliques;// to reduce memory footprint (set to zero if no limit)
+
+  clique_stop_criteria<T>* criteria;
+
+  bool cancel_search;
+  const map<T, std::vector<double> > &hex_ranks;
+  graph &G;
+  graph_data_no_hash Q;// the current clique
 };
 
 
 
 template<class T>
 class cliques_losses_graph : public cliques_compatibility_graph<T> {
-//    typedef  set<T> graph_data;
-//    typedef  map<T, graph_data > graph;
-//    typedef tr1::unordered_set<T> graph_data;
-//    typedef tr1::unordered_map<T, graph_data > graph;
-  public:
-    typedef unsigned long long hash_key;
-    typedef multimap<hash_key, T> graph_data;
-    typedef multimap<hash_key, pair<T, graph_data > > graph;
-//    typedef tr1::unordered_multimap<hash_key, T> graph_data;
-    typedef void (*ptrfunction_export)(cliques_compatibility_graph<T>&, int, string);
-
-    cliques_losses_graph(graph &_g, const map<T, std::vector<double> > &_hex_ranks, unsigned int _max_nb_cliques, unsigned int _nb_hex_potentiels, clique_stop_criteria<T> *csc, ptrfunction_export fct);
-    ~cliques_losses_graph();
-
-  protected:
-    // returns false if two nodes are connected in the losses graph (i.e. true if connected in compatibility graph)
-    virtual bool compatibility(const T &u, const hash_key &u_key, const T &v, const hash_key &v_key);
-    graph &G;
+  //    typedef  set<T> graph_data;
+  //    typedef  map<T, graph_data > graph;
+  //    typedef tr1::unordered_set<T> graph_data;
+  //    typedef tr1::unordered_map<T, graph_data > graph;
+public:
+  typedef unsigned long long hash_key;
+  typedef multimap<hash_key, T> graph_data;
+  typedef multimap<hash_key, pair<T, graph_data > > graph;
+  //    typedef tr1::unordered_multimap<hash_key, T> graph_data;
+  typedef void (*ptrfunction_export)(cliques_compatibility_graph<T>&, int, string);
+
+  cliques_losses_graph(graph &_g, const map<T, std::vector<double> > &_hex_ranks, unsigned int _max_nb_cliques, unsigned int _nb_hex_potentiels, clique_stop_criteria<T> *csc, ptrfunction_export fct);
+  ~cliques_losses_graph();
+
+protected:
+  // returns false if two nodes are connected in the losses graph (i.e. true if connected in compatibility graph)
+  virtual bool compatibility(const T &u, const hash_key &u_key, const T &v, const hash_key &v_key);
+  graph &G;
 };
 
 
 
 class Hex{
-  private:
-    double quality;
-    unsigned long long hash;
-    MVertex *a,*b,*c,*d,*e,*f,*g,*h;
-    void set_hash();
-  public:
-    Hex();
-    Hex(MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*);
-    ~Hex();
-    double get_quality();
-    void set_quality(double);
-    MVertex* get_a();
-    MVertex* get_b();
-    MVertex* get_c();
-    MVertex* get_d();
-    MVertex* get_e();
-    MVertex* get_f();
-    MVertex* get_g();
-    MVertex* get_h();
-    MVertex* getVertex(int n);
-    bool hasVertex(const MVertex *v);
-    bool same_vertices(Hex *h);
-    void set_vertices(MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*);
-    unsigned long long get_hash();
-    bool operator<( Hex&) ;
+private:
+  double quality;
+  unsigned long long hash;
+  MVertex *a,*b,*c,*d,*e,*f,*g,*h;
+  void set_hash();
+public:
+  Hex();
+  Hex(MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*);
+  ~Hex();
+  double get_quality();
+  void set_quality(double);
+  MVertex* get_a();
+  MVertex* get_b();
+  MVertex* get_c();
+  MVertex* get_d();
+  MVertex* get_e();
+  MVertex* get_f();
+  MVertex* get_g();
+  MVertex* get_h();
+  MVertex* getVertex(int n);
+  bool hasVertex(const MVertex *v);
+  bool same_vertices(Hex *h);
+  void set_vertices(MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*);
+  unsigned long long get_hash();
+  bool operator<( Hex&) ;
 };
 
 class Facet{
-  private:
-    MVertex *a,*b,*c;
-    unsigned long long hash;
-  public:
-    Facet();
-    Facet(MVertex*,MVertex*,MVertex*);
-    ~Facet();
-    MVertex* get_a();
-    MVertex* get_b();
-    MVertex* get_c();
-    void set_vertices(MVertex*,MVertex*,MVertex*);
-    bool same_vertices(Facet);
-    void compute_hash();
-    unsigned long long get_hash() const;
-    bool operator<(const Facet&) const;
+private:
+  MVertex *a,*b,*c;
+  unsigned long long hash;
+public:
+  Facet();
+  Facet(MVertex*,MVertex*,MVertex*);
+  ~Facet();
+  MVertex* get_a();
+  MVertex* get_b();
+  MVertex* get_c();
+  void set_vertices(MVertex*,MVertex*,MVertex*);
+  bool same_vertices(Facet);
+  void compute_hash();
+  unsigned long long get_hash() const;
+  bool operator<(const Facet&) const;
 };
 
 class Diagonal{
-  private:
-    MVertex *a,*b;
-    unsigned long long hash;
-  public:
-    Diagonal();
-    Diagonal(MVertex*,MVertex*);
-    ~Diagonal();
-    MVertex* get_a();
-    MVertex* get_b();
-    void set_vertices(MVertex*,MVertex*);
-    bool same_vertices(Diagonal);
-    void compute_hash();
-    unsigned long long get_hash() const;
-    bool operator<(const Diagonal&) const;
+private:
+  MVertex *a,*b;
+  unsigned long long hash;
+public:
+  Diagonal();
+  Diagonal(MVertex*,MVertex*);
+  ~Diagonal();
+  MVertex* get_a();
+  MVertex* get_b();
+  void set_vertices(MVertex*,MVertex*);
+  bool same_vertices(Diagonal);
+  void compute_hash();
+  unsigned long long get_hash() const;
+  bool operator<(const Diagonal&) const;
 };
 
 class Tuple{
-  private:
-    MVertex *v1,*v2,*v3;
-    MElement* element;
-    GFace* gf;
-    unsigned long long hash;
-  public:
-    Tuple();
-    Tuple(MVertex*,MVertex*,MVertex*,MElement*,GFace*);
-    Tuple(MVertex*,MVertex*,MVertex*);
-    ~Tuple();
-    MVertex* get_v1();
-    MVertex* get_v2();
-    MVertex* get_v3();
-    MElement* get_element() const;
-    GFace* get_gf() const;
-    bool same_vertices(Tuple);
-    unsigned long long get_hash() const;
-    bool operator<(const Tuple&) const;
+private:
+  MVertex *v1,*v2,*v3;
+  MElement* element;
+  GFace* gf;
+  unsigned long long hash;
+public:
+  Tuple();
+  Tuple(MVertex*,MVertex*,MVertex*,MElement*,GFace*);
+  Tuple(MVertex*,MVertex*,MVertex*);
+  ~Tuple();
+  MVertex* get_v1();
+  MVertex* get_v2();
+  MVertex* get_v3();
+  MElement* get_element() const;
+  GFace* get_gf() const;
+  bool same_vertices(Tuple);
+  unsigned long long get_hash() const;
+  bool operator<(const Tuple&) const;
 };
 
 //inline std::ostream& operator<<(std::ostream& s, const PETriangle& t){
@@ -265,116 +265,116 @@ class Tuple{
 //};
 
 class Recombinator{
-  protected:
-    std::vector<Hex*> potential;
-    std::map<MElement*,bool> markings;
-    std::multiset<Facet> hash_tableA;
-    std::multiset<Diagonal> hash_tableB;
-    std::multiset<Diagonal> hash_tableC;
-    std::multiset<Tuple> tuples;
-    std::set<MElement*> triangles;
-  public:
-    Recombinator();
-    ~Recombinator();
-
-    std::map<MVertex*,std::set<MVertex*> > vertex_to_vertices;
-    std::map<MVertex*,std::set<MElement*> > vertex_to_elements;
-
-    virtual void execute();
-    virtual void execute(GRegion*);
-
-    void init_markings(GRegion*);
-    virtual void pattern1(GRegion*);
-    virtual void pattern2(GRegion*);
-    virtual void pattern3(GRegion*);
-    virtual void merge(GRegion*);
-    void improved_merge(GRegion*);
-    void rearrange(GRegion*);
-    void statistics(GRegion*);
-    // tuples are triangles on geometrical faces (region boundaries...)
-    void build_tuples(GRegion*);
-    void modify_surfaces(GRegion*);
-    void modify_surfaces(MVertex*,MVertex*,MVertex*,MVertex*);
-
-    bool sliver(MElement*,Hex&);
-    double diagonal(MElement*,int&,int&);
-    double distance(MVertex*,MVertex*);
-    double distance(MVertex*,MVertex*,MVertex*);
-    double scalar(MVertex*,MVertex*,MVertex*,MVertex*);
-    void two_others(int,int,int&,int&);
-    // soit une face du cube: abcd
-    // en principe, on doit avoir soit les facets (abc) et (acd), soit les facets (abd) et(bcd) qui sont inclues dans un des tets qui forment l'hex.
-    // si c'est le cas pour toutes les 6 faces de l'hex, return true.
-    // ce test permet probablement de virer les hex "avec des trous" (avec 8 noeuds ok, mais un tet manquant, ce qui peut occasionner un hex à 14 faces, par exemple, si l'on compte les faces à partir des tets inclus)
-    bool valid(Hex&,const std::set<MElement*>&);
-    // renvoie true si le "MQuadrangle::etaShapeMeasure" des 6 faces est plus grand que 0.000001
-    bool valid(Hex&);
-    double eta(MVertex*,MVertex*,MVertex*,MVertex*);
-    bool linked(MVertex*,MVertex*);
-
-    void find(MVertex*,MVertex*,const std::vector<MVertex*>&,std::set<MVertex*>&);
-    void find(MVertex*,MVertex*,MVertex*,const std::vector<MVertex*>&,std::set<MVertex*>&);
-    void find(MVertex*,MVertex*,std::set<MElement*>&);
-    void find(MVertex*,Hex,std::set<MElement*>&);
-    MVertex* find(MVertex*,MVertex*,MVertex*,MVertex*,const std::set<MElement*>&);
-
-    void intersection(const std::set<MVertex*>&,const std::set<MVertex*>&,const std::vector<MVertex*>&,std::set<MVertex*>&);
-    void intersection(const std::set<MVertex*>&,const std::set<MVertex*>&,const std::set<MVertex*>&,const std::vector<MVertex*>&,std::set<MVertex*>&);
-    void intersection(const std::set<MElement*>&,const std::set<MElement*>&,std::set<MElement*>&);
-
-    // return true if vertex belong to hex
-    bool inclusion(MVertex*,Hex);
-    // renvoie true si vertex se trouve dans [a,b,c]
-    bool inclusion(MVertex*,MVertex*,MVertex*,MVertex*,MVertex*);
-    // return true if all three vertices v1,v2 and v3 belong to one tet
-    bool inclusion(MVertex*,MVertex*,MVertex*,const std::set<MElement*>&);
-    // return true si la facet existe dans la table A
-    bool inclusion(Facet);
-    // return true si la diagonal existe dans la table B
-    bool inclusion(Diagonal);
-    // return true si la diagonal existe dans la table C !!!!!!!!!   Sinon, c'est exactement la même fonction !!!!! avec un nom différent ! 
-    bool duplicate(Diagonal);
-
-    // return true si un hex est "conforme A"
-    // est "conforme A" un hex dont les 6 faces sont "conforme A"
-    bool conformityA(Hex&);
-    // est "conforme A" une face si ses 4 facets existent dans tableA, ou bien si aucune des ses facets ne se trouve dans table A
-    bool conformityA(MVertex*,MVertex*,MVertex*,MVertex*);
-// return false si:
-//- une des 12 arrêtes de l'hex se trouve dans tableB !!! (pas C !!!), càd si une arrete a été utilisée comme diagonale d'un autre hex 
-//- (ou bien) si, pour chaque face de l'hex, on a une diagonale dans tableB et pas l'autre 
-    bool conformityB(Hex&);
-// return false si une des 12 diagonales du cube se trouve dans tableC, càd a été utilisée comme arrête
-    bool conformityC(Hex&);
-
-    // return true si les 6 faces de l'hex sont "faces_statuquo"
-    bool faces_statuquo(Hex&);
-    // return false si, parmis les deux paires de facets de la face, il existe un couple de facet qui soient toutes les deux des tuples, mais correspondant à des geometric faces différentes. Bref, une arrête géométrique confondue avec une diagonale de la face.
-    bool faces_statuquo(MVertex*,MVertex*,MVertex*,MVertex*);
-
-    void build_vertex_to_vertices(GRegion*);
-    void build_vertex_to_elements(GRegion*);
-    void build_hash_tableA(Hex);
-    void build_hash_tableA(MVertex*,MVertex*,MVertex*,MVertex*);
-    void build_hash_tableA(Facet);
-    void build_hash_tableB(Hex);
-    void build_hash_tableB(MVertex*,MVertex*,MVertex*,MVertex*);
-    void build_hash_tableB(Diagonal);
-    void build_hash_tableC(Hex);
-    void build_hash_tableC(Diagonal);
-
-    void print_vertex_to_vertices(GRegion*);
-    void print_vertex_to_elements(GRegion*);
-    void print_hash_tableA();
-    void print_segment(SPoint3,SPoint3,std::ofstream&);
-
-    double scaled_jacobian(MVertex*,MVertex*,MVertex*,MVertex*);
-    double max_scaled_jacobian(MElement*,int&);
-    double min_scaled_jacobian(Hex&);
+protected:
+  std::vector<Hex*> potential;
+  std::map<MElement*,bool> markings;
+  std::multiset<Facet> hash_tableA;
+  std::multiset<Diagonal> hash_tableB;
+  std::multiset<Diagonal> hash_tableC;
+  std::multiset<Tuple> tuples;
+  std::set<MElement*> triangles;
+public:
+  Recombinator();
+  ~Recombinator();
+
+  std::map<MVertex*,std::set<MVertex*> > vertex_to_vertices;
+  std::map<MVertex*,std::set<MElement*> > vertex_to_elements;
+
+  virtual void execute();
+  virtual void execute(GRegion*);
+
+  void init_markings(GRegion*);
+  virtual void pattern1(GRegion*);
+  virtual void pattern2(GRegion*);
+  virtual void pattern3(GRegion*);
+  virtual void merge(GRegion*);
+  void improved_merge(GRegion*);
+  void rearrange(GRegion*);
+  void statistics(GRegion*);
+  // tuples are triangles on geometrical faces (region boundaries...)
+  void build_tuples(GRegion*);
+  void modify_surfaces(GRegion*);
+  void modify_surfaces(MVertex*,MVertex*,MVertex*,MVertex*);
+
+  bool sliver(MElement*,Hex&);
+  double diagonal(MElement*,int&,int&);
+  double distance(MVertex*,MVertex*);
+  double distance(MVertex*,MVertex*,MVertex*);
+  double scalar(MVertex*,MVertex*,MVertex*,MVertex*);
+  void two_others(int,int,int&,int&);
+  // soit une face du cube: abcd
+  // en principe, on doit avoir soit les facets (abc) et (acd), soit les facets (abd) et(bcd) qui sont inclues dans un des tets qui forment l'hex.
+  // si c'est le cas pour toutes les 6 faces de l'hex, return true.
+  // ce test permet probablement de virer les hex "avec des trous" (avec 8 noeuds ok, mais un tet manquant, ce qui peut occasionner un hex à 14 faces, par exemple, si l'on compte les faces à partir des tets inclus)
+  bool valid(Hex&,const std::set<MElement*>&);
+  // renvoie true si le "MQuadrangle::etaShapeMeasure" des 6 faces est plus grand que 0.000001
+  bool valid(Hex&);
+  double eta(MVertex*,MVertex*,MVertex*,MVertex*);
+  bool linked(MVertex*,MVertex*);
+
+  void find(MVertex*,MVertex*,const std::vector<MVertex*>&,std::set<MVertex*>&);
+  void find(MVertex*,MVertex*,MVertex*,const std::vector<MVertex*>&,std::set<MVertex*>&);
+  void find(MVertex*,MVertex*,std::set<MElement*>&);
+  void find(MVertex*,Hex,std::set<MElement*>&);
+  MVertex* find(MVertex*,MVertex*,MVertex*,MVertex*,const std::set<MElement*>&);
+
+  void intersection(const std::set<MVertex*>&,const std::set<MVertex*>&,const std::vector<MVertex*>&,std::set<MVertex*>&);
+  void intersection(const std::set<MVertex*>&,const std::set<MVertex*>&,const std::set<MVertex*>&,const std::vector<MVertex*>&,std::set<MVertex*>&);
+  void intersection(const std::set<MElement*>&,const std::set<MElement*>&,std::set<MElement*>&);
+
+  // return true if vertex belong to hex
+  bool inclusion(MVertex*,Hex);
+  // renvoie true si vertex se trouve dans [a,b,c]
+  bool inclusion(MVertex*,MVertex*,MVertex*,MVertex*,MVertex*);
+  // return true if all three vertices v1,v2 and v3 belong to one tet
+  bool inclusion(MVertex*,MVertex*,MVertex*,const std::set<MElement*>&);
+  // return true si la facet existe dans la table A
+  bool inclusion(Facet);
+  // return true si la diagonal existe dans la table B
+  bool inclusion(Diagonal);
+  // return true si la diagonal existe dans la table C !!!!!!!!!   Sinon, c'est exactement la même fonction !!!!! avec un nom différent !
+  bool duplicate(Diagonal);
+
+  // return true si un hex est "conforme A"
+  // est "conforme A" un hex dont les 6 faces sont "conforme A"
+  bool conformityA(Hex&);
+  // est "conforme A" une face si ses 4 facets existent dans tableA, ou bien si aucune des ses facets ne se trouve dans table A
+  bool conformityA(MVertex*,MVertex*,MVertex*,MVertex*);
+  // return false si:
+  //- une des 12 arrêtes de l'hex se trouve dans tableB !!! (pas C !!!), càd si une arrete a été utilisée comme diagonale d'un autre hex
+  //- (ou bien) si, pour chaque face de l'hex, on a une diagonale dans tableB et pas l'autre
+  bool conformityB(Hex&);
+  // return false si une des 12 diagonales du cube se trouve dans tableC, càd a été utilisée comme arrête
+  bool conformityC(Hex&);
+
+  // return true si les 6 faces de l'hex sont "faces_statuquo"
+  bool faces_statuquo(Hex&);
+  // return false si, parmis les deux paires de facets de la face, il existe un couple de facet qui soient toutes les deux des tuples, mais correspondant à des geometric faces différentes. Bref, une arrête géométrique confondue avec une diagonale de la face.
+  bool faces_statuquo(MVertex*,MVertex*,MVertex*,MVertex*);
+
+  void build_vertex_to_vertices(GRegion*);
+  void build_vertex_to_elements(GRegion*);
+  void build_hash_tableA(Hex);
+  void build_hash_tableA(MVertex*,MVertex*,MVertex*,MVertex*);
+  void build_hash_tableA(Facet);
+  void build_hash_tableB(Hex);
+  void build_hash_tableB(MVertex*,MVertex*,MVertex*,MVertex*);
+  void build_hash_tableB(Diagonal);
+  void build_hash_tableC(Hex);
+  void build_hash_tableC(Diagonal);
+
+  void print_vertex_to_vertices(GRegion*);
+  void print_vertex_to_elements(GRegion*);
+  void print_hash_tableA();
+  void print_segment(SPoint3,SPoint3,std::ofstream&);
+
+  double scaled_jacobian(MVertex*,MVertex*,MVertex*,MVertex*);
+  double max_scaled_jacobian(MElement*,int&);
+  double min_scaled_jacobian(Hex&);
 };
 
 class Recombinator_Graph : public Recombinator{
-  public:
+public:
   typedef size_t my_hash_key;
   typedef multimap<my_hash_key,PETriangle*> trimap;
   typedef map<PETriangle*, PETriangle*> tripair;
@@ -383,7 +383,7 @@ class Recombinator_Graph : public Recombinator{
   //typedef tr1::unordered_multimap<my_hash_key,PETriangle*> trimap;
   typedef trimap::iterator iter;
   typedef trimap::const_iterator citer;
-  
+
   typedef multimap<my_hash_key,PELine*> linemap;
   //typedef tr1::unordered_multimap<my_hash_key,PELine*> linemap;
 
@@ -392,224 +392,224 @@ class Recombinator_Graph : public Recombinator{
 
   set<Hex*>& getHexInGraph(){return set_of_all_hex_in_graph;};
 
-  protected:
-    bool debug,debug_graph;
-    std::map<Hex*, std::set<MElement*> > hex_to_tet;
-    std::map<Hex*, std::set<PELine*> > hex_to_edges;
-    std::map<PELine*, std::set<Hex*> > edges_to_hex;
-    std::map<Hex*, std::set<PETriangle*> > hex_to_faces;
-    std::map<PETriangle*, std::set<Hex*> > faces_to_hex;
-    std::map<PETriangle*, unsigned int > faces_connectivity;// # of adjacent tets (1 or 2) 
-    std::map<MElement*, std::set<Hex*> >tet_to_hex;
-    std::map<Hex*, std::vector<double> > hex_ranks;
-
-
-    typedef unsigned long long hash_key;
-//    typedef tr1::unordered_multimap<hash_key, Hex*> graph_data;
-//    typedef tr1::unordered_multimap<hash_key, pair<Hex*, graph_data > > graph;
-    typedef multimap<hash_key, Hex*> graph_data;
-    typedef multimap<hash_key, pair<Hex*, graph_data > > graph;
-
-    graph incompatibility_graph;
-    set<Hex*> set_of_all_hex_in_graph;
-
-    std::multimap<unsigned long long, Hex*>created_potential_hex;
-   
-    void create_faces_connectivity();
-    void add_face_connectivity(MElement *tet, int i, int j, int k);
-
-    void add_edges(Hex *hex);
-    void fill_edges_table(const MVertex *a, const MVertex *b, const MVertex *c, const MVertex *d, Hex *hex);
-    void add_face(const MVertex *a,const MVertex* b,const MVertex *c,Hex *hex);
-    void add_face(const MVertex *a,const MVertex* b,const MVertex *c,std::multimap<unsigned long long, pair<PETriangle*,int> > &f);
-    std::multimap<double,Hex*> degree;// degree = the final ranking of hexahedra
-    std::multimap<int,Hex*> idegree;// idegree = number of connected hex in indirect neighbors graph
-    std::multimap<int,Hex*> ndegree;// ndegree = number of direct neighbors !!! not chosen yet !!!
-    std::map<Hex*,int> reverse_idegree;
-    std::map<Hex*,int> reverse_ndegree;
-    // each tet has at least one neighbor, at most four. For all not chosen hex, check this data to find how many direct neighbors...
-//    std::map<MElement*,set<PETriangle*> > tet_to_triangle;
-    std::map<PETriangle*,set<MElement*> > triangle_to_tet;
-    std::map<MElement*,int> tet_degree;
-
-    bool find_face_in_blossom_info(MVertex *a, MVertex *b, MVertex *c, MVertex *d);
-    void compute_hex_ranks_blossom();
-    PETriangle* get_triangle(MVertex*a, MVertex* b, MVertex *c);
-    bool is_blossom_pair(PETriangle *t1, PETriangle *t2);
-
-    tripair blossom_info;
-    trimap triangular_faces;
-    linemap edges_and_diagonals;
-
-    map<PETriangle*, GFace*> tri_to_gface_info;
-
-    vector<Hex*> chosen_hex;
-    vector<MElement*> chosen_tet;
-
-    citer find_the_triangle(PETriangle *t, const trimap &list);
-    linemap::const_iterator  find_the_line(PELine *t, const linemap &list);
-    std::multimap<unsigned long long, pair<PETriangle*,int> >::iterator  find_the_triangle(PETriangle *t, std::multimap<unsigned long long, pair<PETriangle*, int> > &list);
-    std::multimap<unsigned long long, Hex* >::const_iterator  find_the_created_potential_hex(Hex *t, const std::multimap<unsigned long long, Hex*>  &list);
-
-    int nbhex_in_losses_graph;
-    double average_connectivity;
-    bool post_check_validation(Hex* current_hex);
-
-    PETriangle* get_triangle(MElement *element, int i, int j, int k);
-
-    void compute_hex_ranks();
-
-    // check if the hex is good enough to be put into the graph. If not in the graph, it cannot be chosen... 
-    bool is_not_good_enough(Hex* hex);
-
-    // fills incompatibility_graph if two hex share a common (non-sliver!) tet 
-    void create_indirect_neighbors_graph();
-
-    graph::iterator find_hex_in_graph(Hex* hex);
-    graph_data::iterator find_hex_in_graphrow(Hex* hex, graph_data &row);
-    bool find_hex_couple_in_graph(Hex* hex, Hex* other_hex);
-    void add_graph_entry(Hex* hex, Hex* other_hex);
-
-    // fills incompatibility_graph if two hex are incompatible direct neighbors, 
-    // i.e. they have one (or more) common face or common edge and are not compatible 
-    void create_direct_neighbors_incompatibility_graph();
-    void evaluate_hex_couple(Hex* hex, Hex* other_hex);
-
-    // if two hex are not connected in the incompatibility_graph, they are compatible
-    void create_losses_graph(GRegion *gr);
-
-    void merge_clique(GRegion* gr, cliques_losses_graph<Hex*> &cl,int clique_number=0);
-
-    void fill_tet_to_hex_table(Hex *hex);
-
-    virtual void pattern1(GRegion*);
-    virtual void pattern2(GRegion*);
-    virtual void pattern3(GRegion*);
-
-    void merge(GRegion*);
-
-    // ------- exports --------
-    void export_tets(set<MElement*> &tetset, Hex* hex, string s);
-    void export_single_hex_all(Hex* hex,string s);
-    void export_single_hex(Hex* hex,string s);
-    void export_single_hex_faces(Hex* hex,string s);
-    void export_single_hex_tet(Hex* hex,string s);
-    void export_all_hex(int &file,GRegion *gr);
-    void export_hexmesh_so_far(int &file);
-    void export_direct_neighbor_table(int max);
-    void export_hex_init_degree(GRegion *gr, const std::map<Hex*,int> &init_degree, const vector<Hex*> &chosen_hex);
-
-    int max_nb_cliques;
-    string graphfilename;
-
-  public:
-    Recombinator_Graph(unsigned int max_nb_cliques, string filename=string());
-    ~Recombinator_Graph();
-    virtual void execute();
-    virtual void execute(GRegion*);
-    virtual void buildGraphOnly(unsigned int max_nb_cliques, string filename=string());
-    virtual void buildGraphOnly(GRegion*, unsigned int max_nb_cliques, string filename=string());
-    virtual void execute_blossom(unsigned int max_nb_cliques, string filename=string());
-    virtual void execute_blossom(GRegion*, unsigned int max_nb_cliques, string filename=string());
-    virtual void createBlossomInfo();
-    void createBlossomInfo(GRegion *gr);
+protected:
+  bool debug,debug_graph;
+  std::map<Hex*, std::set<MElement*> > hex_to_tet;
+  std::map<Hex*, std::set<PELine*> > hex_to_edges;
+  std::map<PELine*, std::set<Hex*> > edges_to_hex;
+  std::map<Hex*, std::set<PETriangle*> > hex_to_faces;
+  std::map<PETriangle*, std::set<Hex*> > faces_to_hex;
+  std::map<PETriangle*, unsigned int > faces_connectivity;// # of adjacent tets (1 or 2)
+  std::map<MElement*, std::set<Hex*> >tet_to_hex;
+  std::map<Hex*, std::vector<double> > hex_ranks;
+
+
+  typedef unsigned long long hash_key;
+  //    typedef tr1::unordered_multimap<hash_key, Hex*> graph_data;
+  //    typedef tr1::unordered_multimap<hash_key, pair<Hex*, graph_data > > graph;
+  typedef multimap<hash_key, Hex*> graph_data;
+  typedef multimap<hash_key, pair<Hex*, graph_data > > graph;
+
+  graph incompatibility_graph;
+  set<Hex*> set_of_all_hex_in_graph;
+
+  std::multimap<unsigned long long, Hex*>created_potential_hex;
+
+  void create_faces_connectivity();
+  void add_face_connectivity(MElement *tet, int i, int j, int k);
+
+  void add_edges(Hex *hex);
+  void fill_edges_table(const MVertex *a, const MVertex *b, const MVertex *c, const MVertex *d, Hex *hex);
+  void add_face(const MVertex *a,const MVertex* b,const MVertex *c,Hex *hex);
+  void add_face(const MVertex *a,const MVertex* b,const MVertex *c,std::multimap<unsigned long long, pair<PETriangle*,int> > &f);
+  std::multimap<double,Hex*> degree;// degree = the final ranking of hexahedra
+  std::multimap<int,Hex*> idegree;// idegree = number of connected hex in indirect neighbors graph
+  std::multimap<int,Hex*> ndegree;// ndegree = number of direct neighbors !!! not chosen yet !!!
+  std::map<Hex*,int> reverse_idegree;
+  std::map<Hex*,int> reverse_ndegree;
+  // each tet has at least one neighbor, at most four. For all not chosen hex, check this data to find how many direct neighbors...
+  //    std::map<MElement*,set<PETriangle*> > tet_to_triangle;
+  std::map<PETriangle*,set<MElement*> > triangle_to_tet;
+  std::map<MElement*,int> tet_degree;
+
+  bool find_face_in_blossom_info(MVertex *a, MVertex *b, MVertex *c, MVertex *d);
+  void compute_hex_ranks_blossom();
+  PETriangle* get_triangle(MVertex*a, MVertex* b, MVertex *c);
+  bool is_blossom_pair(PETriangle *t1, PETriangle *t2);
+
+  tripair blossom_info;
+  trimap triangular_faces;
+  linemap edges_and_diagonals;
+
+  map<PETriangle*, GFace*> tri_to_gface_info;
+
+  vector<Hex*> chosen_hex;
+  vector<MElement*> chosen_tet;
+
+  citer find_the_triangle(PETriangle *t, const trimap &list);
+  linemap::const_iterator  find_the_line(PELine *t, const linemap &list);
+  std::multimap<unsigned long long, pair<PETriangle*,int> >::iterator  find_the_triangle(PETriangle *t, std::multimap<unsigned long long, pair<PETriangle*, int> > &list);
+  std::multimap<unsigned long long, Hex* >::const_iterator  find_the_created_potential_hex(Hex *t, const std::multimap<unsigned long long, Hex*>  &list);
+
+  int nbhex_in_losses_graph;
+  double average_connectivity;
+  bool post_check_validation(Hex* current_hex);
+
+  PETriangle* get_triangle(MElement *element, int i, int j, int k);
+
+  void compute_hex_ranks();
+
+  // check if the hex is good enough to be put into the graph. If not in the graph, it cannot be chosen...
+  bool is_not_good_enough(Hex* hex);
+
+  // fills incompatibility_graph if two hex share a common (non-sliver!) tet
+  void create_indirect_neighbors_graph();
+
+  graph::iterator find_hex_in_graph(Hex* hex);
+  graph_data::iterator find_hex_in_graphrow(Hex* hex, graph_data &row);
+  bool find_hex_couple_in_graph(Hex* hex, Hex* other_hex);
+  void add_graph_entry(Hex* hex, Hex* other_hex);
+
+  // fills incompatibility_graph if two hex are incompatible direct neighbors,
+  // i.e. they have one (or more) common face or common edge and are not compatible
+  void create_direct_neighbors_incompatibility_graph();
+  void evaluate_hex_couple(Hex* hex, Hex* other_hex);
+
+  // if two hex are not connected in the incompatibility_graph, they are compatible
+  void create_losses_graph(GRegion *gr);
+
+  void merge_clique(GRegion* gr, cliques_losses_graph<Hex*> &cl,int clique_number=0);
+
+  void fill_tet_to_hex_table(Hex *hex);
+
+  virtual void pattern1(GRegion*);
+  virtual void pattern2(GRegion*);
+  virtual void pattern3(GRegion*);
+
+  void merge(GRegion*);
+
+  // ------- exports --------
+  void export_tets(set<MElement*> &tetset, Hex* hex, string s);
+  void export_single_hex_all(Hex* hex,string s);
+  void export_single_hex(Hex* hex,string s);
+  void export_single_hex_faces(Hex* hex,string s);
+  void export_single_hex_tet(Hex* hex,string s);
+  void export_all_hex(int &file,GRegion *gr);
+  void export_hexmesh_so_far(int &file);
+  void export_direct_neighbor_table(int max);
+  void export_hex_init_degree(GRegion *gr, const std::map<Hex*,int> &init_degree, const vector<Hex*> &chosen_hex);
+
+  int max_nb_cliques;
+  string graphfilename;
+
+public:
+  Recombinator_Graph(unsigned int max_nb_cliques, string filename=string());
+  ~Recombinator_Graph();
+  virtual void execute();
+  virtual void execute(GRegion*);
+  virtual void buildGraphOnly(unsigned int max_nb_cliques, string filename=string());
+  virtual void buildGraphOnly(GRegion*, unsigned int max_nb_cliques, string filename=string());
+  virtual void execute_blossom(unsigned int max_nb_cliques, string filename=string());
+  virtual void execute_blossom(GRegion*, unsigned int max_nb_cliques, string filename=string());
+  virtual void createBlossomInfo();
+  void createBlossomInfo(GRegion *gr);
 };
 
 class Prism{
-  private:
-    double quality;
-    MVertex *a,*b,*c,*d,*e,*f;
-  public:
-    Prism();
-    Prism(MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*);
-    ~Prism();
-    double get_quality() const;
-    void set_quality(double);
-    MVertex* get_a();
-    MVertex* get_b();
-    MVertex* get_c();
-    MVertex* get_d();
-    MVertex* get_e();
-    MVertex* get_f();
-    void set_vertices(MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*);
-    bool operator<(const Prism&) const;
+private:
+  double quality;
+  MVertex *a,*b,*c,*d,*e,*f;
+public:
+  Prism();
+  Prism(MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*);
+  ~Prism();
+  double get_quality() const;
+  void set_quality(double);
+  MVertex* get_a();
+  MVertex* get_b();
+  MVertex* get_c();
+  MVertex* get_d();
+  MVertex* get_e();
+  MVertex* get_f();
+  void set_vertices(MVertex*,MVertex*,MVertex*,MVertex*,MVertex*,MVertex*);
+  bool operator<(const Prism&) const;
 };
 
 class Supplementary{
-  private:
-    std::vector<Prism> potential;
-    std::map<MElement*,bool> markings;
-    std::map<MVertex*,std::set<MVertex*> > vertex_to_vertices;
-    std::map<MVertex*,std::set<MElement*> > vertex_to_tetrahedra;
-    std::multiset<Facet> hash_tableA;
-    std::multiset<Diagonal> hash_tableB;
-    std::multiset<Diagonal> hash_tableC;
-    std::multiset<Tuple> tuples;
-    std::set<MElement*> triangles;
-  public:
-    Supplementary();
-    ~Supplementary();
-
-    void execute();
-    void execute(GRegion*);
-
-    void init_markings(GRegion*);
-    void pattern(GRegion*);
-    void merge(GRegion*);
-    void rearrange(GRegion*);
-    void statistics(GRegion*);
-    void build_tuples(GRegion*);
-    void modify_surfaces(GRegion*);
-    void modify_surfaces(MVertex*,MVertex*,MVertex*,MVertex*);
-
-    bool four(MElement*);
-    bool five(MElement*);
-    bool six(MElement*);
-    bool eight(MElement*);
-
-    bool sliver(MElement*,Prism);
-    bool valid(Prism,const std::set<MElement*>&);
-    bool valid(Prism);
-    double eta(MVertex*,MVertex*,MVertex*,MVertex*);
-    bool linked(MVertex*,MVertex*);
-
-    void find(MVertex*,MVertex*,const std::vector<MVertex*>&,std::set<MVertex*>&);
-    void find(MVertex*,Prism,std::set<MElement*>&);
-
-    void intersection(const std::set<MVertex*>&,const std::set<MVertex*>&,const std::vector<MVertex*>&,std::set<MVertex*>&);
-
-    bool inclusion(MVertex*,Prism);
-    bool inclusion(MVertex*,MVertex*,MVertex*,MVertex*,MVertex*);
-    bool inclusion(MVertex*,MVertex*,MVertex*,const std::set<MElement*>&);
-    bool inclusion(Facet);
-    bool inclusion(Diagonal);
-    bool duplicate(Diagonal);
-
-    bool conformityA(Prism);
-    bool conformityA(MVertex*,MVertex*,MVertex*,MVertex*);
-    bool conformityB(Prism);
-    bool conformityC(Prism);
-
-    bool faces_statuquo(Prism);
-    bool faces_statuquo(MVertex*,MVertex*,MVertex*,MVertex*);	
-
-    void build_vertex_to_vertices(GRegion*);
-    void build_vertex_to_tetrahedra(GRegion*);
-    void build_hash_tableA(Prism);
-    void build_hash_tableA(MVertex*,MVertex*,MVertex*,MVertex*);
-    void build_hash_tableA(Facet);
-    void build_hash_tableB(Prism);
-    void build_hash_tableB(MVertex*,MVertex*,MVertex*,MVertex*);
-    void build_hash_tableB(Diagonal);
-    void build_hash_tableC(Prism);
-    void build_hash_tableC(Diagonal);
-
-    double scaled_jacobian(MVertex*,MVertex*,MVertex*,MVertex*);
-    double min_scaled_jacobian(Prism);
+private:
+  std::vector<Prism> potential;
+  std::map<MElement*,bool> markings;
+  std::map<MVertex*,std::set<MVertex*> > vertex_to_vertices;
+  std::map<MVertex*,std::set<MElement*> > vertex_to_tetrahedra;
+  std::multiset<Facet> hash_tableA;
+  std::multiset<Diagonal> hash_tableB;
+  std::multiset<Diagonal> hash_tableC;
+  std::multiset<Tuple> tuples;
+  std::set<MElement*> triangles;
+public:
+  Supplementary();
+  ~Supplementary();
+
+  void execute();
+  void execute(GRegion*);
+
+  void init_markings(GRegion*);
+  void pattern(GRegion*);
+  void merge(GRegion*);
+  void rearrange(GRegion*);
+  void statistics(GRegion*);
+  void build_tuples(GRegion*);
+  void modify_surfaces(GRegion*);
+  void modify_surfaces(MVertex*,MVertex*,MVertex*,MVertex*);
+
+  bool four(MElement*);
+  bool five(MElement*);
+  bool six(MElement*);
+  bool eight(MElement*);
+
+  bool sliver(MElement*,Prism);
+  bool valid(Prism,const std::set<MElement*>&);
+  bool valid(Prism);
+  double eta(MVertex*,MVertex*,MVertex*,MVertex*);
+  bool linked(MVertex*,MVertex*);
+
+  void find(MVertex*,MVertex*,const std::vector<MVertex*>&,std::set<MVertex*>&);
+  void find(MVertex*,Prism,std::set<MElement*>&);
+
+  void intersection(const std::set<MVertex*>&,const std::set<MVertex*>&,const std::vector<MVertex*>&,std::set<MVertex*>&);
+
+  bool inclusion(MVertex*,Prism);
+  bool inclusion(MVertex*,MVertex*,MVertex*,MVertex*,MVertex*);
+  bool inclusion(MVertex*,MVertex*,MVertex*,const std::set<MElement*>&);
+  bool inclusion(Facet);
+  bool inclusion(Diagonal);
+  bool duplicate(Diagonal);
+
+  bool conformityA(Prism);
+  bool conformityA(MVertex*,MVertex*,MVertex*,MVertex*);
+  bool conformityB(Prism);
+  bool conformityC(Prism);
+
+  bool faces_statuquo(Prism);
+  bool faces_statuquo(MVertex*,MVertex*,MVertex*,MVertex*);
+
+  void build_vertex_to_vertices(GRegion*);
+  void build_vertex_to_tetrahedra(GRegion*);
+  void build_hash_tableA(Prism);
+  void build_hash_tableA(MVertex*,MVertex*,MVertex*,MVertex*);
+  void build_hash_tableA(Facet);
+  void build_hash_tableB(Prism);
+  void build_hash_tableB(MVertex*,MVertex*,MVertex*,MVertex*);
+  void build_hash_tableB(Diagonal);
+  void build_hash_tableC(Prism);
+  void build_hash_tableC(Diagonal);
+
+  double scaled_jacobian(MVertex*,MVertex*,MVertex*,MVertex*);
+  double min_scaled_jacobian(Prism);
 };
 
 class PostOp{
- private:
+private:
   int nbr,nbr8,nbr6,nbr5,nbr4;
   double vol,vol8,vol6,vol5,vol4;
   int estimate1;
@@ -620,54 +620,54 @@ class PostOp{
   std::map<MVertex*,std::set<MElement*> > vertex_to_pyramids;
   std::multiset<Tuple> tuples;
   std::set<MElement*> triangles;
- public:
+public:
   PostOp();
   ~PostOp();
 
-    void execute(bool);
-    void execute(GRegion*,bool);
-
-    inline int get_nb_hexahedra()const{return nbr8;};
-    inline double get_vol_hexahedra()const{return vol8;};
-    inline int get_nb_elements()const{return nbr;};
-    inline double get_vol_elements()const{return vol;};
-
-    void init_markings(GRegion*);
-    void pyramids1(GRegion*);
-    void pyramids2(GRegion*);
-    void pyramids1(MVertex*,MVertex*,MVertex*,MVertex*,GRegion*);
-    void pyramids2(MVertex*,MVertex*,MVertex*,MVertex*,GRegion*);
-    void rearrange(GRegion*);
-    void statistics(GRegion*);
-    void build_tuples(GRegion*);
-    void modify_surfaces(GRegion*);
-    void modify_surfaces(MVertex*,MVertex*,MVertex*,MVertex*);
-
-    bool four(MElement*);
-    bool five(MElement*);
-    bool six(MElement*);
-    bool eight(MElement*);
-
-    bool equal(MVertex*,MVertex*,MVertex*,MVertex*);
-    bool different(MVertex*,MVertex*,MVertex*,MVertex*);
-    MVertex* other(MElement*,MVertex*,MVertex*);
-    MVertex* other(MElement*,MVertex*,MVertex*,MVertex*);
-    void mean(const std::set<MVertex*>&,MVertex*,const std::vector<MElement*>&);
-    double workaround(MElement*);
-
-    MVertex* find(MVertex*,MVertex*,MVertex*,MVertex*,MElement*);
-    void find_tetrahedra(MVertex*,MVertex*,std::set<MElement*>&);
-    void find_pyramids(MVertex*,MVertex*,std::set<MElement*>&);
-
-    void intersection(const std::set<MElement*>&,const std::set<MElement*>&,std::set<MElement*>&);
-
-    void build_vertex_to_tetrahedra(GRegion*);
-    void build_vertex_to_tetrahedra(MElement*);
-    void erase_vertex_to_tetrahedra(MElement*);
-
-    void build_vertex_to_pyramids(GRegion*);
-    void build_vertex_to_pyramids(MElement*);
-    void erase_vertex_to_pyramids(MElement*);
+  void execute(bool);
+  void execute(GRegion*,bool);
+
+  inline int get_nb_hexahedra()const{return nbr8;};
+  inline double get_vol_hexahedra()const{return vol8;};
+  inline int get_nb_elements()const{return nbr;};
+  inline double get_vol_elements()const{return vol;};
+
+  void init_markings(GRegion*);
+  void pyramids1(GRegion*);
+  void pyramids2(GRegion*);
+  void pyramids1(MVertex*,MVertex*,MVertex*,MVertex*,GRegion*);
+  void pyramids2(MVertex*,MVertex*,MVertex*,MVertex*,GRegion*);
+  void rearrange(GRegion*);
+  void statistics(GRegion*);
+  void build_tuples(GRegion*);
+  void modify_surfaces(GRegion*);
+  void modify_surfaces(MVertex*,MVertex*,MVertex*,MVertex*);
+
+  bool four(MElement*);
+  bool five(MElement*);
+  bool six(MElement*);
+  bool eight(MElement*);
+
+  bool equal(MVertex*,MVertex*,MVertex*,MVertex*);
+  bool different(MVertex*,MVertex*,MVertex*,MVertex*);
+  MVertex* other(MElement*,MVertex*,MVertex*);
+  MVertex* other(MElement*,MVertex*,MVertex*,MVertex*);
+  void mean(const std::set<MVertex*>&,MVertex*,const std::vector<MElement*>&);
+  double workaround(MElement*);
+
+  MVertex* find(MVertex*,MVertex*,MVertex*,MVertex*,MElement*);
+  void find_tetrahedra(MVertex*,MVertex*,std::set<MElement*>&);
+  void find_pyramids(MVertex*,MVertex*,std::set<MElement*>&);
+
+  void intersection(const std::set<MElement*>&,const std::set<MElement*>&,std::set<MElement*>&);
+
+  void build_vertex_to_tetrahedra(GRegion*);
+  void build_vertex_to_tetrahedra(MElement*);
+  void erase_vertex_to_tetrahedra(MElement*);
+
+  void build_vertex_to_pyramids(GRegion*);
+  void build_vertex_to_pyramids(MElement*);
+  void erase_vertex_to_pyramids(MElement*);
 };
 
 #endif