Skip to content
Snippets Groups Projects
Commit 2a69e34f authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

spheres are ok

parent e57a860a
Branches
Tags
No related merge requests found
...@@ -483,14 +483,6 @@ BDS_Face *BDS_Mesh::add_triangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3) ...@@ -483,14 +483,6 @@ BDS_Face *BDS_Mesh::add_triangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3)
return t; 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) void BDS_Mesh::del_face(BDS_Face *t)
{ {
t->e1->del(t); t->e1->del(t);
...@@ -624,46 +616,6 @@ BDS_Mesh::~BDS_Mesh() ...@@ -624,46 +616,6 @@ BDS_Mesh::~BDS_Mesh()
DESTROOOY(triangles.begin(), triangles.end()); 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) bool BDS_Mesh::split_edge(BDS_Edge *e, BDS_Point *mid)
{ {
......
...@@ -403,12 +403,9 @@ struct EdgeToRecover ...@@ -403,12 +403,9 @@ struct EdgeToRecover
class BDS_Mesh class BDS_Mesh
{ {
public: public:
int MAXPOINTNUMBER, SNAP_SUCCESS, SNAP_FAILURE; int MAXPOINTNUMBER;
double Min[3], Max[3], LC; double Min[3], Max[3], LC;
BDS_Mesh(int _MAXX = 0) : MAXPOINTNUMBER(_MAXX){} 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(); virtual ~BDS_Mesh();
BDS_Mesh(const BDS_Mesh &other); BDS_Mesh(const BDS_Mesh &other);
std::set<BDS_GeomEntity*, GeomLessThan> geom; std::set<BDS_GeomEntity*, GeomLessThan> geom;
...@@ -429,9 +426,9 @@ class BDS_Mesh ...@@ -429,9 +426,9 @@ class BDS_Mesh
BDS_Edge *find_edge(BDS_Point *p1, BDS_Point *p2, BDS_Face *t) const; BDS_Edge *find_edge(BDS_Point *p1, BDS_Point *p2, BDS_Face *t) const;
// Triangles & Quadrangles // Triangles & Quadrangles
BDS_Face *add_triangle(int p1, int p2, int p3); 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_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); void del_face(BDS_Face *t);
BDS_Face *find_triangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3); 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); BDS_Face *find_quadrangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3, BDS_Edge *e4);
...@@ -444,13 +441,9 @@ class BDS_Mesh ...@@ -444,13 +441,9 @@ class BDS_Mesh
BDS_Edge *recover_edge_fast(BDS_Point *p1, BDS_Point *p2); BDS_Edge *recover_edge_fast(BDS_Point *p1, BDS_Point *p2);
bool swap_edge(BDS_Edge*, const BDS_SwapEdgeTest &theTest); bool swap_edge(BDS_Edge*, const BDS_SwapEdgeTest &theTest);
bool collapse_edge_parametric(BDS_Edge*, BDS_Point*); 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_parametric(BDS_Point *p, GFace *gf);
bool smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality=false); 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_edge(BDS_Edge *, BDS_Point *);
bool split_face(BDS_Face *, BDS_Point *);
bool edge_constraint(BDS_Point *p1, BDS_Point *p2); bool edge_constraint(BDS_Point *p1, BDS_Point *p2);
// Global operators // Global operators
void cleanup(); void cleanup();
......
...@@ -14,7 +14,6 @@ set(SRC ...@@ -14,7 +14,6 @@ set(SRC
meshGFaceTransfinite.cpp meshGFaceExtruded.cpp meshGFaceTransfinite.cpp meshGFaceExtruded.cpp
meshGFaceBamg.cpp meshGFaceBDS.cpp meshGFaceDelaunayInsertion.cpp meshGFaceBamg.cpp meshGFaceBDS.cpp meshGFaceDelaunayInsertion.cpp
meshGFaceLloyd.cpp meshGFaceOptimize.cpp meshGFaceLloyd.cpp meshGFaceOptimize.cpp
meshGFaceQuadrilateralize.cpp
meshGRegion.cpp meshGRegion.cpp
meshDiscreteRegion.cpp meshDiscreteRegion.cpp
meshGRegionDelaunayInsertion.cpp meshGRegionTransfinite.cpp meshGRegionDelaunayInsertion.cpp meshGRegionTransfinite.cpp
......
...@@ -2246,6 +2246,13 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) ...@@ -2246,6 +2246,13 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
gf->meshStatistics.nbGoodLength);*/ gf->meshStatistics.nbGoodLength);*/
gf->meshStatistics.status = GFace::DONE; 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) { // if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine || 1) {
// backgroundMesh::unset(); // backgroundMesh::unset();
// } // }
...@@ -2293,7 +2300,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) ...@@ -2293,7 +2300,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
} }
// recoverMap.insert(new_relations.begin(), new_relations.end()); // 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 // fill the small gmsh structures
{ {
......
...@@ -54,65 +54,6 @@ inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f) ...@@ -54,65 +54,6 @@ inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f)
return l1 + l2; 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) inline double computeEdgeLinearLength(BDS_Edge *e, GFace *f)
{ {
// FIXME !!! // FIXME !!!
...@@ -444,9 +385,9 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split) ...@@ -444,9 +385,9 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
std::sort(edges.begin(), edges.end(), edges_sort); std::sort(edges.begin(), edges.end(), edges_sort);
bool isPeriodic = gf->periodic(0) || gf->periodic(1) ; bool isPeriodic = gf->periodic(0) || gf->periodic(1) ;
SPoint3 center; // SPoint3 center;
double radius; // double radius;
bool isSphere = gf->isSphere (radius, center); // bool isSphere = gf->isSphere (radius, center);
for (unsigned int i = 0; i < edges.size(); ++i){ for (unsigned int i = 0; i < edges.size(); ++i){
BDS_Edge *e = edges[i].second; BDS_Edge *e = edges[i].second;
...@@ -500,11 +441,43 @@ double getMaxLcWhenCollapsingEdge(GFace *gf, BDS_Mesh &m, BDS_Edge *e, BDS_Point ...@@ -500,11 +441,43 @@ double getMaxLcWhenCollapsingEdge(GFace *gf, BDS_Mesh &m, BDS_Edge *e, BDS_Point
return maxLc; return maxLc;
} }
static bool revertTriangleSphere (SPoint3 &center, 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) void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_collaps)
{ {
std::list<BDS_Edge*>::iterator it = m.edges.begin(); std::list<BDS_Edge*>::iterator it = m.edges.begin();
std::vector<std::pair<double, BDS_Edge*> > edges; std::vector<std::pair<double, BDS_Edge*> > edges;
double radius;
SPoint3 center;
bool isSphere = gf->isSphere (radius, center);
while (it != m.edges.end()){ while (it != m.edges.end()){
if(!(*it)->deleted && (*it)->numfaces() == 2 && (*it)->g->classif_degree == 2){ 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 ...@@ -547,8 +520,11 @@ void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_c
else if (collapseP1Allowed && !collapseP2Allowed) else if (collapseP1Allowed && !collapseP2Allowed)
p = e->p2; p = e->p2;
bool flip = false;
if (p && isSphere)flip = revertTriangleSphere(center, p, e->othervertex(p));
bool res = false; bool res = false;
if(p) if(!flip && p)
res = m.collapse_edge_parametric(e, p); res = m.collapse_edge_parametric(e, p);
if(res) if(res)
nb_collaps++; nb_collaps++;
...@@ -556,31 +532,6 @@ void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_c ...@@ -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) 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, ...@@ -735,10 +686,6 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
Msg::Debug(" ---------------------------------------"); Msg::Debug(" ---------------------------------------");
} }
void allowAppearanceofEdge (BDS_Point *p1, BDS_Point *p2)
{
}
void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*, MVertex*,PointLessThan> *recoverMap, void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*, MVertex*,PointLessThan> *recoverMap,
std::set<BDS_Edge*, EdgeLessThan> &toSplit) std::set<BDS_Edge*, EdgeLessThan> &toSplit)
{ {
...@@ -831,9 +778,69 @@ int solveInvalidPeriodic(GFace *gf, BDS_Mesh &m, ...@@ -831,9 +778,69 @@ int solveInvalidPeriodic(GFace *gf, BDS_Mesh &m,
return toSplit.size(); 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, void optimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
std::map<BDS_Point*,MVertex*,PointLessThan> *recoverMap=0) std::map<BDS_Point*,MVertex*,PointLessThan> *recoverMap=0)
{ {
// return;
int nb_swap; int nb_swap;
delaunayizeBDS(gf, m, nb_swap); delaunayizeBDS(gf, m, nb_swap);
...@@ -859,17 +866,11 @@ void optimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT, ...@@ -859,17 +866,11 @@ void optimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
while(solveInvalidPeriodic(gf, m, recoverMap)){ while(solveInvalidPeriodic(gf, m, recoverMap)){
} }
} }
TRYTOFIXSPHERES(gf,m,recoverMap);
} }
// DELAUNAY BDS // 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 // build the BDS from a list of GFace
// This is a TRUE copy // This is a TRUE copy
BDS_Mesh *gmsh2BDS(std::list<GFace*> &l) BDS_Mesh *gmsh2BDS(std::list<GFace*> &l)
...@@ -904,20 +905,3 @@ BDS_Mesh *gmsh2BDS(std::list<GFace*> &l) ...@@ -904,20 +905,3 @@ BDS_Mesh *gmsh2BDS(std::list<GFace*> &l)
} }
return m; 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;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment