Skip to content
Snippets Groups Projects
Commit bac66ad2 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

fix pyramids

parent 54cfe24f
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/gmshEdge.cpp
+ 3
− 2
View file @ bac66ad2
...@@ -288,9 +288,10 @@ SPoint2 gmshEdge::reparamOnFace(const GFace *face, double epar,int dir) const ...@@ -288,9 +288,10 @@ SPoint2 gmshEdge::reparamOnFace(const GFace *face, double epar,int dir) const
} }
return SPoint2(U, V); return SPoint2(U, V);
} }
else else{
return GEdge::reparamOnFace(face, epar, dir); return GEdge::reparamOnFace(face, epar, dir);
} }
}
void gmshEdge::writeGEO(FILE *fp) void gmshEdge::writeGEO(FILE *fp)
{ {
......
Geo/gmshFace.cpp
+ 18
− 18
View file @ bac66ad2
Mesh/meshGFace.cpp
+ 9
− 9
View file @ bac66ad2
Mesh/meshGRegion.cpp
+ 97
− 88
View file @ bac66ad2
...@@ -34,60 +34,73 @@ ...@@ -34,60 +34,73 @@
// hybrid mesh recovery structure // hybrid mesh recovery structure
class splitQuadRecovery { class splitQuadRecovery {
std::multimap<GEntity*, std::pair<MVertex*,MFace> >_data; std::multimap<GEntity*, std::pair<MVertex*,MFace> >_data;
bool _empty;
public : public :
void add (const MFace &f, MVertex *v, GEntity*ge) ; splitQuadRecovery() : _empty(true) {}
int buildPyramids (GModel *gm); bool empty(){ return _empty; }
}; void setEmpty(bool val){ _empty = val; }
void add (const MFace &f, MVertex *v, GEntity*ge)
void splitQuadRecovery::add (const MFace &f, MVertex *v, GEntity*ge) { {
_data.insert(std::make_pair(ge, std::make_pair(v,f))); _data.insert(std::make_pair(ge, std::make_pair(v,f)));
} }
int buildPyramids(GModel *gm)
// re-create quads and {
int splitQuadRecovery::buildPyramids (GModel *gm) { if(empty()) return 0;
int NBPY = 0; int NBPY = 0;
std::set<MFace,Less_Face> allFaces;
for (GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it){ 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++){ for (unsigned int i = 0; i < (*it)->triangles.size(); i++){
allFaces.insert ((*it)->triangles[i]->getFace(0)); allFaces.insert ((*it)->triangles[i]->getFace(0));
delete (*it)->triangles[i]; delete (*it)->triangles[i];
} }
(*it)->triangles.clear(); (*it)->triangles.clear();
for ( std::multimap<GEntity*, std::pair<MVertex*,MFace> >::iterator it2 = _data.lower_bound(*it) ; it2 != _data.upper_bound(*it) ; ++it2) { 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; MVertex *v = it2->second.first;
v->onWhat()->mesh_vertices.erase v->onWhat()->mesh_vertices.erase(std::find(v->onWhat()->mesh_vertices.begin(),
(std::find(v->onWhat()->mesh_vertices.begin(), v->onWhat()->mesh_vertices.end(), v));
v->onWhat()->mesh_vertices.end(),
v));
const MFace &f = it2->second.second; 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 itf0 = allFaces.find(MFace(f.getVertex(0),
std::set<MFace,Less_Face> :: iterator itf1 = allFaces.find(MFace(f.getVertex(1),f.getVertex(2),v)); f.getVertex(1),v));
std::set<MFace,Less_Face> :: iterator itf2 = allFaces.find(MFace(f.getVertex(2),f.getVertex(3),v)); std::set<MFace, Less_Face>::iterator itf1 = allFaces.find(MFace(f.getVertex(1),
std::set<MFace,Less_Face> :: iterator itf3 = allFaces.find(MFace(f.getVertex(3),f.getVertex(0),v)); f.getVertex(2),v));
if (itf0 != allFaces.end() && itf1 != allFaces.end() && itf2 != allFaces.end() && itf3 != allFaces.end() ){ std::set<MFace, Less_Face>::iterator itf2 = allFaces.find(MFace(f.getVertex(2),
(*it)->quadrangles.push_back(new MQuadrangle(f.getVertex(0),f.getVertex(1),f.getVertex(2),f.getVertex(3))); 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(*itf0);
allFaces.erase(*itf1); allFaces.erase(*itf1);
allFaces.erase(*itf2); allFaces.erase(*itf2);
allFaces.erase(*itf3); 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++){ for (int iReg = 0; iReg < (*it)->numRegions(); iReg++){
// printf("one pyramid is created\n"); 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)); } if (iReg == 1) {
else v->setEntity ((*it)->getRegion(iReg)); Msg::Error("cannot build pyramids on non manifold faces ...");
(*it)->getRegion(iReg)->pyramids.push_back(new MPyramid(f.getVertex(0),f.getVertex(1),f.getVertex(2),f.getVertex(3),v)); 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); (*it)->getRegion(iReg)->mesh_vertices.push_back(v);
NBPY++; NBPY++;
} }
} }
} }
for (std::set<MFace,Less_Face> :: iterator itf = allFaces.begin() ; itf != allFaces.end(); ++itf){ for (std::set<MFace, Less_Face>::iterator itf = allFaces.begin();
(*it)->triangles.push_back(new MTriangle(itf->getVertex(0),itf->getVertex(1),itf->getVertex(2))); 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) void printVoronoi(GRegion *gr, std::set<SPoint3> &candidates)
{ {
...@@ -194,7 +207,6 @@ 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) void skeletonFromVoronoi(GRegion *gr, std::set<SPoint3> &voronoiPoles)
{ {
std::vector<MTetrahedron*> elements = gr->tetrahedra; std::vector<MTetrahedron*> elements = gr->tetrahedra;
std::list<GFace*> allFaces = gr->faces(); std::list<GFace*> allFaces = gr->faces();
...@@ -316,8 +328,8 @@ void skeletonFromVoronoi(GRegion *gr, std::set<SPoint3> &voronoiPoles) ...@@ -316,8 +328,8 @@ void skeletonFromVoronoi(GRegion *gr, std::set<SPoint3> &voronoiPoles)
void getBoundingInfoAndSplitQuads(GRegion *gr, void getBoundingInfoAndSplitQuads(GRegion *gr,
std::map<MFace,GEntity*,Less_Face> &allBoundingFaces, std::map<MFace,GEntity*,Less_Face> &allBoundingFaces,
std::set<MVertex*> &allBoundingVertices, std::set<MVertex*> &allBoundingVertices,
splitQuadRecovery &sqr){ splitQuadRecovery &sqr)
{
std::map<MFace, GEntity*, Less_Face> allBoundingFaces_temp; std::map<MFace, GEntity*, Less_Face> allBoundingFaces_temp;
// Get all the faces that are on the boundary // Get all the faces that are on the boundary
...@@ -331,8 +343,8 @@ void getBoundingInfoAndSplitQuads (GRegion *gr, ...@@ -331,8 +343,8 @@ void getBoundingInfoAndSplitQuads (GRegion *gr,
++it; ++it;
} }
// if some elements pre-exist in the mesh, then use the // if some elements pre-exist in the mesh, then use the internal faces of
// internal faces of those // those
for (unsigned int i=0;i<gr->getNumMeshElements();i++){ for (unsigned int i=0;i<gr->getNumMeshElements();i++){
MElement *e = gr->getMeshElement(i); MElement *e = gr->getMeshElement(i);
...@@ -348,11 +360,11 @@ void getBoundingInfoAndSplitQuads (GRegion *gr, ...@@ -348,11 +360,11 @@ void getBoundingInfoAndSplitQuads (GRegion *gr,
const MFace &f = itx->first; const MFace &f = itx->first;
// SPLIT that face into 4 triangular faces // SPLIT that face into 4 triangular faces
if (f.getNumVertices() == 4){ if (f.getNumVertices() == 4){
sqr.setEmpty(false);
MVertex *v0 = f.getVertex(0); MVertex *v0 = f.getVertex(0);
MVertex *v1 = f.getVertex(1); MVertex *v1 = f.getVertex(1);
MVertex *v2 = f.getVertex(2); MVertex *v2 = f.getVertex(2);
MVertex *v3 = f.getVertex(3); MVertex *v3 = f.getVertex(3);
MVertex *newv = new MVertex ((v0->x() + v1->x() + v2->x() + v3->x())*0.25, 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->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);
...@@ -378,7 +390,6 @@ void getBoundingInfoAndSplitQuads (GRegion *gr, ...@@ -378,7 +390,6 @@ void getBoundingInfoAndSplitQuads (GRegion *gr,
} }
} }
void getAllBoundingVertices(GRegion *gr, std::set<MVertex*> &allBoundingVertices) void getAllBoundingVertices(GRegion *gr, std::set<MVertex*> &allBoundingVertices)
{ {
std::list<GFace*> faces = gr->faces(); std::list<GFace*> faces = gr->faces();
...@@ -404,8 +415,6 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb ...@@ -404,8 +415,6 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb
std::map<MFace,GEntity*,Less_Face> allBoundingFaces; std::map<MFace,GEntity*,Less_Face> allBoundingFaces;
getBoundingInfoAndSplitQuads(gr, allBoundingFaces, allBoundingVertices, sqr); getBoundingInfoAndSplitQuads(gr, allBoundingFaces, allBoundingVertices, sqr);
//getAllBoundingVertices(gr, allBoundingVertices);
in.mesh_dim = 3; in.mesh_dim = 3;
in.firstnumber = 1; in.firstnumber = 1;
in.numberofpoints = allBoundingVertices.size() + Filler::get_nbr_new_vertices() + in.numberofpoints = allBoundingVertices.size() + Filler::get_nbr_new_vertices() +
...@@ -772,7 +781,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> &regions) ...@@ -772,7 +781,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> &regions)
buildVertexToElement(regions[k]->pyramids, adj); buildVertexToElement(regions[k]->pyramids, adj);
v2t_cont::iterator it = adj.begin(); v2t_cont::iterator it = adj.begin();
while (it != adj.end()){ while (it != adj.end()){
//_relocateVertex( it->first, it->second); _relocateVertex( it->first, it->second);
++it; ++it;
} }
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment