diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 86309a383197ef79a292f08af560c11c41386778..8bbc7cef22b5ab8bb6aa11a4843c3ff4368b7563 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -582,8 +582,8 @@ void GFaceCompound::one2OneMap() const
bool badCavity = closedCavity(v,vTri) ? checkCavity(vTri, vCoord) : false;
if(badCavity){
- Msg::Debug("Wrong cavity around vertex %d (onwhat=%d).",
- v->getNum(), v->onWhat()->dim());
+ Msg::Debug("Wrong cavity around vertex %d .",
+ v->getNum());
Msg::Debug("--> Place vertex at center of gravity of %d-Polygon kernel." ,
vTri.size());
@@ -638,6 +638,10 @@ bool GFaceCompound::parametrize() const
Msg::Debug("Parametrizing surface %d with 'conformal map'", tag());
fillNeumannBCS();
bool hasOverlap = parametrize_conformal_spectral();
+ if (hasOverlap || !checkOrientation(0) ){
+ Msg::Warning("!!! Overlap or Flipping: parametrization switched to 'FE conformal' map");
+ hasOverlap = parametrize_conformal(0, NULL, NULL);
+ }
if (hasOverlap || !checkOrientation(0) ){
printStuff(33);
Msg::Warning("$$$ Overlap or Flipping: parametrization switched to 'harmonic' map");
@@ -1798,9 +1802,6 @@ void GFaceCompound::buildOct() const
bool GFaceCompound::checkTopology() const
{
if (_mapping == RBF) return true;
-
- //fixme tristan
- //return true;
bool correctTopo = true;
if(allNodes.empty()) buildAllNodes();
@@ -1832,7 +1833,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: Genus=%d, Nb boundaries=%d, AR=%d ", G, Nb, AR);
if (_allowPartition == 1){
nbSplit = -2;
Msg::Info("-----------------------------------------------------------");
@@ -1862,8 +1863,6 @@ bool GFaceCompound::checkTopology() const
double GFaceCompound::checkAspectRatio() const
{
- //if ((*(_compound.begin()))->geomType() != GEntity::DiscreteSurface)
- // return true;
if(allNodes.empty()) buildAllNodes();
@@ -1898,19 +1897,13 @@ double GFaceCompound::checkAspectRatio() const
std::list<GEdge*>::const_iterator it0 = _U0.begin();
double tot_length = 0;
- for( ; it0 != _U0.end(); ++it0 ){
- for(unsigned int i = 0; i < (*it0)->lines.size(); i++ ){
- MVertex *v0 = (*it0)->lines[i]->getVertex(0);
- MVertex *v1 = (*it0)->lines[i]->getVertex(1);
- const double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
- (v0->y() - v1->y()) * (v0->y() - v1->y()) +
- (v0->z() - v1->z()) * (v0->z() - v1->z()));
- tot_length += length;
- }
- }
- double AR = 2*3.14*area3D/(tot_length*tot_length);
+ for( ; it0 != _U0.end(); ++it0 )
+ for(unsigned int i = 0; i < (*it0)->lines.size(); i++ )
+ tot_length += (*it0)->lines[i]->getInnerRadius();
+
+ double AR = M_PI*area3D/(tot_length*tot_length);
- if (areaMin > 0 && areaMin < limit && nb > 2) {
+ if (areaMin > 0 && areaMin < limit && nb > 3) {
Msg::Warning("Too small triangles in mapping (a_2D=%g)", areaMin);
}
else {
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 4c1a79ec2a410669a7673bfcb8a990cbf5908760..c69d79b201a785a5d99dfe77819b1831c8ce45bb 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -27,6 +27,10 @@
#include "GFace.h"
#include "discreteEdge.h"
#include "discreteFace.h"
+#include "GEdgeCompound.h"
+#include "GFaceCompound.h"
+#include "meshGFace.h"
+#include "meshGEdge.h"
#if defined(HAVE_ANN)
#include <ANN/ANN.h>
@@ -644,6 +648,144 @@ void Centerline::buildKdTree(){
}
+void Centerline::remeshSplitMesh(){
+
+ // Remesh new faces (Compound Lines and Compound Surfaces)
+ Msg::Info("*** Starting parametrize compounds:");
+ double t0 = Cpu();
+
+ //Parametrize Compound Lines
+ int NE = split->getMaxElementaryNumber(1);
+ for (int i=0; i < NE; i++){
+ std::vector<GEdge*>e_compound;
+ GEdge *pe = split->getEdgeByTag(i+1);//split edge
+ e_compound.push_back(pe);
+ int num_gec = NE+i+1;
+ Msg::Info("Parametrize Compound Line (%d) = %d discrete edge",
+ num_gec, pe->tag());
+ GEdgeCompound *gec = new GEdgeCompound(split, num_gec, e_compound);
+ split->add(gec);
+ gec->parametrize();
+ }
+
+ // Parametrize Compound surfaces
+ std::set<MVertex*> allNod;
+ std::list<GEdge*> U0;
+ int NF = split->getMaxElementaryNumber(2);
+ for (int i=0; i < NF; i++){
+ std::list<GFace*> f_compound;
+ GFace *pf = split->getFaceByTag(i+1);//split face
+ f_compound.push_back(pf);
+ int num_gfc = NF+i+1;
+ Msg::Info("Parametrize Compound Surface (%d) = %d discrete face",
+ num_gfc, pf->tag());
+ GFaceCompound::typeOfMapping typ = GFaceCompound::CONFORMAL;
+ GFaceCompound *gfc = new GFaceCompound(split, num_gfc, f_compound, U0,
+ typ, 0);
+ int recombine = 0;//EMI TODO RECUPERER LE RECOMBINE
+ gfc->meshAttributes.recombine = recombine;
+ split->add(gfc);
+ gfc->parametrize();
+ }
+ double t1 = Cpu();
+ Msg::Info("*** Parametrize compounds done (%g s)", t1-t0);
+
+ //set centerline field
+ FieldManager *fields = split->getFields();
+ fields->reset();
+ int id = fields->newId();
+ (*fields)[id] = this;
+ fields->background_field = id;
+
+ Msg::Info("*** Starting meshing 1D edges ...:");
+ for (int i = 0; i < NE; i++){
+ GEdge *gec = split->getEdgeByTag(NE+i+1);
+ meshGEdge mge;
+ mge(gec);
+ }
+ double t2 = Cpu();
+ Msg::Info("*** Meshing 1D edges done (%gs)", t2-t1);
+
+ Msg::Info("*** Starting mesh of surface ");
+ for (int i=0; i < NF; i++){
+ GFace *gfc = split->getFaceByTag(NF+i+1);
+ meshGFace mgf;
+ mgf(gfc);
+ }
+
+ // // Removing discrete Vertices - Edges - Faces
+ // int NV = split->getMaxElementaryNumber(0) - numv + 1;
+ // for (int i=0; i < NV; i++){
+ // GVertex *pv = gf->model()->getVertexByTag(numv+i);
+ // gf->model()->remove(pv);
+ // }
+ // for (int i=0; i < NE; i++){
+ // GEdge *gec = split->getEdgeByTag(nume+NE+i);
+ // GEdge *pe = split->getEdgeByTag(nume+i);
+ // gf->model()->remove(pe);
+ // gf->model()->remove(gec);
+ // }
+ // for (int i=0; i < NF; i++){
+ // GFace *gfc = split->getFaceByTag(numf+NF+i);
+ // GFace *pf = split->getFaceByTag(numf+i);
+ // gf->model()->remove(pf);
+ // gf->model()->remove(gfc);
+ // }
+
+ // // Put new mesh in a new discreteFace
+ // for(std::set<MVertex*>::iterator it = allNod.begin(); it != allNod.end(); ++it){
+ // gf->mesh_vertices.push_back(*it);
+ // }
+
+ // // 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);
+ // }
+ // }
+
+ // double t3 = Cpu();
+ // Msg::Info("*** Mesh of surface %d done by assembly %d remeshed faces (%g s)",
+ // gf->tag(), NF, t3-t2);
+ // Msg::Info("-----------------------------------------------------------");
+
+ // gf->coherenceNormals();
+ // gf->meshStatistics.status = GFace::DONE;
+
+}
void Centerline::createFaces(){
std::vector<std::vector<MTriangle*> > faces;
@@ -715,11 +857,8 @@ void Centerline::createFaces(){
myVertices.begin(), myVertices.end());
}
- // printf("%d discrete Faces (bef =%d) (after =%d)\n", numAfter-numBef, numBef-1, numAfter-1);
-
}
-
void Centerline::splitMesh(){
Msg::Info("Splitting surface mesh (%d tris) with centerline field ", triangles.size());
@@ -733,7 +872,7 @@ void Centerline::splitMesh(){
double AR = L/R;
printf("*** branch =%d \n", i, AR);
if( AR > 8.){
- int nbSplit = (int)ceil(AR/4);
+ int nbSplit = (int)ceil(AR/3.5);
double li = L/nbSplit;
double lc = 0.0;
for (int j= 0; j < lines.size(); j++){
@@ -771,6 +910,10 @@ void Centerline::splitMesh(){
Msg::Info("Writing splitted mesh 'myPARTS.msh'");
split->writeMSH("myPARTS.msh", 2.2, false, true);
+ //remesh splitted mesh
+ remeshSplitMesh();
+ split->writeMSH("myMESH.msh", 2.2, false, true);
+
//print
// FILE * f2 = fopen("myCUTLINES.pos","w");
// fprintf(f2, "View \"\"{\n");
@@ -801,7 +944,7 @@ 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);
- const double EPS = 0.005;
+ 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++)
@@ -809,7 +952,7 @@ void Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
}
bool closedCut = false;
int step = 0;
- while (!closedCut && step < 10){
+ while (!closedCut && step < 20){
double rad = 1.3*maxRad+0.15*step*maxRad;
std::map<MEdge,MVertex*,Less_Edge> cutEdges;
std::set<MVertex*> cutVertices;
@@ -850,7 +993,7 @@ void Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
break;
}
else {
- if (step ==9) {printf("no closed cut %d \n", (int)newCut.size()); };
+ if (step ==19) {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 5b5e3c86add3794c5ecacaa0a61ea20f31f817b3..e8c51669f800a4b97bebb64af66a53cfc3b21cd4 100644
--- a/Mesh/CenterlineField.h
+++ b/Mesh/CenterlineField.h
@@ -133,6 +133,7 @@ class Centerline : public Field{
//create discrete faces
void createFaces();
+ void remeshSplitMesh();
//Print for debugging
void printSplit() const;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 8a4f82a191a85a5959ad93e103f7c1ed8e7f1f08..8cb6b54a1dbe74306fe3cbea633037398b6cc655 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1958,7 +1958,6 @@ void partitionAndRemesh(GFaceCompound *gf)
double t1 = Cpu();
Msg::Info("*** Parametrize compounds done (%g s)", t1-t0);
-
Msg::Info("*** Starting meshing 1D edges ...:");
for (int i = 0; i < NE; i++){
GEdge *gec = gf->model()->getEdgeByTag(nume + NE + i);
@@ -1970,7 +1969,6 @@ void partitionAndRemesh(GFaceCompound *gf)
Msg::Info("*** Starting Mesh of surface %d ...", gf->tag());
- // lloydAlgorithm
for (int i=0; i < NF; i++){
GFace *gfc = gf->model()->getFaceByTag(numf + NF + i );
meshGFace mgf;