diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 41dac287f5deac19d81462afe135a4f73049e1b4..e5b892c86f0bb6b9738d9de0ab45b4f8da3a13fe 100644
--- a/Mesh/CenterlineField.cpp
+++ b/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());
@@ -308,7 +446,7 @@ void Centerline::createBranches(int maxN){
}
}
printf("in/outlets =%d branches =%d \n", color_edges.size(), edges.size());
-
+
//create children
for(unsigned int i = 0; i < edges.size(); ++i) {
MVertex *vE = edges[i].vE;
@@ -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
-
- 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()));
- }
- }
+ //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::set<MVertex*> allVS;
+ for(int j = 0; j < triangles.size(); j++)
+ for(int k = 0; k<3; k++) allVS.insert(triangles[j]->getVertex(k));
+
+ 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++;
}
- 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);}
-
+ 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 ");
-
- for(unsigned int i = 0; i < edges.size(); ++i){
- std::vector<MLine*> lines = edges[i].lines;
- MVertex *v1 = edges[i].vE;
- 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);
- 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);
- }
+ Msg::Info("Splitting surface mesh (%d tris) with centerline field ", triangles.size());
+
+ for(unsigned int i = 0; i < edges.size(); i++){
+ printf("*** branch =%d \n", i);
+ std::vector<MLine*> lines = edges[i].lines;
+ 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;
+ 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());
+ 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(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);
- }
- 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);
- }
+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 triangles.push_back(&T);
-
- for(int i=0;i<nbV;i++) delete nodeLs[i];
- delete []nodeLs;
+ 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;
+ }
+ 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);
}
}
- printf("split mesh done \n");
+ // //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;
}
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
index e773cde8e87771b0bf36e3a26300d5f2315ace30..e7d0450fe038160a906ecf1f7f6f50f032f5f9ac 100644
--- a/Mesh/CenterlineField.h
+++ b/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;
diff --git a/contrib/DiscreteIntegration/Integration3D.h b/contrib/DiscreteIntegration/Integration3D.h
index 3e5f22156cf2bce5505498a0ee5988c9f85039cd..4e7f1cf3aa5ffb3a5c4adc61f353fec5bf9a777e 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/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