splitMesh OK for centerlines (next step is classification of compounds)

Mesh/CenterlineField.cpp

@@ -19,12 +19,12 @@
 #include "OS.h"
 #include "GModel.h"
 #include "MElement.h"
+#include "MTriangle.h"
 #include "MVertex.h"
 #include "MLine.h"
 #include "GEntity.h"
 #include "Field.h"
 #include "GFace.h"
-#include "Integration3D.h"
 
 #if defined(HAVE_ANN)
 #include <ANN/ANN.h>
@@ -44,7 +44,48 @@ double computeLength(std::vector<MLine*> lines){
   return length;
 }
 
-void orderMLines(std::vector<MLine*> &lines) { //, std::set<MVertex*> & junctions, MVertex *vB, MVertex *vE){
+bool isClosed(std::set<MEdge, Less_Edge> &theCut){
+
+  std::multiset<MVertex*> boundV;
+  std::set<MEdge,Less_Edge>::iterator it = theCut.begin();
+  for (; it!= theCut.end(); it++){
+    boundV.insert(it->getVertex(0));
+    boundV.insert(it->getVertex(1));
+  }
+  std::multiset<MVertex*>::iterator itb = boundV.begin();
+  for ( ; itb != boundV.end(); ++itb){
+    if (boundV.count(*itb) != 2) {
+      return false;
+    }
+  }
+  return true;
+
+  // std::list<MEdge> segments;
+  // std::map<MVertex*, MEdge> boundv;
+  // std::set<MEdge,Less_Edge>::iterator it = theCut.begin();
+  // for (; it!= theCut.end(); it++){
+  //   segments.push_back(*it);
+  // }
+
+  // // find a lonely MEdge
+  // for (std::list<MEdge>::iterator it = segments.begin(); 
+  //      it != segments.end(); ++it){
+  //   MVertex *vL = it->getVertex(0);
+  //   MVertex *vR = it->getVertex(1);
+  //   std::map<MVertex*,MEdge>::iterator it1 = boundv.find(vL);
+  //   if (it1 == boundv.end()) boundv.insert(std::make_pair(vL,*it));
+  //   else  boundv.erase(it1);
+  //   std::map<MVertex*,MEdge>::iterator it2 = boundv.find(vR);
+  //   if (it2 == boundv.end()) boundv.insert(std::make_pair(vR,*it));
+  //   else boundv.erase(it2);
+  // }
+
+  // if (boundv.size() == 0) return true; 
+  // else return false;
+
+}
+
+void orderMLines(std::vector<MLine*> &lines, MVertex *vB, MVertex *vE){
 
   std::vector<MLine*> _m;
   std::list<MLine*> segments;
@@ -72,7 +113,13 @@ void orderMLines(std::vector<MLine*> &lines) { //, std::set<MVertex*> & junction
   // find the first MLine and erase it from the list segments
   MLine *firstLine;
   if (boundv.size() == 2){   // non periodic
-    firstLine = (boundv.begin())->second;
+    MVertex *v = (boundv.begin())->first;
+    if ( v == vB) firstLine = (boundv.begin())->second;
+    else{
+      MVertex *v = (boundv.rbegin())->first;
+      if (v == vB) firstLine = (boundv.rbegin())->second;
+      else{ Msg::Error("begin vertex not found for branch"); exit(1);}
+    }
     for (std::list<MLine*>::iterator it = segments.begin(); 
          it != segments.end(); ++it){
       if (*it == firstLine){
@@ -81,13 +128,8 @@ void orderMLines(std::vector<MLine*> &lines) { //, std::set<MVertex*> & junction
       }
     }
   }
-  else if (boundv.size() == 0){ // periodic
-    firstLine = *(segments.begin());
-    segments.erase(segments.begin());
-  }
   else{
-    Msg::Error("line is wrong (it has %d end points)", 
-               boundv.size());
+    Msg::Error("line is wrong (it has %d end points)",  boundv.size());
   }
 
   // loop over all segments to order segments and store it in the list _m
@@ -152,6 +194,92 @@ void orderMLines(std::vector<MLine*> &lines) { //, std::set<MVertex*> & junction
   // printf("in order vB= %d =%d \n", vB->getNum(), vE->getNum());
 }
 
+void cutTriangle(MTriangle *tri, 
+		 std::map<MEdge,MVertex*,Less_Edge> &cutEdges, 
+		 std::set<MVertex*> &cutVertices,
+		 std::vector<MTriangle*> &newTris, 
+		 std::set<MEdge,Less_Edge> &newCut){
+
+  MVertex *c[3] = {0,0,0};
+  for (int j=0;j<3;j++){
+    MEdge ed = tri->getEdge(j); 
+    std::map<MEdge,MVertex*,Less_Edge> :: iterator it = cutEdges.find(ed);
+    if (it != cutEdges.end()){
+      c[j] = it->second;
+    }
+  }
+  if (c[0] && c[1]){
+    newTris.push_back(new MTriangle (c[0],tri->getVertex(1),c[1]));
+    newTris.push_back(new MTriangle (tri->getVertex(0),c[0],tri->getVertex(2)));
+    newTris.push_back(new MTriangle (tri->getVertex(2),c[0],c[1]));
+    newCut.insert(MEdge(c[0],c[1]));
+  }
+  else if (c[0] && c[2]){
+    newTris.push_back(new MTriangle (tri->getVertex(0),c[0],c[2]));
+    newTris.push_back(new MTriangle (c[0],tri->getVertex(1),tri->getVertex(2)));
+    newTris.push_back(new MTriangle (tri->getVertex(2),c[2],c[0]));
+    newCut.insert(MEdge(c[0],c[2]));
+  }
+  else if (c[1] && c[2]){
+    newTris.push_back(new MTriangle (tri->getVertex(2),c[2],c[1]));
+    newTris.push_back(new MTriangle (tri->getVertex(0),tri->getVertex(1),c[2]));
+    newTris.push_back(new MTriangle (c[2],tri->getVertex(1),c[1]));
+    newCut.insert(MEdge(c[1],c[2]));
+  }
+  else if (c[0]){
+    newTris.push_back(new MTriangle (tri->getVertex(0),c[0],tri->getVertex(2)));
+    newTris.push_back(new MTriangle (tri->getVertex(2),c[0],tri->getVertex(1)));
+    if (cutVertices.find (tri->getVertex(0)) != cutVertices.end()){
+      newCut.insert(MEdge(c[0],tri->getVertex(0)));
+    }
+    else if (cutVertices.find (tri->getVertex(1)) != cutVertices.end()) {
+      newCut.insert(MEdge(c[0],tri->getVertex(1)));
+    }
+    else{
+      newCut.insert(MEdge(c[0],tri->getVertex(2)));
+    }
+  }
+  else if (c[1]){
+    newTris.push_back(new MTriangle (tri->getVertex(1),c[1],tri->getVertex(0)));
+    newTris.push_back(new MTriangle (tri->getVertex(0),c[1],tri->getVertex(2)));
+    if (cutVertices.find (tri->getVertex(0)) != cutVertices.end()){
+      newCut.insert(MEdge(c[1],tri->getVertex(0)));
+    }
+    else if (cutVertices.find (tri->getVertex(1)) != cutVertices.end()) {
+      newCut.insert(MEdge(c[1],tri->getVertex(1)));
+    }
+    else{
+      newCut.insert(MEdge(c[1],tri->getVertex(2)));
+    }
+  }
+  else if (c[2]){
+    newTris.push_back(new MTriangle (tri->getVertex(0),tri->getVertex(1), c[2]));
+    newTris.push_back(new MTriangle (tri->getVertex(1),tri->getVertex(2), c[2]));
+    if (cutVertices.find (tri->getVertex(0)) != cutVertices.end()){
+      newCut.insert(MEdge(c[2],tri->getVertex(0)));
+    }
+    else if (cutVertices.find (tri->getVertex(1)) != cutVertices.end()) {
+      newCut.insert(MEdge(c[2],tri->getVertex(1)));
+    }
+    else{
+      newCut.insert(MEdge(c[2],tri->getVertex(2)));
+    }
+  }
+  else {
+    newTris.push_back(tri);
+    if (cutVertices.find (tri->getVertex(0)) != cutVertices.end() &&
+	cutVertices.find (tri->getVertex(1)) != cutVertices.end())
+      newCut.insert(MEdge(tri->getVertex(0),tri->getVertex(1)));
+    if (cutVertices.find (tri->getVertex(0)) != cutVertices.end() &&
+	cutVertices.find (tri->getVertex(1)) != cutVertices.end())
+      newCut.insert(MEdge(tri->getVertex(0),tri->getVertex(2)));
+    if (cutVertices.find (tri->getVertex(2)) != cutVertices.end() &&
+	cutVertices.find (tri->getVertex(1)) != cutVertices.end())
+      newCut.insert(MEdge(tri->getVertex(2),tri->getVertex(1)));
+  }
+
+}
+
 Centerline::Centerline(std::string fileName): kdtree(0), nodes(0){
 
   index = new ANNidx[1];
@@ -186,7 +314,17 @@ void Centerline::importFile(std::string fileName){
   std::vector<GEntity*> entities ;
   current->getEntities(entities) ; 
   for(unsigned int i = 0; i < entities.size(); i++){
-    if(entities[i]->dim() == 2) surfaces.push_back((GFace*)entities[i]);
+    if(entities[i]->dim() != 2) continue; 
+    for(int j = 0; j < entities[i]->getNumMeshElements(); j++){ 
+      MElement *e = entities[i]->getMeshElement(j);
+      if (e->getType() != TYPE_TRI){
+  	Msg::Error("Centerline split only implemented for tri meshes so far ..."); 
+  	exit(1);
+      }
+      else{
+	triangles.push_back((MTriangle*)e);
+      }	
+    }
   }
   entities.clear();
 
@@ -300,7 +438,7 @@ void Centerline::createBranches(int maxN){
 	else itl++;
       }
       if (vB == vE) { Msg::Error("Begin and end points branch are the same \n");break;}
-      orderMLines(myLines); 
+      orderMLines(myLines, vB, vE); 
       std::vector<Branch> children;
       Branch newBranch ={ tag++, myLines, computeLength(myLines), vB, vE, children, 1.e6, 0.0};
       //printf("branch = %d VB = %d VE %d \n", myLines.size(), vB->getNum(), vE->getNum());
@@ -342,48 +480,83 @@ void Centerline::createBranches(int maxN){
 void Centerline::distanceToLines(){
   Msg::Info("Centerline: computing distance to lines ");
 
-  //EMI TO DO DO THIS WITH ANN !
-  //MUCH FASTER
+  //SOLUTION 1: REVERSE ANN TREE (SURFACE POINTS IN TREE)
+  ANNkd_tree *kdtreeR; 
+  ANNpointArray nodesR;
+  ANNidxArray indexR = new ANNidx[1];
+  ANNdistArray distR = new ANNdist[1];
 
-  std::vector<SPoint3> pts;
-  std::set<SPoint3> pts_set;
-  std::vector<double> distances;
-  std::vector<double> distancesE;
-  std::vector<SPoint3> closePts;
+  std::set<MVertex*> allVS; 
+  for(int j = 0; j < triangles.size(); j++)
+    for(int k = 0; k<3; k++) allVS.insert(triangles[j]->getVertex(k));
   
-  for(unsigned int i = 0; i < surfaces.size(); i++){ 
-    for(int j = 0; j < surfaces[i]->getNumMeshElements(); j++){ 
-      MElement *e = surfaces[i]->getMeshElement(j);
-      for (int k = 0; k < e->getNumVertices(); k++){
-	MVertex *v = e->getVertex(k);
-	pts_set.insert(SPoint3(v->x(), v->y(),v->z()));
-      }
-    }
-  }
-  std::set<SPoint3>::iterator its =  pts_set.begin();
-  for (; its != pts_set.end(); its++){
-    pts.push_back(*its);
+  int nbSNodes = allVS.size();
+  nodesR = annAllocPts(nbSNodes, 3);
+  int ind = 0;
+  std::set<MVertex*>::iterator itp = allVS.begin();
+  while (itp != allVS.end()){
+    MVertex *v = *itp;
+    nodesR[ind][0] = v->x(); 
+    nodesR[ind][1] = v->y(); 
+    nodesR[ind][2] = v->z(); 
+    itp++; ind++;
   }
-  if (pts.size() == 0) {Msg::Error("Centerline needs a GModel with a surface \n"); exit(1);}
+  kdtreeR = new ANNkd_tree(nodesR, nbSNodes, 3);
 
   for(unsigned int i = 0; i < lines.size(); i++){
-    std::vector<double> iDistances;
-    std::vector<SPoint3> iClosePts;
-    std::vector<double> iDistancesE;
     MLine *l = lines[i];
     MVertex *v1 = l->getVertex(0);
     MVertex *v2 = l->getVertex(1);
-    SPoint3 p1(v1->x(), v1->y(), v1->z());
-    SPoint3 p2(v2->x(), v2->y(), v2->z());
-    signedDistancesPointsLine(iDistances, iClosePts, pts, p1,p2);
-
-    double minRad = std::abs(iDistances[0]);
-    for (unsigned int kk = 1; kk< pts.size(); kk++) {
-      if (std::abs(iDistances[kk]) < minRad)
-	minRad = std::abs(iDistances[kk]);
-    }
+    double midp[3] = {0.5*(v1->x()+v2->x()), 0.5*(v1->y()+v1->y()),0.5*(v1->z()+v2->z())};
+    kdtreeR->annkSearch(midp, 1, indexR, distR);
+    double minRad = sqrt(distR[0]);
     radiusl.insert(std::make_pair(lines[i], minRad)); 
   }
+  delete kdtreeR;
+  annDeallocPts(nodesR);
+  delete[]indexR;
+  delete[]distR;
+  
+  // //SOLUTION 2: DISTANCE TO LINES (QUITE SLOW)
+  // std::vector<SPoint3> pts;
+  // std::set<SPoint3> pts_set;
+  // std::vector<double> distances;
+  // std::vector<double> distancesE;
+  // std::vector<SPoint3> closePts;
+
+  // for(unsigned int i = 0; i < surfaces.size(); i++){ 
+  //   for(int j = 0; j < surfaces[i]->getNumMeshElements(); j++){ 
+  //     MElement *e = surfaces[i]->getMeshElement(j);
+  //     for (int k = 0; k < e->getNumVertices(); k++){
+  // 	MVertex *v = e->getVertex(k);
+  // 	pts_set.insert(SPoint3(v->x(), v->y(),v->z()));
+  //     }
+  //   }
+  // }
+  // std::set<SPoint3>::iterator its =  pts_set.begin();
+  // for (; its != pts_set.end(); its++){
+  //   pts.push_back(*its);
+  // } 
+  // if (pts.size() == 0) {Msg::Error("Centerline needs a GModel with a surface \n"); exit(1);}
+ 
+  // for(unsigned int i = 0; i < lines.size(); i++){
+  //   std::vector<double> iDistances;
+  //   std::vector<SPoint3> iClosePts;
+  //   std::vector<double> iDistancesE;
+  //   MLine *l = lines[i];
+  //   MVertex *v1 = l->getVertex(0);
+  //   MVertex *v2 = l->getVertex(1);
+  //   SPoint3 p1(v1->x(), v1->y(), v1->z());
+  //   SPoint3 p2(v2->x(), v2->y(), v2->z());
+  //   signedDistancesPointsLine(iDistances, iClosePts, pts, p1,p2);
+
+  //   double minRad = std::abs(iDistances[0]);
+  //   for (unsigned int kk = 1; kk< pts.size(); kk++) {
+  //     if (std::abs(iDistances[kk]) < minRad)
+  // 	minRad = std::abs(iDistances[kk]);
+  //   }
+  //   radiusl.insert(std::make_pair(lines[i], minRad)); 
+  // }
      
 }
 void Centerline::computeRadii(){
@@ -448,96 +621,178 @@ void Centerline::buildKdTree(){
 
 void Centerline::splitMesh(){
 
-  Msg::Info("Splitting surface mesh with centerline field ");
+  Msg::Info("Splitting surface mesh (%d tris) with centerline field ", triangles.size());
 
- for(unsigned int i = 0; i < edges.size(); ++i){
+  for(unsigned int i = 0; i < edges.size(); i++){
+    printf("*** branch =%d \n", i);
     std::vector<MLine*> lines = edges[i].lines;
-   MVertex *v1 = edges[i].vE;
+    double L = edges[i].length;
+    double R = edges[i].maxRad;
+    if(edges[i].children.size() > 0.0 && L/R > 1.5){
+      MVertex *v1 = lines[lines.size()-1]->getVertex(1);//end vertex
       MVertex *v2;
-   SPoint3 pt(v1->x(), v1->y(), v1->z());
-   if (lines[0]->getVertex(0) == v1) v2 = lines[0]->getVertex(1);
-   else if (lines[0]->getVertex(1) == v1) v2 = lines[0]->getVertex(0);
-   else if (lines.back()->getVertex(0) == v1) v2 = lines.back()->getVertex(1);
-   else if (lines.back()->getVertex(1) == v1) v2 = lines.back()->getVertex(0);
+      if(L/R < 2.0) v2 = lines[0]->getVertex(0);
+      else if (lines.size() > 4) v2 = lines[lines.size()-4]->getVertex(0);
+      else v2 = lines[lines.size()-1]->getVertex(0);
+      SVector3 pt(v1->x(), v1->y(), v1->z());
       SVector3 dir(v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
-   if(edges[i].children.size() > 0.0){
-     cutByDisk(pt, dir, 1.2*edges[i].maxRad);
-   }
- }
-
-}
-
-void Centerline::cutByDisk(SPoint3 pt, SVector3 norm, double rad){
-  
-  printf("Cutting disk centered at pt  %g %g %g dir %g %g %g \n", pt.x(), pt.y(), pt.z(), norm.x(), norm.y(), norm.z());
-  double a = norm.x();
-  double b = norm.y();
-  double c = norm.z();
-  double d = -a * pt[0] - b * pt[1] - c * pt[2];
-  printf("Plane = %g %g %g %g \n", a, b, c, d);
-  //ls = a * x + b * y + c * z + d;
-
-  //variables for using the cutting tools
-  std::vector<DI_Line *> lines;
-  std::vector<DI_Triangle *> triangles;
-  std::vector<DI_Quad *> quads;
-  DI_Point::Container cp;//cut points
-  std::vector<DI_IntegrationPoint *> ipV;//integration points vol
-  std::vector<DI_IntegrationPoint *> ipS;//integration points surf
-  bool isCut = false;
-  int integOrder = 1;
-  int recur = 0;
-  gLevelset *GLS = new gLevelsetPlane(0., 0., 0., 0.);
-  std::vector<gLevelset *> RPN;
-  RPN.push_back(GLS);
-
-  //loop over all surface mesh elements
-  for(unsigned int i = 0; i < surfaces.size(); i++){ 
-    for(int j = 0; j < surfaces[i]->getNumMeshElements(); j++){ 
-      MElement *e = surfaces[i]->getMeshElement(j);
-      if (e->getNumVertices() != 3){
-  	Msg::Error("Centerline split only implemented for tri meshes so far ..."); 
-  	exit(1);
+      cutByDisk(pt, dir, edges[i].maxRad);
+    }
+ }
+
+  FILE * f2 = fopen("myCUT.pos","w");
+  fprintf(f2, "View \"\"{\n");
+  std::set<MEdge,Less_Edge>::iterator itp =  theCut.begin();
+   while (itp != theCut.end()){
+     MEdge l = *itp;
+     fprintf(f2, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n",
+              l.getVertex(0)->x(), l.getVertex(0)->y(), l.getVertex(0)->z(),
+              l.getVertex(1)->x(), l.getVertex(1)->y(), l.getVertex(1)->z(),
+              1.0,1.0);
+     itp++; 
+  }
+  //  for (int i= 0; i< triangles.size(); i++){
+  //    MTriangle *t = triangles[i];
+  //    fprintf(f2, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+  //             t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
+  //             t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
+  //             t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
+  //             1.0,1.0,1.0);
+  // }
+  fprintf(f2,"};\n");
+  fclose(f2);
+
+  Msg::Info("Splitting mesh by centerlines done ");
+
+}
+
+void Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
+   
+  double a = NORM.x();
+  double b = NORM.y();
+  double c = NORM.z();
+  double d = -a * PT.x() - b * PT.y() - c * PT.z();
+  printf("cut disk (R=%g)= %g %g %g %g \n", maxRad, a, b, c, d);
+ 
+  //EMI STUFF
+  const double EPS = 0.001;
+  std::set<MEdge,Less_Edge> allEdges;
+  for(unsigned int i = 0; i < triangles.size(); i++){ 
+    for ( unsigned int j= 0; j <  3; j++)
+      allEdges.insert(triangles[i]->getEdge(j)); 
+  }
+  bool closedCut = false;
+  int step = 0;
+  while (!closedCut && step < 10){
+    double rad = 1.3*maxRad+0.15*step*maxRad;
+    //printf("radius(%d) =%g \n", step, rad);
+    std::map<MEdge,MVertex*,Less_Edge> cutEdges;
+    std::set<MVertex*> cutVertices;
+    std::vector<MTriangle*> newTris; 
+    std::set<MEdge,Less_Edge> newCut;
+    cutEdges.clear();
+    cutVertices.clear();
+    newTris.clear();
+    newCut.clear();
+    for (std::set<MEdge,Less_Edge>::iterator it = allEdges.begin();
+	 it != allEdges.end() ; ++it){
+      MEdge me = *it;
+      SVector3 P1(me.getVertex(0)->x(),me.getVertex(0)->y(), me.getVertex(0)->z());
+      SVector3 P2(me.getVertex(1)->x(),me.getVertex(1)->y(), me.getVertex(1)->z());
+      double V1 = a * P1.x() + b * P1.y() + c * P1.z() + d;
+      double V2 = a * P2.x() + b * P2.y() + c * P2.z() + d;
+      bool inters = (V1*V2<0.0) ? true: false;
+      double rdist = -V1/(V2-V1);
+      if (inters && rdist > EPS && rdist < 1.-EPS){
+	SVector3 PZ = P1+rdist*(P2-P1);
+	MVertex *newv = new MVertex (PZ.x(), PZ.y(), PZ.z());
+	if ( norm(PZ-PT) < rad ) cutEdges.insert(std::make_pair(*it,newv));
+      }
+      else if (inters && rdist <= EPS && norm(P1-PT) < rad ) 
+	cutVertices.insert(me.getVertex(0));
+      else if (inters && rdist >= 1.-EPS && norm(P2-PT) < rad) 
+	cutVertices.insert(me.getVertex(1));
+    }
+    for(unsigned int i = 0; i < triangles.size(); i++){ 
+      cutTriangle(triangles[i], cutEdges,cutVertices, newTris, newCut);
+    }
+    if (isClosed(newCut)) {
+      printf("closed cut \n");
+      triangles = newTris;
+      theCut.insert(newCut.begin(),newCut.end());
+      //if (step > 1) printf("YOUPIIIIIIIIIIIIIIIIIIIIIII \n");
+      break;
     }
-      int nbV = e->getNumVertices();
-
-      double **nodeLs= new double*[nbV];
-      DI_Triangle T(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
-  		    e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
-  		    e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z());
-      SPoint3 cg(0.333*(e->getVertex(0)->x()+e->getVertex(1)->x()+e->getVertex(2)->x()), 
-  		 0.333*(e->getVertex(0)->y()+e->getVertex(1)->y()+e->getVertex(2)->y()),
-  		 0.333*(e->getVertex(0)->z()+e->getVertex(1)->z()+e->getVertex(2)->z()));
-      T.setPolynomialOrder(1);
-      double x0 = e->getVertex(0)->x(), y0= e->getVertex(0)->y(), z0=e->getVertex(0)->z();
-      double val0 = a * x0 + b * y0 + c * z0 + d;
-      double sign0 = val0;
-      bool zeroElem = false;
-      nodeLs[0] = new double[1];//one ls value
-      nodeLs[0][0]= val0;
-      for(int i=1;i<nbV;i++){
-	double x = e->getVertex(i)->x(), y= e->getVertex(i)->y(), z=e->getVertex(i)->z();
-  	double val = a * x + b * y + c * z + d;
-    	double sign = sign0*val;
-  	if (sign < 0.0 && !zeroElem) zeroElem = true;
-  	nodeLs[i] = new double[1];//one ls value
-	nodeLs[i][0]= val;
-      }
-      if (zeroElem){
-  	double radius = sqrt((cg.x()-pt.x())*(cg.x()-pt.x())+(cg.y()-pt.y())*(cg.y()-pt.y())+(cg.z()-pt.z())*(cg.z()-pt.z()));
-  	if (radius <  rad){
-      	  isCut = T.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
-      			quads, triangles, lines, recur, nodeLs);
-  	}
-      }
-      else triangles.push_back(&T);
-
-      for(int i=0;i<nbV;i++) delete nodeLs[i];
-      delete []nodeLs;
-    }
-  }
-
-  printf("split mesh done \n");
+    else {
+      if (step ==9) printf("no closed cut %d \n", newCut.size());
+      step++;
+      // // triangles = newTris;
+      // // theCut.insert(newCut.begin(),newCut.end());
+      // char name[256];
+      // sprintf(name, "myCUT-%d.pos", step);
+      // FILE * f2 = fopen(name,"w");
+      // fprintf(f2, "View \"\"{\n");
+      // std::set<MEdge,Less_Edge>::iterator itp =  newCut.begin();
+      // while (itp != newCut.end()){
+      // 	MEdge l = *itp;
+      // 	fprintf(f2, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n",
+      // 		  l.getVertex(0)->x(), l.getVertex(0)->y(), l.getVertex(0)->z(),
+      // 		  l.getVertex(1)->x(), l.getVertex(1)->y(), l.getVertex(1)->z(),
+      // 		  1.0,1.0);
+      // 	  itp++; 
+      // 	}
+      // 	fprintf(f2,"};\n");
+      // 	fclose(f2);
+    }
+  }
+
+  // //GAETAN STUFF
+  // //variables for using the cutting tools
+  // DI_Point::Container cp;//cut points
+  // std::vector<DI_IntegrationPoint *> ipV;//integration points vol
+  // std::vector<DI_IntegrationPoint *> ipS;//integration points surf
+  // bool isCut = false;
+  // gLevelset *GLS = new gLevelsetPlane(0., 0., 0., 0.);
+  // std::vector<gLevelset *> RPN;
+  // RPN.push_back(GLS);
+
+  // //loop over all surface mesh elements
+  // for(unsigned int i = 0; i < triangles.size(); i++){ 
+  //   MTriangle *e = triangles[i];
+  //   double **nodeLs= new double*[3];
+  //   DI_Triangle *T = new DI_Triangle(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
+  // 				     e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
+  // 				     e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z());
+  //   SPoint3 cg(0.333*(e->getVertex(0)->x()+e->getVertex(1)->x()+e->getVertex(2)->x()), 
+  // 	       0.333*(e->getVertex(0)->y()+e->getVertex(1)->y()+e->getVertex(2)->y()),
+  // 	       0.333*(e->getVertex(0)->z()+e->getVertex(1)->z()+e->getVertex(2)->z()));
+  //   double x0 = e->getVertex(0)->x(), y0= e->getVertex(0)->y(), z0=e->getVertex(0)->z();
+  //   double val0 = a * x0 + b * y0 + c * z0 + d;
+  //   double sign0 = val0;
+  //   bool zeroElem = false;
+  //   nodeLs[0] = new double[1];//one ls value
+  //   nodeLs[0][0]= val0;
+  //   for(int i=1;i<3;i++){
+  //     double x = e->getVertex(i)->x(), y= e->getVertex(i)->y(), z=e->getVertex(i)->z();
+  //     double val = a * x + b * y + c * z + d;
+  //     double sign = sign0*val;
+  //     if (sign < 0.0 && !zeroElem) zeroElem = true;
+  //     nodeLs[i] = new double[1];//one ls value
+  //     nodeLs[i][0]= val;
+  //   }
+  //   if (zeroElem){
+  //     double radius = sqrt((cg.x()-PT.x())*(cg.x()-PT.x())+(cg.y()-PT.y())*(cg.y()-PT.y())+(cg.z()-PT.z())*(cg.z()-PT.z()));
+  //     if (radius <  rad){
+  // 	isCut = T->cut(RPN, ipV, ipS, cp, 1, 1, 1,
+  // 		       cutQUAD, cutTRI, cutLIN, 0, nodeLs);
+  //     }
+  //     else cutTRI.push_back(T);
+  //   }
+  //   else  cutTRI.push_back(T);
+  //   for(int i=0;i<3;i++) delete nodeLs[i];
+  //   delete []nodeLs;
+  // }
+
   return;
 
 }

Mesh/CenterlineField.h

@@ -14,11 +14,13 @@
 #include <set>
 #include <string>
 #include "Field.h"
+#include "MEdge.h"
 class GModel;
 class GFace;
 class MLine;
 class MVertex;
 class GEntity;
+class MTriangle;   
 
 #if defined(HAVE_ANN)
 #include <ANN/ANN.h>
@@ -66,7 +68,9 @@ class Centerline : public Field{
   std::map<MLine*,int> colorl;
 
   //the tubular surface mesh
-  std::vector<GFace*> surfaces; 
+  std::vector<MTriangle*> triangles;
+  //the lines cut of the tubular mesh by planes
+  std::set<MEdge,Less_Edge> theCut;
 
  public:
   Centerline(std::string fileName);
@@ -116,7 +120,7 @@ class Centerline : public Field{
 
   // Cut the tubular structure with a disk
   // perpendicular to the tubular structure
-  void cutByDisk(SPoint3 pt, SVector3 dir, double rad);
+  void cutByDisk(SVector3 &pt, SVector3 &dir, double &maxRad);
 
   //Print for debugging
   void printSplit() const;

contrib/DiscreteIntegration/Integration3D.h

@@ -324,13 +324,13 @@ class DI_Element
   void getCuttingPoints (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
                          std::vector<DI_CuttingPoint*> &cp) const;
   // return the ith point
-  virtual DI_Point* pt(int i) const {return (i < nbVert()) ? &(pts_[i]) : &(mid_[i - nbVert()]);}
+  DI_Point* pt(int i) const {return (i < nbVert()) ? &(pts_[i]) : &(mid_[i - nbVert()]);}
   // return the ith middle point
   inline DI_Point* mid(int i) const {return mid_ ? &(mid_[i]) : NULL;}
   // return the coordinates of the ith point
-  virtual double x(int i) const {return pt(i)->x();}
-  virtual double y(int i) const {return pt(i)->y();}
-  virtual double z(int i) const {return pt(i)->z();}
+  double x(int i) const {return pt(i)->x(); }
+  double y(int i) const {return pt(i)->y();}
+  double z(int i) const {return pt(i)->z();}
   // return the last levelset value of the ith point
   virtual double ls(int i) const {return pt(i)->ls();}
   // return the jth levelset value of the ith point