Skip to content
Snippets Groups Projects
Commit 3f6beb30 authored by Emilie Marchandise's avatar Emilie Marchandise
Browse files

Centerline with close and create volume works fine, next step is to add some...

Centerline with close and create volume works fine, next step is to add some boundary layers in the volume
parent 84b0dc98
No related branches found
No related tags found
No related merge requests found
...@@ -106,7 +106,10 @@ void GFace::delFreeEdge(GEdge *e) ...@@ -106,7 +106,10 @@ void GFace::delFreeEdge(GEdge *e)
void GFace::deleteMesh() void GFace::deleteMesh()
{ {
for(unsigned int i = 0; i < mesh_vertices.size(); i++) delete mesh_vertices[i]; for(unsigned int i = 0; i < mesh_vertices.size(); i++) {
printf("delete mesh i=%d num =%d \n", i, mesh_vertices[i]->getNum());
delete mesh_vertices[i];
}
mesh_vertices.clear(); mesh_vertices.clear();
transfinite_vertices.clear(); transfinite_vertices.clear();
for(unsigned int i = 0; i < triangles.size(); i++) delete triangles[i]; for(unsigned int i = 0; i < triangles.size(); i++) delete triangles[i];
......
...@@ -71,7 +71,7 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> > ...@@ -71,7 +71,7 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
for(int j = 0; j < ne; j++){ for(int j = 0; j < ne; j++){
GEdge *ge = edges[i][j]; GEdge *ge = edges[i][j];
int numEdge = ge->tag(); int numEdge = ge->tag();
//create curve //create curve if it does not exist
Curve *c = FindCurve(numEdge); Curve *c = FindCurve(numEdge);
if(!c){ if(!c){
GVertex *gvb = ge->getBeginVertex(); GVertex *gvb = ge->getBeginVertex();
...@@ -119,7 +119,10 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> > ...@@ -119,7 +119,10 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
} }
int num = gm->getMaxElementaryNumber(2) + 1+i; int num = gm->getMaxElementaryNumber(2) + 1+i;
if(FindEdgeLoop(num)) num++; while (FindSurfaceLoop(num)){
num++;
if (!FindSurfaceLoop(num)) break;
}
sortEdgesInLoop(num, temp); sortEdgesInLoop(num, temp);
EdgeLoop *l = Create_EdgeLoop(num, temp); EdgeLoop *l = Create_EdgeLoop(num, temp);
vecLoops.push_back(l); vecLoops.push_back(l);
...@@ -127,6 +130,7 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> > ...@@ -127,6 +130,7 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
List_Delete(temp); List_Delete(temp);
} }
//create surface
int numf = gm->getMaxElementaryNumber(2)+1; int numf = gm->getMaxElementaryNumber(2)+1;
Surface *s = Create_Surface(numf, MSH_SURF_PLAN); Surface *s = Create_Surface(numf, MSH_SURF_PLAN);
List_T *temp = List_Create(nLoops, nLoops, sizeof(int)); List_T *temp = List_Create(nLoops, nLoops, sizeof(int));
...@@ -148,37 +152,39 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> > ...@@ -148,37 +152,39 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
GRegion* GeoFactory::addVolume (GModel *gm, std::vector<std::vector<GFace *> > faces) GRegion* GeoFactory::addVolume (GModel *gm, std::vector<std::vector<GFace *> > faces)
{ {
//surface loop printf("add volume \n");
std::vector<SurfaceLoop *> vecLoops;
//create surface loop
int nLoops = faces.size(); int nLoops = faces.size();
std::vector<SurfaceLoop *> vecLoops;
for (int i=0; i< nLoops; i++){ for (int i=0; i< nLoops; i++){
int nl=(int)faces[i].size();
List_T *temp = List_Create(nl, nl, sizeof(int));
for(int j = 0; j < nl; j++){
int numFace = faces[i][j]->tag();
List_Add(temp, &numFace);
}
int numfl = gm->getMaxElementaryNumber(2) + 1; int numfl = gm->getMaxElementaryNumber(2) + 1;
while (FindSurfaceLoop(numfl)){ while (FindSurfaceLoop(numfl)){
numfl++; numfl++;
if (!FindSurfaceLoop(numfl)) break; if (!FindSurfaceLoop(numfl)) break;
} }
int nl=(int)faces[i].size(); SurfaceLoop *l = Create_SurfaceLoop(numfl, temp);
List_T *iListl = List_Create(nl, nl, sizeof(int));
for(int j = 0; j < nl; j++){
int numFace = faces[i][j]->tag();
List_Add(iListl, &numFace);
}
SurfaceLoop *l = Create_SurfaceLoop(numfl, iListl);
vecLoops.push_back(l); vecLoops.push_back(l);
Tree_Add(gm->getGEOInternals()->SurfaceLoops, &l); Tree_Add(gm->getGEOInternals()->SurfaceLoops, &l);
List_Delete(iListl); List_Delete(temp);
} }
//volume //create volume
int numv = gm->getMaxElementaryNumber(3) + 1; int numv = gm->getMaxElementaryNumber(3) + 1;
Volume *v = Create_Volume(numv, MSH_VOLUME); Volume *v = Create_Volume(numv, MSH_VOLUME);
List_T *iList = List_Create(nLoops, nLoops, sizeof(int)); List_T *temp = List_Create(nLoops, nLoops, sizeof(int));
for (unsigned int i = 0; i < vecLoops.size(); i++){ for (unsigned int i = 0; i < vecLoops.size(); i++){
int numl = vecLoops[i]->Num; int numl = vecLoops[i]->Num;
List_Add(iList, &numl); List_Add(temp, &numl);
} }
setVolumeSurfaces(v, iList); setVolumeSurfaces(v, temp);
List_Delete(iList); List_Delete(temp);
Tree_Add(gm->getGEOInternals()->Volumes, &v); Tree_Add(gm->getGEOInternals()->Volumes, &v);
v->Typ = MSH_VOLUME; v->Typ = MSH_VOLUME;
v->Num = numv; v->Num = numv;
...@@ -228,22 +234,21 @@ std::vector<GFace *> GeoFactory::addRuledFaces(GModel *gm, ...@@ -228,22 +234,21 @@ std::vector<GFace *> GeoFactory::addRuledFaces(GModel *gm,
if (!FindEdgeLoop(numl)) break; if (!FindEdgeLoop(numl)) break;
} }
int nl=(int)edges[i].size(); int nl=(int)edges[i].size();
List_T *iListl = List_Create(nl, nl, sizeof(int)); List_T *temp = List_Create(nl, nl, sizeof(int));
for(int j = 0; j < nl; j++){ for(int j = 0; j < nl; j++){
int numEdge = edges[i][j]->tag(); int numEdge = edges[i][j]->tag();
List_Add(iListl, &numEdge); List_Add(temp, &numEdge);
} }
int type=ENT_LINE; int type=ENT_LINE;
if(select_contour(type, edges[i][0]->tag(), iListl)) if(select_contour(type, edges[i][0]->tag(), temp))
{ {
sortEdgesInLoop(numl, iListl); sortEdgesInLoop(numl, temp);
EdgeLoop *l = Create_EdgeLoop(numl, iListl); EdgeLoop *l = Create_EdgeLoop(numl, temp);
vecLoops.push_back(l); vecLoops.push_back(l);
Tree_Add(gm->getGEOInternals()->EdgeLoops, &l); Tree_Add(gm->getGEOInternals()->EdgeLoops, &l);
l->Num = numl; l->Num = numl;
} }
List_Delete(temp);
List_Delete(iListl);
} }
//create plane surfaces //create plane surfaces
......
...@@ -53,6 +53,48 @@ double computeLength(std::vector<MLine*> lines){ ...@@ -53,6 +53,48 @@ double computeLength(std::vector<MLine*> lines){
return length; return length;
} }
void removeBoundVertices(GFace *gf){
// Remove mesh_vertices that belong to l_edges
std::list<GEdge*> l_edges = gf->edges();
for(std::list<GEdge*>::iterator it = l_edges.begin(); it != l_edges.end(); it++){
std::vector<MVertex*> edge_vertices = (*it)->mesh_vertices;
std::vector<MVertex*>::const_iterator itv = edge_vertices.begin();
for(; itv != edge_vertices.end(); itv++){
std::vector<MVertex*>::iterator itve = std::find
(gf->mesh_vertices.begin(), gf->mesh_vertices.end(), *itv);
if (itve != gf->mesh_vertices.end()) gf->mesh_vertices.erase(itve);
}
MVertex *vB = (*it)->getBeginVertex()->mesh_vertices[0];
std::vector<MVertex*>::iterator itvB = std::find
(gf->mesh_vertices.begin(), gf->mesh_vertices.end(), vB);
if (itvB != gf->mesh_vertices.end()) gf->mesh_vertices.erase(itvB);
MVertex *vE = (*it)->getEndVertex()->mesh_vertices[0];
std::vector<MVertex*>::iterator itvE = std::find
(gf->mesh_vertices.begin(), gf->mesh_vertices.end(), vE);
if (itvE != gf->mesh_vertices.end()) gf->mesh_vertices.erase(itvE);
//if l_edge is a compond
if((*it)->getCompound()){
GEdgeCompound *gec = (*it)->getCompound();
std::vector<MVertex*> edge_vertices = gec->mesh_vertices;
std::vector<MVertex*>::const_iterator itv = edge_vertices.begin();
for(; itv != edge_vertices.end(); itv++){
std::vector<MVertex*>::iterator itve = std::find
(gf->mesh_vertices.begin(), gf->mesh_vertices.end(), *itv);
if (itve != gf->mesh_vertices.end()) gf->mesh_vertices.erase(itve);
}
MVertex *vB = (*it)->getBeginVertex()->mesh_vertices[0];
std::vector<MVertex*>::iterator itvB = std::find
(gf->mesh_vertices.begin(), gf->mesh_vertices.end(), vB);
if (itvB != gf->mesh_vertices.end()) gf->mesh_vertices.erase(itvB);
MVertex *vE = (*it)->getEndVertex()->mesh_vertices[0];
std::vector<MVertex*>::iterator itvE = std::find
(gf->mesh_vertices.begin(), gf->mesh_vertices.end(), vE);
if (itvE != gf->mesh_vertices.end()) gf->mesh_vertices.erase(itvE);
}
}
}
bool isClosed(std::set<MEdge, Less_Edge> &theCut){ bool isClosed(std::set<MEdge, Less_Edge> &theCut){
std::multiset<MVertex*> boundV; std::multiset<MVertex*> boundV;
...@@ -315,7 +357,7 @@ void cutTriangle(MTriangle *tri, ...@@ -315,7 +357,7 @@ void cutTriangle(MTriangle *tri,
} }
} }
Centerline::Centerline(std::string fileName): kdtree(0), nodes(0){ Centerline::Centerline(std::string fileName): kdtree(0), kdtreeR(0), nodes(0), nodesR(0){
recombine = CTX::instance()->mesh.recombineAll; recombine = CTX::instance()->mesh.recombineAll;
...@@ -332,7 +374,7 @@ Centerline::Centerline(std::string fileName): kdtree(0), nodes(0){ ...@@ -332,7 +374,7 @@ Centerline::Centerline(std::string fileName): kdtree(0), nodes(0){
is_closed = false; is_closed = false;
} }
Centerline::Centerline(): kdtree(0), nodes(0){ Centerline::Centerline(): kdtree(0), kdtreeR(0), nodes(0), nodesR(0){
index = new ANNidx[1]; index = new ANNidx[1];
dist = new ANNdist[1]; dist = new ANNdist[1];
...@@ -353,7 +395,9 @@ Centerline::Centerline(): kdtree(0), nodes(0){ ...@@ -353,7 +395,9 @@ Centerline::Centerline(): kdtree(0), nodes(0){
Centerline::~Centerline(){ Centerline::~Centerline(){
if (mod) delete mod; if (mod) delete mod;
if(kdtree) delete kdtree; if(kdtree) delete kdtree;
if(kdtreeR) delete kdtreeR;
if(nodes) annDeallocPts(nodes); if(nodes) annDeallocPts(nodes);
if(nodesR) annDeallocPts(nodesR);
delete[]index; delete[]index;
delete[]dist; delete[]dist;
} }
...@@ -367,11 +411,13 @@ void Centerline::importFile(std::string fileName){ ...@@ -367,11 +411,13 @@ void Centerline::importFile(std::string fileName){
if (gf->geomType() == GEntity::DiscreteSurface){ if (gf->geomType() == GEntity::DiscreteSurface){
for(unsigned int j = 0; j < gf->triangles.size(); j++) for(unsigned int j = 0; j < gf->triangles.size(); j++)
triangles.push_back(gf->triangles[j]); triangles.push_back(gf->triangles[j]);
if (is_cut){
gf->triangles.clear(); gf->triangles.clear();
gf->deleteVertexArrays(); gf->deleteVertexArrays();
current->remove(gf); current->remove(gf);
} }
} }
}
if(triangles.empty()){ if(triangles.empty()){
Msg::Error("Current GModel has no triangles ..."); Msg::Error("Current GModel has no triangles ...");
...@@ -431,24 +477,16 @@ void Centerline::createBranches(int maxN){ ...@@ -431,24 +477,16 @@ void Centerline::createBranches(int maxN){
junctions.insert(*it); junctions.insert(*it);
} }
} }
// printf("junctions =%d \n", junctions.size());
// std::set<MVertex*>::iterator itt = junctions.begin();
// for ( ; itt != junctions.end(); ++itt){
// MVertex *v = *itt;
// printf("JUNC = %d \n", v->getNum());
// }
//split edges //split edges
int tag = 0; int tag = 0;
for(unsigned int i = 0; i < color_edges.size(); ++i){ for(unsigned int i = 0; i < color_edges.size(); ++i){
std::vector<MLine*> lines = color_edges[i]; std::vector<MLine*> lines = color_edges[i];
//printf("EDGE %d line = %d \n", lines.size());
while (!lines.empty()) { while (!lines.empty()) {
std::vector<MLine*> myLines; std::vector<MLine*> myLines;
std::vector<MLine*>::iterator itl = lines.begin(); std::vector<MLine*>::iterator itl = lines.begin();
MVertex *vB = (*itl)->getVertex(0); MVertex *vB = (*itl)->getVertex(0);
MVertex *vE = (*itl)->getVertex(1); MVertex *vE = (*itl)->getVertex(1);
//printf("VB =%d vE = %d \n", vB->getNum(), vE->getNum());
myLines.push_back(*itl); myLines.push_back(*itl);
erase(lines, *itl); erase(lines, *itl);
itl = lines.begin(); itl = lines.begin();
...@@ -456,7 +494,6 @@ void Centerline::createBranches(int maxN){ ...@@ -456,7 +494,6 @@ void Centerline::createBranches(int maxN){
junctions.find(vB) != junctions.end()) ) { junctions.find(vB) != junctions.end()) ) {
MVertex *v1 = (*itl)->getVertex(0); MVertex *v1 = (*itl)->getVertex(0);
MVertex *v2 = (*itl)->getVertex(1); MVertex *v2 = (*itl)->getVertex(1);
//printf("line %d, VB = %d vE = %d v1=%d v2=%d \n", (*itl)->getNum(), vB->getNum(), vE->getNum(), v1->getNum(), v2->getNum());
bool goVE = (junctions.find(vE) == junctions.end()) ? true : false ; bool goVE = (junctions.find(vE) == junctions.end()) ? true : false ;
bool goVB = (junctions.find(vB) == junctions.end()) ? true : false; bool goVB = (junctions.find(vB) == junctions.end()) ? true : false;
if (v1 == vE && goVE){ if (v1 == vE && goVE){
...@@ -489,7 +526,6 @@ void Centerline::createBranches(int maxN){ ...@@ -489,7 +526,6 @@ void Centerline::createBranches(int maxN){
orderMLines(myLines, vB, vE); orderMLines(myLines, vB, vE);
std::vector<Branch> children; std::vector<Branch> children;
Branch newBranch ={ tag++, myLines, computeLength(myLines), vB, vE, children, 1.e6, 0.0}; 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());
edges.push_back(newBranch) ; edges.push_back(newBranch) ;
} }
} }
...@@ -510,32 +546,15 @@ void Centerline::createBranches(int maxN){ ...@@ -510,32 +546,15 @@ void Centerline::createBranches(int maxN){
distanceToSurface(); distanceToSurface();
computeRadii(); computeRadii();
//print for debug
printSplit(); printSplit();
//print info
// for(unsigned int i = 0; i < edges.size(); ++i) {
// printf("EDGE =%d tag=%d length = %g childs = %d \n", i, edges[i].tag, edges[i].length, edges[i].children.size());
// for(unsigned int j = 0; j < edges[i].children.size(); ++j) {
// printf("children (%d) =%d \n", edges[i].children.size(), edges[i].children[j].tag);
// }
// }
// std::set<MVertex*>::iterator itj = junctions.begin();
// for ( ; itj != junctions.end(); ++itj){
// MVertex *v = *itj;
// printf("JUNC = %d \n", v->getNum());
// }
} }
void Centerline::distanceToSurface(){ void Centerline::distanceToSurface(){
Msg::Info("Centerline: computing distance to surface mesh "); Msg::Info("Centerline: computing distance to surface mesh ");
//COMPUTE WITH 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];
ANNdistArray distR = new ANNdist[1];
std::set<MVertex*> allVS; std::set<MVertex*> allVS;
for(int j = 0; j < triangles.size(); j++) for(int j = 0; j < triangles.size(); j++)
for(int k = 0; k<3; k++) allVS.insert(triangles[j]->getVertex(k)); for(int k = 0; k<3; k++) allVS.insert(triangles[j]->getVertex(k));
...@@ -557,14 +576,10 @@ void Centerline::distanceToSurface(){ ...@@ -557,14 +576,10 @@ void Centerline::distanceToSurface(){
MVertex *v1 = l->getVertex(0); MVertex *v1 = l->getVertex(0);
MVertex *v2 = l->getVertex(1); MVertex *v2 = l->getVertex(1);
double midp[3] = {0.5*(v1->x()+v2->x()), 0.5*(v1->y()+v1->y()),0.5*(v1->z()+v2->z())}; 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); kdtreeR->annkSearch(midp, 1, index, dist);
double minRad = sqrt(distR[0]); double minRad = sqrt(dist[0]);
radiusl.insert(std::make_pair(lines[i], minRad)); radiusl.insert(std::make_pair(lines[i], minRad));
} }
delete kdtreeR;
annDeallocPts(nodesR);
delete[]indexR;
delete[]distR;
} }
void Centerline::computeRadii(){ void Centerline::computeRadii(){
...@@ -633,6 +648,7 @@ void Centerline::createSplitCompounds(){ ...@@ -633,6 +648,7 @@ void Centerline::createSplitCompounds(){
NV = current->getMaxElementaryNumber(0); NV = current->getMaxElementaryNumber(0);
NE = current->getMaxElementaryNumber(1); NE = current->getMaxElementaryNumber(1);
NF = current->getMaxElementaryNumber(2); NF = current->getMaxElementaryNumber(2);
NR = current->getMaxElementaryNumber(3);
// Remesh new faces (Compound Lines and Compound Surfaces) // Remesh new faces (Compound Lines and Compound Surfaces)
Msg::Info("*** Starting parametrize compounds:"); Msg::Info("*** Starting parametrize compounds:");
...@@ -665,15 +681,24 @@ void Centerline::createSplitCompounds(){ ...@@ -665,15 +681,24 @@ void Centerline::createSplitCompounds(){
GFaceCompound *gfc = new GFaceCompound(current, num_gfc, f_compound, U0, GFaceCompound *gfc = new GFaceCompound(current, num_gfc, f_compound, U0,
typ, 0); typ, 0);
gfc->meshAttributes.recombine = recombine; gfc->meshAttributes.recombine = recombine;
gfc->addPhysicalEntity(i+1); if (i < discFaces.size()){
gfc->addPhysicalEntity(100);
current->setPhysicalName("newsurf", 2, 100);
}
else{
gfc->addPhysicalEntity(200);
current->setPhysicalName("in/out", 2, 200);
}
current->add(gfc); current->add(gfc);
} }
} }
void Centerline::cleanMesh(){ void Centerline::cleanMesh(){
return;
if (!is_cut) return;
for (int i=0; i < NF; i++){ for (int i=0; i < NF; i++){
std::ostringstream oss; std::ostringstream oss;
oss << "new" << NF+i+1 ; oss << "new" << NF+i+1 ;
...@@ -683,15 +708,16 @@ void Centerline::cleanMesh(){ ...@@ -683,15 +708,16 @@ void Centerline::cleanMesh(){
Msg::Info("Writing new splitted mesh mySPLITMESH.msh"); Msg::Info("Writing new splitted mesh mySPLITMESH.msh");
current->writeMSH("mySPLITMESH.msh", 2.2, false, false); current->writeMSH("mySPLITMESH.msh", 2.2, false, false);
if (!is_cut) return;
std::set<MVertex*> allNod; std::set<MVertex*> allNod;
std::list<GEdge*> U0;
discreteFace * mySplitMesh; discreteFace * mySplitMesh;
std::vector<std::set<MVertex*> > inOutNod;
std::vector<discreteFace* > inOutMesh;
mySplitMesh= new discreteFace(current, 2*NF+1); mySplitMesh= new discreteFace(current, 2*NF+1);
mySplitMesh->addPhysicalEntity(2*NF+1);
current->setPhysicalName("surface mesh", 2, 2*NF+1);
current->add(mySplitMesh); current->add(mySplitMesh);
for (int i=0; i < NF; i++){ for (int i=0; i < discFaces.size(); i++){
GFace *gfc = current->getFaceByTag(NF+i+1); GFace *gfc = current->getFaceByTag(NF+i+1);
for(unsigned int j = 0; j < gfc->triangles.size(); ++j){ for(unsigned int j = 0; j < gfc->triangles.size(); ++j){
MTriangle *t = gfc->triangles[j]; MTriangle *t = gfc->triangles[j];
...@@ -713,6 +739,33 @@ void Centerline::cleanMesh(){ ...@@ -713,6 +739,33 @@ void Centerline::cleanMesh(){
} }
} }
// int nbInOut = NF - discFaces.size();
// inOutMesh.resize(nbInOut);
// inOutNod.resize(nbInOut);
// for (int i=0; i < nbInOut; i++){
// inOutMesh[i]= new discreteFace(current, 2*NF+2+i);
// current->add(inOutMesh[i]);
// GFace *gfc = current->getFaceByTag(NF+discFaces.size()+i+1);
// for(unsigned int j = 0; j < gfc->triangles.size(); ++j){
// MTriangle *t = gfc->triangles[j];
// std::vector<MVertex *> v(3);
// for(int k = 0; k < 3; k++){
// v[k] = t->getVertex(k);
// inOutNod[i].insert(v[k]);
// }
// inOutMesh[i]->triangles.push_back(new MTriangle(v[0], v[1], v[2]));
// }
// for(unsigned int j = 0; j < gfc->quadrangles.size(); ++j){
// MQuadrangle *q = gfc->quadrangles[j];
// std::vector<MVertex *> v(4);
// for(int k = 0; k < 4; k++){
// v[k] = q->getVertex(k);
// inOutNod[i].insert(v[k]);
// }
// inOutMesh[i]->quadrangles.push_back(new MQuadrangle(v[0], v[1], v[2], v[3]));
// }
// }
//Removing discrete Vertices - Edges - Faces //Removing discrete Vertices - Edges - Faces
for (int i=0; i < NV; i++){ for (int i=0; i < NV; i++){
GVertex *gv = current->getVertexByTag(i+1); GVertex *gv = current->getVertexByTag(i+1);
...@@ -726,17 +779,30 @@ void Centerline::cleanMesh(){ ...@@ -726,17 +779,30 @@ void Centerline::cleanMesh(){
} }
for (int i=0; i < NF; i++){ for (int i=0; i < NF; i++){
GFace *gf = current->getFaceByTag(i+1); GFace *gf = current->getFaceByTag(i+1);
GFace *gfc = current->getFaceByTag(NF+i+1);
current->remove(gf); current->remove(gf);
}
for (int i=0; i < discFaces.size(); i++){
GFace *gfc = current->getFaceByTag(NF+i+1);
current->remove(gfc); current->remove(gfc);
} }
//Put new mesh in a new discreteFace //Put new mesh in a new discreteFace
for(std::set<MVertex*>::iterator it = allNod.begin(); it != allNod.end(); ++it) for(std::set<MVertex*>::iterator it = allNod.begin(); it != allNod.end(); ++it)
mySplitMesh->mesh_vertices.push_back(*it); mySplitMesh->mesh_vertices.push_back(*it);
removeBoundVertices(mySplitMesh);
mySplitMesh->meshStatistics.status = GFace::DONE; mySplitMesh->meshStatistics.status = GFace::DONE;
current->createTopologyFromMesh();
// for (int i = 0; i< nbInOut; i++){
// for(std::set<MVertex*>::iterator it = inOutNod[i].begin(); it != inOutNod[i].end(); ++it){
// printf("mesh vertex =%d \n", (*it)->getNum());
// inOutMesh[i]->mesh_vertices.push_back(*it);
// }
// printf("inOutMesh[%d]->mesh vertices =%d \n", inOutMesh[i]->tag(), inOutMesh[i]->mesh_vertices.size());
// inOutMesh[i]->meshStatistics.status = GFace::DONE;
// }
current->createTopologyFromMesh();
} }
void Centerline::createFaces(){ void Centerline::createFaces(){
...@@ -811,13 +877,35 @@ void Centerline::createInOutFaces(){ ...@@ -811,13 +877,35 @@ void Centerline::createInOutFaces(){
//create the inlet and outlet planar face //create the inlet and outlet planar face
current->setFactory("Gmsh"); current->setFactory("Gmsh");
//std::vector<std::vector<GFace *> > myFaceLoops;
//std::vector<GFace *> myFaces;
for (int i = 0; i< boundEdges.size(); i++){ for (int i = 0; i< boundEdges.size(); i++){
std::vector<std::vector<GEdge *> > myEdges; std::vector<std::vector<GEdge *> > myEdgeLoops;
std::vector<GEdge *> myEdge; std::vector<GEdge *> myEdges;
myEdge.push_back(boundEdges[i]); myEdges.push_back(boundEdges[i]);
myEdges.push_back(myEdge); myEdgeLoops.push_back(myEdges);
GFace *newFace = current->addPlanarFace(myEdges); GFace *newFace = current->addPlanarFace(myEdgeLoops);
//myFaces.push_back(newFace);
}
//for (int i= 0; i< discFaces.size(); i++)
// myFaces.push_back((GFace*)discFaces[i]);
//myFaceLoops.push_back(myFaces);
//GRegion *reg = current->addVolume(myFaceLoops);
}
void Centerline::createClosedVolume(){
std::vector<std::vector<GFace *> > myFaceLoops;
std::vector<GFace *> myFaces;
for (int i = 0; i < NF; i++){
GFace * gf = current->getFaceByTag(NF+i+1);
myFaces.push_back(gf);
} }
myFaceLoops.push_back(myFaces);
GRegion *reg = current->addVolume(myFaceLoops);
reg->addPhysicalEntity(reg->tag());
current->setPhysicalName("newvol", 3, reg->tag());
} }
...@@ -906,6 +994,7 @@ void Centerline::cutMesh(){ ...@@ -906,6 +994,7 @@ void Centerline::cutMesh(){
//create compounds //create compounds
createSplitCompounds(); createSplitCompounds();
if (is_closed) createClosedVolume();
Msg::Info("Splitting mesh by centerlines done "); Msg::Info("Splitting mesh by centerlines done ");
...@@ -967,7 +1056,7 @@ void Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){ ...@@ -967,7 +1056,7 @@ void Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
break; break;
} }
else { else {
if (step == 9) {printf("no closed cut %d \n", (int)newCut.size()); }; //if (step == 9) {printf("no closed cut %d \n", (int)newCut.size()); };
step++; step++;
// // triangles = newTris; // // triangles = newTris;
// // theCut.insert(newCut.begin(),newCut.end()); // // theCut.insert(newCut.begin(),newCut.end());
...@@ -1007,6 +1096,21 @@ double Centerline::operator() (double x, double y, double z, GEntity *ge){ ...@@ -1007,6 +1096,21 @@ double Centerline::operator() (double x, double y, double z, GEntity *ge){
} }
double xyz[3] = {x,y,z}; double xyz[3] = {x,y,z};
bool isCompound = (ge->dim() == 2 && ge->geomType() == GEntity::CompoundSurface) ? true: false;
std::list<GFace*> cFaces;
if (isCompound) cFaces = ((GFaceCompound*)ge)->getCompounds();
//take xyz = closest point on boundary in case we are on the planar in/out faces
if ( ge->dim() == 3 || (ge->dim() == 2 && ge->geomType() == GEntity::Plane) ||
(isCompound && (*cFaces.begin())->geomType() == GEntity::Plane) ){
int num_neighbours = 1;
kdtreeR->annkSearch(xyz, num_neighbours, index, dist);
xyz[0] = nodesR[index[0]][0];
xyz[1] = nodesR[index[0]][1];
xyz[2] = nodesR[index[0]][2];
}
int num_neighbours = 1; int num_neighbours = 1;
kdtree->annkSearch(xyz, num_neighbours, index, dist); kdtree->annkSearch(xyz, num_neighbours, index, dist);
double d = sqrt(dist[0]); double d = sqrt(dist[0]);
......
...@@ -55,14 +55,14 @@ class Centerline : public Field{ ...@@ -55,14 +55,14 @@ class Centerline : public Field{
GModel *current; //current GModel GModel *current; //current GModel
GModel *mod; //centerline GModel GModel *mod; //centerline GModel
GModel *split; //split GModel GModel *split; //split GModel
ANNkd_tree *kdtree; ANNkd_tree *kdtree, *kdtreeR;
ANNpointArray nodes; ANNpointArray nodes, nodesR;
ANNidxArray index; ANNidxArray index;
ANNdistArray dist; ANNdistArray dist;
std::string fileName; std::string fileName;
int nbPoints; int nbPoints;
double recombine; double recombine;
int NF, NV, NE; int NF, NV, NE, NR;
bool is_cut; bool is_cut;
bool is_closed; bool is_closed;
...@@ -147,6 +147,7 @@ class Centerline : public Field{ ...@@ -147,6 +147,7 @@ class Centerline : public Field{
//create discrete faces //create discrete faces
void createFaces(); void createFaces();
void createInOutFaces(); void createInOutFaces();
void createClosedVolume();
void createSplitCompounds(); void createSplitCompounds();
//Print for debugging //Print for debugging
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment