x

Mesh/BackgroundMesh.cpp

@@ -484,6 +484,9 @@ void backgroundMesh::propagateCrossField(GFace *_gf, simpleFunction<double> *ONE
   propagateValuesOnFace(_gf,_cosines4,ONE,false);
   propagateValuesOnFace(_gf,_sines4,ONE,false);
 
+  //    print("cos4.pos",0,_cosines4,0);
+  //    print("sin4.pos",0,_sines4,0);
+  
   std::map<MVertex*,MVertex*>::iterator itv2 = _2Dto3D.begin();
   for ( ; itv2 != _2Dto3D.end(); ++itv2){
     MVertex *v_2D = itv2->first;

Mesh/delaunay3d.cpp

@@ -687,6 +687,7 @@ static bool fixDelaunayCavity (Vert *v,
   starShapeness (v, bndK, _negatives);
 
   if (_negatives.empty())return false;
+  
   //  return true;
   // unset all tets of the cavity
   for (unsigned int i=0; i< cavity.size(); i++)cavity[i]->unset(myThread,K);

Mesh/delaunay3d_private.h

@@ -146,16 +146,22 @@ public :
 struct edgeContainer
 {
   std::vector<std::vector<Edge> > _hash;
-  size_t _size;
+  size_t _size, _size_obj;
   edgeContainer(unsigned int N = 1000000)
   {
     _size = 0;
     _hash.resize(N);
+    _size_obj = sizeof(Edge); 
   }
 
+  inline size_t H (const Edge &e) const {
+    const size_t h = ((size_t)e.first) ;
+    //    printf("%lu %lu %lu %lu\n",h,(h/2)%_hash.size(),h/64,h>>6);
+    return (h/_size_obj) %_hash.size();  
+  }
+  
   inline bool find (const Edge &e) const {
-    size_t h = ((size_t) e.first >> 8) ;
-    const std::vector<Edge> &v = _hash[h %_hash.size()];
+    const std::vector<Edge> &v = _hash[H(e)];
     for (unsigned int i=0; i< v.size();i++)if (e == v[i]) {return true;}
     return false;
   }
@@ -164,8 +170,7 @@ struct edgeContainer
   
   bool addNewEdge (const Edge &e)
   {
-    size_t h = ((size_t) e.first >> 8) ;
-    std::vector<Edge> &v = _hash[h %_hash.size()];
+    std::vector<Edge> &v = _hash[H(e)];
     for (unsigned int i=0; i< v.size();i++)if (e == v[i]) {return false;}
     v.push_back(e);
     _size++;

Mesh/delaunay_refinement.cpp

@@ -142,7 +142,7 @@ void saturateEdge(Edge &e, std::vector<Vert*> &S, std::stack<IPT> &temp)
     bgm = fields->get(fields->getBackgroundField());
   }  
 
-  if (1){
+  if (0){
     const double length = p1.distance(p2);
     const double lc_avg = (e.first->lc() + e.second->lc())*.5;
     const double l_adim_approx =   length / lc_avg;
@@ -221,7 +221,7 @@ void saturateEdges(edgeContainer &ec,
   const int N = T.size(0);
   for (int i=0;i<N;i++){
     Tet *t = T(0,i);
-   if (t->V[0] && t->_modified){
+    if (t->V[0] && t->_modified){
       t->_modified = false;
       for (int iEdge=0;iEdge<6;iEdge++){
 	Edge ed = t->getEdge(iEdge);

Mesh/meshGFaceTransfinite.cpp

@@ -422,7 +422,7 @@ int MeshTransfiniteSurface(GFace *gf)
         MVertex *v4 = tab[i][j + 1];
         if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine)
           gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
-        else if(gf->meshAttributes.transfiniteArrangement == 1 ||
+	else if(gf->meshAttributes.transfiniteArrangement == 1 ||
                 (gf->meshAttributes.transfiniteArrangement == 2 &&
                  ((i % 2 == 0 && j % 2 == 1) ||
                   (i % 2 == 1 && j % 2 == 0))) ||
@@ -430,9 +430,10 @@ int MeshTransfiniteSurface(GFace *gf)
                  ((i % 2 == 0 && j % 2 == 0) ||
                   (i % 2 == 1 && j % 2 == 1)))
 		){
-          gf->triangles.push_back(new MTriangle(v1, v2, v3));
-          gf->triangles.push_back(new MTriangle(v3, v4, v1));
-        }
+	  //        else if(rand() % 2 == 0){
+	  gf->triangles.push_back(new MTriangle(v1, v2, v3));
+	  gf->triangles.push_back(new MTriangle(v3, v4, v1));
+	}
         else{
           gf->triangles.push_back(new MTriangle(v1, v2, v4));
           gf->triangles.push_back(new MTriangle(v4, v2, v3));

Mesh/meshGRegionBoundaryRecovery.cpp

@@ -227,7 +227,9 @@ bool tetgenmesh::reconstructmesh(void *p)
       ver2tetarray[i] = NULL;
     }
 
+    
     // Create the tetrahedra and connect those that share a common face.
+    printf("DOING for %d tets\n",tets.size());
     for (unsigned int i = 0; i < tets.size(); i++) {
       // Get the four vertices.
       for (int j = 0; j < 4; j++) {
@@ -756,6 +758,7 @@ bool tetgenmesh::reconstructmesh(void *p)
           }
         }
         assert(ge != NULL);
+	Msg::Info("Steiner points exist on GEdge %d", ge->tag());
         // Delete the old triangles.
         for(unsigned int i = 0; i < ge->lines.size(); i++)
           delete ge->lines[i];

Mesh/meshGRegionDelaunayInsertion.cpp

@@ -23,6 +23,42 @@
 
 int MTet4::radiusNorm = 2;
 
+struct edgeContainerB
+{
+  std::vector<std::vector<MEdge> > _hash;
+  size_t _size, _size_obj;
+  edgeContainerB(unsigned int N = 1000000)
+  {
+    _size = 0;
+    _hash.resize(N);
+    _size_obj = sizeof(MEdge); 
+  }
+
+  inline size_t H (const MEdge &e) const {
+    const size_t h = ((size_t)e.getSortedVertex(0)) ;
+    return (h/_size_obj) %_hash.size();  
+  }
+  
+  inline bool find (const MEdge &e) const {
+    const std::vector<MEdge> &v = _hash[H(e)];
+    for (unsigned int i=0; i< v.size();i++)if (e == v[i]) {return true;}
+    return false;
+  }
+
+  bool empty () const {return _size == 0;}
+  
+  bool addNewEdge (const MEdge &e)
+  {
+    std::vector<MEdge> &v = _hash[H(e)];
+    for (unsigned int i=0; i< v.size();i++)if (e == v[i]) {return false;}
+    v.push_back(e);
+    _size++;
+    return true;
+  }
+};
+
+
+
 static void createAllEmbeddedEdges (GRegion *gr, std::set<MEdge, Less_Edge> &allEmbeddedEdges)
 {
   std::list<GEdge*> e = gr->embeddedEdges();
@@ -34,6 +70,16 @@ static void createAllEmbeddedEdges (GRegion *gr, std::set<MEdge, Less_Edge> &all
   }
 }
 
+static void createAllEmbeddedEdges (GRegion *gr, edgeContainerB &embedded){
+  std::list<GEdge*> e = gr->embeddedEdges();
+  for (std::list<GEdge*>::iterator it = e.begin() ; it != e.end(); ++it){
+    for (unsigned int i = 0; i < (*it)->lines.size(); i++){
+      embedded.addNewEdge(MEdge((*it)->lines[i]->getVertex(0),(*it)->lines[i]->getVertex(1)));
+    }
+  }
+}
+
+
 static void createAllEmbeddedFaces (GRegion *gr, std::set<MFace, Less_Face> &allEmbeddedFaces)
 {
   std::list<GFace*> f = gr->embeddedFaces();
@@ -77,11 +123,44 @@ struct faceXtet{
     MVertex *v1 = t1->tet()->getVertex(faces[iFac][1]);
     MVertex *v2 = t1->tet()->getVertex(faces[iFac][2]);
 
-    unsorted[0] = v[0] = v0;
-    unsorted[1] = v[1] = v1;
-    unsorted[2] = v[2] = v2;
-    MVertexLessThanNum lll;
-    std::sort(v, v + 3, lll);
+    unsorted[0] = v0;
+    unsorted[1] = v1;
+    unsorted[2] = v2;
+
+    if (v0->getNum() < v1->getNum() && v0->getNum() <  v2->getNum()) {
+      v[0] = v0;
+      if (v1->getNum() < v2->getNum()){
+	v[1] = v1;
+	v[2] = v2;
+      }
+      else {
+	v[1] = v2;
+	v[2] = v1;
+      }
+    }
+    else if (v1->getNum() < v0->getNum() && v1->getNum() <  v2->getNum()) {
+      v[0] = v1;
+      if (v0->getNum() < v2->getNum()){
+	v[1] = v0;
+	v[2] = v2;
+      }      
+      else {
+	v[1] = v2;
+	v[2] = v0;
+      }
+    }
+    else {
+      v[0] = v2;
+      if (v0->getNum() < v1->getNum()){
+	v[1] = v0;
+	v[2] = v1;
+      }      
+      else {
+	v[1] = v1;
+	v[2] = v0;
+      }
+    }
+    //  std::sort(v, v + 3, lll);
   }
 
   inline MVertex * getVertex (int i) const { return unsorted[i];}
@@ -117,7 +196,9 @@ struct faceXtet{
 void connectTets_vector2(std::vector<MTet4*> &t, std::vector<faceXtet> &conn)
 {
   conn.clear();
+  conn.reserve(4*t.size());
   //  unsigned int k = 0;
+  //  printf("COUCOU1 %d tets\n",t.size());
   for (unsigned int i=0;i<t.size();i++){
     if (!t[i]->isDeleted()){
       for (int j = 0; j < 4; j++){
@@ -126,8 +207,10 @@ void connectTets_vector2(std::vector<MTet4*> &t, std::vector<faceXtet> &conn)
     }
   }
   if (!conn.size())return;
+  //  printf("COUCOU2 %d faces\n",conn.size());
   std::sort(conn.begin(), conn.end());
 
+  //  printf("COUCOU3\n");
   for (unsigned int i=0;i<conn.size()-1;i++){
     faceXtet &f1  = conn[i];
     faceXtet &f2  = conn[i+1];
@@ -137,6 +220,7 @@ void connectTets_vector2(std::vector<MTet4*> &t, std::vector<faceXtet> &conn)
       ++i;
     }
   }
+  //  printf("COUCOU4\n");  
 }
 
 
@@ -349,16 +433,10 @@ bool insertVertexB(std::list<faceXtet> &shell,
     double d3 = distance (it->v[2],v);
     
     double lc = Extend2dMeshIn3dVolumes() ? std::min(lc1, lc2) : lc2;
-    if (d1 < lc * .05 || d2 < lc * .05 || d3 < lc * .05) {
-      //      printf("coucou %g %g %g %g\n",lc,d1,d2,d3);
-      onePointIsTooClose = true;
-    }
+    if (d1 < lc * .05 || d2 < lc * .05 || d3 < lc * .05) onePointIsTooClose = true;
 
     newTets[k++] = t4;
-    // all new tets are pushed front in order to ba able to destroy
-    // them if the cavity is not star shaped around the new vertex.
-    // here, we better use robust perdicates that implies to orient
-    // all faces and ensure that the tets are all oriented the same.
+
     new_cavity.push_back(t4);
     MTet4 *otherSide = it->t1->getNeigh(it->i1);
 
@@ -381,6 +459,7 @@ bool insertVertexB(std::list<faceXtet> &shell,
       (*ittet)->setDeleted(false);
       ++ittet;
     }
+    delete [] newTets;
     return false;
   }
 }
@@ -532,6 +611,7 @@ void adaptMeshGRegion::operator () (GRegion *gr)
   createAllEmbeddedFaces (gr, allEmbeddedFaces);
   std::set<MEdge, Less_Edge> allEmbeddedEdges;
   createAllEmbeddedEdges (gr, allEmbeddedEdges);
+
   connectTets(allTets.begin(),allTets.end(),&allEmbeddedFaces);
 
   double t1 = Cpu();
@@ -715,7 +795,7 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
   if (!gr->tetrahedra.size())return;
 
 
-  typedef std::list<MTet4 *> CONTAINER ;
+  typedef std::vector<MTet4 *> CONTAINER ;
   CONTAINER allTets;
   for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
     MTet4 * t = new MTet4(gr->tetrahedra[i], qm);
@@ -726,10 +806,18 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
 
   std::set<MFace, Less_Face> allEmbeddedFaces;
   createAllEmbeddedFaces (gr, allEmbeddedFaces);
-  connectTets(allTets.begin(),allTets.end(),&allEmbeddedFaces);
+
+
   std::set<MEdge, Less_Edge> allEmbeddedEdges;
   createAllEmbeddedEdges (gr, allEmbeddedEdges);
 
+  // daaaaaaamn slow !!!
+  //  connectTets(allTets.begin(),allTets.end(),allEmbeddedFaces.empty() ? NULL : &allEmbeddedFaces);
+  {
+    std::vector<faceXtet> conn;
+    connectTets_vector2(allTets, conn);
+  }
+  
   double t1 = Cpu();
   std::vector<MTet4*> illegals;
   const int nbRanges = 10;
@@ -755,7 +843,7 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
         }
       }
     }
-    Msg::Info("Opti : START with %12.5E QBAD %12.5E QAVG %12.5E",
+    Msg::Info("Opti : STARTS with %12.5E QBAD %12.5E QAVG %12.5E",
               totalVolumeb, worst, avg / count);
     for (int i = 0; i < nbRanges; i++){
       double low  = (double)i / nbRanges;
@@ -765,11 +853,13 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
     }
   }
 
-  double qMin = 0.5;
-  double sliverLimit = 0.1;
+  double qMin = 0.3;
+  double sliverLimit = 0.04;
 
   int nbESwap = 0, nbFSwap = 0, nbReloc = 0;
 
+  double worstA = 0.0;
+  
   while (1){
     //    printf("coucou\n");
     std::vector<MTet4*> newTets;
@@ -778,7 +868,7 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
         double qq = (*it)->getQuality();
         if (qq < qMin){
           for (int i = 0; i < 4; i++){
-            if (0 && faceSwap(newTets, *it, i, qm)){
+            if (faceSwap(newTets, *it, i, qm)){
               nbFSwap++;
               break;
             }
@@ -869,6 +959,8 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
     double t2 = Cpu();
     Msg::Info("Opti : (%d,%d,%d) = %12.5E QBAD %12.5E QAVG %12.5E (%8.3f sec)",
               nbESwap, nbFSwap, nbReloc, totalVolumeb, worst, avg / count, t2 - t1);
+    if (worstA != 0.0 && worst - worstA < 1.e-6)break;
+    worstA = worst;
   }
 
   if (illegals.size()){
@@ -989,6 +1081,31 @@ static void memoryCleanup(MTet4Factory &myFactory, std::set<MTet4*, compareTet4P
   //  Msg::Info("cleaning up the memory %d -> %d", n1, allTets.size());
 }
 
+
+int isCavityCompatibleWithEmbeddedEdges(std::list<MTet4*> &cavity,
+					std::list<faceXtet> &shell,   
+					edgeContainerB &allEmbeddedEdges){
+  
+
+  std::vector<MEdge> ed;
+  for (std::list<faceXtet>::iterator it = shell.begin(); it != shell.end();it++){
+    ed.push_back(MEdge(it->v[0],it->v[1]));
+    ed.push_back(MEdge(it->v[1],it->v[2]));
+    ed.push_back(MEdge(it->v[2],it->v[0]));
+  }
+  
+  for (std::list<MTet4*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc){
+    for (int j=0;j<6;j++){
+      MEdge e = (*itc)->tet()->getEdge(j);
+      if (std::find(ed.begin(), ed.end(), e) == ed.end() && allEmbeddedEdges.find(e)){
+	return 0;
+      }
+    }
+  }
+  return 1;
+} 
+
+
 int isCavityCompatibleWithEmbeddedEdges(std::list<MTet4*> &cavity,
 					std::list<faceXtet> &shell,
 					std::set<MEdge, Less_Edge> &allEmbeddedEdges){
@@ -1025,30 +1142,29 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
   std::set<MTet4*, compareTet4Ptr> &allTets = myFactory.getAllTets();
   int NUM = 0;
 
-	//	printTets ("before.pos", cavity, true);
-  std::set<MEdge, Less_Edge> allEmbeddedEdges;
-  createAllEmbeddedEdges (gr, allEmbeddedEdges);
 
   { // leave this in a block so the map gets deallocated directly
     std::map<MVertex*, double,MVertexLessThanNum> vSizesMap;
     std::set<MVertex*,MVertexLessThanNum> bndVertices;
 
-    std::set<MEdge, Less_Edge>::iterator it = allEmbeddedEdges.begin();
-    while (it != allEmbeddedEdges.end()){
-      MVertex *vi = it->getVertex(0);
-      MVertex *vj = it->getVertex(1);
-      double dx = vi->x()-vj->x();
-      double dy = vi->y()-vj->y();
-      double dz = vi->z()-vj->z();
-      double l = sqrt(dx * dx + dy * dy + dz * dz);
-      std::map<MVertex*,double,MVertexLessThanNum>::iterator iti = vSizesMap.find(vi);
-      std::map<MVertex*,double,MVertexLessThanNum>::iterator itj = vSizesMap.find(vj);
-      // smallest tet edge
-      if (iti == vSizesMap.end() || iti->second > l) vSizesMap[vi] = l;
-      if (itj == vSizesMap.end() || itj->second > l) vSizesMap[vj] = l;
+    std::list<GEdge*> e = gr->embeddedEdges();
+    for (std::list<GEdge*>::iterator it = e.begin() ; it != e.end(); ++it){
+      for (unsigned int i = 0; i < (*it)->lines.size(); i++){
+	MVertex *vi = (*it)->lines[i]->getVertex(0);
+	MVertex *vj = (*it)->lines[i]->getVertex(1);
+	double dx = vi->x()-vj->x();
+	double dy = vi->y()-vj->y();
+	double dz = vi->z()-vj->z();
+	double l = sqrt(dx * dx + dy * dy + dz * dz);
+	std::map<MVertex*,double,MVertexLessThanNum>::iterator iti = vSizesMap.find(vi);
+	std::map<MVertex*,double,MVertexLessThanNum>::iterator itj = vSizesMap.find(vj);
+	// smallest tet edge
+	if (iti == vSizesMap.end() || iti->second > l) vSizesMap[vi] = l;
+	if (itj == vSizesMap.end() || itj->second > l) vSizesMap[vj] = l;
+      }
       ++it;
     }
-
+  
     for(GModel::fiter it = gr->model()->firstFace(); it != gr->model()->lastFace(); ++it){
       GFace *gf = *it;
       for(unsigned int i = 0; i < gf->triangles.size(); i++){
@@ -1057,6 +1173,7 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
     }
     for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
       setLcs(gr->tetrahedra[i], vSizesMap, bndVertices);
+
     for(std::map<MVertex*, double, MVertexLessThanNum>::iterator it = vSizesMap.begin();
         it != vSizesMap.end(); ++it){
       it->first->setIndex(NUM++);
@@ -1123,6 +1240,8 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
   // store all embedded faces
   std::set<MFace, Less_Face> allEmbeddedFaces;
   createAllEmbeddedFaces (gr, allEmbeddedFaces);
+  edgeContainerB allEmbeddedEdges;
+  createAllEmbeddedEdges (gr, allEmbeddedEdges);
   connectTets(allTets.begin(), allTets.end(),&allEmbeddedFaces);
   Msg::Debug("All %d tets were connected", allTets.size());
 
@@ -1153,7 +1272,7 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
       allTets.erase(allTets.begin());
     }
     else{
-      if(ITER++ % 5000 == 0)
+      if(ITER++ % 500 == 0)
         Msg::Info("%d points created - Worst tet radius is %g (PTS removed %d %d)",
                  REALCOUNT, worst->getRadius(), COUNT_MISS_1,COUNT_MISS_2);
       if(worst->getRadius() < 1) break;
@@ -1214,6 +1333,7 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
 	  if (k == -1){starShaped = false ; break;}
 	  else if (k == 0) break;
 	  else if (k == 1) correctCavity = true;
+	  //	  printf("coucou\n");
 	}
 	if (correctCavity && starShaped) {
 	  NB_CORRECTION_OF_CAVITY ++;

Mesh/meshGRegionLocalMeshMod.cpp

@@ -106,7 +106,7 @@ bool buildEdgeCavity(MTet4 *t, int iLocalEdge, MVertex **v1, MVertex **v2,
     }
     // FIXME when hybrid mesh, this loops for ever
     if (cavity.size() > 1000) {
-      printf("cavity size gets laaaarge\n");
+      //      printf("cavity size gets laaaarge\n");
       return false;
     }
     //    printf("%d %d\n",ITER++, cavity.size());
@@ -210,8 +210,7 @@ bool edgeSwap(std::vector<MTet4 *> &newTets,
   //static int edges[6][2] =    {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
   int permut[6] = {0,3,1,2,5,4};
   iLocalEdge = permut[iLocalEdge];
-  
-  
+ 
   std::vector<MTet4*> cavity;
   std::vector<MTet4*> outside;
   std::vector<MVertex*> ring;

Mesh/qualityMeasures.cpp

@@ -3,6 +3,7 @@
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to the public mailing list <gmsh@onelab.info>.
 
+#include "robustPredicates.h"
 #include "qualityMeasures.h"
 #include "BDS.h"
 #include "MVertex.h"
@@ -636,17 +637,13 @@ double qmTetrahedron::eta(const double &x1, const double &y1, const double &z1,
                           const double &x4, const double &y4, const double &z4,
                           double *volume)
 {
-  double mat[3][3];
-  mat[0][0] = x2 - x1;
-  mat[0][1] = x3 - x1;
-  mat[0][2] = x4 - x1;
-  mat[1][0] = y2 - y1;
-  mat[1][1] = y3 - y1;
-  mat[1][2] = y4 - y1;
-  mat[2][0] = z2 - z1;
-  mat[2][1] = z3 - z1;
-  mat[2][2] = z4 - z1;
-  *volume = fabs(det3x3(mat)) / 6.;
+  double p0[3] = {x1, y1, z1};
+  double p1[3] = {x2, y2, z2};
+  double p2[3] = {x3, y3, z3};
+  double p3[3] = {x4, y4, z4};
+
+  *volume =fabs (robustPredicates::orient3d (p0,p1,p2,p3)) / 6.0;
+
   double l = ((x2 - x1) * (x2 - x1) +
               (y2 - y1) * (y2 - y1) +
               (z2 - z1) * (z2 - z1));
@@ -665,40 +662,46 @@ double qmTetrahedron::gamma(const double &x1, const double &y1, const double &z1
                             const double &x4, const double &y4, const double &z4,
                             double *volume)
 {
-  double mat[3][3];
-  mat[0][0] = x2 - x1;
-  mat[0][1] = x3 - x1;
-  mat[0][2] = x4 - x1;
-  mat[1][0] = y2 - y1;
-  mat[1][1] = y3 - y1;
-  mat[1][2] = y4 - y1;
-  mat[2][0] = z2 - z1;
-  mat[2][1] = z3 - z1;
-  mat[2][2] = z4 - z1;
-  *volume = fabs(det3x3(mat)) / 6.;
   double p0[3] = {x1, y1, z1};
   double p1[3] = {x2, y2, z2};
   double p2[3] = {x3, y3, z3};
   double p3[3] = {x4, y4, z4};
+
+  *volume =fabs (robustPredicates::orient3d (p0,p1,p2,p3)) / 6.0;
+  
+  // double mat[3][3];
+  // mat[0][0] = x2 - x1;
+  // mat[0][1] = x3 - x1;
+  // mat[0][2] = x4 - x1;
+  // mat[1][0] = y2 - y1;
+  // mat[1][1] = y3 - y1;
+  // mat[1][2] = y4 - y1;
+  // mat[2][0] = z2 - z1;
+  // mat[2][1] = z3 - z1;
+  // mat[2][2] = z4 - z1;
+  // *volume = fabs(det3x3(mat)) / 6.;
+
+  // printf("vv %g volume %g\n",vv,*volume);
+
   double s1 = fabs(triangle_area(p0, p1, p2));
   double s2 = fabs(triangle_area(p0, p2, p3));
   double s3 = fabs(triangle_area(p0, p1, p3));
   double s4 = fabs(triangle_area(p1, p2, p3));
   double rhoin = 3. * fabs(*volume) / (s1 + s2 + s3 + s4);
-  double l = sqrt((x2 - x1) * (x2 - x1) +
-                  (y2 - y1) * (y2 - y1) +
-                  (z2 - z1) * (z2 - z1));
-  l = std::max(l, sqrt((x3 - x1) * (x3 - x1) + (y3 - y1) * (y3 - y1) +
-                       (z3 - z1) * (z3 - z1)));
-  l = std::max(l, sqrt((x4 - x1) * (x4 - x1) + (y4 - y1) * (y4 - y1) +
-                       (z4 - z1) * (z4 - z1)));
-  l = std::max(l, sqrt((x3 - x2) * (x3 - x2) + (y3 - y2) * (y3 - y2) +
-                       (z3 - z2) * (z3 - z2)));
-  l = std::max(l, sqrt((x4 - x2) * (x4 - x2) + (y4 - y2) * (y4 - y2) +
-                       (z4 - z2) * (z4 - z2)));
-  l = std::max(l, sqrt((x3 - x4) * (x3 - x4) + (y3 - y4) * (y3 - y4) +
-                       (z3 - z4) * (z3 - z4)));
-  return 2. * sqrt(6.) * rhoin / l;
+  double l = (x2 - x1) * (x2 - x1) +
+    (y2 - y1) * (y2 - y1) +
+    (z2 - z1) * (z2 - z1);
+  l = std::max(l, ((x3 - x1) * (x3 - x1) + (y3 - y1) * (y3 - y1) +
+		   (z3 - z1) * (z3 - z1)));
+  l = std::max(l, ((x4 - x1) * (x4 - x1) + (y4 - y1) * (y4 - y1) +
+		   (z4 - z1) * (z4 - z1)));
+  l = std::max(l, ((x3 - x2) * (x3 - x2) + (y3 - y2) * (y3 - y2) +
+		   (z3 - z2) * (z3 - z2)));
+  l = std::max(l, ((x4 - x2) * (x4 - x2) + (y4 - y2) * (y4 - y2) +
+		   (z4 - z2) * (z4 - z2)));
+  l = std::max(l, ((x3 - x4) * (x3 - x4) + (y3 - y4) * (y3 - y4) +
+		   (z3 - z4) * (z3 - z4)));
+  return sqrt(24./l) * rhoin ;
 }
 
 

benchmarks/boolean/aneurysmBL.py

+from gmshpy import *
+
+GmshSetOption('General', 'Verbosity', 99.0) 
+GmshSetOption('Mesh', 'SaveAll', 1.0) 
+GmshSetOption('Mesh', 'Algorithm', 1.0)
+
+g =  GModel()
+g.setFactory('Gmsh')
+
+g.load("aneurysm.stl");
+g.createTopologyFromMesh();
+
+face = g.getFaceByTag(1);
+newent  = g.extrudeBoundaryLayer(face, 4, 0.5, 0, -1)
+newent2 = g.extrudeBoundaryLayer(face, 4, -0.5, 1, -1)
+
+print "*** face = %d " % face.tag()
+print "*** new face = %d newface2 = %d " % (newent[0].tag(), newent2[0].tag())
+
+g.mesh(2)
+g.save("aneurysmBL.msh")
+

benchmarks/boolean/neuron.geo

+SetFactory("OpenCASCADE");
+
+Mesh.Algorithm = 6;
+Mesh.CharacteristicLengthMin = 1;
+Mesh.CharacteristicLengthMax = 1;
+
+Macro dendrite
+  For i In {1:5}
+    z = -2+7*i;
+    r = 1 + 0.6*Sin(2*Pi*i/5.);
+    Point(nump+1) = {x,0,z};
+    Point(nump+2) = {x+r,0,z};
+    Point(nump+3) = {x,r,z};
+    Point(nump+4) = {x-r,0,z};
+    Point(nump+5) = {x,-r,z};
+    Circle(numc+1) = {nump+2,nump+1,nump+3};
+    Circle(numc+2) = {nump+3,nump+1,nump+4};
+    Circle(numc+3) = {nump+4,nump+1,nump+5};
+    Circle(numc+4) = {nump+5,nump+1,nump+2};
+    edges~{i}() = {numc+1:numc+4};
+    nump += 5;
+    numc += 4;
+  EndFor
+  numr += 1;
+  ThruSections(numr) = {edges~{1}(),edges~{2}(),edges~{3}(),edges~{4}(),edges~{5}()};
+  reg() += numr;
+Return
+
+Sphere(1) = {0, 0, 0, 8};
+
+reg() = {};
+nump = 0; numc = 0; numr = 2;
+For x In{-4:4:4}
+  Call dendrite;
+EndFor
+
+DefineConstant[
+  op = {0, Choices{0="None", 1="Union", 2="Intersection", 3="Subtraction"},
+    Name "Boolean operation" }
+];
+
+// boolean operations can explicitly create an entity with the form
+// "op(tag)={}{};", or let Gmsh decide with the form "op{}{}". The first form
+// can only be used if the result of the boolean operation is a single
+// shape. Only the second form returns the list of created entities.
+If(op == 1)
+  BooleanUnion(100) = { Volume{1}; Delete; }{ Volume{reg(0)}; Delete; };
+  BooleanUnion(101) = { Volume{100}; Delete; }{ Volume{reg(1)}; Delete; };
+  BooleanUnion(102) = { Volume{101}; Delete; }{ Volume{reg(2)}; Delete; };
+ElseIf(op == 2)
+  BooleanIntersection(100) = { Volume{1}; }{ Volume{reg(0)}; Delete; };
+  BooleanIntersection(101) = { Volume{1}; }{ Volume{reg(1)}; Delete; };
+  BooleanIntersection(102) = { Volume{1}; Delete; }{ Volume{reg(2)}; Delete; };
+ElseIf(op == 3)
+  BooleanSubtraction(100) = { Volume{1}; Delete; }{ Volume{reg(0)}; Delete; };
+  BooleanSubtraction(101) = { Volume{100}; Delete; }{ Volume{reg(1)}; Delete; };
+  BooleanSubtraction(102) = { Volume{101}; Delete; }{ Volume{reg(2)}; Delete; };
+EndIf

benchmarks/boolean/sphere.py

+from gmshpy import *
+
+GmshSetOption('Mesh', 'CharacteristicLengthFactor', 0.9)
+
+R = 0.3;
+myModel = GModel();
+myModel.addSphere(0.5,0.5,0.5,R);
+
+myModel.setAsCurrent();
+
+myModel.mesh(2);
+myModel.save("sphere.stl");

benchmarks/boolean/square1.py

+from gmshpy import *
+lc = 0.5
+GmshSetOption('Mesh', 'CharacteristicLengthFactor', lc)
+ 
+g = GModel()
+
+v1 = g.addVertex(0, 0, 0, lc)
+v2 = g.addVertex(1, 0, 0, lc)
+v3 = g.addVertex(1, 1, 0, lc)
+v4 = g.addVertex(0, 1, 0, lc)
+e1 = g.addLine(v2, v1)
+e2 = g.addLine(v3, v2)
+e3 = g.addLine(v4, v3)
+e4 = g.addLine(v4, v1)
+v11 = g.addVertex(.4, .4, 0, lc)
+v12 = g.addVertex(.6, .4, 0, lc)
+v13 = g.addVertex(.6, .5, 0, lc)
+v14 = g.addVertex(.4, .6, 0, lc)
+e11 = g.addLine(v11, v12)
+e12 = g.addLine(v12, v13)
+e13 = g.addLine(v13, v14)
+e14 = g.addLine(v14, v11)
+
+f = g.addPlanarFace ([[e1,e2,e3,e4],[e11,e12,e13,e14]])
+
+g.mesh(2)
+g.save("square1.msh")
+
+#v100 = g.addVertex(0, 0, 0, .1)
+#v200 = g.addVertex(0, 0, 1, .1)
+#v300 = g.addVertex(0, .1, .33, .1)
+#v400 = g.addVertex(.1, .1, .66, .1)
+#line = g.addBezier(v100,v200,{{v300:x(),v300:y(),v300:z()},{v400:x(),v400:y(),v400:z()}});
+#g.addPipe (f, {line})
+#g.glue(1.e-9);
+
+#myTool = GModel();
+#myTool:addSphere(0.2,0.2,0.1,.52012);
+
+#g.addSphere(1,1.3,1,.3);
+#g.computeDifference(myTool,0);
+
+#g.setAsCurrent();

benchmarks/boolean/wikipedia.geo

+SetFactory("OpenCASCADE");
+
+Mesh.Algorithm = 6;
+Mesh.CharacteristicLengthMin = 0.4;
+Mesh.CharacteristicLengthMax = 0.4;
+
+R = DefineNumber[ 1.4 , Name "Parameters/Dimension of the Block" ];
+s = DefineNumber[  .7 , Name "Parameters/Cylinder of radius s *R" ];
+t = DefineNumber[ 1.25, Name "Parameters/Sphere   of radius t *R" ];
+
+Block(1) = {-R,-R,-R, R,R,R};
+
+Sphere(2) = {0,0,0,R*t};
+
+BooleanIntersection(3) = { Volume{1}; Delete; }{ Volume{2}; Delete; };
+
+Cylinder(4) = {-2*R,0,0, 2*R,0,0, R*s};
+Cylinder(5) = {0,-2*R,0, 0,2*R,0, R*s};
+Cylinder(6) = {0,0,-2*R, 0,0,2*R, R*s};
+
+BooleanUnion(7) = { Volume{4}; Delete; }{ Volume{5}; Delete; };
+BooleanUnion(8) = { Volume{6}; Delete; }{ Volume{7}; Delete; };
+BooleanSubtraction(9) = { Volume{3}; Delete; }{ Volume{8}; Delete; };
+
+

benchmarks/boolean/wikipedia.py

+#!/usr/bin/env python
+
+from gmshpy import *
+
+# from http://en.wikipedia.org/wiki/Constructive_solid_geometry
+
+R = 1.4;
+s = .7;
+t = 1.25;
+myModel = GModel();
+myModel.addBlock([-R,-R,-R],[R,R,R]);
+
+myTool = GModel();
+myTool.addSphere(0,0,0,R*t);
+
+myModel.computeBooleanIntersection(myTool);
+
+myTool2 = GModel();
+myTool2.addCylinder([-2*R,0,0],[2*R,0,0],R*s);
+
+myTool3 = GModel();
+myTool3.addCylinder([0,-2*R,0],[0,2*R,0],R*s);
+
+myModel2 = GModel();
+myModel2.addCylinder([0,0,-2*R],[0,0,2*R],R*s);
+myModel2.computeBooleanUnion(myTool2);
+myModel2.computeBooleanUnion(myTool3);
+
+myModel.computeBooleanDifference(myModel2);
+
+myModel2.setVisibility(0);
+myModel.setAsCurrent();
+myModel.setVisibility(1);
+
+GmshSetOption("Mesh", "CharacteristicLengthFactor", 0.4);
+
+#myModel.mesh(3);
+#myModel.save("wikipedia.msh");
+#myModel.save("wikipedia.brep");
+
+FlGui.instance();
+FlGui.run();