From 84947cb6d8a74063d1c65287afad38b7fe67707d Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 25 Apr 2010 10:48:49 +0000
Subject: [PATCH] fix mesh workflow
---
Geo/GEdge.cpp | 17 ++++++-------
Geo/GEdge.h | 4 ++-
Geo/GFace.cpp | 10 ++++----
Geo/GFace.h | 11 ++++----
Geo/GRegion.cpp | 18 ++++++++------
Geo/GRegion.h | 5 ++--
Mesh/Generator.cpp | 43 ++++++++++++++------------------
Mesh/meshGEdge.cpp | 15 ++++++-----
Mesh/meshGFace.cpp | 51 +++++++++++++++-----------------------
Post/PViewDataGModelIO.cpp | 4 +++
10 files changed, 85 insertions(+), 93 deletions(-)
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 0ed4d7de19..a98a9ede1e 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -21,7 +21,7 @@ GEdge::GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1)
{
if(v0) v0->addEdge(this);
if(v1 && v1 != v0) v1->addEdge(this);
- _mStatus = GEdge::PENDING;
+ meshStatistics.status = GEdge::PENDING;
resetMeshAttributes();
}
@@ -45,12 +45,10 @@ void GEdge::deleteMesh()
void GEdge::reverse()
{
- GVertex* tmp=v0;
- v0=v1;
- v1=tmp;
- for (std::vector<MLine*>::iterator line = lines.begin();
- line != lines.end();
- line++)
+ GVertex* tmp = v0;
+ v0 = v1;
+ v1 = tmp;
+ for(std::vector<MLine*>::iterator line = lines.begin(); line != lines.end(); line++)
(*line)->revert();
}
@@ -347,8 +345,9 @@ bool GEdge::XYZToU(const double X, const double Y, const double Z,
return false;
}
-void GEdge::replaceEndingPoints (GVertex *replOfv0, GVertex *replOfv1){
- replaceEndingPointsInternals (replOfv0,replOfv1);
+void GEdge::replaceEndingPoints (GVertex *replOfv0, GVertex *replOfv1)
+{
+ replaceEndingPointsInternals (replOfv0, replOfv1);
if (replOfv0 != v0){
v0->delEdge(this);
replOfv0->addEdge(this);
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 9123c5b580..f1acad2c81 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -191,7 +191,9 @@ class GEdge : public GEntity {
} meshAttributes ;
typedef enum {PENDING, DONE, FAILED} meshGenerationStatus;
- mutable meshGenerationStatus _mStatus;
+ struct {
+ mutable meshGenerationStatus status;
+ } meshStatistics;
std::vector<MLine*> lines;
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index a10fcdd155..6e2174bcce 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1049,19 +1049,20 @@ bool GFace::fillPointCloud(double maxDist, std::vector<SPoint3> *points,
return true;
}
-void GFace::lloyd(int nbiter, int infn){
+void GFace::lloyd(int nbiter, int infn)
+{
lloydAlgorithm algo(nbiter, infn);
algo(this);
}
- // replace edges (gor gluing)
-void GFace::replaceEdges (std::list<GEdge*> & new_edges) {
+void GFace::replaceEdges (std::list<GEdge*> &new_edges)
+{
replaceEdgesInternal (new_edges);
std::list<GEdge*>::iterator it = l_edges.begin();
std::list<GEdge*>::iterator it2 = new_edges.begin();
std::list<int>::iterator it3 = l_dirs.begin();
std::list<int> newdirs;
- for ( ; it != l_edges.end() ; ++it,++it2,++it3){
+ for ( ; it != l_edges.end(); ++it, ++it2, ++it3){
(*it)->delFace(this);
(*it2)->addFace(this);
if ((*it2)->getBeginVertex() == (*it)->getBeginVertex())
@@ -1073,7 +1074,6 @@ void GFace::replaceEdges (std::list<GEdge*> & new_edges) {
l_dirs = newdirs;
}
-
#include "Bindings.h"
void GFace::registerBindings(binding *b)
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 0bc5ee2ad4..23447b4b2d 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -50,8 +50,9 @@ class GFace : public GEntity
std::list<GEdge *> embedded_edges;
std::list<GVertex *> embedded_vertices;
GFaceCompound *compound; // this model edge belongs to a compound
- // replace edges (gor gluing) for specific modelers, we have to re-create
- // internal data ...
+
+ // replace edges (gor gluing) for specific modelers, we have to
+ // re-create internal data ...
virtual void replaceEdgesInternal (std::list<GEdge*> &){}
public: // this will become protected or private
@@ -161,7 +162,8 @@ class GFace : public GEntity
virtual bool containsParam(const SPoint2 &pt) const;
// return the point on the face closest to the given point
- virtual GPoint closestPoint(const SPoint3 & queryPoint, const double initialGuess[2]) const;
+ virtual GPoint closestPoint(const SPoint3 & queryPoint,
+ const double initialGuess[2]) const;
// return the normal to the face at the given parameter location
virtual SVector3 normal(const SPoint2 ¶m) const;
@@ -175,8 +177,7 @@ class GFace : public GEntity
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const = 0;
// return the curvature computed as the divergence of the normal
- inline double curvature(const SPoint2 ¶m) const
- {return curvatureMax(param);}
+ inline double curvature(const SPoint2 ¶m) const { return curvatureMax(param); }
virtual double curvatureDiv(const SPoint2 ¶m) const;
// return the maximum curvature at a point
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
index 10b67155ba..d139779d3e 100644
--- a/Geo/GRegion.cpp
+++ b/Geo/GRegion.cpp
@@ -96,10 +96,13 @@ MElement *GRegion::getMeshElement(unsigned int index)
return hexahedra[index - tetrahedra.size()];
else if(index < tetrahedra.size() + hexahedra.size() + prisms.size())
return prisms[index - tetrahedra.size() - hexahedra.size()];
- else if(index < tetrahedra.size() + hexahedra.size() + prisms.size() + pyramids.size())
+ else if(index < tetrahedra.size() + hexahedra.size() + prisms.size() +
+ pyramids.size())
return pyramids[index - tetrahedra.size() - hexahedra.size() - prisms.size()];
- else if(index < tetrahedra.size() + hexahedra.size() + prisms.size() + pyramids.size() + polyhedra.size())
- return polyhedra[index - tetrahedra.size() - hexahedra.size() - prisms.size() - pyramids.size()];
+ else if(index < tetrahedra.size() + hexahedra.size() + prisms.size() +
+ pyramids.size() + polyhedra.size())
+ return polyhedra[index - tetrahedra.size() - hexahedra.size() - prisms.size() -
+ pyramids.size()];
return 0;
}
@@ -260,7 +263,8 @@ bool GRegion::edgeConnected(GRegion *r) const
}
// replace faces (gor gluing)
-void GRegion::replaceFaces (std::list<GFace*> & new_faces) {
+void GRegion::replaceFaces (std::list<GFace*> & new_faces)
+{
replaceFacesInternal (new_faces);
std::list<GFace*>::iterator it = l_faces.begin();
std::list<GFace*>::iterator it2 = new_faces.begin();
@@ -269,10 +273,10 @@ void GRegion::replaceFaces (std::list<GFace*> & new_faces) {
for ( ; it != l_faces.end() ; ++it,++it2,++it3){
(*it)->delRegion(this);
(*it2)->addRegion(this);
- // if ((*it2)->getBeginVertex() == (*it)->getBeginVertex())
+ // if ((*it2)->getBeginVertex() == (*it)->getBeginVertex())
newdirs.push_back(*it3);
- // else
- // newdirs.push_back(-(*it3));
+ // else
+ // newdirs.push_back(-(*it3));
}
l_faces = new_faces;
l_dirs = newdirs;
diff --git a/Geo/GRegion.h b/Geo/GRegion.h
index 5105ad9acd..c72e768ae9 100644
--- a/Geo/GRegion.h
+++ b/Geo/GRegion.h
@@ -26,8 +26,9 @@ class GRegion : public GEntity {
protected:
std::list<GFace*> l_faces;
std::list<int> l_dirs;
- // replace faces (for gluing) for specific modelers, we have to re-create
- // internal data ...
+
+ // replace faces (for gluing) for specific modelers, we have to
+ // re-create internal data ...
virtual void replaceFacesInternal (std::list<GFace*> &) {}
public:
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index f9759546a9..4085019b47 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -339,12 +339,15 @@ static void Mesh1D(GModel *m)
Msg::StatusBar(1, true, "Meshing 1D...");
double t1 = Cpu();
+ for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
+ (*it)->meshStatistics.status = GEdge::PENDING;
+
int nIter = 0;
while(1){
meshGEdge mesher;
int nbPending = 0;
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it){
- if ((*it)->_mStatus == GEdge::PENDING){
+ if ((*it)->meshStatistics.status == GEdge::PENDING){
mesher(*it);
nbPending++;
}
@@ -353,8 +356,6 @@ static void Mesh1D(GModel *m)
if(nIter++ > 10) break;
}
- // std::for_each(m->firstEdge(), m->lastEdge(), meshGEdge());
-
double t2 = Cpu();
CTX::instance()->meshTimer[0] = t2 - t1;
Msg::Info("Mesh 1D complete (%g s)", CTX::instance()->meshTimer[0]);
@@ -414,6 +415,9 @@ static void Mesh2D(GModel *m)
Msg::StatusBar(1, true, "Meshing 2D...");
double t1 = Cpu();
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
+ (*it)->meshStatistics.status = GFace::PENDING;
+
// skip short mesh edges
//geomTresholdVertexEquivalence inst(m);
@@ -421,41 +425,30 @@ static void Mesh2D(GModel *m)
// and curve meshes) is global as it depends on a smooth normal
// field generated from the surface mesh of the source surfaces
if(!Mesh2DWithBoundaryLayers(m)){
-
- std::set<GFace*> classFaces;
- std::set<GFace*> compFaces;
+
+ std::set<GFace*> faces, compoundFaces;
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
- if ((*it)->geomType() != GEntity::CompoundSurface)
- classFaces.insert(*it);
+ if ((*it)->geomType() == GEntity::CompoundSurface)
+ compoundFaces.insert(*it);
else
- compFaces.insert(*it);
+ faces.insert(*it);
}
- std::for_each(classFaces.begin(), classFaces.end(), meshGFace());
- std::for_each(compFaces.begin(), compFaces.end(), meshGFace());
+ std::for_each(faces.begin(), faces.end(), meshGFace());
+ std::for_each(compoundFaces.begin(), compoundFaces.end(), meshGFace());
// lloyd optimization
if (CTX::instance()->mesh.optimizeLloyd){
- Msg::Info("------------------------------");
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
GFace *gf = *it;
if(gf->geomType() == GEntity::DiscreteSurface) continue;
- if(gf->geomType() == GEntity::CompoundSurface ) {
+ if(gf->geomType() == GEntity::CompoundSurface) {
GFaceCompound *gfc = (GFaceCompound*) gf;
- if (gfc->getNbSplit() != 0) continue;
+ if(gfc->getNbSplit() != 0) continue;
}
int recombine = gf->meshAttributes.recombine;
-
Msg::StatusBar(2, true, "Lloyd optimization for face %d", gf->tag());
- gf->lloyd(25,recombine);
-
- if(recombine) recombineIntoQuads(gf);
-
- // computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
-// gf->meshStatistics.average_element_shape,
-// gf->meshStatistics.best_element_shape,
-// gf->meshStatistics.nbTriangle,
-// gf->meshStatistics.nbGoodQuality);
-
+ gf->lloyd(25, recombine);
+ if(recombine) recombineIntoQuads(gf);
}
}
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 63625ce1e3..06cea7868f 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -256,6 +256,7 @@ void deMeshGEdge::operator() (GEdge *ge)
{
if(ge->geomType() == GEntity::DiscreteCurve) return;
ge->deleteMesh();
+ ge->meshStatistics.status = GEdge::PENDING;
}
void meshGEdge::operator() (GEdge *ge)
@@ -278,14 +279,14 @@ void meshGEdge::operator() (GEdge *ge)
if (ge->meshMaster() != ge->tag()){
GEdge *gef = ge->model()->getEdgeByTag(abs(ge->meshMaster()));
- if (gef->_mStatus == GEdge::PENDING)return;
- Msg::Info("Meshing curve %d (%s) as a copy of %d",ge->tag(),ge->getTypeString().c_str(),ge->meshMaster());
- copyMesh(gef,ge,ge->meshMaster());
- ge->_mStatus = GEdge::DONE;
+ if (gef->meshStatistics.status == GEdge::PENDING) return;
+ Msg::Info("Meshing curve %d (%s) as a copy of %d", ge->tag(),
+ ge->getTypeString().c_str(), ge->meshMaster());
+ copyMesh(gef, ge, ge->meshMaster());
+ ge->meshStatistics.status = GEdge::DONE;
return;
}
-
Msg::Info("Meshing curve %d (%s)", ge->tag(), ge->getTypeString().c_str());
// compute bounds
@@ -401,7 +402,5 @@ void meshGEdge::operator() (GEdge *ge)
v0->y() = beg_p.y();
v0->z() = beg_p.z();
}
- ge->_mStatus = GEdge::DONE;
+ ge->meshStatistics.status = GEdge::DONE;
}
-
-
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index b5bec6d868..d189bd41ba 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -38,8 +38,8 @@
#include "multiscalePartition.h"
#include "meshGFaceLloyd.h"
-static void copyMesh (GFace *source, GFace *target){
-
+static void copyMesh (GFace *source, GFace *target)
+{
std::map<MVertex*,MVertex*> vs2vt;
std::list<GEdge*> edges = target->edges();
@@ -47,17 +47,16 @@ static void copyMesh (GFace *source, GFace *target){
GEdge *te = *it;
int master = te->meshMaster();
if (master == te->tag()){
- Msg::Error ("Periodic face %d does not have periodic edges (master %d -- edge %d)",target->tag(),master,te->tag());
+ Msg::Error("Periodic face %d does not have periodic edges (master %d -- edge %d)",
+ target->tag(), master, te->tag());
}
GEdge *se = source->model()->getEdgeByTag(abs(master));
- // printf("%d %d\n",se->tag(),te->tag());
if (master > 0){
vs2vt[se->getBeginVertex()->mesh_vertices[0]] = te->getBeginVertex()->mesh_vertices[0];
vs2vt[se->getEndVertex()->mesh_vertices[0]] = te->getEndVertex()->mesh_vertices[0];
for (int i=0;i<se->mesh_vertices.size();i++){
MVertex *vs = se->mesh_vertices[i];
MVertex *vt = te->mesh_vertices[i];
- // printf("D %g %g %g vs %g %g %g\n",vs->x(),vs->y(),vs->z(),vt->x(),vt->y(),vt->z());
vs2vt[vs] = vt;
}
}
@@ -67,7 +66,6 @@ static void copyMesh (GFace *source, GFace *target){
for (int i=0;i<se->mesh_vertices.size();i++){
MVertex *vs = se->mesh_vertices[i];
MVertex *vt = te->mesh_vertices[se->mesh_vertices.size()-i-1];
- // printf("R %g %g %g vs %g %g %g\n",vs->x(),vs->y(),vs->z(),vt->x(),vt->y(),vt->z());
vs2vt[vs] = vt;
}
}
@@ -83,12 +81,11 @@ static void copyMesh (GFace *source, GFace *target){
if (vs->onWhat()->dim() == 1){
bool success1 = reparamMeshVertexOnFace(vs, source, param_source[count]);
bool success2 = reparamMeshVertexOnFace(vt, target, param_target[count++]);
- // printf("%g %g %g vs %g %g %g\n",vs->x(),vs->y(),vs->z(),vt->x(),vt->y(),vt->z());
- if (count == 2)break;
+ if (count == 2) break;
}
}
- if (count < 2)return;
+ if (count < 2) return;
const double t1u = param_target[0].x();
const double t1v = param_target[0].y();
@@ -104,21 +101,10 @@ static void copyMesh (GFace *source, GFace *target){
SVector3 _c = crossprod(_a,_b);
double sinA = _c.z();
double cosA = dot(_a,_b);
- // printf("%g %g-- %g %g\n",_a.x(),_a.y(),_b.x(),_b.y());
const double theta = atan2(sinA, cosA);
const double c = cos(theta);
const double s = sin(theta);
- // printf("s1 %g %g s2 %g %g\n",s1u,s1v,s2u,s2v);
- // printf("t1 %g %g t2 %g %g\n",t1u,t1v,t2u,t2v);
- // printf("theta = %g\n",theta*180/M_PI);
- // {
- // double u = param_source[1].x();
- // double v = param_source[1].y();
- // printf("transfo of s2 = %g %g\n",c * (u-s1u) + s * (v-s1v) + t1u,
- // -s * (u-s1u) + c * (v-s1v) + t1v);
- // }
-
for(unsigned int i = 0; i < source->mesh_vertices.size(); i++){
MVertex *vs = source->mesh_vertices[i];
double u,v;
@@ -146,10 +132,8 @@ static void copyMesh (GFace *source, GFace *target){
MVertex *v4 = vs2vt[source->quadrangles[i]->getVertex(3)];
target->quadrangles.push_back(new MQuadrangle(v1,v2,v3,v4));
}
-
}
-
void fourthPoint(double *p1, double *p2, double *p3, double *p4)
{
double c[3];
@@ -1372,14 +1356,13 @@ void meshGFace::operator() (GFace *gf)
if(MeshTransfiniteSurface(gf)) return;
if(MeshExtrudedSurface(gf)) return;
if(gf->meshMaster() != gf->tag()){
- printf("AAAAAAAAARGH %d %d\n",gf->meshMaster(),gf->tag());
GFace *gff = gf->model()->getFaceByTag(abs(gf->meshMaster()));
if (gff->meshStatistics.status != GFace::DONE){
- // Msg::Info("Meshing face %d (%s) as a copy of %d",gf->tag(),gf->getTypeString().c_str(),gf->meshMaster());
gf->meshStatistics.status = GFace::PENDING;
return;
}
- Msg::Info("Meshing face %d (%s) as a copy of %d",gf->tag(),gf->getTypeString().c_str(),gf->meshMaster());
+ Msg::Info("Meshing face %d (%s) as a copy of %d", gf->tag(),
+ gf->getTypeString().c_str(), gf->meshMaster());
copyMesh(gff,gf);
gf->meshStatistics.status = GFace::DONE;
return;
@@ -1602,14 +1585,17 @@ void partitionAndRemesh(GFaceCompound *gf)
std::vector<MVertex*> edge_vertices = (*it)->mesh_vertices;
std::vector<MVertex*>::const_iterator itv = edge_vertices.begin();
for(; itv != edge_vertices.end(); itv++){
- std::vector<MVertex*>::iterator itve = std::find(gf->mesh_vertices.begin(), gf->mesh_vertices.end(), *itv);
+ std::vector<MVertex*>::iterator itve = std::find
+ (gf->mesh_vertices.begin(), gf->mesh_vertices.end(), *itv);
if (itve != gf->mesh_vertices.end()) gf->mesh_vertices.erase(itve);
}
MVertex *vB = (*it)->getBeginVertex()->mesh_vertices[0];
- std::vector<MVertex*>::iterator itvB = std::find(gf->mesh_vertices.begin(), gf->mesh_vertices.end(), vB);
+ std::vector<MVertex*>::iterator itvB = std::find
+ (gf->mesh_vertices.begin(), gf->mesh_vertices.end(), vB);
if (itvB != gf->mesh_vertices.end()) gf->mesh_vertices.erase(itvB);
MVertex *vE = (*it)->getEndVertex()->mesh_vertices[0];
- std::vector<MVertex*>::iterator itvE = std::find(gf->mesh_vertices.begin(), gf->mesh_vertices.end(), vE);
+ std::vector<MVertex*>::iterator itvE = std::find
+ (gf->mesh_vertices.begin(), gf->mesh_vertices.end(), vE);
if (itvE != gf->mesh_vertices.end()) gf->mesh_vertices.erase(itvE);
//if l_edge is a compond
@@ -1618,14 +1604,17 @@ void partitionAndRemesh(GFaceCompound *gf)
std::vector<MVertex*> edge_vertices = gec->mesh_vertices;
std::vector<MVertex*>::const_iterator itv = edge_vertices.begin();
for(; itv != edge_vertices.end(); itv++){
- std::vector<MVertex*>::iterator itve = std::find(gf->mesh_vertices.begin(), gf->mesh_vertices.end(), *itv);
+ std::vector<MVertex*>::iterator itve = std::find
+ (gf->mesh_vertices.begin(), gf->mesh_vertices.end(), *itv);
if (itve != gf->mesh_vertices.end()) gf->mesh_vertices.erase(itve);
}
MVertex *vB = (*it)->getBeginVertex()->mesh_vertices[0];
- std::vector<MVertex*>::iterator itvB = std::find(gf->mesh_vertices.begin(), gf->mesh_vertices.end(), vB);
+ std::vector<MVertex*>::iterator itvB = std::find
+ (gf->mesh_vertices.begin(), gf->mesh_vertices.end(), vB);
if (itvB != gf->mesh_vertices.end()) gf->mesh_vertices.erase(itvB);
MVertex *vE = (*it)->getEndVertex()->mesh_vertices[0];
- std::vector<MVertex*>::iterator itvE = std::find(gf->mesh_vertices.begin(), gf->mesh_vertices.end(), vE);
+ std::vector<MVertex*>::iterator itvE = std::find
+ (gf->mesh_vertices.begin(), gf->mesh_vertices.end(), vE);
if (itvE != gf->mesh_vertices.end()) gf->mesh_vertices.erase(itvE);
}
}
diff --git a/Post/PViewDataGModelIO.cpp b/Post/PViewDataGModelIO.cpp
index 69ed6c2ca8..97a266e439 100644
--- a/Post/PViewDataGModelIO.cpp
+++ b/Post/PViewDataGModelIO.cpp
@@ -102,6 +102,10 @@ bool PViewDataGModel::readMSH(std::string fileName, int fileIndex, FILE *fp,
_steps[step]->getPartitions().insert(partition);
+ // FIXME: we should do this at a higher-level, since this will be
+ // very slow for large multi-step, multi-partition datasets (we
+ // recompute the min/max for all the previously loaded
+ // steps/partitions -> loop over all elements many many times...)
finalize();
return true;
}
--
GitLab