From c933defa4ceb3b8ce30c8018076638f829d77760 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Tue, 3 Aug 2010 10:47:38 +0000
Subject: [PATCH] make sure we mesh extruded faces *after* other faces

---
 Common/Octree.cpp                                   | 11 ++++++-----
 Geo/GFace.cpp                                       |  8 ++++++++
 Geo/GFace.h                                         |  3 +++
 Geo/MVertexOctree.cpp                               | 10 +++++++---
 Geo/MVertexOctree.h                                 |  4 +++-
 Mesh/Generator.cpp                                  | 12 ++++++++----
 Mesh/meshGRegion.cpp                                |  1 -
 benchmarks/extrude/t1_boundary_layer_connection.geo |  3 ++-
 8 files changed, 37 insertions(+), 15 deletions(-)

diff --git a/Common/Octree.cpp b/Common/Octree.cpp
index 055a8ccb34..4bcdc092ef 100644
--- a/Common/Octree.cpp
+++ b/Common/Octree.cpp
@@ -14,8 +14,8 @@ Octree* Octree_Create(int maxElements, double origin[3], double size[3],
                       int   (*InEle)(void *, double *))
 {
   Octree *myOctree = new Octree;
-  initializeOctantBuckets (origin, size, maxElements,
-                           &(myOctree->root), &(myOctree->info));       
+  initializeOctantBuckets(origin, size, maxElements,
+                          &(myOctree->root), &(myOctree->info));       
   myOctree->function_BB = BB;
   myOctree->function_centroid = Centroid;
   myOctree->function_inElement = InEle;
@@ -53,7 +53,7 @@ void Octree_Delete(Octree *myOctree)
   delete myOctree;
 }
 
-void Octree_Insert(void * element, Octree *myOctree)
+void Octree_Insert(void *element, Octree *myOctree)
 {
   if(!myOctree) return;
   double minBB[3], maxBB[3], centroid[3];
@@ -61,8 +61,9 @@ void Octree_Insert(void * element, Octree *myOctree)
   (*(myOctree->function_BB))(element, minBB, maxBB);
   (*(myOctree->function_centroid))(element, centroid);
   bucket = findElementBucket(myOctree->root, centroid);
-  addElement2Bucket(bucket, element, minBB, maxBB,
-                    centroid, myOctree->info);  
+  if(bucket)
+    addElement2Bucket(bucket, element, minBB, maxBB,
+                      centroid, myOctree->info);  
 }
 
 void Octree_Arrange(Octree *myOctree)
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index ec953370dd..760ae369a7 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -16,6 +16,7 @@
 #include "fullMatrix.h"
 #include "Numeric.h"
 #include "GaussLegendre1D.h"
+#include "ExtrudeParams.h"
 #include "Context.h"
 #include "meshGFaceLloyd.h"
 #include "Bindings.h"
@@ -137,6 +138,13 @@ void GFace::resetMeshAttributes()
   meshAttributes.extrude = 0;
 }
 
+bool GFace::isMeshExtruded()
+{
+  ExtrudeParams *ep = meshAttributes.extrude;
+  if(ep && ep->mesh.ExtrudeMesh) return true;
+  return false;
+}
+
 SBoundingBox3d GFace::bounds() const
 {
   SBoundingBox3d res;
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 00aee5a100..f9f3fef914 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -235,6 +235,9 @@ class GFace : public GEntity
   // reset the mesh attributes to default values
   virtual void resetMeshAttributes();
 
+  // tells if the mesh is obtained by extrusion
+  bool isMeshExtruded();
+
   // for periodic faces, move parameters into the range chosen
   // for that face
   void moveToValidRange(SPoint2 &pt) const;
diff --git a/Geo/MVertexOctree.cpp b/Geo/MVertexOctree.cpp
index 0dedfcd98f..9dcc936a81 100644
--- a/Geo/MVertexOctree.cpp
+++ b/Geo/MVertexOctree.cpp
@@ -9,10 +9,12 @@
 #include "MVertexOctree.h"
 #include "SBoundingBox3d.h"
 
+double MVertexOctree::tolerance = 1.e-6;
+
 static void MVertexBB(void *a, double *min, double *max)
 {
   MVertex *v = (MVertex*)a;
-  double tol = MVertexLessThanLexicographic::tolerance;
+  double tol = MVertexOctree::tolerance;
   min[0] = v->x() - tol;
   max[0] = v->x() + tol;
   min[1] = v->y() - tol;
@@ -35,12 +37,14 @@ static int MVertexInEle(void *a, double *x)
   double d = sqrt((x[0] - v->x()) * (x[0] - v->x()) +
                   (x[1] - v->y()) * (x[1] - v->y()) +
                   (x[2] - v->z()) * (x[2] - v->z()));
-  return (d < MVertexLessThanLexicographic::tolerance);
+  return (d < MVertexOctree::tolerance);
 }
 
-MVertexOctree::MVertexOctree(GModel *m)
+MVertexOctree::MVertexOctree(GModel *m, double tol)
 {
+  tolerance = tol;
   SBoundingBox3d bb = m->bounds();
+  bb *= 1.2;
   double min[3] = {bb.min().x(), bb.min().y(), bb.min().z()};
   double size[3] = {bb.max().x() - bb.min().x(),
                     bb.max().y() - bb.min().y(),
diff --git a/Geo/MVertexOctree.h b/Geo/MVertexOctree.h
index e94e5b9104..eb4b5d53c0 100644
--- a/Geo/MVertexOctree.h
+++ b/Geo/MVertexOctree.h
@@ -9,7 +9,9 @@ class MVertexOctree{
  private:
   Octree *_octree;
  public:
-  MVertexOctree(GModel *);
+  static double tolerance;
+ public:
+  MVertexOctree(GModel *, double);
   ~MVertexOctree();
   void insert(MVertex *);
   void finalize();
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index e8fd043438..6ab042e0eb 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -439,14 +439,18 @@ static void Mesh2D(GModel *m)
   // field generated from the surface mesh of the source surfaces
   if(!Mesh2DWithBoundaryLayers(m)){
   
-    std::set<GFace*> faces, compoundFaces;
+    std::set<GFace*> compoundFaces, extrudedFaces, otherFaces;
     for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
-      if ((*it)->geomType() == GEntity::CompoundSurface)
+      GFace *gf = *it;
+      if (gf->geomType() == GEntity::CompoundSurface)
         compoundFaces.insert(*it);
+      else if(gf->isMeshExtruded())
+        extrudedFaces.insert(gf);
       else
-        faces.insert(*it);
+        otherFaces.insert(gf);
     }
-    std::for_each(faces.begin(), faces.end(), meshGFace());
+    std::for_each(otherFaces.begin(), otherFaces.end(), meshGFace());
+    std::for_each(extrudedFaces.begin(), extrudedFaces.end(), meshGFace());
     std::for_each(compoundFaces.begin(), compoundFaces.end(), meshGFace());
 
     // lloyd optimization
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 689d63208c..008eb286df 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -583,7 +583,6 @@ void meshGRegion::operator() (GRegion *gr)
   if(CTX::instance()->mesh.meshOnlyVisible && !gr->getVisibility()) return;
 
   ExtrudeParams *ep = gr->meshAttributes.extrude;
-
   if(ep && ep->mesh.ExtrudeMesh) return;
 
   // destroy the mesh if it exists
diff --git a/benchmarks/extrude/t1_boundary_layer_connection.geo b/benchmarks/extrude/t1_boundary_layer_connection.geo
index 41a3643f92..21f04bb8da 100644
--- a/benchmarks/extrude/t1_boundary_layer_connection.geo
+++ b/benchmarks/extrude/t1_boundary_layer_connection.geo
@@ -16,9 +16,10 @@ Extrude { Surface{6}; Layers{5, 0.01}; Recombine; }
 
 Point(100) = {0,  0, 0.15, lc} ;
 Point(101) = {.1, 0, 0.15, lc} ;
-
+/*
 Line(100) = {8, 100};
 Line(101) = {12, 101};
 Line(102) = {100, 101};
 Line Loop(100) = {9,101,-102,-100};
 Plane Surface(100) = 100;
+*/
-- 
GitLab