From da1194ac61d099bb9faaa786d8d53aa2b6435c72 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Wed, 29 Feb 2012 12:45:40 +0000
Subject: [PATCH] Centerline with closeVolume action
---
Geo/GModelFactory.cpp | 120 +++++++++++++++++++++++++--------------
Geo/GModelIO_Geo.cpp | 2 +
Mesh/CenterlineField.cpp | 101 ++++++++++++++++----------------
Mesh/CenterlineField.h | 7 ++-
Mesh/meshGEdge.cpp | 3 +-
benchmarks/2d/square.geo | 6 +-
6 files changed, 141 insertions(+), 98 deletions(-)
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 4780eb6743..22903d89d5 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -26,7 +26,7 @@ GVertex *GeoFactory::addVertex(GModel *gm, double x, double y, double z, double
if(lc == 0.) lc = MAX_LC; // no mesh size given at the point
Vertex *p;
p = Create_Vertex(num, x, y, z, lc, 1.0);
- Tree_Add(GModel::current()->getGEOInternals()->Points, &p);
+ Tree_Add(gm->getGEOInternals()->Points, &p);
p->Typ = MSH_POINT;
p->Num = num;
@@ -47,7 +47,7 @@ GEdge *GeoFactory::addLine(GModel *gm, GVertex *start, GVertex *end)
Curve *c = Create_Curve(num, MSH_SEGM_LINE, 1, iList, NULL,
-1, -1, 0., 1.);
- Tree_Add(GModel::current()->getGEOInternals()->Curves, &c);
+ Tree_Add(gm->getGEOInternals()->Curves, &c);
CreateReversedCurve(c);
List_Delete(iList);
c->Typ = MSH_SEGM_LINE;
@@ -61,49 +61,83 @@ GEdge *GeoFactory::addLine(GModel *gm, GVertex *start, GVertex *end)
GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> > edges)
{
- //create line loops
- std::vector<EdgeLoop *> vecLoops;
+
+ //create line loops
int nLoops = edges.size();
+ std::vector<EdgeLoop *> vecLoops;
for (int i=0; i< nLoops; i++){
- int numl = gm->getMaxElementaryNumber(1) + i;
- while (FindEdgeLoop(numl)){
- numl++;
- if (!FindEdgeLoop(numl)) break;
- }
- int nl=(int)edges[i].size();
- List_T *iListl = List_Create(nl, nl, sizeof(int));
- for(int j = 0; j < nl; j++){
- int numEdge = edges[i][j]->tag();
- List_Add(iListl, &numEdge);
- }
- int type=ENT_LINE;
-
- select_contour(type, edges[i][0]->tag(), iListl);
+ int ne=(int)edges[i].size();
+ List_T *temp = List_Create(ne, ne, sizeof(int));
+ for(int j = 0; j < ne; j++){
+ GEdge *ge = edges[i][j];
+ int numEdge = ge->tag();
+ //create curve
+ Curve *c = FindCurve(numEdge);
+ if(!c){
+ GVertex *gvb = ge->getBeginVertex();
+ GVertex *gve = ge->getEndVertex();
+ Vertex *vertb = FindPoint((int)fabs(gvb->tag()));
+ Vertex *verte = FindPoint((int)fabs(gve->tag()));
+ if (!vertb){
+ vertb = Create_Vertex(gvb->tag(), gvb->x(), gvb->y(), gvb->z(),
+ gvb->prescribedMeshSizeAtVertex(), 1.0);
+ Tree_Add(gm->getGEOInternals()->Points, &vertb);
+ }
+ if (!verte){
+ verte = Create_Vertex(gve->tag(), gve->x(), gve->y(), gve->z(),
+ gve->prescribedMeshSizeAtVertex(), 1.0);
+ Tree_Add(gm->getGEOInternals()->Points, &verte);
+ }
+
+ List_T *temp = List_Create(2, 2, sizeof(int));
+ List_Add(temp, &vertb);
+ List_Add(temp, &verte);
+
+ if (ge->geomType() == GEntity::Line){
+ printf("creating line \n");
+ c = Create_Curve(numEdge, MSH_SEGM_LINE, 1, temp, NULL, -1, -1, 0., 1.);
+ }
+ else if (ge->geomType() == GEntity::DiscreteCurve){
+ printf("creating discrete line \n");
+ c = Create_Curve(numEdge, MSH_SEGM_DISCRETE, 1, NULL, NULL, -1, -1, 0., 1.);
+ }
+ else {
+ printf("type not implemented \n"); exit(1);
+ }
+ c->Control_Points = List_Create(2, 1, sizeof(Vertex *));
+ List_Add(c->Control_Points, &vertb);
+ List_Add(c->Control_Points, &verte);
+ c->beg = vertb;
+ c->end = verte;
+ End_Curve(c);
+
+ Tree_Add(gm->getGEOInternals()->Curves, &c);
+ CreateReversedCurve(c);
+ List_Delete(temp);
+ }
+ List_Add(temp, &numEdge);
+ }
- sortEdgesInLoop(numl, iListl);
- EdgeLoop *l = Create_EdgeLoop(numl, iListl);
+ int num = gm->getMaxElementaryNumber(2)+1+i;
+ if(FindEdgeLoop(num)) num++;
+ sortEdgesInLoop(num, temp);
+ EdgeLoop *l = Create_EdgeLoop(num, temp);
vecLoops.push_back(l);
- Tree_Add(GModel::current()->getGEOInternals()->EdgeLoops, &l);
- l->Num = numl;
-
-
- List_Delete(iListl);
- }
-
- //create plane surfaces
- int numf = gm->getMaxElementaryNumber(2) + 1;
+ Tree_Add(gm->getGEOInternals()->EdgeLoops, &l);
+ List_Delete(temp);
+ }
+
+ int numf = gm->getMaxElementaryNumber(2)+1;
Surface *s = Create_Surface(numf, MSH_SURF_PLAN);
- List_T *iList = List_Create(nLoops, nLoops, sizeof(int));
- for (unsigned int i=0; i< vecLoops.size(); i++){
+ List_T *temp = List_Create(nLoops, nLoops, sizeof(int));
+ for (unsigned int i=0; i< nLoops; i++){
int numl = vecLoops[i]->Num;
- List_Add(iList, &numl);
+ List_Add(temp, &numl);
}
- setSurfaceGeneratrices(s, iList);
+ setSurfaceGeneratrices(s, temp);
+ List_Delete(temp);
End_Surface(s);
- Tree_Add(GModel::current()->getGEOInternals()->Surfaces, &s);
- s->Typ= MSH_SURF_PLAN;
- s->Num = numf;
- List_Delete(iList);
+ Tree_Add(gm->getGEOInternals()->Surfaces, &s);
//gmsh surface
GFace *gf = new gmshFace(gm,s);
@@ -131,7 +165,7 @@ GRegion* GeoFactory::addVolume (GModel *gm, std::vector<std::vector<GFace *> > f
}
SurfaceLoop *l = Create_SurfaceLoop(numfl, iListl);
vecLoops.push_back(l);
- Tree_Add(GModel::current()->getGEOInternals()->SurfaceLoops, &l);
+ Tree_Add(gm->getGEOInternals()->SurfaceLoops, &l);
List_Delete(iListl);
}
@@ -145,7 +179,7 @@ GRegion* GeoFactory::addVolume (GModel *gm, std::vector<std::vector<GFace *> > f
}
setVolumeSurfaces(v, iList);
List_Delete(iList);
- Tree_Add(GModel::current()->getGEOInternals()->Volumes, &v);
+ Tree_Add(gm->getGEOInternals()->Volumes, &v);
v->Typ = MSH_VOLUME;
v->Num = numv;
@@ -169,7 +203,7 @@ GEdge* GeoFactory::addCircleArc(GModel *gm,GVertex *begin, GVertex *center, GVer
Curve *c = Create_Curve(num, MSH_SEGM_CIRC, 1, iList, NULL,
-1, -1, 0., 1.);
- Tree_Add(GModel::current()->getGEOInternals()->Curves, &c);
+ Tree_Add(gm->getGEOInternals()->Curves, &c);
CreateReversedCurve(c);
List_Delete(iList);
c->Typ = MSH_SEGM_CIRC;
@@ -200,12 +234,12 @@ std::vector<GFace *> GeoFactory::addRuledFaces(GModel *gm,
List_Add(iListl, &numEdge);
}
int type=ENT_LINE;
- if(select_contour(type, edges[0][0]->tag(), iListl))
+ if(select_contour(type, edges[i][0]->tag(), iListl))
{
sortEdgesInLoop(numl, iListl);
EdgeLoop *l = Create_EdgeLoop(numl, iListl);
vecLoops.push_back(l);
- Tree_Add(GModel::current()->getGEOInternals()->EdgeLoops, &l);
+ Tree_Add(gm->getGEOInternals()->EdgeLoops, &l);
l->Num = numl;
}
@@ -222,7 +256,7 @@ std::vector<GFace *> GeoFactory::addRuledFaces(GModel *gm,
}
setSurfaceGeneratrices(s, iList);
End_Surface(s);
- Tree_Add(GModel::current()->getGEOInternals()->Surfaces, &s);
+ Tree_Add(gm->getGEOInternals()->Surfaces, &s);
s->Typ= MSH_SURF_TRIC;
s->Num = numf;
List_Delete(iList);
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index a0f921a258..49b21febce 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -118,6 +118,7 @@ int GModel::exportDiscreteGEOInternals()
int GModel::importGEOInternals()
{
+
if(Tree_Nbr(_geo_internals->Points)) {
List_T *points = Tree2List(_geo_internals->Points);
for(int i = 0; i < List_Nbr(points); i++){
@@ -163,6 +164,7 @@ int GModel::importGEOInternals()
if(!c->Visible) e->setVisibility(0);
if(c->Color.type) e->setColor(c->Color.mesh);
if(c->degenerated) {
+ printf("setting edge too small c-degenerated \n");
e->setTooSmall(true);
}
}
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index d3c827efea..6ee94d46bc 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -329,6 +329,7 @@ Centerline::Centerline(std::string fileName): kdtree(0), nodes(0){
update_needed = false;
is_cut = false;
+ is_closed = false;
}
Centerline::Centerline(): kdtree(0), nodes(0){
@@ -342,8 +343,10 @@ Centerline::Centerline(): kdtree(0), nodes(0){
options["FileName"] = new FieldOptionString (fileName, "File name for the centerlines", &update_needed);
options["nbPoints"] = new FieldOptionInt(nbPoints, "Number of mesh elements in a circle");
+ callbacks["closeVolume"] = new FieldCallbackGeneric<Centerline>(this, &Centerline::closeVolume, "Create In/Outlet planar faces \n");
callbacks["cutMesh"] = new FieldCallbackGeneric<Centerline>(this, &Centerline::cutMesh, "Cut the initial mesh in different mesh partitions using the centerlines \n");
is_cut = false;
+ is_closed = false;
}
@@ -490,7 +493,7 @@ void Centerline::createBranches(int maxN){
edges.push_back(newBranch) ;
}
}
- printf("in/outlets =%d branches =%d \n", (int)color_edges.size(), (int)edges.size());
+ printf("in/outlets =%d branches =%d \n", (int)color_edges.size()+1, (int)edges.size());
//create children
for(unsigned int i = 0; i < edges.size(); ++i) {
@@ -504,7 +507,7 @@ void Centerline::createBranches(int maxN){
}
//compute radius
- distanceToLines();
+ distanceToSurface();
computeRadii();
printSplit();
@@ -524,10 +527,10 @@ void Centerline::createBranches(int maxN){
}
-void Centerline::distanceToLines(){
- Msg::Info("Centerline: computing distance to lines ");
+void Centerline::distanceToSurface(){
+ Msg::Info("Centerline: computing distance to surface mesh ");
- //SOLUTION 1: REVERSE ANN TREE (SURFACE POINTS IN TREE)
+ //COMPUTE WITH REVERSE ANN TREE (SURFACE POINTS IN TREE)
ANNkd_tree *kdtreeR;
ANNpointArray nodesR;
ANNidxArray indexR = new ANNidx[1];
@@ -562,47 +565,6 @@ void Centerline::distanceToLines(){
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(){
@@ -774,6 +736,7 @@ void Centerline::cleanMesh(){
mySplitMesh->mesh_vertices.push_back(*it);
mySplitMesh->meshStatistics.status = GFace::DONE;
current->createTopologyFromMesh();
+
}
void Centerline::createFaces(){
@@ -830,6 +793,42 @@ void Centerline::createFaces(){
}
+void Centerline::createInOutFaces(){
+
+ //identify the boundary edges by looping over all discreteFaces
+ std::vector<GEdge*> boundEdges;
+ int nbFaces = current->getMaxElementaryNumber(2);
+ for (int i= 0; i< nbFaces; i++){
+ GFace *gf = current->getFaceByTag(i+1);
+ std::list<GEdge*> l_edges = gf->edges();
+ for(std::list<GEdge*>::iterator it = l_edges.begin(); it != l_edges.end(); it++){
+ std::vector<GEdge*>::iterator ite = std::find(boundEdges.begin(), boundEdges.end(), *it);
+ if (ite != boundEdges.end()) boundEdges.erase(ite);
+ else boundEdges.push_back(*it);
+ }
+ }
+ printf("boundEdges size =%d \n", boundEdges.size());
+
+ //create the inlet and outlet planar face
+ current->setFactory("Gmsh");
+ for (int i = 0; i< boundEdges.size(); i++){
+ std::vector<std::vector<GEdge *> > myEdges;
+ std::vector<GEdge *> myEdge;
+ myEdge.push_back(boundEdges[i]);
+ myEdges.push_back(myEdge);
+ GFace *newFace = current->addPlanarFace(myEdges);
+ }
+
+}
+
+void Centerline::closeVolume(){
+
+ is_closed = true;
+ //printf("calling closed volume \n");
+ //exit(1);
+
+}
+
void Centerline::cutMesh(){
is_cut = true;
@@ -861,11 +860,11 @@ void Centerline::cutMesh(){
for(unsigned int i = 0; i < edges.size(); i++){
std::vector<MLine*> lines = edges[i].lines;
double L = edges[i].length;
- double R = 0.5*(edges[i].minRad+edges[i].maxRad);
+ double R = edges[i].minRad; //0.5*(edges[i].minRad+edges[i].maxRad);
double AR = L/R;
- printf("*** branch =%d \n", i, AR);
+ printf("*** branch =%d AR=%g \n", i, AR);
if( AR > 4.5){
- int nbSplit = (int)ceil(AR/4.0);
+ int nbSplit = (int)ceil(AR/3.0);
double li = L/nbSplit;
double lc = 0.0;
for (int j= 0; j < lines.size(); j++){
@@ -898,6 +897,8 @@ void Centerline::cutMesh(){
//create discreteFaces
createFaces();
current->createTopologyFromFaces(discFaces);
+ current->exportDiscreteGEOInternals();
+ if (is_closed) createInOutFaces();
//write
Msg::Info("Writing splitted mesh 'myPARTS.msh'");
@@ -966,7 +967,7 @@ void Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
break;
}
else {
- if (step ==19) {printf("no closed cut %d \n", (int)newCut.size()); };
+ if (step == 9) {printf("no closed cut %d \n", (int)newCut.size()); };
step++;
// // triangles = newTris;
// // theCut.insert(newCut.begin(),newCut.end());
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
index cac7103cea..73b1947f3c 100644
--- a/Mesh/CenterlineField.h
+++ b/Mesh/CenterlineField.h
@@ -64,6 +64,7 @@ class Centerline : public Field{
double recombine;
int NF, NV, NE;
bool is_cut;
+ bool is_closed;
//all (unique) lines of centerlines
std::vector<MLine*> lines;
@@ -131,17 +132,21 @@ class Centerline : public Field{
void computeRadii();
//Computes for each MLine the minRadius
- void distanceToLines();
+ void distanceToSurface();
// Cut the mesh in different parts of small aspect ratio
void cutMesh();
+ //Create In and Outlet Planar Faces
+ void closeVolume();
+
// Cut the tubular structure with a disk
// perpendicular to the tubular structure
void cutByDisk(SVector3 &pt, SVector3 &dir, double &maxRad);
//create discrete faces
void createFaces();
+ void createInOutFaces();
void createSplitCompounds();
//Print for debugging
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 57df16408b..35b22c060b 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -344,8 +344,9 @@ void meshGEdge::operator() (GEdge *ge)
ge->setLength(length);
Points.clear();
- if(length < CTX::instance()->mesh.toleranceEdgeLength)
+ if(length < CTX::instance()->mesh.toleranceEdgeLength){
ge->setTooSmall(true);
+ }
// Integrate detJ/lc du
double a;
diff --git a/benchmarks/2d/square.geo b/benchmarks/2d/square.geo
index 366bae5b86..75d3f2c59a 100644
--- a/benchmarks/2d/square.geo
+++ b/benchmarks/2d/square.geo
@@ -25,10 +25,10 @@ Plane Surface(10) = {5};
//----------------------
-Compound Line(10)={1,2,3,4};
-Compound Surface(100)={10};
+//Compound Line(10)={1,2,3,4};
+//Compound Surface(100)={10};
-Line {1,2,3,4} In Surface{100};
+//Line {1,2,3,4} In Surface{100};
//Physical Surface(100)={10};
//Physical Line(200)={1,2,3,4};
--
GitLab