- put back necessary ifnedf in GModel.h for linking on Windows with python bindings

 - new options to create all possible tets for a Gregion with createTetrahedralMesh
 - changes in yamaka.h and cpp - put const and inline functions

Geo/GModel.h

@@ -164,10 +164,11 @@ class GModel
  public:
   GModel(std::string name="");
   virtual ~GModel();
-
-
+// Required for python bindings on Windows
+#ifndef SWIG
   // the static list of all loaded models
   static std::vector<GModel*> list;
+# endif
 
   // return the current model, and sets the current model index if
   // index >= 0

Mesh/Generator.h

@@ -19,7 +19,7 @@ void SmoothMesh(GModel *m);
 void RefineMesh(GModel *m, bool linear, bool splitIntoQuads=false,
                 bool splitIntoHexas=false);
 void RecombineMesh(GModel *m);
-GRegion * createTetrahedralMesh ( GModel *gm, fullMatrix<double> & pts, fullMatrix<int> &triangles ) ;
+GRegion * createTetrahedralMesh ( GModel *gm, fullMatrix<double> & pts, fullMatrix<int> &triangles, bool all_tets=false ) ;
   //GRegion * createTetrahedralMesh ( GModel *gm, unsigned int nbPts , double *pts, unsigned int nbTriangles, int *triangles );
 
 #endif

Mesh/meshDiscreteRegion.cpp

@@ -32,11 +32,49 @@
 #include "MPrism.h"
 #include "MHexahedron.h"
 
+// Recursive function to generate all combinations of 4 indices between start 
+// and end indices included 
+// Jeanne - HEXTREME
+void combination_of_4( int combination[4], int start, int end, int index,
+  const std::vector<MVertex*>& vertices, std::vector<MTetrahedron*>& tets)
+{
+  if (index == 4 ) {
+    // Create the tet and get out
+    MVertex* v1 = vertices[combination[0]];
+    MVertex* v2 = vertices[combination[1]];
+    MVertex* v3 = vertices[combination[2]];
+    MVertex* v4 = vertices[combination[3]];
+    MTetrahedron* tet = new MTetrahedron(v1, v2, v3, v4);
+    tets.push_back(tet);
+    return;
+  }
+  for (int i = start; i <= end+1 -(4-index); i++) {
+    combination[index] = i;
+    combination_of_4(combination, i+1, end, index+1, vertices, tets);
+  }
+}
+
+// Fill a region with all possible tets genereated from the
+// combination of points assigned to it 
+// Jeanne - HEXTREME
+void create_all_possible_tets(GRegion* region, const std::vector<MVertex*>& vertices) {
+  unsigned int nb_points = vertices.size();
+
+  int combinaison[4];
+  std::vector<MTetrahedron*> tets;
+  combination_of_4(combinaison, 0, nb_points - 1, 0, vertices, tets);
+  std::cout << " Number of tets created - all possible combinations - " << tets.size() << std::endl;
+
+  for (int i = 0; i < tets.size(); ++i) {
+    region->addTetrahedron(tets[i]);
+  }
+}
+
 // triangles are defining the boundary
 // internal points are allowed
 // This has been done for HEXTREME
 GRegion * createDiscreteRegionFromRawData(GModel *gm, fullMatrix<double> &pts,
-                                          fullMatrix<int> &triangles)
+                                          fullMatrix<int> &triangles, bool all_tets)
 {
   GRegion *gr = new discreteRegion(gm, NEWREG());
   GFace *gf = new discreteFace(gm, NEWREG());
@@ -61,8 +99,8 @@ GRegion * createDiscreteRegionFromRawData(GModel *gm, fullMatrix<double> &pts,
       vs[i] = v;
     }
     else {
-      MVertex *v = new MFaceVertex(pts(i,0), pts(i,1), pts(i,2), gr, 0, 0);
-      gr->mesh_vertices.push_back(v);
+      MVertex *v = new MFaceVertex(pts(i,0), pts(i,1), pts(i,2), gf, 0, 0);
+      gf->mesh_vertices.push_back(v);
       vs[i] = v;
     }
   }
@@ -75,24 +113,32 @@ GRegion * createDiscreteRegionFromRawData(GModel *gm, fullMatrix<double> &pts,
     gf->triangles.push_back(t);
   }
 
-  delaunayTriangulation(1, 1, vs, gr->tetrahedra);
+//  delaunayTriangulation(1, 1, vs, gr->tetrahedra); // tetrahedralization is done when recovering the boundary
+
+  // create all tets
+  if (all_tets) {
+    create_all_possible_tets(gr, vs);
+  }
 
   return gr;
 }
 
 GRegion *createTetrahedralMesh(GModel *gm, fullMatrix<double> &pts,
-                               fullMatrix<int> &triangles)
+                               fullMatrix<int> &triangles, bool all_tets)
 {
-  GRegion *gr = createDiscreteRegionFromRawData(gm, pts, triangles);
-  try{
-    meshGRegionBoundaryRecovery(gr);
-  }
-  catch(int err){
-    if(err == 3){
-      Msg::Warning("Self-intersecting surface mesh: TODO!");
+  GRegion *gr = createDiscreteRegionFromRawData(gm, pts, triangles, all_tets);
+
+  if (!all_tets){
+    try {
+      meshGRegionBoundaryRecovery(gr);
     }
-    else{
-      Msg::Error("Could not recover boundary: error %d", err);
+    catch (int err) {
+      if (err == 3) {
+        Msg::Warning("Self-intersecting surface mesh: TODO!");
+      }
+      else {
+        Msg::Error("Could not recover boundary: error %d", err);
+      }
     }
   }
   return gr;

Mesh/meshGRegionBoundaryRecovery.cpp

@@ -135,6 +135,8 @@ bool tetgenmesh::reconstructmesh(void *p)
         all.insert(gv->points[i]->getVertex(0));
       }
     }
+    all.insert(_gr->mesh_vertices.begin(), _gr->mesh_vertices.end());
+
     _vertices.insert(_vertices.begin(), all.begin(), all.end());
   }
 

Mesh/yamakawa.h

@@ -170,68 +170,164 @@ protected:
 
 
 
-class Hex{
+class Hex {
 private:
   double quality;
   unsigned long long hash;
-  MVertex *a,*b,*c,*d,*e,*f,*g,*h;
-  void set_hash();
+  MVertex* vertices_[8];
+  //MVertex *a,*b,*c,*d,*e,*f,*g,*h;
+private:
+  void set_hash(){
+    hash = 0.;
+    for (int i = 0; i < 8; ++i) {
+      hash += vertices_[i]->getNum();
+    }    
+  }
+
 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&) ;
+  Hex() : quality(0.), hash(0.), vertices_{NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL} {};
+
+  Hex(MVertex* a2, MVertex* b2, MVertex* c2, MVertex* d2, MVertex* e2, MVertex* f2, MVertex* g2, MVertex* h2) :
+    quality(0.), vertices_{a2,b2,c2, d2, e2, f2, g2, h2}
+  {
+    set_hash();
+  }
+  ~Hex() {};
+  double get_quality() const { return quality; }
+  void set_quality(double new_quality) { quality = new_quality; }
+  MVertex* get_a() const { return vertices_[0]; }
+  MVertex* get_b() const { return vertices_[1]; }
+  MVertex* get_c() const { return vertices_[2]; }
+  MVertex* get_d() const { return vertices_[3]; }
+  MVertex* get_e() const { return vertices_[4]; }
+  MVertex* get_f() const { return vertices_[5]; }
+  MVertex* get_g() const { return vertices_[6]; }
+  MVertex* get_h() const { return vertices_[7]; }
+  MVertex* getVertex(unsigned int n) const {
+    if (n < 8) {
+      return vertices_[n];
+    }
+    else {
+      cout << "Hex: unknown vertex number " << n << endl;
+      throw;
+      return NULL;
+    }
+  }
+
+  bool hasVertex(const MVertex *v) const {
+    for (int i = 0; i < 8; i++) {
+      if (getVertex(i) == v) {
+        return true;
+      }
+    }
+    return false;
+  }
+
+  bool same_vertices(Hex *h) const {
+    for (int i = 0; i < 8; i++) {
+      if (!(h->hasVertex(getVertex(i)))) {
+        return false;
+      }
+    }
+    return true;
+  }
+
+  void set_vertices(MVertex* a2, MVertex* b2, MVertex* c2, MVertex* d2, MVertex* e2, MVertex* f2, MVertex* g2, MVertex* h2) {
+    vertices_[0] = a2;
+    vertices_[1] = b2;
+    vertices_[2] = c2;
+    vertices_[3] = d2;
+    vertices_[4] = e2;
+    vertices_[5] = f2;
+    vertices_[6] = g2;
+    vertices_[7] = h2;
+  }
+
+  unsigned long long get_hash() {
+    if (hash == 0. && vertices_[0]!=NULL) {
+      set_hash();
+    }
+    return hash;
+  }
+  bool operator<(const Hex& hex) const {
+    return quality > hex.get_quality(); // Why > ??? Shouldn't it be < ?? Jeanne.
+  }
+
 };
 
-class Facet{
+class Facet {
 private:
-  MVertex *a,*b,*c;
+  MVertex *a, *b, *c;
   int num[3];
   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;
+  Facet() : a(NULL), b(NULL), c(NULL), num{ -1, -1, -1 }, hash(0.){}
+  Facet(MVertex* a2, MVertex* b2, MVertex* c2) :
+    a(a2), b(b2), c(c2), num {-1,-1,-1}, hash(0.){
+    compute_hash();
+  }
+  ~Facet() {};
+  MVertex* get_a() const { return a; }
+  MVertex* get_b() const { return b; }
+  MVertex* get_c() const { return c; }
+  void set_vertices(MVertex*a2, MVertex*b2, MVertex*c2) {
+    a = a2;
+    b = b2;
+    c = c2;
+    compute_hash();
+  }
+  bool same_vertices(const Facet& facet) const {    
+    bool c1 = (a == facet.get_a()) || (a == facet.get_b()) || (a == facet.get_c());
+    bool c2 = (b == facet.get_a()) || (b == facet.get_b()) || (b == facet.get_c());
+    bool c3 = (c == facet.get_a()) || (c == facet.get_b()) || (c == facet.get_c());
+    return c1 && c2 && c3;
+  }
+  void compute_hash() {
+    num[0] = a->getNum();
+    num[1] = b->getNum();
+    num[2] = c->getNum();
+    std::sort(num, num + 3);
+    hash = num[2] + 1e4*num[1] + 1e8*num[0];
+  }
+  unsigned long long get_hash() const { return hash; }
+  bool operator<(const Facet& rhs) const {
+    return hash<rhs.get_hash();
+  }
 };
 
 class Diagonal{
 private:
   MVertex *a,*b;
   unsigned long long hash;
+private:
+  void compute_hash() {
+    hash = a->getNum() + b->getNum();
+  }
+
 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;
+  Diagonal() :a(NULL), b(NULL), hash() {};
+  Diagonal(MVertex*a2, MVertex*b2) :a(a2), b(b2){
+    compute_hash();
+  }
+  ~Diagonal() {};
+  MVertex* get_a() const {return a;}
+  MVertex* get_b() const {return b;}
+  void set_vertices(MVertex*a2, MVertex*b2) {
+    a = a2;
+    b = b2;
+    compute_hash();
+  }
+  bool same_vertices(Diagonal diagonal) const {
+    bool c1 = (a == diagonal.get_a()) || (a == diagonal.get_b());
+    bool c2 = (b == diagonal.get_a()) || (b == diagonal.get_b());
+    return c1 && c2;
+  }
+  unsigned long long get_hash() const {
+    return hash;
+  }
+  bool operator<(const Diagonal& rhs) const {
+    return hash<rhs.get_hash();
+  }
 };
 
 class Tuple{
@@ -241,20 +337,44 @@ private:
   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;
+  Tuple() : v1(NULL), v2(NULL), v3(NULL), element(NULL), gf(NULL), hash(0.) {}
+  Tuple(MVertex* a, MVertex* b, MVertex* c, MElement* element2, GFace* gf2)
+    : Tuple(a,b,c)
+  {    
+    element = element2;
+    gf = gf2;  
+  }
+  Tuple(MVertex* a, MVertex* b, MVertex* c)
+    :element(NULL), gf(NULL) 
+  {
+    MVertex* tmp[3] = { a,b,c };
+    std::sort(tmp, tmp + 3);
+    v1 = tmp[0];
+    v2 = tmp[1];
+    v2 = tmp[2];
+    hash = a->getNum() + b->getNum() + c->getNum();
+  }
+  ~Tuple() {};
+  MVertex* get_v1() const {return v1;}
+  MVertex* get_v2() const {return v2;}
+  MVertex* get_v3() const {return v3;}
+  MElement* get_element() const { return element; }
+  GFace* get_gf() const { return gf; }
+  bool same_vertices(const Tuple& rhs) const {
+    if (v1 == rhs.get_v1() && v2 == rhs.get_v2() && v3 == rhs.get_v3()) {
+      return true;
+    } else {
+      return false;
+    }
+  }
+  unsigned long long get_hash() const { return hash; }
+  bool operator<(const Tuple& rhs) const {
+    return hash< rhs.get_hash();
+  }
 };
 
+
+
 //inline std::ostream& operator<<(std::ostream& s, const PETriangle& t){
 //  const MVertex *v;
 //  for (int i=0;i<3;i++){
@@ -310,26 +430,26 @@ public:
   // 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 valid(Hex&) const;
+  double eta(MVertex*,MVertex*,MVertex*,MVertex*) const;
   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 find(MVertex*,MVertex*,const std::vector<MVertex*>&,std::set<MVertex*>&) const;
+  void find(MVertex*,MVertex*,MVertex*,const std::vector<MVertex*>&,std::set<MVertex*>&) const;
+  void find(MVertex*,MVertex*,std::set<MElement*>&) const;
+  void find(MVertex*,Hex,std::set<MElement*>&) const;
+  MVertex* find(MVertex*,MVertex*,MVertex*,MVertex*,const std::set<MElement*>&) const;
 
-  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*>&);
+  void intersection(const std::set<MVertex*>&,const std::set<MVertex*>&,const std::vector<MVertex*>&,std::set<MVertex*>&) const;
+  void intersection(const std::set<MVertex*>&,const std::set<MVertex*>&,const std::set<MVertex*>&,const std::vector<MVertex*>&,std::set<MVertex*>&) const;
+  void intersection(const std::set<MElement*>&,const std::set<MElement*>&,std::set<MElement*>&) const;
 
   // return true if vertex belong to hex
-  bool inclusion(MVertex*,Hex);
+  bool inclusion(MVertex*,Hex) const;
   // renvoie true si vertex se trouve dans [a,b,c]
-  bool inclusion(MVertex*,MVertex*,MVertex*,MVertex*,MVertex*);
+  bool inclusion(MVertex*,MVertex*,MVertex*,MVertex*,MVertex*) const ;
   // return true if all three vertices v1,v2 and v3 belong to one tet
-  bool inclusion(MVertex*,MVertex*,MVertex*,const std::set<MElement*>&);
+  bool inclusion(MVertex*,MVertex*,MVertex*,const std::set<MElement*>&) const;
   // return true si la facet existe dans la table A
   bool inclusion(Facet);
   // return true si la diagonal existe dans la table B