hopla

Mesh/meshGRegionDelaunayInsertion.cpp

@@ -75,8 +75,38 @@ struct faceXtet{
     if (v[2] < other.v[2]) return true;
     return false;
   }
+  inline bool operator == (const faceXtet & other) const
+  {
+    return (v[0] == other.v[0] &&
+	    v[1] == other.v[1] &&
+	    v[2] == other.v[2] );
+  }
 };
 
+template <class ITER>
+void connectTets_vector(ITER beg, ITER end)
+{
+  //  std::set<faceXtet> conn;
+  std::vector<faceXtet> conn;
+  while (beg != end){
+    if (!(*beg)->isDeleted()){
+      for (int i = 0; i < 4; i++){
+        faceXtet fxt(*beg, i);
+	std::vector<faceXtet>::iterator found  = std::find(conn.begin(), conn.end(), fxt); 
+	//        std::set<faceXtet>::iterator found = conn.find(fxt);
+        if (found == conn.end())
+	  conn.push_back(fxt);
+	// conn.insert(fxt);
+        else if (found->t1 != *beg){
+          found->t1->setNeigh(found->i1, *beg);
+          (*beg)->setNeigh(i, found->t1);
+        }
+      }
+    }
+    ++beg;
+  }
+}
+
 template <class ITER>
 void connectTets(ITER beg, ITER end)
 {
@@ -98,6 +128,7 @@ void connectTets(ITER beg, ITER end)
   }
 }
 
+
 static void updateActiveFaces(MTet4 *t, double limit_, std::set<MFace,Less_Face> &front)
 {
   if (t->isDeleted()) return;
@@ -160,8 +191,50 @@ void recurFindCavity(std::list<faceXtet> & shell,
         shell.push_back(faceXtet(t, i));
     }
   }
+  //  printf("cavity size %d\n",cavity.size());
 }
 
+void nonrecurFindCavity(std::list<faceXtet> & shell,
+                     std::list<MTet4*> & cavity,
+                     MVertex *v ,
+                     MTet4 *t)
+{
+  // Msg::Info("tet %d %d %d %d",t->tet()->getVertex(0)->getNum(),
+  //     t->tet()->getVertex(1)->getNum(),
+  //     t->tet()->getVertex(2)->getNum(),
+  //     t->tet()->getVertex(3)->getNum());
+
+  // invariant : this one has to be inserted in the cavity
+  // consider this tet deleted
+  // remove its reference to its neighbors
+
+  std::stack<MTet4*> _stack;
+  _stack.push(t);
+  while(!_stack.empty()){
+    t = _stack.top();
+    _stack.pop();
+    t->setDeleted(true);
+    // the cavity that has to be removed
+    // because it violates delaunay criterion
+    cavity.push_back(t);
+    
+    for (int i = 0; i < 4; i++){
+      MTet4 *neigh = t->getNeigh(i) ;
+      if (!neigh)
+	shell.push_back(faceXtet(t, i));
+      else  if (!neigh->isDeleted()){
+	int circ = neigh->inCircumSphere(v);
+	if (circ && (neigh->onWhat() == t->onWhat()))
+	  _stack.push(neigh);
+	else
+	  shell.push_back(faceXtet(t, i));
+      }
+    }
+  }
+  //  printf("cavity size %d\n",cavity.size());
+}
+
+
 bool insertVertex(MVertex *v,
                   MTet4 *t,
                   MTet4Factory &myFactory,
@@ -237,7 +310,7 @@ bool insertVertex(MVertex *v,
   // OK, the cavity is star shaped
   if (fabs(oldVolume - newVolume) < 1.e-10 * oldVolume &&
       !onePointIsTooClose){
-    connectTets(new_cavity.begin(), new_cavity.end());
+    connectTets_vector(new_cavity.begin(), new_cavity.end());
     allTets.insert(newTets, newTets + shell.size());
 
     if (activeTets){
@@ -344,20 +417,56 @@ GRegion *getRegionFromBoundingFaces(GModel *model,
   return 0;
 }
 
-void recur_classify(MTet4 *t, std::list<MTet4*> &theRegion,
+// void recur_classify(MTet4 *t, std::list<MTet4*> &theRegion,
+//                     std::set<GFace*> &faces_bound, GRegion *bidon,
+//                     GModel *model, const fs_cont &search)
+// {
+//   if (!t) Msg::Error("a tet is not connected by a boundary face");
+//   if (t->onWhat()) {
+//     return; // should never return here...
+//   }
+//   theRegion.push_back(t);
+//   t->setOnWhat(bidon);
+
+//   bool FF[4] = {0,0,0,0};
+
+//   for (int i = 0; i < 4; i++){
+//     //      if (!t->getNeigh(i) || !t->getNeigh(i)->onWhat())
+//     {
+//       GFace* gfound = findInFaceSearchStructure (t->tet()->getVertex(faces[i][0]),
+//                                                  t->tet()->getVertex(faces[i][1]),
+//                                                  t->tet()->getVertex(faces[i][2]),
+//                                                  search);
+//       if (gfound){
+//         FF[i] = true;
+//         if (faces_bound.find(gfound) == faces_bound.end())
+//           faces_bound.insert(gfound);
+//       }
+//     }
+//   }
+//   for (int i = 0; i < 4; i++){
+//     if (!FF[i])
+//       recur_classify(t->getNeigh(i), theRegion, faces_bound, bidon, model, search );
+//   }
+// }
+
+void non_recursive_classify(MTet4 *t, std::list<MTet4*> &theRegion,
 			    std::set<GFace*> &faces_bound, GRegion *bidon,
 			    GModel *model, const fs_cont &search)
 {
+
+  std::stack<MTet4*> _stackounette;
+  _stackounette.push(t);
+
+  while(!_stackounette.empty()){
+    t = _stackounette.top();
+    _stackounette.pop();
     if (!t) Msg::Error("a tet is not connected by a boundary face");
-  if (t->onWhat()) return; // should never return here...
+    if (!t->onWhat()) {
       theRegion.push_back(t);
       t->setOnWhat(bidon);
-
       bool FF[4] = {0,0,0,0};
-
       for (int i = 0; i < 4; i++){
-    //      if (!t->getNeigh(i) || !t->getNeigh(i)->onWhat())
-    {
 	GFace* gfound = findInFaceSearchStructure (t->tet()->getVertex(faces[i][0]),
 						   t->tet()->getVertex(faces[i][1]),
 						   t->tet()->getVertex(faces[i][2]),
@@ -367,23 +476,17 @@ void recur_classify(MTet4 *t, std::list<MTet4*> &theRegion,
 	  if (faces_bound.find(gfound) == faces_bound.end())
 	    faces_bound.insert(gfound);
 	}
-
-//        MTriangle tri (t->tet()->getVertex (faces[i][0]),
-//                       t->tet()->getVertex (faces[i][1]),
-//                       t->tet()->getVertex (faces[i][2]));
-//        GFace *gfound;
-//        if (FF[i] = find_triangle_in_model(model, &tri, &gfound, false)){
-//          if (faces_bound.find(gfound) == faces_bound.end())
-//            faces_bound.insert(gfound);
-//        }
-    }
       }
       for (int i = 0; i < 4; i++){
 	if (!FF[i])
-      recur_classify(t->getNeigh(i), theRegion, faces_bound, bidon, model, search );
+	  _stackounette.push(t->getNeigh(i));
+      }
+    }
   }
 }
 
+
+
 void adaptMeshGRegion::operator () (GRegion *gr)
 {
   const qualityMeasure4Tet qm = QMTET_2;
@@ -744,7 +847,7 @@ void optimizeMesh(GRegion *gr, const qualityMeasure4Tet &qm)
 }
 
 
-static double tetcircumcenter(double a[3], double b[3], double c[3], double d[3],
+double tetcircumcenter(double a[3], double b[3], double c[3], double d[3],
 			      double circumcenter[3], double *xi, double *eta, double *zeta){
   double xba, yba, zba, xca, yca, zca, xda, yda, zda;
   double balength, calength, dalength;
@@ -865,7 +968,7 @@ void insertVerticesInRegion (GRegion *gr)
       GRegion *bidon = (GRegion*)123;
       double _t1 = Cpu();
       Msg::Debug("start with a non classified tet");
-      recur_classify(*it, theRegion, faces_bound, bidon, gr->model(), search);
+      non_recursive_classify(*it, theRegion, faces_bound, bidon, gr->model(), search);
       double _t2 = Cpu();
       Msg::Debug("found %d tets with %d faces (%g sec for the classification)",
                  theRegion.size(), faces_bound.size(), _t2 - _t1);
@@ -894,6 +997,8 @@ void insertVerticesInRegion (GRegion *gr)
 
   int ITER = 0;
 
+  int COUNT_MISS_1 = 0;
+  int COUNT_MISS_2 = 0;
   while(1){
     if(allTets.empty()){
       Msg::Error("No tetrahedra in region %d %d", gr->tag(), allTets.size());
@@ -908,8 +1013,8 @@ void insertVerticesInRegion (GRegion *gr)
     }
     else{
       if(ITER++ %5000 == 0)
-        Msg::Info("%d points created -- Worst tet radius is %g",
-                  vSizes.size(), worst->getRadius());
+        Msg::Info("%d points created -- Worst tet radius is %g (MISSES %d %d)",
+                  vSizes.size(), worst->getRadius(), COUNT_MISS_1,COUNT_MISS_2);
       if(worst->getRadius() < 1) break;
       double center[3];
       double uvw[3];
@@ -956,6 +1061,7 @@ void insertVerticesInRegion (GRegion *gr)
         vSizesBGM.push_back(lc);
         // compute mesh spacing there
         if(!insertVertex(v, worst, myFactory, allTets, vSizes,vSizesBGM)){
+	  COUNT_MISS_1++;
           myFactory.changeTetRadius(allTets.begin(), 0.);
           delete v;
         }
@@ -965,6 +1071,7 @@ void insertVerticesInRegion (GRegion *gr)
       else{
 	//	printf("point outside %12.5E %12.5E %12.5E %12.5E %12.5E\n",VV,uvw[0], uvw[1], uvw[2],1-uvw[0]-uvw[1]-uvw[2]);
         myFactory.changeTetRadius(allTets.begin(), 0.0);
+	COUNT_MISS_2++;
       }
     }
 
@@ -986,6 +1093,7 @@ void insertVerticesInRegion (GRegion *gr)
     }
   }
   
+
     // relocate vertices
   int nbReloc = 0;
   for (int SM=0;SM<CTX::instance()->mesh.nbSmoothing;SM++){
@@ -1094,7 +1202,7 @@ void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
       GRegion *bidon = (GRegion*)123;
       double _t1 = Cpu();
       Msg::Debug("start with a non classified tet");
-      recur_classify(*it, theRegion, faces_bound, bidon, gr->model(), search);
+      non_recursive_classify(*it, theRegion, faces_bound, bidon, gr->model(), search);
       double _t2 = Cpu();
       Msg::Debug("found %d tets with %d faces (%g sec for the classification)",
                  theRegion.size(), faces_bound.size(), _t2 - _t1);