diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 256032c68aaec3fa13f7597abdd455bb22d3ec8f..9e066aad1daf07698de1af13a9054ba298ec6dfc 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2199,6 +2199,13 @@ GEntity *GModel::extrude(GEntity *e, std::vector<double> p1, std::vector<double>
return 0;
}
+GEntity *GModel::extrudeBoundaryLayer(GEntity *e, int nbLayers, double hLayers, int dir, int view)
+{
+ if(_factory)
+ return _factory->extrudeBoundaryLayer(this, e, nbLayers,hLayers, dir, view);
+ return 0;
+}
+
GEntity *GModel::addPipe(GEntity *e, std::vector<GEdge *> edges)
{
if(_factory)
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 18d062f313ad2c41aa820050fb3e20def52a8ebe..f2812d1359f682eede8796114563fa0fc4e5cb1e 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -450,6 +450,7 @@ class GModel
GEntity *revolve(GEntity *e, std::vector<double> p1, std::vector<double> p2,
double angle);
GEntity *extrude(GEntity *e, std::vector<double> p1, std::vector<double> p2);
+ GEntity *extrudeBoundaryLayer(GEntity *e, int nbLayers, double hLayers, int dir=1, int view=-1);
GEntity *addPipe(GEntity *e, std::vector<GEdge *> edges);
void addRuledFaces(std::vector<std::vector<GEdge *> > edges);
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 9b2d65edcafc871f3d774b5b42607bb635fdfb4d..c88833d1118f597f311c74b3ac208b6f1b899b01 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -16,6 +16,8 @@
#include "MLine.h"
#include "GModel.h"
#include "Numeric.h"
+#include "ExtrudeParams.h"
+#include "Geo.h"
GVertex *GeoFactory::addVertex(GModel *gm, double x, double y, double z, double lc)
{
@@ -283,6 +285,101 @@ std::vector<GFace *> GeoFactory::addRuledFaces(GModel *gm,
return faces;
}
+GEntity* GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers, double hLayer, int dir, int view)
+{
+
+ printf("coucou dans geo extrude normals \n");
+ Msg::Debug("Gmsh model (GModel) before extrude has:");
+ Msg::Debug("%d Vertices", gm->getNumVertices());
+ Msg::Debug("%d Edges", gm->getNumEdges());
+ Msg::Debug("%d Faces", gm->getNumFaces());
+ Msg::Debug("%d Regions", gm->getNumRegions());
+
+ ExtrudeParams ep;
+ printf("***** BL INDEX =%d \n", dir);
+ ep.mesh.BoundaryLayerIndex = dir;
+ ep.mesh.ViewIndex = view;
+ ep.mesh.NbLayer = 1; //this may be more general for defining different layers
+ ep.mesh.hLayer.clear();
+ ep.mesh.hLayer.push_back(hLayer);
+ ep.mesh.NbElmLayer.clear();
+ ep.mesh.NbElmLayer.push_back(nbLayers);
+ ep.mesh.ExtrudeMesh = true;
+
+ int type = BOUNDARY_LAYER; //TRANSLATE;
+ double T0=0., T1=0., T2=0.;
+ double A0=0., A1=0., A2=0.;
+ double X0=0., X1=0., X2=0.,alpha=0.;
+
+ //extrude shape dans geo.cpp
+ Shape shape;
+ if(e->dim() == 0){
+ Vertex *v = FindPoint(e->tag());
+ if(!v) printf("vertex %d not found \n", e->tag());
+ shape.Num = v->Num;
+ shape.Type = v->Typ;
+ }
+ else if (e->dim() == 1){
+ ((GEdge*)e)->meshAttributes.extrude = &ep;
+ Curve *c = FindCurve(e->tag());
+ if(!c) printf("curve %d not found \n", e->tag());
+ shape.Num = c->Num;
+ shape.Type = c->Typ;
+ }
+ else if (e->dim() == 2){
+ ((GFace*)e)->meshAttributes.extrude = &ep;
+ Surface *s = FindSurface(e->tag());
+ if(!s) printf("surface %d not found \n", e->tag());
+ else printf("surface %d found type =%d\n", e->tag(), s->Typ);
+ shape.Num = s->Num;
+ shape.Type = s->Typ;
+ }
+
+ //extrude shape
+ List_T *list_out= List_Create(2, 1, sizeof(Shape));
+ List_T *tmp = List_Create(1, 1, sizeof(Shape));
+ List_Add(tmp, &shape);
+ ExtrudeShapes(type, tmp,
+ T0, T1, T2,
+ A0, A1, A2,
+ X0, X1, X2, alpha,
+ &ep,
+ list_out);
+
+ //create GEntities
+ gm->importGEOInternals();
+
+ //return the new created entity
+ GEntity *newEnt=0;
+ if (e->dim() == 1){
+ for(GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); it++){
+ GEdge *ge = *it;
+ if(ge->getNativeType() == GEntity::GmshModel &&
+ ge->geomType() == GEntity::BoundaryLayerCurve){
+ ExtrudeParams *ep = ge->meshAttributes.extrude;
+ if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == COPIED_ENTITY &&
+ std::abs(ep->geo.Source) ==e->tag() )
+ newEnt = ge;
+ }
+ }
+ }
+ else if (e->dim() ==2){
+ for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); it++){
+ GFace *gf = *it;
+ if(gf->getNativeType() == GEntity::GmshModel &&
+ gf->geomType() == GEntity::BoundaryLayerSurface){
+ ExtrudeParams *ep = gf->meshAttributes.extrude;
+ if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == COPIED_ENTITY
+ && std::abs(ep->geo.Source) == e->tag())
+ newEnt = gf;
+ }
+ }
+ }
+
+ return newEnt;
+
+};
+
#if defined(HAVE_OCC)
#include "OCCIncludes.h"
#include "GModelIO_OCC.h"
@@ -557,8 +654,6 @@ GEntity *OCCFactory::extrude(GModel *gm, GEntity* base,
const double z2 = p2[2];
gp_Vec direction(gp_Pnt(x1, y1, z1), gp_Pnt(x2, y2, z2));
- gp_Ax1 axisOfRevolution(gp_Pnt(x1, y1, z1), direction);
-
BRepPrimAPI_MakePrism MP(*(TopoDS_Shape*)base->getNativePtr(), direction,
Standard_False);
GEntity *ret = 0;
diff --git a/Geo/GModelFactory.h b/Geo/GModelFactory.h
index 1e9f363141a46968f4c7bb3f0cd1926a3a26e940..076c805c1b9cdf35ac949ae6995facd4d2ac9477 100644
--- a/Geo/GModelFactory.h
+++ b/Geo/GModelFactory.h
@@ -90,6 +90,11 @@ class GModelFactory {
Msg::Error("extrude not implemented yet");
return 0;
}
+ virtual GEntity* extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers, double hLayers, int dir, int view)
+ {
+ Msg::Error("extrude normals not implemented yet");
+ return 0;
+ }
virtual GEntity *addPipe(GModel *gm, GEntity *base, std::vector<GEdge *> wire)
{
Msg::Error("addPipe not implemented yet");
@@ -177,6 +182,7 @@ class GeoFactory : public GModelFactory {
GRegion *addVolume(GModel *gm, std::vector<std::vector<GFace *> > faces);
GEdge *addCircleArc(GModel *gm,GVertex *begin, GVertex *center, GVertex *end);
std::vector<GFace *> addRuledFaces(GModel *gm, std::vector<std::vector<GEdge *> > edges);
+ GEntity* extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers, double hLayers, int dir, int view);
};
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index a8b0ab1165277221cc8f0e3ed53c2f1ff4550630..34f20c910803f5ac7c23b7a42bc3b1a7af2110ae 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2134,12 +2134,12 @@ void ProtudeXYZ(double &x, double &y, double &z, ExtrudeParams *e)
List_Reset(ListOfTransformedPoints);
}
-static int Extrude_ProtudePoint(int type, int ip,
- double T0, double T1, double T2,
- double A0, double A1, double A2,
- double X0, double X1, double X2, double alpha,
- Curve **pc, Curve **prc, int final,
- ExtrudeParams *e)
+int Extrude_ProtudePoint(int type, int ip,
+ double T0, double T1, double T2,
+ double A0, double A1, double A2,
+ double X0, double X1, double X2, double alpha,
+ Curve **pc, Curve **prc, int final,
+ ExtrudeParams *e)
{
double matrix[4][4], T[3], Ax[3], d;
Vertex V, *pv, *newp, *chapeau;
@@ -2298,12 +2298,12 @@ static int Extrude_ProtudePoint(int type, int ip,
return chapeau->Num;
}
-static int Extrude_ProtudeCurve(int type, int ic,
- double T0, double T1, double T2,
- double A0, double A1, double A2,
- double X0, double X1, double X2, double alpha,
- Surface **ps, int final,
- ExtrudeParams *e)
+int Extrude_ProtudeCurve(int type, int ic,
+ double T0, double T1, double T2,
+ double A0, double A1, double A2,
+ double X0, double X1, double X2, double alpha,
+ Surface **ps, int final,
+ ExtrudeParams *e)
{
double matrix[4][4], T[3], Ax[3];
Curve *CurveBeg, *CurveEnd;
@@ -2469,11 +2469,11 @@ static int Extrude_ProtudeCurve(int type, int ic,
return chapeau->Num;
}
-static int Extrude_ProtudeSurface(int type, int is,
- double T0, double T1, double T2,
- double A0, double A1, double A2,
- double X0, double X1, double X2, double alpha,
- Volume **pv, ExtrudeParams *e)
+int Extrude_ProtudeSurface(int type, int is,
+ double T0, double T1, double T2,
+ double A0, double A1, double A2,
+ double X0, double X1, double X2, double alpha,
+ Volume **pv, ExtrudeParams *e)
{
double matrix[4][4], T[3], Ax[3];
Curve *c, *c2;
@@ -2684,6 +2684,7 @@ void ExtrudeShapes(int type, List_T *list_in,
ExtrudeParams *e,
List_T *list_out)
{
+
for(int i = 0; i < List_Nbr(list_in); i++){
Shape shape;
List_Read(list_in, i, &shape);
@@ -2753,6 +2754,7 @@ void ExtrudeShapes(int type, List_T *list_in,
top.Num = Extrude_ProtudeSurface(type, shape.Num, T0, T1, T2,
A0, A1, A2, X0, X1, X2, alpha,
&pv, e);
+ printf("extruded surface =%d \n", top.Num);
Surface *ps = FindSurface(top.Num);
top.Type = ps ? ps->Typ : 0;
@@ -2761,6 +2763,7 @@ void ExtrudeShapes(int type, List_T *list_in,
Shape body;
body.Num = pv->Num;
body.Type = pv->Typ;
+ printf("extruded volume =%d \n", pv->Num);
List_Add(list_out, &body);
if(CTX::instance()->geom.extrudeReturnLateral){
for(int j = 0; j < List_Nbr(pv->Surfaces); j++){
@@ -2783,6 +2786,7 @@ void ExtrudeShapes(int type, List_T *list_in,
break;
}
}
+
}
// Duplicate removal
diff --git a/Geo/Geo.h b/Geo/Geo.h
index e670cd0f623d048fc5330ab7f8492f25343c2990..371da95bf959d56ff7e84ccf90460fc9ab1623dd 100644
--- a/Geo/Geo.h
+++ b/Geo/Geo.h
@@ -297,6 +297,24 @@ void ExtrudeShapes(int extrude_type, List_T *in,
ExtrudeParams *e,
List_T *out);
+int Extrude_ProtudePoint(int type, int ip,
+ double T0, double T1, double T2,
+ double A0, double A1, double A2,
+ double X0, double X1, double X2, double alpha,
+ Curve **pc, Curve **prc, int final,
+ ExtrudeParams *e);
+int Extrude_ProtudeCurve(int type, int ic,
+ double T0, double T1, double T2,
+ double A0, double A1, double A2,
+ double X0, double X1, double X2, double alpha,
+ Surface **ps, int final,
+ ExtrudeParams *e);
+int Extrude_ProtudeSurface(int type, int is,
+ double T0, double T1, double T2,
+ double A0, double A1, double A2,
+ double X0, double X1, double X2, double alpha,
+ Volume **pv, ExtrudeParams *e);
+
void ProtudeXYZ(double &x, double &y, double &z, ExtrudeParams *e);
void ReplaceAllDuplicates();
diff --git a/Mesh/BoundaryLayers.cpp b/Mesh/BoundaryLayers.cpp
index f12bcfaeb22a26e9e5fed7601cdb8aca82019c94..844a4ff60fa4f38884930d4067f6869c4d3a48c2 100644
--- a/Mesh/BoundaryLayers.cpp
+++ b/Mesh/BoundaryLayers.cpp
@@ -29,7 +29,7 @@ static void addExtrudeNormals(std::vector<T*> &elements, int invert,
Msg::Error("Boundary layer index should be 0 or 1");
return;
}
-
+
if(octree && !gouraud){ // get extrusion direction from post-processing view
std::set<MVertex*> verts;
for(unsigned int i = 0; i < elements.size(); i++)
@@ -45,6 +45,7 @@ static void addExtrudeNormals(std::vector<T*> &elements, int invert,
}
}
else{ // get extrusion direction from Gouraud-shaded element normals
+ printf("get extrusion from normals \n");
for(unsigned int i = 0; i < elements.size(); i++){
MElement *ele = elements[i];
SVector3 n(0, 0, 0);
@@ -82,6 +83,7 @@ static void addExtrudeNormals(std::set<T*> &entities,
OctreePost *octree = 0;
#if defined(HAVE_POST)
if(view != -1){
+ printf("extrude normals with view %d \n", view);
if(view >= 0 && view < (int)PView::list.size()){
octree = new OctreePost(PView::list[view]);
if(PView::list[view]->getData()->getNumVectors())
@@ -186,6 +188,7 @@ int Mesh2DWithBoundaryLayers(GModel *m)
Msg::Error("Unknown source curve %d for boundary layer", ep->geo.Source);
return 0;
}
+ //printf("index edge (%d) =%d \n", ge->tag(), ep->mesh.BoundaryLayerIndex);
std::pair<bool, std::pair<int, int> > tags
(ep->geo.Source < 0, std::pair<int, int>
(ep->mesh.BoundaryLayerIndex, ep->mesh.ViewIndex));
@@ -207,6 +210,7 @@ int Mesh2DWithBoundaryLayers(GModel *m)
Msg::Error("Unknown source face %d for boundary layer", ep->geo.Source);
return 0;
}
+ //printf("index surface (%d) =%d \n", gf->tag(), ep->mesh.BoundaryLayerIndex);
std::pair<bool, std::pair<int, int> > tags
(ep->geo.Source < 0, std::pair<int, int>
(ep->mesh.BoundaryLayerIndex, ep->mesh.ViewIndex));
@@ -218,8 +222,7 @@ int Mesh2DWithBoundaryLayers(GModel *m)
}
}
- if(sourceEdges.empty() && sourceFaces.empty())
- return 0;
+ if(sourceEdges.empty() && sourceFaces.empty()) return 0;
// compute mesh dependencies in source faces (so we can e.g. create
// a boundary layer on an extruded mesh)
@@ -280,6 +283,7 @@ int Mesh2DWithBoundaryLayers(GModel *m)
vdest = ge->getEndVertex();
}
GPoint p = vsrc->point();
+
ep->Extrude(ep->mesh.NbLayer - 1, ep->mesh.NbElmLayer[ep->mesh.NbLayer - 1],
p.x(), p.y(), p.z());
vdest->setPosition(p);
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index cbfb2b4c076ab11cea53924386e14b4ff03fac51..6b68c987063a10f6bac339802d986c68cd0b5037 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -326,10 +326,13 @@ Centerline::Centerline(std::string fileName): kdtree(0), kdtreeR(0), nodes(0), n
importFile(fileName);
buildKdTree();
nbPoints = 25;
+ hLayer = 0.3;
+ nbElemLayer = 3;
update_needed = false;
is_cut = false;
is_closed = false;
+ is_extruded = false;
}
Centerline::Centerline(): kdtree(0), kdtreeR(0), nodes(0), nodesR(0){
@@ -340,13 +343,21 @@ Centerline::Centerline(): kdtree(0), kdtreeR(0), nodes(0), nodesR(0){
recombine = CTX::instance()->mesh.recombineAll;
fileName = "centerlines.vtk";//default
nbPoints = 25;
+ hLayer = 0.3;
+ nbElemLayer = 3;
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");
+
+ callbacks["extrudeWall"] = new FieldCallbackGeneric<Centerline>(this, &Centerline::extrudeWall, "Extrude wall \n");
+ options["nbElemLayer"] = new FieldOptionInt(nbElemLayer, "Number of mesh elements the extruded layer");
+ options["hLayer"] = new FieldOptionDouble(hLayer, "Thickness (% of radius) of the extruded layer");
+
is_cut = false;
is_closed = false;
+ is_extruded = false;
}
@@ -827,6 +838,14 @@ void Centerline::closeVolume(){
}
+void Centerline::extrudeWall(){
+
+ is_extruded = true;
+ //printf("calling extrude wall \n");
+ //exit(1);
+
+}
+
void Centerline::cutMesh(){
is_cut = true;
@@ -860,9 +879,13 @@ void Centerline::cutMesh(){
double L = edges[i].length;
double D = (edges[i].minRad+edges[i].maxRad);
double AR = L/D;
- printf("*** Centerline branch %d (AR=%d) \n", i, (int)floor(AR + 0.5));
- if( AR > 2.5){
- int nbSplit = (int)floor(AR/3. + 0.5);
+ printf("*** Centerline branch %d (AR=%g): ", edges[i].tag, AR);
+ printf("children (%d) = ", edges[i].children.size());
+ for (int k= 0; k< edges[i].children.size() ; k++) printf("%d ", edges[i].children[k].tag);
+ printf("\n");
+
+ int nbSplit = (int)floor(AR/2. + 0.5);
+ if( AR > 3.0){
double li = L/nbSplit;
double lc = 0.0;
for (unsigned int j= 0; j < lines.size(); j++){
@@ -873,7 +896,8 @@ void Centerline::cutMesh(){
SVector3 pt(v1->x(), v1->y(), v1->z());
SVector3 dir(v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
printf("->> cut length %d split \n", nbSplit);
- cutByDisk(pt, dir, edges[i].maxRad);
+ std::map<MLine*,double>::iterator itr = radiusl.find(lines[j]);
+ bool cutted = cutByDisk(pt, dir, itr->second);
nbSplit--;
lc = 0.0;
}
@@ -888,7 +912,18 @@ void Centerline::cutMesh(){
SVector3 pt(v1->x(), v1->y(), v1->z());
SVector3 dir(v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
printf("-->> cut bifurcation \n");
- cutByDisk(pt, dir, edges[i].maxRad);
+ std::map<MLine*,double>::iterator itr = radiusl.find(lines[lines.size()-1]);
+ bool cutted = cutByDisk(pt, dir, itr->second);
+ if(!cutted){
+ int l = lines.size()-1-lines.size()/(4*nbSplit);
+ v1 = lines[l]->getVertex(1);
+ v2 = lines[l]->getVertex(0);
+ pt = SVector3(v1->x(), v1->y(), v1->z());
+ dir = SVector3(v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
+ printf("-->> cut bifurcation NEW \n");
+ itr = radiusl.find(lines[l]);
+ cutted = cutByDisk(pt, dir, itr->second);
+ }
}
}
@@ -905,11 +940,14 @@ void Centerline::cutMesh(){
createSplitCompounds();
if (is_closed) createClosedVolume();
+ //extrude wall
+ //if(is_extruded) extrudeNormalWall();
+
Msg::Info("Splitting mesh by centerlines done ");
}
-void Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
+bool Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
double a = NORM.x();
double b = NORM.y();
@@ -917,15 +955,17 @@ void Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
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);
+ int maxStep = 40;
const double EPS = 0.007;
+
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 < 20){
- double rad = 1.2*maxRad+0.1*step*maxRad;
+ while (!closedCut && step < maxStep){
+ double rad = 1.1*maxRad+0.05*step*maxRad;
std::map<MEdge,MVertex*,Less_Edge> cutEdges;
std::set<MVertex*> cutVertices;
std::vector<MTriangle*> newTris;
@@ -961,12 +1001,11 @@ void Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
triangles.clear();
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", (int)newCut.size()); };
step++;
+ //if (step == maxStep) {printf("no closed cut %d \n", (int)newCut.size()); };
// // triangles = newTris;
// // theCut.insert(newCut.begin(),newCut.end());
// char name[256];
@@ -987,9 +1026,15 @@ void Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
}
}
-
-
- return;
+
+ if (step < maxStep){
+ //printf("cutByDisk OK step =%d \n", step);
+ return true;
+ }
+ else {
+ //printf("cutByDisk not succeeded \n");
+ return false;
+ }
}
@@ -1072,14 +1117,14 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
buildOrthoBasis2(dir_a, dir_n, dir_t);
double lc = 2*M_PI*radMax/nbPoints;
- double lc_a = 3.*lc;
+ double lc_a = 3.5*lc;
double lc_n, lc_t;
if (onTubularSurface){
lc_n = lc_t = lc;
}
else{
- double e = radMax/3.;
+ double e = radMax/4.;
double hn = e/10.;
double rm = std::max(radMax-ds, radMax-e);
lc_t = 2*M_PI*rm/nbPoints;
@@ -1193,7 +1238,7 @@ void Centerline::printSplit() const{
fprintf(f, "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(),
- (double)i, (double)i);
+ (double)edges[i].tag, (double)edges[i].tag);
}
}
fprintf(f,"};\n");
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
index 6e82d6361fbd709ab12d7a23de85d4760add4fa3..6d1711817749d123fa49cda73819cda2ef98c52d 100644
--- a/Mesh/CenterlineField.h
+++ b/Mesh/CenterlineField.h
@@ -65,7 +65,9 @@ class Centerline : public Field{
int NF, NV, NE, NR;
bool is_cut;
bool is_closed;
-
+ bool is_extruded;
+ double hLayer;
+ int nbElemLayer;
//all (unique) lines of centerlines
std::vector<MLine*> lines;
@@ -144,9 +146,12 @@ class Centerline : public Field{
//Create In and Outlet Planar Faces
void closeVolume();
+ //Create extruded wall
+ void extrudeWall();
+
// Cut the tubular structure with a disk
// perpendicular to the tubular structure
- void cutByDisk(SVector3 &pt, SVector3 &dir, double &maxRad);
+ bool cutByDisk(SVector3 &pt, SVector3 &dir, double &maxRad);
//create discrete faces
void createFaces();