From c5ccbdb583302ab28420de57fa5e917191b8bd09 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Wed, 21 Dec 2011 11:29:19 +0000
Subject: [PATCH] fixed bug meshBamg with compound edges
---
Geo/Curvature.cpp | 6 +-
Geo/GFaceCompound.cpp | 6 +-
Mesh/meshGFace.cpp | 5 -
Mesh/meshGFace.h | 5 +
Mesh/meshGFaceBamg.cpp | 120 ++++++++++++++----
benchmarks/2d/recombine.geo | 4 +-
benchmarks/stl/BifurcationRemeshCurvature.geo | 23 ++++
contrib/bamg/bamg-gmsh.cpp | 3 +-
8 files changed, 130 insertions(+), 42 deletions(-)
create mode 100644 benchmarks/stl/BifurcationRemeshCurvature.geo
diff --git a/Geo/Curvature.cpp b/Geo/Curvature.cpp
index 9a41dee7f4..6266b5aed9 100644
--- a/Geo/Curvature.cpp
+++ b/Geo/Curvature.cpp
@@ -963,9 +963,9 @@ void Curvature::buildEdgeList()
NbEdges += _VertexToEdgeList[i].size();
}
- std::cout << "Euler formula:" << std::endl;
- std::cout << "Edges + 2 = " << NbEdges + 2 << std::endl;
- std::cout << "Elements + Nodes = " << _VertexToInt.size() + _ElementToInt.size() << std::endl;
+ //std::cout << "Euler formula:" << std::endl;
+ //std::cout << "Edges + 2 = " << NbEdges + 2 << std::endl;
+ //std::cout << "Elements + Nodes = " << _VertexToInt.size() + _ElementToInt.size() << std::endl;
}
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 4d33e5f1e2..36bb476cd6 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -607,8 +607,8 @@ bool GFaceCompound::parametrize() const
if(oct) return paramOK;
if(trivial()) return paramOK;
- // if (_mapping != RBF)
- // coordinates.clear();
+ if (_mapping != RBF)
+ coordinates.clear();
computeNormals();
@@ -1681,7 +1681,7 @@ void GFaceCompound::getTriangle(double u, double v,
// }
// }
- *lt = (GFaceCompoundTriangle*)Octree_Search(uv, oct);
+ *lt = (GFaceCompoundTriangle*)Octree_Search(uv, oct);
// if(!(*lt)) {
// for(int i=0;i<nbT;i++){
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index ba6203e8c3..a33434abe0 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -433,11 +433,6 @@ void BDS2GMSH(BDS_Mesh *m, GFace *gf, std::map<BDS_Point*, MVertex*> &recoverMap
}
}
-bool meshGenerator(GFace *gf, int RECUR_ITER,
- bool repairSelfIntersecting1dMesh,
- bool onlyInitialMesh,
- bool debug = true,
- std::list<GEdge*> *replacement_edges = 0);
static void addOrRemove(MVertex *v1, MVertex *v2, std::set<MEdge,Less_Edge> & bedges)
{
diff --git a/Mesh/meshGFace.h b/Mesh/meshGFace.h
index 62fa4a61d0..2cb1d2338a 100644
--- a/Mesh/meshGFace.h
+++ b/Mesh/meshGFace.h
@@ -54,5 +54,10 @@ int MeshExtrudedSurface(GFace *gf, std::set<std::pair<MVertex*, MVertex*> >
*constrainedEdges=0);
void partitionAndRemesh(GFaceCompound *gf);
bool checkMeshCompound(GFaceCompound *gf, std::list<GEdge*> &edges);
+bool meshGenerator(GFace *gf, int RECUR_ITER,
+ bool repairSelfIntersecting1dMesh,
+ bool onlyInitialMesh,
+ bool debug = true,
+ std::list<GEdge*> *replacement_edges = 0);
#endif
diff --git a/Mesh/meshGFaceBamg.cpp b/Mesh/meshGFaceBamg.cpp
index a4463d3807..caff64c033 100644
--- a/Mesh/meshGFaceBamg.cpp
+++ b/Mesh/meshGFaceBamg.cpp
@@ -18,6 +18,10 @@
#include <map>
#include "BackgroundMesh.h"
#include "meshGFaceDelaunayInsertion.h"
+#include "Options.h"
+#include "meshGFace.h"
+#include "GFaceCompound.h"
+#include "MElementOctree.h"
#if defined(HAVE_BAMG)
@@ -69,11 +73,11 @@ static void computeMeshMetricsForBamg(GFace *gf, int numV,
void meshGFaceBamg(GFace *gf){
//Replace edges by their compounds
- std::list<GEdge*> edges = gf->edges();
- std::set<GEdge*> mySet;
- std::list<GEdge*>::iterator it = edges.begin();
- bool hasCompounds = false;
- while(it != edges.end()){
+ std::list<GEdge*> edges = gf->edges();
+ std::set<GEdge*> mySet;
+ std::list<GEdge*>::iterator it = edges.begin();
+ bool hasCompounds = false;
+ while(it != edges.end()){
if((*it)->getCompound()){
hasCompounds = true;
GEdge *gec = (GEdge*)(*it)->getCompound();
@@ -83,22 +87,16 @@ void meshGFaceBamg(GFace *gf){
mySet.insert(*it);
}
++it;
- }
- edges.clear();
- edges.insert(edges.begin(), mySet.begin(), mySet.end());
- std::set<MVertex*> bcVertex;
- for (std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); it++){
+ }
+ edges.clear();
+ edges.insert(edges.begin(), mySet.begin(), mySet.end());
+ std::set<MVertex*> bcVertex;
+ for (std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); it++){
for (unsigned int i = 0; i < (*it)->lines.size(); i++){
bcVertex.insert((*it)->lines[i]->getVertex(0));
bcVertex.insert((*it)->lines[i]->getVertex(1));
}
}
- //if has compound edges use a frontal mesher to remesh compound surface
- if (hasCompounds){
- printf("tris compound = %d \n", gf->triangles.size());
- bowyerWatsonFrontal(gf);
- printf("after tris compound = %d \n", gf->triangles.size());
- }
//fill mesh data fo bamg (bamgVertices, bamgTriangles, bamgBoundary)
std::set<MVertex*> all;
@@ -108,13 +106,11 @@ void meshGFaceBamg(GFace *gf){
all.insert(gf->triangles[i]->getVertex(j));
}
- //FIXME EMI: how should we define bamgVertices when we have remeshed compound edges not with lines not belonging to triangles of the compound face
- int numVertices = all.size()+bcVertex.size();
- Vertex2 *bamgVertices = new Vertex2[numVertices];
+ Vertex2 *bamgVertices = new Vertex2[all.size()];
int index = 0;
- // for(std::set<MVertex*>::iterator it = all.begin(); it!=all.end(); ++it){
- // if ((*it)->onWhat()->dim() <= 1){
- for(std::set<MVertex*>::iterator it = bcVertex.begin(); it!=bcVertex.end(); ++it){
+ for(std::set<MVertex*>::iterator it = all.begin(); it!=all.end(); ++it){
+ if ((*it)->onWhat()->dim() <= 1){
+ //for(std::set<MVertex*>::iterator it = bcVertex.begin(); it!=bcVertex.end(); ++it){
SPoint2 p;
bool success = reparamMeshVertexOnFace(*it, gf, p);
bamgVertices[index][0] = p.x();
@@ -122,7 +118,7 @@ void meshGFaceBamg(GFace *gf){
bamgVertices[index].lab = index;
recover[index] = *it;
(*it)->setIndex(index++);
- //}
+ }
}
int nbFixedVertices = index;
for(std::set<MVertex*>::iterator it = all.begin(); it!=all.end(); ++it){
@@ -136,6 +132,8 @@ void meshGFaceBamg(GFace *gf){
(*it)->setIndex(index++);
}
}
+
+ std::vector<MElement*> myParamElems;
Triangle2 *bamgTriangles = new Triangle2[gf->triangles.size()];
for (unsigned int i = 0; i < gf->triangles.size(); i++){
int nodes [3] = {gf->triangles[i]->getVertex(0)->getIndex(),
@@ -147,6 +145,13 @@ void meshGFaceBamg(GFace *gf){
double v1(bamgVertices[nodes[0]][1]);
double v2(bamgVertices[nodes[1]][1]);
double v3(bamgVertices[nodes[2]][1]);
+ if (hasCompounds){
+ MVertex *vv1 = new MVertex(u1,v1,0.0);
+ MVertex *vv2 = new MVertex(u2,v2,0.0);
+ MVertex *vv3 = new MVertex(u3,v3,0.0);
+ MTriangle *tri = new MTriangle(vv1,vv2,vv3, i);
+ myParamElems.push_back(tri);
+ }
double sign = (u2-u1)*(v3-v1) - (u3-u1)*(v2-v1);
if (sign < 0){
int temp = nodes[0];
@@ -155,7 +160,7 @@ void meshGFaceBamg(GFace *gf){
}
bamgTriangles[i].init(bamgVertices, nodes, gf->tag());
}
-
+
int numEdges = 0;
for (std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); ++it){
numEdges += (*it)->lines.size();
@@ -167,7 +172,7 @@ void meshGFaceBamg(GFace *gf){
for (unsigned int i = 0; i < (*it)->lines.size(); ++i){
int nodes [2] = {(*it)->lines[i]->getVertex(0)->getIndex(),
(*it)->lines[i]->getVertex(1)->getIndex()};
- bamgBoundary[count].init (bamgVertices, nodes, (*it)->tag());
+ bamgBoundary[count].init(bamgVertices, nodes, (*it)->tag());
bamgBoundary[count++].lab = count;
}
}
@@ -175,8 +180,36 @@ void meshGFaceBamg(GFace *gf){
Mesh2 *bamgMesh = new Mesh2 (all.size(), gf->triangles.size(), numEdges,
bamgVertices, bamgTriangles, bamgBoundary);
+ MElementOctree *_octree ;
+ if (hasCompounds){
+ _octree = new MElementOctree(myParamElems);
+ // //EMI PRINT TRIS
+ // FILE * fi = fopen("TRIS.pos","w");
+ // fprintf(fi, "View \"\"{\n");
+ // for( int i =0; i< gf->triangles.size(); i++){
+ // int nodes [3] = {gf->triangles[i]->getVertex(0)->getIndex(),
+ // gf->triangles[i]->getVertex(1)->getIndex(),
+ // gf->triangles[i]->getVertex(2)->getIndex()};
+ // double u1(bamgVertices[nodes[0]][0]);
+ // double u2(bamgVertices[nodes[1]][0]);
+ // double u3(bamgVertices[nodes[2]][0]);
+ // double v1(bamgVertices[nodes[0]][1]);
+ // double v2(bamgVertices[nodes[1]][1]);
+ // double v3(bamgVertices[nodes[2]][1]);
+ // fprintf(fi, "ST(%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,"
+ // "%22.15E,%22.15E){%d,%d,%d};\n",
+ // u1, v1, 0.0,
+ // u2, v2, 0.0,
+ // u3, v3, 0.0,
+ // i, i, i);
+ // }
+ // fprintf(fi,"};\n");
+ // fclose(fi);
+ // //END EMI PRINT TRIS
+ }
+
Mesh2 *refinedBamgMesh = 0;
- int iterMax = 11;
+ int iterMax = 1;
for (int k= 0; k < iterMax; k++){
int nbVert = bamgMesh->nv;
@@ -217,12 +250,42 @@ void meshGFaceBamg(GFace *gf){
Vertex2 &v = refinedBamgMesh->vertices[i];
if (i >= nbFixedVertices){
GPoint gp = gf->point(SPoint2(v[0], v[1]));
+ // if (gp.x() > 2.){ printf("wrong vertex index=%d %g %g %g (%g %g)\n", i, gp.x(), gp.y(), gp.z(), v[0], v[1]);
+ // }
+ //If point not found because compound edges have been remeshed and boundary triangles have changed
+ //then we call our new octree
+ if ( !gp.succeeded() && hasCompounds){
+ SPoint3 uvw(v[0],v[1], 0.0), UV;
+ double initialTol = MElement::getTolerance();
+ MElement::setTolerance(1.e-2);
+ MElement *e = _octree->find(v[0],v[1], 0.0, -1);
+ MElement::setTolerance(initialTol);
+ if (e){
+ e->xyz2uvw(uvw,UV);
+ double *valX = new double[e->getNumVertices()];
+ double *valY = new double[e->getNumVertices()];
+ double *valZ = new double[e->getNumVertices()];
+ for (int i=0;i<e->getNumVertices();i++){
+ int numTri = e->getNum();
+ valX[i] = gf->triangles[numTri]->getVertex(0)->x();
+ valY[i] = gf->triangles[numTri]->getVertex(0)->y();
+ valZ[i] = gf->triangles[numTri]->getVertex(0)->z();
+ }
+ gp.x() = e->interpolate(valX,UV[0],UV[1],UV[2]);
+ gp.y() = e->interpolate(valY,UV[0],UV[1],UV[2]);
+ gp.z() = e->interpolate(valZ,UV[0],UV[1],UV[2]);
+ delete [] valX;
+ delete [] valY;
+ delete [] valZ;
+ }
+ }
MFaceVertex *x = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, v[0], v[1]);
yetAnother[i] = x;
gf->mesh_vertices.push_back(x);
}
else {
- yetAnother[i] = recover[i];
+ MVertex *v = recover[i];
+ yetAnother[i] =v;
}
}
@@ -241,6 +304,9 @@ void meshGFaceBamg(GFace *gf){
}
if (refinedBamgMesh) delete refinedBamgMesh;
+ if ( _octree) delete _octree;
+ for(std::vector<MElement*>::iterator it = myParamElems.begin(); it != myParamElems.end(); it++)
+ delete *it;
}
#else
diff --git a/benchmarks/2d/recombine.geo b/benchmarks/2d/recombine.geo
index b0afab1e8b..7594af00aa 100644
--- a/benchmarks/2d/recombine.geo
+++ b/benchmarks/2d/recombine.geo
@@ -1,8 +1,8 @@
lc = 0.03;
Point(1) = {0, 0, 0, lc/5};
Point(2) = {.1, 0, 0, lc/5};
-Point(3) = {.1, .3, 0, lc};
-Point(4) = {0, .3, 0, lc};
+Point(3) = {.1, .1, 0, lc};
+Point(4) = {0, .1, 0, lc};
Line(1) = {1,2} ;
Line(2) = {3,2} ;
Line(3) = {3,4} ;
diff --git a/benchmarks/stl/BifurcationRemeshCurvature.geo b/benchmarks/stl/BifurcationRemeshCurvature.geo
new file mode 100644
index 0000000000..b70188a8a3
--- /dev/null
+++ b/benchmarks/stl/BifurcationRemeshCurvature.geo
@@ -0,0 +1,23 @@
+Mesh.Algorithm=7; //1=MeshAdapt, 5=Delaunay, 6=Frontal, 7=BAMG
+Mesh.RemeshParametrization=1; //0) harmonic (1) conformal
+Mesh.RemeshAlgorithm=1; //(0) nosplit (1) automatic
+
+Mesh.CharacteristicLengthFactor=0.35;
+Mesh.CharacteristicLengthFromCurvature=1; //-clcurv
+Mesh.CharacteristicLengthMin = 0.02; //-clmin 0.05
+Mesh.CharacteristicLengthMax = 4.00; //-clmax 4.0
+Mesh.LcIntegrationPrecision=1.e-5; //-epslc1d
+Mesh.MinimumCirclePoints=30; //default=7
+Mesh.CharacteristicLengthFromPoints = 0;
+Mesh.CharacteristicLengthExtendFromBoundary=0;
+
+Merge "bifurcation.stl";
+
+CreateTopology;
+
+Compound Surface(20) = {1};
+
+
+Physical Surface("Wall") = {20};
+
+
diff --git a/contrib/bamg/bamg-gmsh.cpp b/contrib/bamg/bamg-gmsh.cpp
index 999226ceed..4fedb9ab8c 100644
--- a/contrib/bamg/bamg-gmsh.cpp
+++ b/contrib/bamg/bamg-gmsh.cpp
@@ -419,8 +419,7 @@ Mesh2 *Bamg(Mesh2 *Thh, double * args,double *mm11,double *mm12,double *mm22, bo
Triangles &Th(*oTh);
// printf("COUCOUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU\n");
- // Th.Write("toto.mesh",bamg::Triangles::AutoMesh);
-
+ //Th.Write("toto.mesh",bamg::Triangles::AutoMesh);
bool mtx=true;
KN_<double> m11(mm11,Thh->nv);
--
GitLab