diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp index 73c43438641cc3c61aba97449743c429ad2f3c38..ba6f8a2e153d34ea9ac036a6536806505a09c3d0 100644 --- a/Geo/gmshEdge.cpp +++ b/Geo/gmshEdge.cpp @@ -246,7 +246,7 @@ SPoint2 gmshEdge::reparamOnFace(const GFace *face, double epar,int dir) const } return SPoint2(U, V); } - else if (s->Typ == MSH_SURF_TRIC){ + else if (s->Typ == MSH_SURF_TRIC){ Curve *C[3]; for(int i = 0; i < 3; i++) List_Read(s->Generatrices, i, &C[i]); @@ -288,8 +288,9 @@ SPoint2 gmshEdge::reparamOnFace(const GFace *face, double epar,int dir) const } return SPoint2(U, V); } - else + else{ return GEdge::reparamOnFace(face, epar, dir); + } } void gmshEdge::writeGEO(FILE *fp) diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp index 4057afbdbde4226cb68be8218d77922a61fc0528..6ddecf55929f39577209bb42ea9ff3c32deb1c58 100644 --- a/Geo/gmshFace.cpp +++ b/Geo/gmshFace.cpp @@ -45,7 +45,7 @@ gmshFace::gmshFace(GModel *m, Surface *face) e->addFace(this); l_dirs.push_back((c->Num > 0) ? 1 : -1); if (List_Nbr(s->Generatrices) == 2){ - e->meshAttributes.minimumMeshSegments = + e->meshAttributes.minimumMeshSegments = std::max(e->meshAttributes.minimumMeshSegments, 2); } } @@ -116,11 +116,11 @@ void gmshFace::resetMeshAttributes() else Msg::Error("Unknown vertex %d in transfinite attributes", corn->Num); } - } + } } Range<double> gmshFace::parBounds(int i) const -{ +{ return Range<double>(0, 1); } @@ -183,7 +183,7 @@ Pair<SVector3, SVector3> gmshFace::firstDer(const SPoint2 ¶m) const } } -void gmshFace::secondDer(const SPoint2 ¶m, +void gmshFace::secondDer(const SPoint2 ¶m, SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const { if(s->Typ == MSH_SURF_PLAN && !s->geometry){ @@ -247,7 +247,7 @@ GPoint gmshFace::closestPoint(const SPoint3 & qp, const double initialGuess[2]) } } sys2x2(MN, BN, UV); - return GPoint(XP, YP, ZP, this, UV); + return GPoint(XP, YP, ZP, this, UV); } Vertex v; @@ -269,7 +269,7 @@ SPoint2 gmshFace::parFromPoint(const SPoint3 &qp, bool onSurface) const double u, v, vec[3] = {qp.x() - x, qp.y() - y, qp.z() - z}; prosca(vec, VX, &u); prosca(vec, VY, &v); - return SPoint2(u, v); + return SPoint2(u, v); } else{ return GFace::parFromPoint(qp, onSurface); @@ -279,25 +279,25 @@ SPoint2 gmshFace::parFromPoint(const SPoint3 &qp, bool onSurface) const GEntity::GeomType gmshFace::geomType() const { switch(s->Typ){ - case MSH_SURF_PLAN: - if(s->geometry) + case MSH_SURF_PLAN: + if(s->geometry) return ParametricSurface; else return Plane; case MSH_SURF_REGL: - case MSH_SURF_TRIC: + case MSH_SURF_TRIC: return RuledSurface; - case MSH_SURF_DISCRETE: + case MSH_SURF_DISCRETE: return DiscreteSurface; case MSH_SURF_BND_LAYER: return BoundaryLayerSurface; - default: + default: return Unknown; } } bool gmshFace::containsPoint(const SPoint3 &pt) const -{ +{ if(s->Typ == MSH_SURF_PLAN){ // OK to use the normal from the mean plane here: we compensate // for the (possibly wrong) orientation at the end @@ -320,7 +320,7 @@ bool gmshFace::containsPoint(const SPoint3 &pt) const } } // we're inside if angle equals 2 * pi - if(fabs(angle) > 2 * M_PI - 0.5 && fabs(angle) < 2 * M_PI + 0.5) + if(fabs(angle) > 2 * M_PI - 0.5 && fabs(angle) < 2 * M_PI + 0.5) return true; return false; } @@ -338,7 +338,7 @@ bool gmshFace::buildSTLTriangulation(bool force) else return true; } - + stl_vertices.clear(); stl_triangles.clear(); @@ -355,21 +355,21 @@ bool gmshFace::buildSTLTriangulation(bool force) meshGFace mesher (false,true); mesher(this); printf("%d triangles face %d\n",triangles_stl.size(),tag()); - CTX::instance()->mesh = _temp; + CTX::instance()->mesh = _temp; std::map<MVertex*,int> _v; int COUNT =0; for (int j = 0; j < triangles_stl.size(); j++){ int C[3]; for (int i=0;i<3;i++){ - std::map<MVertex*,int>::iterator it = + std::map<MVertex*,int>::iterator it = _v.find(triangles_stl[j]->getVertex(j)); if (it != _v.end()){ stl_triangles.push_back(COUNT); _v[triangles_stl[j]->getVertex(j)] = COUNT++; } else stl_triangles.push_back(it->second); - } + } } std::map<MVertex*,int>::iterator itv = _v.begin(); for ( ; itv != _v.end() ; ++itv){ @@ -378,7 +378,7 @@ bool gmshFace::buildSTLTriangulation(bool force) reparamMeshVertexOnFace(v, this, param); stl_vertices.push_back(param); } - + va_geom_triangles = new VertexArray(3, stl_triangles.size() / 3); unsigned int c = CTX::instance()->color.geom.surface; unsigned int col[4] = {c, c, c, c}; diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp index cd67fe73706cade99a9d91f0db2fe85e8a9a9a6f..6ccc042d49dd386f03512375778c45572a5dbe06 100644 --- a/Mesh/meshGFace.cpp +++ b/Mesh/meshGFace.cpp @@ -1107,13 +1107,13 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, // } } - computeMeshSizeFieldAccuracy(gf, *m, gf->meshStatistics.efficiency_index, - gf->meshStatistics.longest_edge_length, - gf->meshStatistics.smallest_edge_length, - gf->meshStatistics.nbEdge, - gf->meshStatistics.nbGoodLength); + computeMeshSizeFieldAccuracy(gf, *m, gf->meshStatistics.efficiency_index, + gf->meshStatistics.longest_edge_length, + gf->meshStatistics.smallest_edge_length, + gf->meshStatistics.nbEdge, + gf->meshStatistics.nbGoodLength); //printf("=== Efficiency index is tau=%g\n", gf->meshStatistics.efficiency_index); - + gf->meshStatistics.status = GFace::DONE; // fill the small gmsh structures @@ -1163,7 +1163,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, gf->meshStatistics.nbEdge, gf->meshStatistics.nbGoodLength); //printf("----- Efficiency index is tau=%g\n", gf->meshStatistics.efficiency_index); - + // delete the mesh delete m; @@ -1446,7 +1446,7 @@ static bool meshGeneratorElliptic(GFace *gf, bool debug = true) bool recombine = (CTX::instance()->mesh.recombineAll); int nbBoundaries = gf->edges().size(); - + if (center && recombine && nbBoundaries == 2) { printf("--> regular periodic grid generator (elliptic smooth) \n"); //bool success = createRegularTwoCircleGrid(center, gf); @@ -2091,7 +2091,7 @@ void partitionAndRemesh(GFaceCompound *gf) gf->meshStatistics.efficiency_index += gfc->meshStatistics.efficiency_index; gf->meshStatistics.longest_edge_length = std::max(gf->meshStatistics.longest_edge_length, gfc->meshStatistics.longest_edge_length); - gf->meshStatistics.smallest_edge_length= std::min(gf->meshStatistics.smallest_edge_length, + gf->meshStatistics.smallest_edge_length= std::min(gf->meshStatistics.smallest_edge_length, gfc->meshStatistics.smallest_edge_length); gf->meshStatistics.nbGoodLength += gfc->meshStatistics.nbGoodLength; gf->meshStatistics.nbGoodQuality += gfc->meshStatistics.nbGoodQuality; diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp index ed0124c00879c2717bb8223f58c833613c60f6a2..c75fe74f47942ab63beeba9bd30a9d824b85fe7c 100644 --- a/Mesh/meshGRegion.cpp +++ b/Mesh/meshGRegion.cpp @@ -34,60 +34,73 @@ // hybrid mesh recovery structure class splitQuadRecovery { std::multimap<GEntity*, std::pair<MVertex*,MFace> >_data; - public : - void add (const MFace &f, MVertex *v, GEntity*ge) ; - int buildPyramids (GModel *gm); -}; - -void splitQuadRecovery::add (const MFace &f, MVertex *v, GEntity*ge) { - _data.insert(std::make_pair(ge, std::make_pair(v,f))); -} - -// re-create quads and -int splitQuadRecovery::buildPyramids (GModel *gm) { - int NBPY = 0; - std::set<MFace,Less_Face> allFaces; - for (GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it){ - for (unsigned int i=0; i < (*it)->triangles.size();i++){ - allFaces.insert ((*it)->triangles[i]->getFace(0)); - delete (*it)->triangles[i]; - } - (*it)->triangles.clear(); - for ( std::multimap<GEntity*, std::pair<MVertex*,MFace> >::iterator it2 = _data.lower_bound(*it) ; it2 != _data.upper_bound(*it) ; ++it2) { - MVertex *v = it2->second.first; - v->onWhat()->mesh_vertices.erase - (std::find(v->onWhat()->mesh_vertices.begin(), - v->onWhat()->mesh_vertices.end(), - v)); - const MFace &f = it2->second.second; - std::set<MFace,Less_Face> :: iterator itf0 = allFaces.find(MFace(f.getVertex(0),f.getVertex(1),v)); - std::set<MFace,Less_Face> :: iterator itf1 = allFaces.find(MFace(f.getVertex(1),f.getVertex(2),v)); - std::set<MFace,Less_Face> :: iterator itf2 = allFaces.find(MFace(f.getVertex(2),f.getVertex(3),v)); - std::set<MFace,Less_Face> :: iterator itf3 = allFaces.find(MFace(f.getVertex(3),f.getVertex(0),v)); - if (itf0 != allFaces.end() && itf1 != allFaces.end() && itf2 != allFaces.end() && itf3 != allFaces.end() ){ - (*it)->quadrangles.push_back(new MQuadrangle(f.getVertex(0),f.getVertex(1),f.getVertex(2),f.getVertex(3))); + bool _empty; + public : + splitQuadRecovery() : _empty(true) {} + bool empty(){ return _empty; } + void setEmpty(bool val){ _empty = val; } + void add (const MFace &f, MVertex *v, GEntity*ge) + { + _data.insert(std::make_pair(ge, std::make_pair(v,f))); + } + int buildPyramids(GModel *gm) + { + if(empty()) return 0; + int NBPY = 0; + for (GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it){ + std::set<MFace, Less_Face> allFaces; + for (unsigned int i = 0; i < (*it)->triangles.size(); i++){ + allFaces.insert ((*it)->triangles[i]->getFace(0)); + delete (*it)->triangles[i]; + } + (*it)->triangles.clear(); + for (std::multimap<GEntity*, std::pair<MVertex*,MFace> >::iterator it2 = + _data.lower_bound(*it); it2 != _data.upper_bound(*it) ; ++it2){ + MVertex *v = it2->second.first; + v->onWhat()->mesh_vertices.erase(std::find(v->onWhat()->mesh_vertices.begin(), + v->onWhat()->mesh_vertices.end(), v)); + const MFace &f = it2->second.second; + std::set<MFace, Less_Face>::iterator itf0 = allFaces.find(MFace(f.getVertex(0), + f.getVertex(1),v)); + std::set<MFace, Less_Face>::iterator itf1 = allFaces.find(MFace(f.getVertex(1), + f.getVertex(2),v)); + std::set<MFace, Less_Face>::iterator itf2 = allFaces.find(MFace(f.getVertex(2), + f.getVertex(3),v)); + std::set<MFace, Less_Face>::iterator itf3 = allFaces.find(MFace(f.getVertex(3), + f.getVertex(0),v)); + if (itf0 != allFaces.end() && itf1 != allFaces.end() && + itf2 != allFaces.end() && itf3 != allFaces.end()){ + (*it)->quadrangles.push_back(new MQuadrangle(f.getVertex(0), f.getVertex(1), + f.getVertex(2), f.getVertex(3))); allFaces.erase(*itf0); allFaces.erase(*itf1); allFaces.erase(*itf2); allFaces.erase(*itf3); - - // printf("some pyramids should be created %d regions\n",(*it)->numRegions()); + printf("some pyramids should be created %d regions\n",(*it)->numRegions()); for (int iReg = 0; iReg < (*it)->numRegions(); iReg++){ - // printf("one pyramid is created\n"); - if (iReg == 1) {Msg::Fatal("cannot build pyramids on non manifold faces ..."); v = new MVertex (v->x(),v->y(),v->z(),(*it)->getRegion(iReg)); } - else v->setEntity ((*it)->getRegion(iReg)); - (*it)->getRegion(iReg)->pyramids.push_back(new MPyramid(f.getVertex(0),f.getVertex(1),f.getVertex(2),f.getVertex(3),v)); + printf("one pyramid is created\n"); + if (iReg == 1) { + Msg::Error("cannot build pyramids on non manifold faces ..."); + v = new MVertex(v->x(), v->y(), v->z(), (*it)->getRegion(iReg)); + } + else + v->setEntity ((*it)->getRegion(iReg)); + (*it)->getRegion(iReg)->pyramids.push_back + (new MPyramid(f.getVertex(0), f.getVertex(1), f.getVertex(2), f.getVertex(3), v)); (*it)->getRegion(iReg)->mesh_vertices.push_back(v); NBPY++; } + } + } + for (std::set<MFace, Less_Face>::iterator itf = allFaces.begin(); + itf != allFaces.end(); ++itf){ + (*it)->triangles.push_back + (new MTriangle(itf->getVertex(0), itf->getVertex(1), itf->getVertex(2))); } } - for (std::set<MFace,Less_Face> :: iterator itf = allFaces.begin() ; itf != allFaces.end(); ++itf){ - (*it)->triangles.push_back(new MTriangle(itf->getVertex(0),itf->getVertex(1),itf->getVertex(2))); - } + return NBPY; } - return NBPY; -} +}; void printVoronoi(GRegion *gr, std::set<SPoint3> &candidates) { @@ -194,7 +207,6 @@ void printVoronoi(GRegion *gr, std::set<SPoint3> &candidates) void skeletonFromVoronoi(GRegion *gr, std::set<SPoint3> &voronoiPoles) { - std::vector<MTetrahedron*> elements = gr->tetrahedra; std::list<GFace*> allFaces = gr->faces(); @@ -313,12 +325,12 @@ void skeletonFromVoronoi(GRegion *gr, std::set<SPoint3> &voronoiPoles) #endif } -void getBoundingInfoAndSplitQuads (GRegion *gr, - std::map<MFace,GEntity*,Less_Face> &allBoundingFaces, - std::set<MVertex*> &allBoundingVertices, - splitQuadRecovery &sqr){ - - std::map<MFace,GEntity*,Less_Face> allBoundingFaces_temp; +void getBoundingInfoAndSplitQuads(GRegion *gr, + std::map<MFace,GEntity*,Less_Face> &allBoundingFaces, + std::set<MVertex*> &allBoundingVertices, + splitQuadRecovery &sqr) +{ + std::map<MFace, GEntity*, Less_Face> allBoundingFaces_temp; // Get all the faces that are on the boundary std::list<GFace*> faces = gr->faces(); @@ -331,38 +343,38 @@ void getBoundingInfoAndSplitQuads (GRegion *gr, ++it; } - // if some elements pre-exist in the mesh, then use the - // internal faces of those + // if some elements pre-exist in the mesh, then use the internal faces of + // those for (unsigned int i=0;i<gr->getNumMeshElements();i++){ MElement *e = gr->getMeshElement(i); - for (int j=0;j<e->getNumFaces();j++){ - std::map<MFace,GEntity*,Less_Face>::iterator it = allBoundingFaces_temp.find(e->getFace(j)); - if (it == allBoundingFaces_temp.end())allBoundingFaces_temp[e->getFace(j)] = gr; + for (int j = 0; j < e->getNumFaces(); j++){ + std::map<MFace, GEntity*, Less_Face>::iterator it = allBoundingFaces_temp.find(e->getFace(j)); + if (it == allBoundingFaces_temp.end()) allBoundingFaces_temp[e->getFace(j)] = gr; else allBoundingFaces_temp.erase(it); } } - - std::map<MFace,GEntity*,Less_Face>::iterator itx = allBoundingFaces_temp.begin(); + + std::map<MFace, GEntity*, Less_Face>::iterator itx = allBoundingFaces_temp.begin(); for (; itx != allBoundingFaces_temp.end();++itx){ const MFace &f = itx->first; - // SPLIT that face into 4 triangular faces + // SPLIT that face into 4 triangular faces if (f.getNumVertices() == 4){ - MVertex *v0 = f.getVertex(0); - MVertex *v1 = f.getVertex(1); - MVertex *v2 = f.getVertex(2); - MVertex *v3 = f.getVertex(3); - + sqr.setEmpty(false); + MVertex *v0 = f.getVertex(0); + MVertex *v1 = f.getVertex(1); + MVertex *v2 = f.getVertex(2); + MVertex *v3 = f.getVertex(3); MVertex *newv = new MVertex ((v0->x() + v1->x() + v2->x() + v3->x())*0.25, (v0->y() + v1->y() + v2->y() + v3->y())*0.25, - (v0->z() + v1->z() + v2->z() + v3->z())*0.25,itx->second); + (v0->z() + v1->z() + v2->z() + v3->z())*0.25,itx->second); sqr.add(f,newv,itx->second); allBoundingFaces[MFace(v0,v1,newv)] = itx->second; allBoundingFaces[MFace(v1,v2,newv)] = itx->second; allBoundingFaces[MFace(v2,v3,newv)] = itx->second; allBoundingFaces[MFace(v3,v0,newv)] = itx->second; itx->second->mesh_vertices.push_back(newv); - // myGr->pyramids.push_back(new MPyramid(v0,v1,v2,v3,newv)); + // myGr->pyramids.push_back(new MPyramid(v0,v1,v2,v3,newv)); allBoundingVertices.insert(v0); allBoundingVertices.insert(v1); allBoundingVertices.insert(v2); @@ -375,10 +387,9 @@ void getBoundingInfoAndSplitQuads (GRegion *gr, allBoundingVertices.insert(f.getVertex(1)); allBoundingVertices.insert(f.getVertex(2)); } - } + } } - void getAllBoundingVertices(GRegion *gr, std::set<MVertex*> &allBoundingVertices) { std::list<GFace*> faces = gr->faces(); @@ -402,9 +413,7 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb { std::set<MVertex*> allBoundingVertices; std::map<MFace,GEntity*,Less_Face> allBoundingFaces; - getBoundingInfoAndSplitQuads (gr, allBoundingFaces,allBoundingVertices,sqr); - - //getAllBoundingVertices(gr, allBoundingVertices); + getBoundingInfoAndSplitQuads(gr, allBoundingFaces, allBoundingVertices, sqr); in.mesh_dim = 3; in.firstnumber = 1; @@ -426,26 +435,26 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb for(int i=0;i<Filler::get_nbr_new_vertices();i++){ MVertex* v; - v = Filler::get_new_vertex(i); - in.pointlist[(I - 1) * 3 + 0] = v->x(); - in.pointlist[(I - 1) * 3 + 1] = v->y(); - in.pointlist[(I - 1) * 3 + 2] = v->z(); - I++; + v = Filler::get_new_vertex(i); + in.pointlist[(I - 1) * 3 + 0] = v->x(); + in.pointlist[(I - 1) * 3 + 1] = v->y(); + in.pointlist[(I - 1) * 3 + 2] = v->z(); + I++; } for(int i=0;i<LpSmoother::get_nbr_interior_vertices();i++){ MVertex* v; - v = LpSmoother::get_interior_vertex(i); - in.pointlist[(I - 1) * 3 + 0] = v->x(); - in.pointlist[(I - 1) * 3 + 1] = v->y(); - in.pointlist[(I - 1) * 3 + 2] = v->z(); - I++; + v = LpSmoother::get_interior_vertex(i); + in.pointlist[(I - 1) * 3 + 0] = v->x(); + in.pointlist[(I - 1) * 3 + 1] = v->y(); + in.pointlist[(I - 1) * 3 + 2] = v->z(); + I++; } in.numberoffacets = allBoundingFaces.size(); in.facetlist = new tetgenio::facet[in.numberoffacets]; in.facetmarkerlist = new int[in.numberoffacets]; - + I = 0; std::map<MFace,GEntity*,Less_Face>::iterator it = allBoundingFaces.begin(); for (; it != allBoundingFaces.end();++it){ @@ -463,8 +472,8 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb p->vertexlist[1] = fac.getVertex(1)->getIndex(); p->vertexlist[2] = fac.getVertex(2)->getIndex(); in.facetmarkerlist[I] = (it->second->dim() == 3) ? -it->second->tag() : it->second->tag(); - ++I; - } + ++I; + } } void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out, @@ -522,7 +531,7 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out, v[2] = numberedV[out.trifacelist[i * 3 + 2] - 1]; // do not recover prismatic faces !!! GFace *gf = gr->model()->getFaceByTag(out.trifacemarkerlist[i]); - + double guess[2] = {0, 0}; if (needParam) { int Count = 0; @@ -562,7 +571,7 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out, guess[1] /= Count; } } - + for(int j = 0; j < 3; j++){ if(v[j]->onWhat()->dim() == 3){ v[j]->onWhat()->mesh_vertices.erase @@ -587,7 +596,7 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out, else{ v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(), gf); } - + gf->mesh_vertices.push_back(v1b); numberedV[out.trifacelist[i * 3 + j] - 1] = v1b; delete v[j]; @@ -761,7 +770,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions) } else if(!Filler::get_nbr_new_vertices() && !LpSmoother::get_nbr_interior_vertices()){ - insertVerticesInRegion(gr); + insertVerticesInRegion(gr); } if (sqr.buildPyramids (gr->model())){ @@ -772,7 +781,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions) buildVertexToElement(regions[k]->pyramids, adj); v2t_cont::iterator it = adj.begin(); while (it != adj.end()){ - //_relocateVertex( it->first, it->second); + _relocateVertex( it->first, it->second); ++it; } }