diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 82a54fa1816e5d0f39bf85cdc262f4dd2110a018..6256adf0dc0fce1c07b24de3e3a7ab2d5ea62771 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -483,14 +483,6 @@ BDS_Face *BDS_Mesh::add_triangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3)
return t;
}
-BDS_Face *BDS_Mesh::add_quadrangle(BDS_Edge *e1, BDS_Edge *e2,
- BDS_Edge *e3, BDS_Edge *e4)
-{
- BDS_Face *t = new BDS_Face(e1, e2, e3, e4);
- triangles.push_back(t);
- return t;
-}
-
void BDS_Mesh::del_face(BDS_Face *t)
{
t->e1->del(t);
@@ -624,46 +616,6 @@ BDS_Mesh::~BDS_Mesh()
DESTROOOY(triangles.begin(), triangles.end());
}
-bool BDS_Mesh::split_face(BDS_Face *f, BDS_Point *mid)
-{
- BDS_Point *p1 = f->e1->commonvertex(f->e2);
- BDS_Point *p2 = f->e3->commonvertex(f->e2);
- BDS_Point *p3 = f->e1->commonvertex(f->e3);
- BDS_Edge *p1_mid = new BDS_Edge(p1, mid);
- edges.push_back(p1_mid);
- BDS_Edge *p2_mid = new BDS_Edge(p2, mid);
- edges.push_back(p2_mid);
- BDS_Edge *p3_mid = new BDS_Edge(p3, mid);
- edges.push_back(p2_mid);
- BDS_Face *t1, *t2, *t3;
- t1 = new BDS_Face(f->e1, p1_mid, p3_mid);
- t2 = new BDS_Face(f->e2, p2_mid, p1_mid);
- t3 = new BDS_Face(f->e3, p3_mid, p2_mid);
-
- t1->g = f->g;
- t2->g = f->g;
- t3->g = f->g;
-
- p1_mid->g = f->g;
- p2_mid->g = f->g;
- p3_mid->g = f->g;
-
- mid->g = f->g;
-
- triangles.push_back(t1);
- triangles.push_back(t2);
- triangles.push_back(t3);
-
- // config has changed
- p1->config_modified = true;
- p2->config_modified = true;
- p3->config_modified = true;
-
- // delete the face
- del_face(f);
-
- return true;
-}
bool BDS_Mesh::split_edge(BDS_Edge *e, BDS_Point *mid)
{
@@ -1273,7 +1225,7 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
if (!gp.succeeded()){
return false;
}
- // if (!gf->containsParam(SPoint2(U,V)))return false;
+ // if (!gf->containsParam(SPoint2(U,V)))return false;
const double oldX = p->X;
const double oldY = p->Y;
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index c6aa237557710a92f97600a4ebbb39f97bfff355..2409a42cc514ef12e7b8fbaa4ab557b4db1fd838 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -403,12 +403,9 @@ struct EdgeToRecover
class BDS_Mesh
{
public:
- int MAXPOINTNUMBER, SNAP_SUCCESS, SNAP_FAILURE;
+ int MAXPOINTNUMBER;
double Min[3], Max[3], LC;
BDS_Mesh(int _MAXX = 0) : MAXPOINTNUMBER(_MAXX){}
- void load(GVertex *gv); // load in BDS all the meshes of the vertex
- void load(GEdge *ge); // load in BDS all the meshes of the edge
- void load(GFace *gf); // load in BDS all the meshes of the surface
virtual ~BDS_Mesh();
BDS_Mesh(const BDS_Mesh &other);
std::set<BDS_GeomEntity*, GeomLessThan> geom;
@@ -429,9 +426,9 @@ class BDS_Mesh
BDS_Edge *find_edge(BDS_Point *p1, BDS_Point *p2, BDS_Face *t) const;
// Triangles & Quadrangles
BDS_Face *add_triangle(int p1, int p2, int p3);
- BDS_Face *add_quadrangle(int p1, int p2, int p3, int p4);
+ // BDS_Face *add_quadrangle(int p1, int p2, int p3, int p4);
BDS_Face *add_triangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3);
- BDS_Face *add_quadrangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3, BDS_Edge *e4);
+ // BDS_Face *add_quadrangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3, BDS_Edge *e4);
void del_face(BDS_Face *t);
BDS_Face *find_triangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3);
BDS_Face *find_quadrangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3, BDS_Edge *e4);
@@ -444,13 +441,9 @@ class BDS_Mesh
BDS_Edge *recover_edge_fast(BDS_Point *p1, BDS_Point *p2);
bool swap_edge(BDS_Edge*, const BDS_SwapEdgeTest &theTest);
bool collapse_edge_parametric(BDS_Edge*, BDS_Point*);
- void snap_point(BDS_Point*, BDS_Mesh *geom = 0);
- bool smooth_point(BDS_Point*, BDS_Mesh *geom = 0);
bool smooth_point_parametric(BDS_Point *p, GFace *gf);
bool smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality=false);
- bool move_point(BDS_Point *p , double X, double Y, double Z);
bool split_edge(BDS_Edge *, BDS_Point *);
- bool split_face(BDS_Face *, BDS_Point *);
bool edge_constraint(BDS_Point *p1, BDS_Point *p2);
// Global operators
void cleanup();
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index ffbe77be989cf9b362d399bd02ad5694f50507ca..a25131ee307a5c4e87540b3a8831a75180fddf36 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -14,7 +14,6 @@ set(SRC
meshGFaceTransfinite.cpp meshGFaceExtruded.cpp
meshGFaceBamg.cpp meshGFaceBDS.cpp meshGFaceDelaunayInsertion.cpp
meshGFaceLloyd.cpp meshGFaceOptimize.cpp
- meshGFaceQuadrilateralize.cpp
meshGRegion.cpp
meshDiscreteRegion.cpp
meshGRegionDelaunayInsertion.cpp meshGRegionTransfinite.cpp
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index f1dc9d3130c3d511991ecf623aa7e35188e21abd..1b6be4db7c89d82b99a2ac5462fc9111ff9af2c2 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -2246,6 +2246,13 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
gf->meshStatistics.nbGoodLength);*/
gf->meshStatistics.status = GFace::DONE;
+
+ if(debug){
+ char name[245];
+ sprintf(name, "surface%d-just-real.pos", gf->tag());
+ outputScalarField(m->triangles, name, 0);
+ }
+
// if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine || 1) {
// backgroundMesh::unset();
// }
@@ -2293,7 +2300,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
}
// recoverMap.insert(new_relations.begin(), new_relations.end());
}
- Msg::Info("%d points that are duplicated for Delaunay meshing", equivalence.size());
+ // Msg::Info("%d points that are duplicated for Delaunay meshing", equivalence.size());
// fill the small gmsh structures
{
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 16ccdeb59cef132730a3fe526d462428f881c485..903bf079158b1693b10f0368a1c2176e366fca43 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -54,65 +54,6 @@ inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f)
return l1 + l2;
}
-inline double computeEdgeLinearLength_new(BDS_Point *p1, BDS_Point *p2, GFace *f)
-{
- const int nbSb = 10;
- GPoint GP[nbSb];
-
- for (int i = 1; i < nbSb; i++){
- double xi = (double)i / nbSb;
- GP[i-1] = f->point(SPoint2(((1-xi) * p1->u + xi * p2->u),
- ((1-xi) * p1->v + xi * p2->v)));
- if (!GP[i-1].succeeded())
- return computeEdgeLinearLength(p1,p2);
- }
- double l = 0;
- for (int i = 0; i < nbSb; i++){
- if (i == 0){
- const double dx1 = p1->X - GP[0].x();
- const double dy1 = p1->Y - GP[0].y();
- const double dz1 = p1->Z - GP[0].z();
- l+= sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
- }
- else if (i == nbSb-1){
- const double dx1 = p2->X - GP[nbSb-1].x();
- const double dy1 = p2->Y - GP[nbSb-1].y();
- const double dz1 = p2->Z - GP[nbSb-1].z();
- l+= sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
- }
- else{
- const double dx1 = GP[i].x() - GP[i-1].x();
- const double dy1 = GP[i].y() - GP[i-1].y();
- const double dz1 = GP[i].z() - GP[i-1].z();
- l+= sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
- }
- }
- return l;
-}
-
-inline double computeEdgeMiddleCoord(BDS_Point *p1, BDS_Point *p2, GFace *f)
-{
- if (f->geomType() == GEntity::Plane)
- return 0.5;
-
- GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u),
- 0.5 * (p1->v + p2->v)));
-
- const double dx1 = p1->X - GP.x();
- const double dy1 = p1->Y - GP.y();
- const double dz1 = p1->Z - GP.z();
- const double l1 = sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
- const double dx2 = p2->X - GP.x();
- const double dy2 = p2->Y - GP.y();
- const double dz2 = p2->Z - GP.z();
- const double l2 = sqrt(dx2 * dx2 + dy2 * dy2 + dz2 * dz2);
-
- if (l1 > l2)
- return 0.25 * (l1 + l2) / l1;
- else
- return 0.25 * (3 * l2 - l1) / l2;
-}
-
inline double computeEdgeLinearLength(BDS_Edge *e, GFace *f)
{
// FIXME !!!
@@ -444,9 +385,9 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
std::sort(edges.begin(), edges.end(), edges_sort);
bool isPeriodic = gf->periodic(0) || gf->periodic(1) ;
- SPoint3 center;
- double radius;
- bool isSphere = gf->isSphere (radius, center);
+ // SPoint3 center;
+ // double radius;
+ // bool isSphere = gf->isSphere (radius, center);
for (unsigned int i = 0; i < edges.size(); ++i){
BDS_Edge *e = edges[i].second;
@@ -500,11 +441,43 @@ double getMaxLcWhenCollapsingEdge(GFace *gf, BDS_Mesh &m, BDS_Edge *e, BDS_Point
return maxLc;
}
+static bool revertTriangleSphere (SPoint3 ¢er, BDS_Point *p, BDS_Point *o){
+ std::list<BDS_Face*> t;
+ p->getTriangles(t);
+ std::list<BDS_Face*>::iterator it = t.begin();
+ std::list<BDS_Face*>::iterator ite = t.end();
+ while(it != ite) {
+ BDS_Face *t = *it;
+ BDS_Point *pts[4];
+ t->getNodes(pts);
+ pts[0] = (pts[0] == p) ? o : pts[0];
+ pts[1] = (pts[1] == p) ? o : pts[1];
+ pts[2] = (pts[2] == p) ? o : pts[2];
+ double norm[3];
+ normal_triangle(pts[0], pts[1], pts[2], norm);
+ double dx = center.x() - o->X;
+ double dy = center.y() - o->Y;
+ double dz = center.z() - o->Z;
+ double ps = dx*norm[0]+dy*norm[1]+dz*norm[2];
+ if (ps < 0){
+ // printf("FLIIIP\n");
+ return true;
+ }
+ ++it;
+ }
+ return false;
+}
+
+
void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_collaps)
{
std::list<BDS_Edge*>::iterator it = m.edges.begin();
std::vector<std::pair<double, BDS_Edge*> > edges;
+ double radius;
+ SPoint3 center;
+ bool isSphere = gf->isSphere (radius, center);
+
while (it != m.edges.end()){
if(!(*it)->deleted && (*it)->numfaces() == 2 && (*it)->g->classif_degree == 2){
@@ -547,8 +520,11 @@ void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_c
else if (collapseP1Allowed && !collapseP2Allowed)
p = e->p2;
+ bool flip = false;
+ if (p && isSphere)flip = revertTriangleSphere(center, p, e->othervertex(p));
+
bool res = false;
- if(p)
+ if(!flip && p)
res = m.collapse_edge_parametric(e, p);
if(res)
nb_collaps++;
@@ -556,31 +532,6 @@ void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_c
}
}
-void collapseEdgePassUnSorted(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP,
- int &nb_collaps)
-{
- int NN1 = m.edges.size();
- int NN2 = 0;
- std::list<BDS_Edge*>::iterator it = m.edges.begin();
-
- while (1){
- if(NN2++ >= NN1) break;
-
- if(!(*it)->deleted){
- double lone = NewGetLc(*it, gf);
- if(!(*it)->deleted && (*it)->numfaces() == 2 && lone < MINE_){
- bool res = false;
- if((*it)->p1->iD > MAXNP)
- res = m.collapse_edge_parametric(*it, (*it)->p1);
- else if((*it)->p2->iD > MAXNP)
- res = m.collapse_edge_parametric(*it, (*it)->p2);
- if(res)
- nb_collaps++;
- }
- }
- ++it;
- }
-}
void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q)
{
@@ -735,10 +686,6 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
Msg::Debug(" ---------------------------------------");
}
-void allowAppearanceofEdge (BDS_Point *p1, BDS_Point *p2)
-{
-}
-
void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*, MVertex*,PointLessThan> *recoverMap,
std::set<BDS_Edge*, EdgeLessThan> &toSplit)
{
@@ -831,9 +778,69 @@ int solveInvalidPeriodic(GFace *gf, BDS_Mesh &m,
return toSplit.size();
}
+void TRYTOFIXSPHERES(GFace *gf, BDS_Mesh &m,
+ std::map<BDS_Point*,MVertex*,PointLessThan> *recoverMap=0)
+{
+ if (!recoverMap)return;
+ double radius;
+ SPoint3 center;
+ bool isSphere = gf->isSphere(radius, center);
+ if (!isSphere)return;
+
+
+ while(1){
+ int count = 0;
+ std::list<BDS_Edge*>::iterator ite = m.edges.begin();
+ while (ite != m.edges.end()){
+ BDS_Edge *e = *ite;
+ if(e->numfaces() == 2){
+ double ps[2] = {1,1};
+ for (int i=0;i<2;i++){
+ BDS_Face *f = e->faces(i);
+ double norm[3];
+ BDS_Point *n[4];
+ f->getNodes(n);
+
+ MVertex *v1 = (recoverMap->find(n[0])==recoverMap->end()) ? NULL : (*recoverMap)[n[0]];
+ MVertex *v2 = (recoverMap->find(n[1])==recoverMap->end()) ? NULL : (*recoverMap)[n[1]];
+ MVertex *v3 = (recoverMap->find(n[2])==recoverMap->end()) ? NULL : (*recoverMap)[n[2]];
+
+ if ((!v1 || (v1 != v2 && v1 != v3)) &&
+ (!v2 || v2 != v3)){
+
+ normal_triangle(n[0], n[1], n[2], norm);
+ double x = (n[0]->X+n[1]->X+n[2]->X)/3.0;
+ double y = (n[0]->Y+n[1]->Y+n[2]->Y)/3.0;
+ double z = (n[0]->Z+n[1]->Z+n[2]->Z)/3.0;
+ double dx = center.x() - x;
+ double dy = center.y() - y;
+ double dz = center.z() - z;
+ ps[i] = dx*norm[0]+dy*norm[1]+dz*norm[2];
+ }
+ }
+ if (ps[0]*ps[1] < 0){
+ printf("Collapsing edge %d %d Because one oof the two triangles is reverted\n",
+ e->p1->iD,e->p2->iD);
+ count++;
+ if (recoverMap->find(e->p1) == recoverMap->end()){
+ m.collapse_edge_parametric(e, e->p1);
+ }
+ else if (recoverMap->find(e->p2) == recoverMap->end()){
+ m.collapse_edge_parametric(e, e->p2);
+ }
+ }
+ }
+ ++ite;
+ }
+ if (!count)break;
+ }
+}
+
+
void optimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
std::map<BDS_Point*,MVertex*,PointLessThan> *recoverMap=0)
{
+ // return;
int nb_swap;
delaunayizeBDS(gf, m, nb_swap);
@@ -859,17 +866,11 @@ void optimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
while(solveInvalidPeriodic(gf, m, recoverMap)){
}
}
+ TRYTOFIXSPHERES(gf,m,recoverMap);
}
// DELAUNAY BDS
-void delaunayPointInsertionBDS(GFace *gf, BDS_Mesh &m, BDS_Point *v, BDS_Face *f)
-{
- m.split_face(f, v);
- int nb_swap = 0;
- delaunayizeBDS(gf, m, nb_swap);
-}
-
// build the BDS from a list of GFace
// This is a TRUE copy
BDS_Mesh *gmsh2BDS(std::list<GFace*> &l)
@@ -904,20 +905,3 @@ BDS_Mesh *gmsh2BDS(std::list<GFace*> &l)
}
return m;
}
-
-void collapseSmallEdges(GModel &gm)
-{
- return;
- // gm.renumberMeshVertices(true);
- std::list<GFace*> faces;
- for (GModel::fiter fit = gm.firstFace(); fit != gm.lastFace(); fit++){
- faces.push_back(*fit);
- }
- BDS_Mesh *pm = gmsh2BDS(faces);
- outputScalarField(pm->triangles, "all.pos", 0);
-
- for (GModel::eiter eit = gm.firstEdge(); eit != gm.lastEdge(); eit++){
- }
-
- delete pm;
-}