diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp index 75bef5ef98e60b3360291aa16d0bf74a45629460..3967e1c14772f5def8508cec46cbe355afddb9dc 100644 --- a/Mesh/BDS.cpp +++ b/Mesh/BDS.cpp @@ -1,4 +1,4 @@ -// $Id: BDS.cpp,v 1.79 2007-10-10 13:59:30 remacle Exp $ +// $Id: BDS.cpp,v 1.80 2007-10-11 08:59:22 remacle Exp $ // // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle // @@ -206,7 +206,7 @@ int Intersect_Edges_2d(double x1, double y1, double x2, double y2, return 0; } -BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2,std::set<EdgeToRecover> *e2r) +BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2r, std::set<EdgeToRecover> *not_recovered ) { BDS_Edge *e = find_edge (num1, num2); @@ -220,6 +220,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2,std::set<EdgeToRecover> *e2r Msg(INFO," edge %d %d has to be recovered",num1,num2); int ix = 0; + int ixMax = 100; while(1) { std::vector<BDS_Edge *> intersected; @@ -238,17 +239,21 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2,std::set<EdgeToRecover> *e2r { std::set<EdgeToRecover>::iterator itr1 = e2r->find(EdgeToRecover(e->p1->iD,e->p2->iD,0)); std::set<EdgeToRecover>::iterator itr2 = e2r->find(EdgeToRecover(num1,num2,0)); - Msg(GERROR," edge %d %d on model edge %d cannot be recovered because it intersects %d %d on model edge %d", + Msg(WARNING," edge %d %d on model edge %d cannot be recovered because it intersects %d %d on model edge %d", num1,num2,itr2->ge->tag(), e->p1->iD,e->p2->iD,itr1->ge->tag()); - return false; + + // now throw a class that contains the diagnostic + not_recovered->insert (EdgeToRecover( num1 , num2, itr2->ge)); + not_recovered->insert (EdgeToRecover( e->p1->iD, e->p2->iD, itr1->ge)); + ixMax = -1; } intersected.push_back(e); } ++it; } - if (!intersected.size() || ix > 100) + if (!intersected.size() || ix > ixMax) { BDS_Edge *eee = find_edge (num1, num2); if (!eee) diff --git a/Mesh/BDS.h b/Mesh/BDS.h index 7a7c584ddc6d405ce7892a6ed911d6993ca2083f..eabd18d6f87937e6142a74ba62592a2b7624fd12 100644 --- a/Mesh/BDS.h +++ b/Mesh/BDS.h @@ -47,6 +47,7 @@ double surface_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3); double surface_triangle_param(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3); double quality_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3); + class BDS_GeomEntity { public: @@ -447,7 +448,7 @@ public: void add_geom(int degree, int tag); BDS_GeomEntity *get_geom(int p1, int p2); // 2D operators - BDS_Edge *recover_edge(int p1, int p2, std::set<EdgeToRecover> *e2r=0); + BDS_Edge *recover_edge(int p1, int p2, std::set<EdgeToRecover> *e2r=0, std::set<EdgeToRecover> *not_recovered = 0); 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); diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp index cccfc7e1f40cefffebacb0701c628b24997de45a..fbd8be9dfb0d40db0d10fc204aa4e7bb4c2a00ec 100644 --- a/Mesh/meshGEdge.cpp +++ b/Mesh/meshGEdge.cpp @@ -1,4 +1,4 @@ -// $Id: meshGEdge.cpp,v 1.44 2007-10-10 08:49:34 remacle Exp $ +// $Id: meshGEdge.cpp,v 1.45 2007-10-11 08:59:22 remacle Exp $ // // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle // @@ -365,12 +365,15 @@ void meshGEdge::operator() (GEdge *ge) const double d = (double)NUMP * b; if((fabs(P2.p) >= fabs(d)) && (fabs(P1.p) < fabs(d))) { double dt = P2.t - P1.t; - double dlc = P2.lc - P1.lc; + double dlc = P2.lc - P1.lc; double dp = P2.p - P1.p; double t = P1.t + dt / dp * (d - P1.p); - double lc = P1.lc + dlc / dp * (d - P1.p); + SVector3 der = ge->firstDer(t); + const double d = norm(der); + double lc = d/(P1.lc + dlc / dp * (d - P1.p)); GPoint V = ge->point(t); ge->mesh_vertices[NUMP - 1] = new MEdgeVertex(V.x(), V.y(), V.z(), ge, t, lc); + // printf("lc = %12.5E %12.5E \n",lc,P1.lc,P2.lc); NUMP++; } else { diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp index 1e584cf0ba1680a7d29610829e11d0564eeca14b..cb57703bb5eae85500d5b85f02d1fb86d31939cf 100644 --- a/Mesh/meshGFace.cpp +++ b/Mesh/meshGFace.cpp @@ -1,4 +1,4 @@ -// $Id: meshGFace.cpp,v 1.91 2007-10-10 13:59:30 remacle Exp $ +// $Id: meshGFace.cpp,v 1.92 2007-10-11 08:59:22 remacle Exp $ // // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle // @@ -485,7 +485,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT) -bool recover_medge ( BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r, int pass_) +bool recover_medge ( BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r, std::set<EdgeToRecover> *not_recovered, int pass_) { BDS_GeomEntity *g=0; if (pass_ == 2) @@ -521,10 +521,10 @@ bool recover_medge ( BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r, int p BDS_GeomEntity *g0 = m->get_geom(vend->getNum(), 0); pend->g = g0; } - BDS_Edge * e = m->recover_edge ( vstart->getNum(), vend->getNum(),e2r); + BDS_Edge * e = m->recover_edge ( vstart->getNum(), vend->getNum(),e2r, not_recovered); if (e)e->g = g; else { - Msg(GERROR,"The unrecoverable edge is on model edge %d",ge->tag()); + // Msg(GERROR,"The unrecoverable edge is on model edge %d",ge->tag()); return false; } return true; @@ -545,11 +545,11 @@ bool recover_medge ( BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r, int p BDS_GeomEntity *g0 = m->get_geom(vstart->getNum(), 0); pstart->g = g0; } - e = m->recover_edge ( vstart->getNum(), vend->getNum(),e2r); + e = m->recover_edge ( vstart->getNum(), vend->getNum(),e2r, not_recovered); if (e)e->g = g; else { - Msg(GERROR,"The unrecoverable edge is on model edge %d",ge->tag()); - return false; + // Msg(GERROR,"The unrecoverable edge is on model edge %d",ge->tag()); + // return false; } } @@ -560,11 +560,11 @@ bool recover_medge ( BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r, int p if (pass_ == 1)e2r->insert (EdgeToRecover (vstart->getNum(), vend->getNum(),ge)); else { - e = m->recover_edge ( vstart->getNum(), vend->getNum(),e2r); + e = m->recover_edge ( vstart->getNum(), vend->getNum(),e2r, not_recovered); if (e)e->g = g; else { - Msg(GERROR,"Unable to recover an edge %g %g && %g %g (%d/%d)",vstart->x(),vstart->y(), vend->x(),vend->y(),i,ge->mesh_vertices.size()); - return false; + // Msg(GERROR,"Unable to recover an edge %g %g && %g %g (%d/%d)",vstart->x(),vstart->y(), vend->x(),vend->y(),i,ge->mesh_vertices.size()); + // return false; } } } @@ -573,11 +573,11 @@ bool recover_medge ( BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r, int p if (pass_ == 1)e2r->insert (EdgeToRecover (vstart->getNum(), vend->getNum(),ge)); else { - e = m->recover_edge ( vstart->getNum(), vend->getNum(),e2r); + e = m->recover_edge ( vstart->getNum(), vend->getNum(),e2r, not_recovered); if (e)e->g = g; else { - Msg(GERROR,"Unable to recover an edge %g %g && %g %g (%d/%d)",vstart->x(),vstart->y(), vend->x(),vend->y(),ge->mesh_vertices.size(),ge->mesh_vertices.size()); - return false; + // Msg(GERROR,"Unable to recover an edge %g %g && %g %g (%d/%d)",vstart->x(),vstart->y(), vend->x(),vend->y(),ge->mesh_vertices.size(),ge->mesh_vertices.size()); + // return false; } BDS_Point *pend = m->find_point(vend->getNum()); if(!pend->g) @@ -769,18 +769,20 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true) BDS_Point *pp = m->add_point ( num, U,V, gf); GEntity *ge = here->onWhat(); - if(ge->dim() == 1) + if(ge->dim() == 0) { - double t; - here->getParameter(0,t); - pp->lcBGM() = BGM_MeshSize(ge,t,-12,here->x(),here->y(),here->z()); + pp->lcBGM() = BGM_MeshSize(ge,0,0,here->x(),here->y(),here->z()); } - else + else if(ge->dim() == 1) { MEdgeVertex *eve = (MEdgeVertex*) here; pp->lcBGM() = eve->getLc(); } + else + pp->lcBGM() = 1.e22; + pp->lc() = pp->lcBGM(); + // printf("dim %d lc = %12.5E\n",ge->dim(),pp->lc()); } Msg(DEBUG1,"Meshing of the convex hull (%d points) done",all_vertices.size()); @@ -821,17 +823,18 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true) // 1D mesh have to be densified std::set<EdgeToRecover> edgesToRecover; + std::set<EdgeToRecover> edgesNotRecovered; it = edges.begin(); while(it != edges.end()) { if(!(*it)->is_mesh_degenerated()) - recover_medge ( m, *it, &edgesToRecover,1); + recover_medge ( m, *it, &edgesToRecover, &edgesNotRecovered, 1); ++it; } it = emb_edges.begin(); while(it != emb_edges.end()) { - recover_medge ( m, *it, &edgesToRecover,1); + recover_medge ( m, *it, &edgesToRecover, &edgesNotRecovered, 1); ++it; } @@ -843,14 +846,17 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true) while(it != edges.end()) { if(!(*it)->is_mesh_degenerated()){ - if (!recover_medge ( m, *it, &edgesToRecover,2)) - { - Msg(GERROR,"Face cannot be not meshed because 1D mesh self intersects !"); - return false; - } + recover_medge ( m, *it, &edgesToRecover, &edgesNotRecovered, 2); } ++it; } + + if (edgesNotRecovered.size()) + { + Msg(GERROR,"%d edges were not recovered on model face %d, gmsh goes back and refines the 1d mesh",edgesNotRecovered.size(),gf->tag()); + return false; + } + // Msg(INFO,"Boundary Edges recovered for surface %d",gf->tag()); // Look for an edge that is on the boundary for which one of the // two neighbors has a negative number node. The other triangle @@ -886,7 +892,7 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true) it = emb_edges.begin(); while(it != emb_edges.end()) { - recover_medge ( m, *it, &edgesToRecover,2); + recover_medge ( m, *it, &edgesToRecover, &edgesNotRecovered, 2); ++it; }