From 522567f5049d01594fec7b54671d55657e00dcea Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Thu, 1 Mar 2012 14:20:54 +0000
Subject: [PATCH] more work on centerline and GFaceCompounds
---
Geo/GFace.cpp | 5 +-
Geo/GFaceCompound.cpp | 99 +++--------------
Geo/GFaceCompound.h | 2 -
Geo/GModelFactory.cpp | 10 +-
Geo/GModelIO_Geo.cpp | 1 -
Geo/GeoInterpolation.cpp | 7 +-
Geo/MLine.cpp | 5 +
Geo/MLine.h | 3 +-
Mesh/CenterlineField.cpp | 162 +++++++---------------------
Mesh/CenterlineField.h | 1 -
Mesh/meshGFaceDelaunayInsertion.cpp | 4 +-
Numeric/fullMatrix.cpp | 4 +-
12 files changed, 73 insertions(+), 230 deletions(-)
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 8904bda21c..16ca3dafb5 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -106,10 +106,7 @@ void GFace::delFreeEdge(GEdge *e)
void GFace::deleteMesh()
{
- 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];
- }
+ for(unsigned int i = 0; i < mesh_vertices.size(); i++) delete mesh_vertices[i];
mesh_vertices.clear();
transfinite_vertices.clear();
for(unsigned int i = 0; i < triangles.size(); i++) delete triangles[i];
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 579a908364..3e4882a637 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -645,13 +645,13 @@ bool GFaceCompound::parametrize() const
// Laplace parametrization
if (_mapping == HARMONIC){
- Msg::Info("Parametrizing surface %d with 'zut harmonic map'", tag());
+ Msg::Info("Parametrizing surface %d with harmonic map'", tag());
fillNeumannBCS();
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
- //parametrize(ITERU,CONVEXCOMBINATION);
- // parametrize(ITERV,CONVEXCOMBINATION);
printStuff(111);
+ //parametrize(ITERU,CONVEXCOMBINATION);
+ //parametrize(ITERV,CONVEXCOMBINATION);
if (_type == MEANPLANE) checkOrientation(0, true);
printStuff(222);
}
@@ -662,24 +662,18 @@ bool GFaceCompound::parametrize() const
std::vector<MVertex *> vert;
bool oriented, hasOverlap;
hasOverlap = parametrize_conformal_spectral();
- printStuff(11);
if (hasOverlap) oriented = checkOrientation(0);
else oriented = checkOrientation(0, true);
- printStuff(22);
hasOverlap = checkOverlap(vert);
if ( !oriented || hasOverlap ){
Msg::Warning("!!! parametrization switched to 'FE conformal' map");
parametrize_conformal(0, NULL, NULL);
- printStuff(33);
checkOrientation(0, true);
- printStuff(44);
}
if (!checkOrientation(0) || checkOverlap(vert)){
- printStuff(55);
Msg::Warning("$$$ parametrization switched to 'harmonic' map");
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
- printStuff(66);
}
}
// Radial-Basis Function parametrization
@@ -715,12 +709,6 @@ bool GFaceCompound::parametrize() const
}
}
- double AR = checkAspectRatio();
- if (floor(AR) > AR_MAX){
- Msg::Warning("Geometrical aspect ratio is high AR=%d ", (int)AR);
- paramOK = true; //false;
- }
-
return paramOK;
}
@@ -790,87 +778,30 @@ void GFaceCompound::getBoundingEdges()
double GFaceCompound::getSizeH() const
{
SBoundingBox3d bb;
- std::vector<SPoint3> vertices;
for(std::set<MVertex *>::iterator itv = allNodes.begin();
itv != allNodes.end() ; ++itv){
SPoint3 pt((*itv)->x(),(*itv)->y(), (*itv)->z());
- vertices.push_back(pt);
bb +=pt;
}
+
double H = norm(SVector3(bb.max(), bb.min()));
-
- //SOrientedBoundingBox obbox = SOrientedBoundingBox::buildOBB(vertices);
- //double H = obbox.getMaxSize();
return H;
}
double GFaceCompound::getSizeBB(const std::list<GEdge* > &elist) const
{
- //SOrientedBoundingBox obboxD = obb_boundEdges(elist);
- //double D = obboxD.getMaxSize();
- SBoundingBox3d bboxD = boundEdges(elist);
+ SBoundingBox3d bboxD;
+ std::list<GEdge*>::const_iterator it = elist.begin();
+ for(; it != elist.end(); it++){
+ bboxD += (*it)->bounds();
+ }
double D = norm(SVector3(bboxD.max(), bboxD.min()));
return D;
}
-SBoundingBox3d GFaceCompound::boundEdges(const std::list<GEdge* > &elist) const
-{
- SBoundingBox3d res;
- std::list<GEdge*>::const_iterator it = elist.begin();
- for(; it != elist.end(); it++)
- res += (*it)->bounds();
-
- return res;
-}
-
-SOrientedBoundingBox GFaceCompound::obb_boundEdges(const std::list<GEdge* > &elist) const
-{
- SOrientedBoundingBox res;
- std::vector<SPoint3> vertices;
-
- std::list<GEdge*>::const_iterator it = elist.begin();
- for(; it != elist.end(); it++) {
-
- if((*it)->getNumMeshVertices() > 0) {
- int N = (*it)->getNumMeshVertices();
- for (int i = 0; i < N; i++) {
- MVertex* mv = (*it)->getMeshVertex(i);
- vertices.push_back(mv->point());
- }
- // Don't forget to add the first and last vertices...
- SPoint3 pt1((*it)->getBeginVertex()->x(),
- (*it)->getBeginVertex()->y(), (*it)->getBeginVertex()->z());
- SPoint3 pt2((*it)->getEndVertex()->x(), (*it)->getEndVertex()->y(),
- (*it)->getEndVertex()->z());
- vertices.push_back(pt1);
- vertices.push_back(pt2);
- }
- else if((*it)->geomType() != DiscreteCurve && (*it)->geomType() !=
- BoundaryLayerCurve){
- Range<double> tr = (*it)->parBounds(0);
- // N can be choosen arbitrarily, but 10 points seems reasonable
- int N = 10;
- for (int i = 0; i < N; i++) {
- double t = tr.low() + (double)i / (double)(N - 1) * (tr.high() -
- tr.low());
- GPoint p = (*it)->point(t);
- SPoint3 pt(p.x(), p.y(), p.z());
- vertices.push_back(pt);
- }
- }
- else {
- SPoint3 dummy(0, 0, 0);
- vertices.push_back(dummy);
- }
-
- }
- res = SOrientedBoundingBox::buildOBB(vertices);
- return res;
-}
-
void GFaceCompound::getUniqueEdges(std::set<GEdge*> &_unique)
{
_unique.clear();
@@ -1937,7 +1868,7 @@ bool GFaceCompound::checkTopology() const
double D = H;
if (_interior_loops.size() > 0) D = getSizeBB(_U0);
int AR1 = (int) checkAspectRatio();
- int AR2 = (int) ceil(H/D);
+ int AR2 = (int) round(H/D);
int AR = std::max(AR1, AR2);
if (G != 0 || Nb < 1){
@@ -1957,7 +1888,7 @@ bool GFaceCompound::checkTopology() const
}
else if (G == 0 && AR > AR_MAX){
correctTopo = false;
- Msg::Warning("Wrong topology: Genus=%d, Nb boundaries=%d, AR=%d ", G, Nb, AR);
+ Msg::Warning("Wrong topology: Aspect ratio is too high AR=%d (AR1=%d AR2=%d)", AR, AR1, AR2);
if (_allowPartition == 1){
nbSplit = -2;
Msg::Info("-----------------------------------------------------------");
@@ -1971,11 +1902,11 @@ bool GFaceCompound::checkTopology() const
tag(), nbSplit);
}
else if (_allowPartition == 0){
- Msg::Warning("The geometrical aspect ratio of your geometry is quite high.\n "
+ Msg::Debug("The geometrical aspect ratio of your geometry is quite high.\n "
"You should enable partitioning of the mesh by activating the\n"
"automatic remeshing algorithm. Add 'Mesh.RemeshAlgorithm=1;'\n "
"in your geo file or through the Fltk window (Options > Mesh >\n "
- "General) \n");
+ "General)");
}
}
else{
@@ -2020,10 +1951,10 @@ double GFaceCompound::checkAspectRatio() const
}
std::list<GEdge*>::const_iterator it0 = _U0.begin();
- double tot_length = 0;
+ double tot_length = 0.0;
for( ; it0 != _U0.end(); ++it0 )
for(unsigned int i = 0; i < (*it0)->lines.size(); i++ )
- tot_length += (*it0)->lines[i]->getInnerRadius();
+ tot_length += (*it0)->lines[i]->getLength();
double AR = M_PI*area3D/(tot_length*tot_length);
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index f3d898e4c3..7684deb6e1 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -103,8 +103,6 @@ class GFaceCompound : public GFace {
bool trivial() const;
double getSizeH() const;
double getSizeBB(const std::list<GEdge* > &elist) const;
- SBoundingBox3d boundEdges(const std::list<GEdge* > &elist) const;
- SOrientedBoundingBox obb_boundEdges(const std::list<GEdge* > &elist) const;
void fillNeumannBCS() const;
/* double sumAngles(std::vector<MVertex*> ordered) const; */
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 6ee3500930..3f3b477c5d 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -11,6 +11,7 @@
#include "gmshEdge.h"
#include "gmshFace.h"
#include "gmshRegion.h"
+#include "GEdgeCompound.h"
#include "MLine.h"
#include "GModel.h"
#include "Numeric.h"
@@ -94,13 +95,17 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
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 if(ge->geomType() == GEntity::CompoundCurve){
+ c = Create_Curve(numEdge, MSH_SEGM_COMPOUND, 1, NULL, NULL, -1, -1, 0., 1.);
+ std::vector<GEdge*> gec = ((GEdgeCompound*)ge)->getCompounds();
+ for(int i = 0; i < gec.size(); i++)
+ c->compound.push_back(gec[i]->tag());
+ }
else {
printf("type not implemented \n"); exit(1);
}
@@ -152,7 +157,6 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
GRegion* GeoFactory::addVolume (GModel *gm, std::vector<std::vector<GFace *> > faces)
{
- printf("add volume \n");
//create surface loop
int nLoops = faces.size();
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index 49b21febce..f3b8fab3fa 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -164,7 +164,6 @@ 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/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index c93a312a8b..f5c06d0c49 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -418,12 +418,13 @@ Vertex InterpolateCurve(Curve *c, double u, int derivee)
break;
case MSH_SEGM_DISCRETE:
-#if defined(HAVE_BFGS)
- // printf("MSH_SEGM_DISCRETE\n");
-#endif
Msg::Debug("Cannot interpolate discrete curve");
break;
+ case MSH_SEGM_COMPOUND:
+ Msg::Debug("Cannot interpolate compound curve");
+ break;
+
default:
Msg::Error("Unknown curve type in interpolation");
break;
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index 51063a5125..5091bc18cc 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -60,3 +60,8 @@ double MLine::getInnerRadius()
{
return _v[0]->distance(_v[1]) * .5;
}
+
+double MLine::getLength()
+{
+ return _v[0]->distance(_v[1]);
+}
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 447c613160..53e77dd72a 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -37,7 +37,8 @@ class MLine : public MElement {
virtual int getDim() const { return 1; }
virtual int getNumVertices() const { return 2; }
virtual MVertex *getVertex(int num){ return _v[num]; }
- virtual double getInnerRadius(); // length of segment line
+ virtual double getInnerRadius(); // half-length of segment line
+ virtual double getLength(); // length of segment line
virtual void getVertexInfo(const MVertex * vertex, int &ithVertex) const
{
ithVertex = _v[0] == vertex ? 0 : 1;
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index cc32bb024d..dc654f23ae 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -48,53 +48,11 @@ double computeLength(std::vector<MLine*> lines){
double length= 0.0;
for (int i = 0; i< lines.size(); i++){
- length += lines[i]->getInnerRadius();
+ length += lines[i]->getLength();
}
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){
std::multiset<MVertex*> boundV;
@@ -645,6 +603,7 @@ void Centerline::buildKdTree(){
void Centerline::createSplitCompounds(){
+ //number of discrete vertices, edges, faces and regions for cut mesh
NV = current->getMaxElementaryNumber(0);
NE = current->getMaxElementaryNumber(1);
NF = current->getMaxElementaryNumber(2);
@@ -664,7 +623,6 @@ void Centerline::createSplitCompounds(){
num_gec, pe->tag());
GEdgeCompound *gec = new GEdgeCompound(current, num_gec, e_compound);
current->add(gec);
- //gec->parametrize();
}
// Parametrize Compound surfaces
@@ -676,19 +634,13 @@ void Centerline::createSplitCompounds(){
int num_gfc = NF+i+1;
Msg::Info("Parametrize Compound Surface (%d) = %d discrete face",
num_gfc, pf->tag());
- //GFaceCompound::typeOfMapping typ = GFaceCompound::HARMONIC;
- GFaceCompound::typeOfMapping typ = GFaceCompound::CONFORMAL;
+ GFaceCompound::typeOfMapping typ = GFaceCompound::HARMONICPLANE;
+ //GFaceCompound::typeOfMapping typ = GFaceCompound::CONFORMAL;
GFaceCompound *gfc = new GFaceCompound(current, num_gfc, f_compound, U0,
typ, 0);
gfc->meshAttributes.recombine = recombine;
- if (i < discFaces.size()){
- gfc->addPhysicalEntity(100);
- current->setPhysicalName("newsurf", 2, 100);
- }
- else{
- gfc->addPhysicalEntity(200);
- current->setPhysicalName("in/out", 2, 200);
- }
+ gfc->addPhysicalEntity(100);
+ current->setPhysicalName("newsurf", 2, 100);
current->add(gfc);
}
@@ -696,8 +648,11 @@ void Centerline::createSplitCompounds(){
void Centerline::cleanMesh(){
- return;
- if (!is_cut) return;
+
+ return; // does not work yet !
+ ////////////////////////////////
+
+ if (!is_cut || !is_closed) return;
for (int i=0; i < NF; i++){
std::ostringstream oss;
@@ -738,33 +693,6 @@ void Centerline::cleanMesh(){
mySplitMesh->quadrangles.push_back(new MQuadrangle(v[0], v[1], v[2], v[3]));
}
}
-
- // 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
for (int i=0; i < NV; i++){
@@ -789,19 +717,8 @@ void Centerline::cleanMesh(){
//Put new mesh in a new discreteFace
for(std::set<MVertex*>::iterator it = allNod.begin(); it != allNod.end(); ++it)
mySplitMesh->mesh_vertices.push_back(*it);
-
- removeBoundVertices(mySplitMesh);
mySplitMesh->meshStatistics.status = GFace::DONE;
- // 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();
}
@@ -833,7 +750,7 @@ void Centerline::createFaces(){
it != touched.end(); ++it)
e2e.erase(*it);
}
- printf("%d faces created \n", (int)faces.size());
+ Msg::Info("Centerline action (cutMesh) has cut surface mesh in %d faces ", (int)faces.size());
//create discFaces
int numBef = current->getMaxElementaryNumber(2) + 1;
@@ -859,12 +776,11 @@ void Centerline::createFaces(){
}
-void Centerline::createInOutFaces(){
+void Centerline::createClosedVolume(){
//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++){
+ for (int i= 0; i< NF; 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++){
@@ -873,40 +789,35 @@ void Centerline::createInOutFaces(){
else boundEdges.push_back(*it);
}
}
- printf("boundEdges size =%d \n", boundEdges.size());
-
- //create the inlet and outlet planar face
+ //printf("boundEdges size =%d \n", boundEdges.size());
current->setFactory("Gmsh");
- //std::vector<std::vector<GFace *> > myFaceLoops;
- //std::vector<GFace *> myFaces;
+ std::vector<std::vector<GFace *> > myFaceLoops;
+ std::vector<GFace *> myFaces;
for (int i = 0; i< boundEdges.size(); i++){
std::vector<std::vector<GEdge *> > myEdgeLoops;
std::vector<GEdge *> myEdges;
- myEdges.push_back(boundEdges[i]);
+ GEdge * gec = current->getEdgeByTag(NE+boundEdges[i]->tag());
+ myEdges.push_back(gec);
myEdgeLoops.push_back(myEdges);
GFace *newFace = current->addPlanarFace(myEdgeLoops);
- //myFaces.push_back(newFace);
+ newFace->addPhysicalEntity(200);
+ current->setPhysicalName("in/out", 2, 200);
+ 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;
+ Msg::Info("Centerline action (closeVolume) has created %d in/out planar faces ", (int)boundEdges.size());
+
for (int i = 0; i < NF; i++){
GFace * gf = current->getFaceByTag(NF+i+1);
- myFaces.push_back(gf);
+ myFaces.push_back(gf);
}
myFaceLoops.push_back(myFaces);
GRegion *reg = current->addVolume(myFaceLoops);
reg->addPhysicalEntity(reg->tag());
current->setPhysicalName("newvol", 3, reg->tag());
+ Msg::Info("Centerline action (closeVolume) has created volume %d ", reg->tag());
+
}
void Centerline::closeVolume(){
@@ -948,15 +859,15 @@ 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 = edges[i].minRad; //0.5*(edges[i].minRad+edges[i].maxRad);
- double AR = L/R;
+ double D = (edges[i].minRad+edges[i].maxRad);
+ double AR = L/D;
printf("*** branch =%d AR=%g \n", i, AR);
- if( AR > 4.5){
- int nbSplit = (int)ceil(AR/3.0);
+ if( AR > 4.0){
+ int nbSplit = (int)round(AR/3.);
double li = L/nbSplit;
double lc = 0.0;
for (int j= 0; j < lines.size(); j++){
- lc += lines[j]->getInnerRadius();
+ lc += lines[j]->getLength();
if (lc > li && nbSplit > 1) {
MVertex *v1 = lines[j]->getVertex(0);
MVertex *v2 = lines[j]->getVertex(1);
@@ -969,10 +880,10 @@ void Centerline::cutMesh(){
}
}
}
- if(edges[i].children.size() > 0.0 && AR > 1.5){
+ if(edges[i].children.size() > 0.0 && AR > 1.0){
MVertex *v1 = lines[lines.size()-1]->getVertex(1);//end vertex
MVertex *v2;
- if(L/R < 2.0) v2 = lines[0]->getVertex(0);
+ if(AR < 1.5) 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());
@@ -986,11 +897,10 @@ void Centerline::cutMesh(){
createFaces();
current->createTopologyFromFaces(discFaces);
current->exportDiscreteGEOInternals();
- if (is_closed) createInOutFaces();
-
+
//write
Msg::Info("Writing splitted mesh 'myPARTS.msh'");
- current->writeMSH("myPARTS.msh", 2.2, false, true);
+ current->writeMSH("myPARTS.msh", 2.2, false, false);
//create compounds
createSplitCompounds();
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
index 2637d899ac..f43e6241fd 100644
--- a/Mesh/CenterlineField.h
+++ b/Mesh/CenterlineField.h
@@ -146,7 +146,6 @@ class Centerline : public Field{
//create discrete faces
void createFaces();
- void createInOutFaces();
void createClosedVolume();
void createSplitCompounds();
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 028c3730ae..ac31bd366b 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1137,11 +1137,9 @@ void optimalPointFrontalB (GFace *gf,
newPoint[1] = uvt(1);
return;
}
- printf("--- NOT CONVERGED ----------\n");
+ //printf("--- NOT CONVERGED ----------\n");
}
-
-
void bowyerWatsonFrontal(GFace *gf)
{
std::set<MTri3*,compareTri3Ptr> AllTris;
diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index 9ed3656841..4b716034d4 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -286,9 +286,9 @@ bool fullMatrix<double>::luSolve(const fullVector<double> &rhs, fullVector<doubl
delete [] ipiv;
if(info == 0) return true;
if(info > 0)
- Msg::Error("U(%d,%d)=0 in LU decomposition", info, info);
+ Msg::Debug("U(%d,%d)=0 in LU decomposition", info, info);
else
- Msg::Error("Wrong %d-th argument in LU decomposition", -info);
+ Msg::Debug("Wrong %d-th argument in LU decomposition", -info);
return false;
}
--
GitLab