From f4215e05f47b4118b4c712b1e9407c4092608fbd Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Fri, 28 Aug 2009 06:13:43 +0000
Subject: [PATCH] pp
---
Mesh/meshGFace.cpp | 48 ++--
Mesh/meshGFaceBDS.cpp | 23 --
Mesh/meshGFaceBDS.h | 2 -
Mesh/meshGFaceDelaunayInsertion.cpp | 8 +-
Mesh/meshGFaceOptimize.cpp | 16 +-
Mesh/meshGFaceQuadrilateralize.cpp | 389 ++++++++++++++--------------
6 files changed, 232 insertions(+), 254 deletions(-)
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index f194b125f9..0d60ef52b3 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -47,7 +47,7 @@ void fourthPoint(double *p1, double *p2, double *p3, double *p4)
p4[2] = c[2] + R * vz[2];
}
-static bool noseam(GFace *gf)
+static bool noSeam(GFace *gf)
{
std::list<GEdge*> edges = gf->edges();
std::list<GEdge*>::iterator it = edges.begin();
@@ -150,16 +150,16 @@ static void remeshUnrecoveredEdges(std::map<MVertex*, BDS_Point*> &recoverMapInv
}
}
-static bool AlgoDelaunay2D(GFace *gf)
+static bool algoDelaunay2D(GFace *gf)
{
- if(noseam(gf) && (CTX::instance()->mesh.algo2d == ALGO_2D_DELAUNAY ||
+ if(noSeam(gf) && (CTX::instance()->mesh.algo2d == ALGO_2D_DELAUNAY ||
CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL))
return true;
return false;
}
-void computeElementShapes(GFace *gf, double &worst, double &avg, double &best,
- int &nT, int &greaterThan)
+static void computeElementShapes(GFace *gf, double &worst, double &avg,
+ double &best, int &nT, int &greaterThan)
{
worst = 1.e22;
avg = 0.0;
@@ -177,13 +177,13 @@ void computeElementShapes(GFace *gf, double &worst, double &avg, double &best,
avg /= nT;
}
-static bool recover_medge(BDS_Mesh *m, GEdge *ge,
- std::map<MVertex*, BDS_Point*> &recoverMapInv,
- std::set<EdgeToRecover> *e2r,
- std::set<EdgeToRecover> *not_recovered, int pass_)
+static bool recoverEdge(BDS_Mesh *m, GEdge *ge,
+ std::map<MVertex*, BDS_Point*> &recoverMapInv,
+ std::set<EdgeToRecover> *e2r,
+ std::set<EdgeToRecover> *notRecovered, int pass)
{
BDS_GeomEntity *g = 0;
- if(pass_ == 2){
+ if(pass == 2){
m->add_geom(ge->tag(), 1);
g = m->get_geom(ge->tag(), 1);
}
@@ -196,10 +196,10 @@ static bool recover_medge(BDS_Mesh *m, GEdge *ge,
if(itpstart != recoverMapInv.end() && itpend != recoverMapInv.end()){
BDS_Point *pstart = itpstart->second;
BDS_Point *pend = itpend->second;
- if(pass_ == 1)
+ if(pass == 1)
e2r->insert(EdgeToRecover(pstart->iD, pend->iD, ge));
else{
- BDS_Edge *e = m->recover_edge(pstart->iD, pend->iD, e2r, not_recovered);
+ BDS_Edge *e = m->recover_edge(pstart->iD, pend->iD, e2r, notRecovered);
if(e) e->g = g;
// else {
// Msg::Error("Unable to recover an edge %g %g && %g %g (%d/%d)",
@@ -211,7 +211,7 @@ static bool recover_medge(BDS_Mesh *m, GEdge *ge,
}
}
- if(pass_ == 2 && ge->getBeginVertex()){
+ if(pass == 2 && ge->getBeginVertex()){
MVertex *vstart = *(ge->getBeginVertex()->mesh_vertices.begin());
MVertex *vend = *(ge->getEndVertex()->mesh_vertices.begin());
std::map<MVertex*, BDS_Point*>::iterator itpstart = recoverMapInv.find(vstart);
@@ -408,14 +408,14 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
it = edges.begin();
while(it != edges.end()){
if(!(*it)->isMeshDegenerated())
- recover_medge
+ recoverEdge
(m, *it, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
++it;
}
it = emb_edges.begin();
while(it != emb_edges.end()){
if(!(*it)->isMeshDegenerated())
- recover_medge
+ recoverEdge
(m, *it, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
++it;
}
@@ -426,7 +426,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
it = edges.begin();
while(it != edges.end()){
if(!(*it)->isMeshDegenerated()){
- recover_medge
+ recoverEdge
(m, *it, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
}
++it;
@@ -509,7 +509,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
it = emb_edges.begin();
while(it != emb_edges.end()){
if(!(*it)->isMeshDegenerated())
- recover_medge
+ recoverEdge
(m, *it, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
++it;
}
@@ -586,7 +586,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
Msg::Debug("Starting to add internal points");
// start mesh generation
- if(!AlgoDelaunay2D(gf)){
+ if(!algoDelaunay2D(gf)){
gmshRefineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, true,
&recoverMapInv);
gmshOptimizeMeshBDS(gf, *m, 2);
@@ -645,7 +645,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
// the delaunay algo is based directly on internal gmsh structures
// BDS mesh is passed in order not to recompute local coordinates of
// vertices
- if(AlgoDelaunay2D(gf)){
+ if(algoDelaunay2D(gf)){
if(CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL)
gmshBowyerWatsonFrontal(gf);
else
@@ -1131,7 +1131,7 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
}
// start mesh generation
- if(!AlgoDelaunay2D(gf)){
+ if(!algoDelaunay2D(gf)){
gmshRefineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, true);
gmshOptimizeMeshBDS(gf, *m, 2);
gmshRefineMeshBDS (gf, *m, -CTX::instance()->mesh.refineSteps, false);
@@ -1147,7 +1147,7 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
// fill the small gmsh structures
{
- std::set<BDS_Point*,PointLessThan>::iterator itp = m->points.begin();
+ std::set<BDS_Point*, PointLessThan>::iterator itp = m->points.begin();
while (itp != m->points.end()){
BDS_Point *p = *itp;
if(recoverMap.find(p) == recoverMap.end()){
@@ -1194,7 +1194,7 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
outputScalarField(m->triangles, name, 1);
}
- if(AlgoDelaunay2D(gf)){
+ if(algoDelaunay2D(gf)){
if(CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL)
gmshBowyerWatsonFrontal(gf);
else
@@ -1254,7 +1254,7 @@ void meshGFace::operator() (GFace *gf)
if(MeshExtrudedSurface(gf)) return;
const char *algo = "Unknown";
- if(AlgoDelaunay2D(gf))
+ if(algoDelaunay2D(gf))
algo = (CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL) ?
"Frontal" : "Delaunay";
else if(CTX::instance()->mesh.algo2d == ALGO_2D_MESHADAPT_OLD)
@@ -1270,7 +1270,7 @@ void meshGFace::operator() (GFace *gf)
Msg::Debug("Computing edge loops");
Msg::Debug("Generating the mesh");
- if(noseam(gf) || gf->getNativeType() == GEntity::GmshModel ||
+ if(noSeam(gf) || gf->getNativeType() == GEntity::GmshModel ||
gf->edgeLoops.empty()){
gmsh2DMeshGenerator(gf, 0, repairSelfIntersecting1dMesh,
debugSurface >= 0 || debugSurface == -100);
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 78a147473d..c932fdda57 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -214,29 +214,6 @@ void computeMeshSizeFieldAccuracy(GFace *gf, BDS_Mesh &m, double &avg,
}
}
-void computeElementShapes(GFace *gf, BDS_Mesh &m, double &worst, double &avg,
- double &best, int &nT, int &nbGQ)
-{
- std::list<BDS_Face*>::iterator it = m.triangles.begin();
- worst = 1.e22;
- avg = 0.0;
- best = 0.0;
- nT = 0;
- nbGQ = 0;
- while (it!= m.triangles.end()){
- if (!(*it)->deleted){
- double q = qmTriangle(*it, QMTRI_RHO);
- if (q > .9) nbGQ++;
- avg += q;
- worst = std::min(worst, q);
- best = std::max(best, q);
- nT++;
- }
- ++it;
- }
- avg /= nT;
-}
-
// SWAP TESTS i.e. tell if swap should be done
bool edgeSwapTestAngle(BDS_Edge *e, double min_cos)
diff --git a/Mesh/meshGFaceBDS.h b/Mesh/meshGFaceBDS.h
index 352f584030..ab408f73fd 100644
--- a/Mesh/meshGFaceBDS.h
+++ b/Mesh/meshGFaceBDS.h
@@ -18,8 +18,6 @@ class BDS_Mesh;
void computeMeshSizeFieldAccuracy(GFace *gf, BDS_Mesh &m, double &avg,
double &max_e, double &min_e, int &nE, int &GS);
-void computeElementShapes(GFace *gf, BDS_Mesh &m, double &worst, double &avg,
- double &best, int &nT, int &nbGQ);
void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
const bool computeNodalSizeField,
std::map<MVertex*, BDS_Point*> *recoverMapInv=0);
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index b5c961b655..28ea0a71f0 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -723,8 +723,8 @@ void gmshBowyerWatson(GFace *gf)
The point isertion is done on the Voronoi Edges
*/
-static double length_metric(const double p[2], const double q[2],
- const double metric[3])
+static double lengthMetric(const double p[2], const double q[2],
+ const double metric[3])
{
return sqrt ( (p[0]-q[0]) * metric [0] * (p[0]-q[0]) +
2 * (p[0]-q[0]) * metric [1] * (p[1]-q[1]) +
@@ -844,8 +844,8 @@ void gmshBowyerWatsonFrontal(GFace *gf)
2 * dir[1] * dir[0] * metric[1] +
dir[1] * dir[1] * metric[2]);
- const double p = 0.5 * length_metric(P, Q, metric); // / RATIO;
- const double q = length_metric (center, midpoint, metric);
+ const double p = 0.5 * lengthMetric(P, Q, metric); // / RATIO;
+ const double q = lengthMetric(center, midpoint, metric);
const double rhoM1 = 0.5 *
(vSizes[base->getVertex(ip1)->getNum()] +
vSizes[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 31f706fadf..88ce8456c6 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -799,12 +799,12 @@ int edgeCollapsePass(double minLC, GFace *gf, std::set<MTri3*,compareTri3Ptr> &a
extern double angle3Points(MVertex *p1, MVertex *p2, MVertex *p3);
-struct recombine_triangle
+struct RecombineTriangle
{
MElement *t1, *t2;
double angle;
MVertex *n1, *n2, *n3, *n4;
- recombine_triangle(const MEdge &me, MElement *_t1, MElement *_t2)
+ RecombineTriangle(const MEdge &me, MElement *_t1, MElement *_t2)
: t1(_t1), t2(_t2)
{
n1 = me.getVertex(0);
@@ -826,7 +826,7 @@ struct recombine_triangle
angle = std::max(fabs(90. - a3),angle);
angle = std::max(fabs(90. - a4),angle);
}
- bool operator < (const recombine_triangle &other) const
+ bool operator < (const RecombineTriangle &other) const
{
return angle < other.angle;
}
@@ -854,20 +854,20 @@ static void _gmshRecombineIntoQuads(GFace *gf)
e2t_cont adj;
buildEdgeToElement(gf->triangles, adj);
- std::set<recombine_triangle> pairs;
+ std::set<RecombineTriangle> pairs;
for (e2t_cont::iterator it = adj.begin(); it!= adj.end(); ++it){
if (it->second.second &&
it->second.first->getNumVertices() == 3 &&
it->second.second->getNumVertices() == 3 &&
(emb_edgeverts.find(it->first.getVertex(0)) == emb_edgeverts.end() ||
emb_edgeverts.find(it->first.getVertex(1)) == emb_edgeverts.end()))
- pairs.insert(recombine_triangle(it->first,
- it->second.first,
- it->second.second));
+ pairs.insert(RecombineTriangle(it->first,
+ it->second.first,
+ it->second.second));
}
std::set<MElement*> touched;
- std::set<recombine_triangle>::iterator itp = pairs.begin();
+ std::set<RecombineTriangle>::iterator itp = pairs.begin();
while(itp != pairs.end()){
// recombine if difference between max quad angle and right
// angle is smaller than tol
diff --git a/Mesh/meshGFaceQuadrilateralize.cpp b/Mesh/meshGFaceQuadrilateralize.cpp
index 4f759a5474..3b28560a4a 100644
--- a/Mesh/meshGFaceQuadrilateralize.cpp
+++ b/Mesh/meshGFaceQuadrilateralize.cpp
@@ -14,70 +14,71 @@
#include "SVector3.h"
class gmshEdgeFront {
-public:
+ public:
typedef std::set<BDS_Edge *>::const_iterator eiter;
-private:
+ private:
BDS_Mesh *m;
GFace *gf;
- void getFrontEdges (BDS_Point *p, eiter & it1, eiter & it2) const;
- void getFrontEdges (BDS_Point *p, std::vector<eiter> & f) const;
-public:
+ void getFrontEdges(BDS_Point *p, eiter & it1, eiter & it2) const;
+ void getFrontEdges(BDS_Point *p, std::vector<eiter> & f) const;
+ public:
std::set<BDS_Edge *> edges;
std::set<BDS_Edge*> stat[5];
- eiter begin() {return edges.begin();}
- eiter end() {return edges.end();}
- gmshEdgeFront (BDS_Mesh *_m, GFace *_gf) : m(_m), gf(_gf){
- }
+ eiter begin() { return edges.begin(); }
+ eiter end() { return edges.end(); }
+ gmshEdgeFront(BDS_Mesh *_m, GFace *_gf) : m(_m), gf(_gf) {}
// initiate edges in the front i.e.
// take all edges that have one neighbor
// and all edges that have a quad and a triangle
// neighbor
- void initiate ();
+ void initiate();
// compute normal vector to an edge that points
// inside the domain
- SVector3 normal (BDS_Edge *e) const;
+ SVector3 normal(BDS_Edge *e) const;
// compute the state of a front edge
// 0 both vertices have edge angles < 60 deg
// 1 e->p1 have and edge angle > 60 deg
// 2 e->p2 have and edge angle > 60 deg
// 3 both vertices have edge angles > 60 deg
- int computeStatus (BDS_Edge *e) const;
- inline bool inFront (BDS_Edge *e) const{return getStatus(e) != -1;}
- inline int getStatus (BDS_Edge *e) const{
- for (int i=0;i<5;i++){
- if (stat[i].find(e) != stat[i].end())return i;
+ int computeStatus(BDS_Edge *e) const;
+ inline bool inFront(BDS_Edge *e) const{ return getStatus(e) != -1; }
+ inline int getStatus(BDS_Edge *e) const{
+ for(int i = 0; i < 5; i++){
+ if(stat[i].find(e) != stat[i].end()) return i;
}
return -1;
}
// get one edge of the front that is tagged "tag"
// and do whatever has to be done to remove it from
// the front/tag. return false if the front is not empty.
- bool emptyFront (int tag);
+ bool emptyFront(int tag);
// update the status of the edge
- void updateStatus (BDS_Edge *e);
+ void updateStatus(BDS_Edge *e);
// form a quad now
- bool formQuad (BDS_Edge *e, BDS_Edge *left, BDS_Edge *right);
+ bool formQuad(BDS_Edge *e, BDS_Edge *left, BDS_Edge *right);
// delete the cavity delimitated by 4 edges
- void emptyCavity (BDS_Edge *bottom, BDS_Edge *top, BDS_Edge *left, BDS_Edge *right);
+ void emptyCavity(BDS_Edge *bottom, BDS_Edge *top, BDS_Edge *left, BDS_Edge *right);
// delete an edge from the front
- void deleteFromFront(BDS_Edge *e){
+ void deleteFromFront(BDS_Edge *e)
+ {
edges.erase(e);
- for (int i=0;i<5;i++){
+ for(int i = 0; i < 5; i++){
std::set<BDS_Edge*>::iterator it = stat[i].find(e);
- if (it !=stat[i].end()){
+ if(it != stat[i].end()){
stat[i].erase(it);
return;
}
}
}
// taking a point from the front, find the optimal edge
- BDS_Edge *findOptimalEdge(BDS_Point *p,BDS_Point *avoid);
+ BDS_Edge *findOptimalEdge(BDS_Point *p, BDS_Point *avoid);
};
-void gmshEdgeFront::updateStatus (BDS_Edge *e){
- for (int i=0;i<5;i++){
+void gmshEdgeFront::updateStatus(BDS_Edge *e)
+{
+ for(int i = 0; i < 5; i++){
std::set<BDS_Edge*>::iterator it = stat[i].find(e);
- if (it !=stat[i].end()){
+ if(it !=stat[i].end()){
stat[i].erase(it);
stat[computeStatus(e)].insert(e);
return;
@@ -86,61 +87,63 @@ void gmshEdgeFront::updateStatus (BDS_Edge *e){
Msg::Error("Something wrong in updateStatus");
}
-SVector3 norm_edge (BDS_Point *p1, BDS_Point *p2){
- SVector3 n (p2->X-p1->X,p2->Y-p1->Y,p2->Z-p1->Z);
+SVector3 norm_edge(BDS_Point *p1, BDS_Point *p2)
+{
+ SVector3 n(p2->X-p1->X,p2->Y-p1->Y,p2->Z-p1->Z);
n.normalize();
return n;
}
-void recur_empty_cavity (BDS_Face *f,
- BDS_Edge *be[4],
- BDS_Point *bv[4],
- std::set<BDS_Face*> & faces,
- std::set<BDS_Edge*> & edges,
- std::set<BDS_Point*> & vertices){
- if (faces.find(f) != faces.end())return;
+void recur_empty_cavity(BDS_Face *f,
+ BDS_Edge *be[4],
+ BDS_Point *bv[4],
+ std::set<BDS_Face*> & faces,
+ std::set<BDS_Edge*> & edges,
+ std::set<BDS_Point*> & vertices)
+{
+ if(faces.find(f) != faces.end()) return;
faces.insert(f);
BDS_Edge *ee[3] = {f->e1,f->e2,f->e3};
- for (int i=0;i<3;i++){
+ for(int i=0;i<3;i++){
BDS_Edge *e = ee[i];
- if (e != be[0] && e != be[1] && e != be[2] && e != be[3])
- {
- edges.insert(e);
- BDS_Face *of = e->otherFace(f);
- recur_empty_cavity (of,be,bv,faces,edges,vertices);
- }
+ if(e != be[0] && e != be[1] && e != be[2] && e != be[3]){
+ edges.insert(e);
+ BDS_Face *of = e->otherFace(f);
+ recur_empty_cavity(of,be,bv,faces,edges,vertices);
+ }
}
}
-
-void gmshEdgeFront::emptyCavity (BDS_Edge *bottom, BDS_Edge *top, BDS_Edge *left, BDS_Edge *right){
+void gmshEdgeFront::emptyCavity(BDS_Edge *bottom, BDS_Edge *top, BDS_Edge *left,
+ BDS_Edge *right)
+{
// not optimal for now, will be improved further away
- BDS_Face *f=0;
- if (bottom->faces(0) && bottom->faces(0)->numEdges() == 3)f=bottom->faces(0);
- else if (bottom->faces(1))f = bottom->faces(1);
+ BDS_Face *f = 0;
+ if(bottom->faces(0) && bottom->faces(0)->numEdges() == 3) f=bottom->faces(0);
+ else if(bottom->faces(1))f = bottom->faces(1);
- std::set<BDS_Face*> m_faces;
- std::set<BDS_Edge*> m_edges;
+ std::set<BDS_Face*> m_faces;
+ std::set<BDS_Edge*> m_edges;
std::set<BDS_Point*> m_vertices;
BDS_Edge *be[4] = {bottom,top,left,right};
BDS_Point *bv[4] = {bottom->commonvertex(left),
- left->commonvertex(top),
- top->commonvertex(right),
- right->commonvertex(bottom)};
+ left->commonvertex(top),
+ top->commonvertex(right),
+ right->commonvertex(bottom)};
-
- recur_empty_cavity (f,be,bv,m_faces,m_edges,m_vertices);
+ recur_empty_cavity(f, be, bv, m_faces, m_edges, m_vertices);
// printf("%d edges %d faces\n",m_edges.size(),m_faces.size());
- for (std::set<BDS_Face*> :: iterator it = m_faces.begin() ; it != m_faces.end() ; ++it)m->del_face(*it);
- for (std::set<BDS_Edge*> :: iterator it = m_edges.begin() ; it != m_edges.end() ; ++it)m->del_edge(*it);
+ for(std::set<BDS_Face*>::iterator it = m_faces.begin();
+ it != m_faces.end(); ++it) m->del_face(*it);
+ for(std::set<BDS_Edge*>::iterator it = m_edges.begin();
+ it != m_edges.end(); ++it) m->del_edge(*it);
}
-SVector3 gmshEdgeFront::normal (BDS_Edge*e) const{
+SVector3 gmshEdgeFront::normal(BDS_Edge*e) const
+{
BDS_Face *t = e->faces(0);
- if (e->numfaces() == 2 && e->faces(1)->numEdges() == 3)t=e->faces(1);
-
- // printf("%d faces e->faces(0)->numEdges() = %d\n",e->numfaces(),e->faces(0)->numEdges());
+ if(e->numfaces() == 2 && e->faces(1)->numEdges() == 3)t=e->faces(1);
/*
points p1, p2 and p3
@@ -153,124 +156,124 @@ SVector3 gmshEdgeFront::normal (BDS_Edge*e) const{
BDS_Point *p1 = e->p1;
BDS_Point *p2 = e->p2;
BDS_Point *p3;
- if (e == t->e1)
+ if(e == t->e1)
p3 = t->e2->commonvertex(t->e3);
- else if (e == t->e2)
+ else if(e == t->e2)
p3 = t->e1->commonvertex(t->e3);
- else if (e == t->e3)
+ else if(e == t->e3)
p3 = t->e2->commonvertex(t->e1);
else {
Msg::Error("Could not compute fron normal");
return SVector3();
}
- SVector3 t1 (p2->X-p1->X,p2->Y-p1->Y,p2->Z-p1->Z);
- SVector3 t2 (p3->X-p1->X,p3->Y-p1->Y,p3->Z-p1->Z);
+ SVector3 t1(p2->X-p1->X,p2->Y-p1->Y,p2->Z-p1->Z);
+ SVector3 t2(p3->X-p1->X,p3->Y-p1->Y,p3->Z-p1->Z);
SVector3 t3 = crossprod(t1,t2);
double m[3][3] = {{t1.x(),t2.x(),t3.x()},
{t1.y(),t2.y(),t3.y()},
{t1.z(),t2.z(),t3.z()}};
double im[3][3];
inv3x3(m,im);
- SVector3 n (im[1][0],im[1][1],im[1][2]);
+ SVector3 n(im[1][0],im[1][1],im[1][2]);
n.normalize();
return n;
}
-void gmshEdgeFront::getFrontEdges (BDS_Point *p, std::vector<eiter> & fe) const{
- for (std::list<BDS_Edge*>::iterator itp = p->edges.begin(); itp != p->edges.end() ; ++ itp){
+void gmshEdgeFront::getFrontEdges(BDS_Point *p, std::vector<eiter> & fe) const
+{
+ for(std::list<BDS_Edge*>::iterator itp = p->edges.begin();
+ itp != p->edges.end(); ++ itp){
eiter it = edges.find(*itp);
- if ( it != edges.end())
+ if(it != edges.end())
fe.push_back(it);
}
}
-
-void gmshEdgeFront::getFrontEdges (BDS_Point *p, eiter & it1, eiter & it2) const{
+void gmshEdgeFront::getFrontEdges(BDS_Point *p, eiter & it1, eiter & it2) const
+{
int count = 0;
- for (std::list<BDS_Edge*>::iterator itp = p->edges.begin(); itp != p->edges.end() ; ++ itp){
- if (count == 0){
+ for(std::list<BDS_Edge*>::iterator itp = p->edges.begin();
+ itp != p->edges.end(); ++ itp){
+ if(count == 0){
it1 = edges.find(*itp);
- if ( it1 != edges.end()) count++;
+ if(it1 != edges.end()) count++;
}
- else if (count == 1){
+ else if(count == 1){
it2 = edges.find(*itp);
- if ( it2 != edges.end()) return;
+ if(it2 != edges.end()) return;
}
}
Msg::Error("point %d is in the front but has only %d edges\n",p->iD,count);
}
-void gmshEdgeFront::initiate () {
+void gmshEdgeFront::initiate()
+{
edges.clear();
- for (int i=0;i<5;i++)stat[i].clear();
+ for(int i = 0; i < 5; i++) stat[i].clear();
std::list<BDS_Edge*>::iterator it = m->edges.begin();
- while (it!= m->edges.end()){
- if (((*it)->numfaces() == 1 && (*it)->faces(0)->e4 == 0) ||
+ while(it != m->edges.end()){
+ if(((*it)->numfaces() == 1 && (*it)->faces(0)->e4 == 0) ||
((*it)->numfaces() == 2 && (*it)->numTriangles() == 1)) {
edges.insert(*it);
}
++it;
}
- for (eiter it2 = begin(); it2 !=end() ; ++it2){
+ for(eiter it2 = begin(); it2 !=end(); ++it2){
int status = computeStatus(*it2);
stat[status].insert(*it2);
}
}
-
-double angle3Points ( BDS_Point *p1, BDS_Point *p2, BDS_Point *p3 ){
- SVector3 a(p1->X-p2->X,p1->Y-p2->Y,p1->Z-p2->Z);
- SVector3 b(p3->X-p2->X,p3->Y-p2->Y,p3->Z-p2->Z);
+double angle3Points(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
+{
+ SVector3 a(p1->X - p2->X, p1->Y - p2->Y, p1->Z - p2->Z);
+ SVector3 b(p3->X - p2->X, p3->Y - p2->Y, p3->Z - p2->Z);
SVector3 c = crossprod(a,b);
double sinA = c.norm();
double cosA = dot(a,b);
// printf("%d %d %d -> %g %g\n",p1->iD,p2->iD,p3->iD,cosA,sinA);
- return atan2 (sinA,cosA);
+ return atan2(sinA, cosA);
}
-int gmshEdgeFront::computeStatus (BDS_Edge *e) const {
+int gmshEdgeFront::computeStatus(BDS_Edge *e) const
+{
eiter it11, it12, it21,it22;
- getFrontEdges (e->p1, it11, it12);
- getFrontEdges (e->p2, it21, it22);
+ getFrontEdges(e->p1, it11, it12);
+ getFrontEdges(e->p2, it21, it22);
BDS_Edge *e1 = (*it11) == e ? *it12 : *it11;
BDS_Edge *e2 = (*it21) == e ? *it22 : *it21;
- double angle1 = angle3Points ( (*it11)->othervertex(e->p1), e->p1, (*it12)->othervertex(e->p1));
- double angle2 = angle3Points ( (*it21)->othervertex(e->p2), e->p2, (*it22)->othervertex(e->p2));
+ double angle1 = angle3Points((*it11)->othervertex(e->p1), e->p1,
+ (*it12)->othervertex(e->p1));
+ double angle2 = angle3Points((*it21)->othervertex(e->p2), e->p2,
+ (*it22)->othervertex(e->p2));
SVector3 n1 = normal(e);
- SVector3 n2 = norm_edge (e->p1,e1->othervertex(e->p1));
- SVector3 n3 = norm_edge (e->p2,e2->othervertex(e->p2));
- if (dot(n1,n2) < 0)angle1 = M_PI;
- if (dot(n1,n3) < 0)angle2 = M_PI;
-
- // printf("edge %d %d angles %g (%d %d) %g (%d %d)\n",e->p1->iD,e->p2->iD,angle1*180/M_PI,e1->p1->iD,e1->p2->iD,angle2*180/M_PI,
- // e2->p1->iD,e2->p2->iD);
+ SVector3 n2 = norm_edge(e->p1, e1->othervertex(e->p1));
+ SVector3 n3 = norm_edge(e->p2, e2->othervertex(e->p2));
+ if(dot(n1,n2) < 0)angle1 = M_PI;
+ if(dot(n1,n3) < 0)angle2 = M_PI;
- const double angleLimit = 3*M_PI/4.;
+ const double angleLimit = 3 * M_PI/4.;
- if (angle1 > angleLimit && angle2 > angleLimit)return 0;
- if (angle1 > angleLimit && angle2 < angleLimit)return 1;
- if (angle1 < angleLimit && angle2 > angleLimit)return 2;
+ if(angle1 > angleLimit && angle2 > angleLimit) return 0;
+ if(angle1 > angleLimit && angle2 < angleLimit) return 1;
+ if(angle1 < angleLimit && angle2 > angleLimit) return 2;
return 3;
}
+bool gmshEdgeFront::formQuad(BDS_Edge *e, BDS_Edge *left, BDS_Edge *right)
+{
-bool gmshEdgeFront::formQuad (BDS_Edge *e,
- BDS_Edge *left,
- BDS_Edge *right){
-
- printf("e (%d,%d), l(%d,%d), r(%d,%d)\n",
+ printf("e (%d,%d), l(%d,%d), r(%d,%d)\n",
e->p1->iD,e->p2->iD,
left->p1->iD,left->p2->iD,
right->p1->iD,right->p2->iD);
-
// outputScalarField(m->triangles, "deb_before.pos", 0);
-
std::vector<BDS_Point*> toUpdate;
BDS_Point *pleft = left->othervertex(e->p1);
@@ -285,18 +288,18 @@ bool gmshEdgeFront::formQuad (BDS_Edge *e,
both left and right parts of the front
*/
- if (!top) {
+ if(!top) {
// printf("recover\n");
- // top = m->recover_edge_fast (pleft,pright);
- // if (!top)
- top = m->recover_edge (pleft->iD,pright->iD);
+ // top = m->recover_edge_fast(pleft,pright);
+ // if(!top)
+ top = m->recover_edge(pleft->iD, pright->iD);
// printf("recover done %p\n",top);
- if (!top)return false;
+ if(!top) return false;
}
// delete internal triangles
- emptyCavity (e,top,left,right);
+ emptyCavity(e, top, left, right);
// add the quad
- m->add_quadrangle(e,left,top,right);
+ m->add_quadrangle(e, left, top, right);
// update the edge status
// bottom edge leaves the front
deleteFromFront(e);
@@ -317,15 +320,15 @@ bool gmshEdgeFront::formQuad (BDS_Edge *e,
// if the left edge was not in the front, then it becomes
// part of it.
// printf("coucou1\n");
- if (inFront(left))deleteFromFront(left);
+ if(inFront(left)) deleteFromFront(left);
else edges.insert(left);
// printf("coucou2\n");
- if (inFront(right))deleteFromFront(right);
+ if(inFront(right)) deleteFromFront(right);
else edges.insert(right);
// printf("coucou3\n");
- if (inFront(top))deleteFromFront(top);
+ if(inFront(top)) deleteFromFront(top);
else edges.insert(top);
toUpdate.push_back(e->p1);
@@ -333,17 +336,18 @@ bool gmshEdgeFront::formQuad (BDS_Edge *e,
toUpdate.push_back(pleft);
toUpdate.push_back(pright);
- for (unsigned int i=0;i<toUpdate.size();i++){
+ for(unsigned int i = 0; i < toUpdate.size(); i++){
toUpdate[i]->config_modified = true;
//bool done =
- m->smooth_point_parametric(toUpdate[i], gf);
+ m->smooth_point_parametric(toUpdate[i], gf);
// printf("smooth done %d (g %d)\n",done,toUpdate[i]->g->classif_degree);
}
- for (unsigned int i=0;i<toUpdate.size();i++){
+ for(unsigned int i = 0; i < toUpdate.size(); i++){
BDS_Point *p = toUpdate[i];
- for (std::list<BDS_Edge*>::iterator itp = p->edges.begin(); itp != p->edges.end() ; ++ itp){
- if (inFront(*itp)){
+ for(std::list<BDS_Edge*>::iterator itp = p->edges.begin();
+ itp != p->edges.end(); ++ itp){
+ if(inFront(*itp)){
updateStatus(*itp);
}
}
@@ -351,67 +355,67 @@ bool gmshEdgeFront::formQuad (BDS_Edge *e,
return true;
}
-BDS_Edge *gmshEdgeFront::findOptimalEdge(BDS_Point *p, BDS_Point *avoid){
-
- eiter it1,it2;
- getFrontEdges (p, it1, it2);
+BDS_Edge *gmshEdgeFront::findOptimalEdge(BDS_Point *p, BDS_Point *avoid)
+{
+ eiter it1, it2;
+ getFrontEdges(p, it1, it2);
// compute bissector of front edges, this is the optimal direction
SVector3 n1 = normal(*it1);
SVector3 n2 = normal(*it2);
- SVector3 n = n1+n2;
+ SVector3 n = n1 + n2;
n.normalize();
// printf("POINT %g %g %g N %g %g %g\n",p->X,p->Y,p->Z,n.x(),n.y(),n.z());
- double lowerBound = cos (M_PI/6.0);
+ double lowerBound = cos(M_PI / 6.0);
BDS_Edge *found = 0;
- for (std::list<BDS_Edge*>::iterator itp = p->edges.begin(); itp != p->edges.end() ; ++ itp){
+ for(std::list<BDS_Edge*>::iterator itp = p->edges.begin();
+ itp != p->edges.end(); ++ itp){
// the edge is not in the front and is not bounded by quads only
- if (*it1 != *itp && *it2 != *itp && (*itp)->numTriangles()){
+ if(*it1 != *itp && *it2 != *itp && (*itp)->numTriangles()){
BDS_Point *q = (*itp)->othervertex(p);
- SVector3 d (q->X-p->X,q->Y-p->Y,q->Z-p->Z);
+ SVector3 d(q->X - p->X, q->Y - p->Y, q->Z - p->Z);
d.normalize();
double COS = dot(n,d);
- if (COS > lowerBound && q != avoid){
+ if(COS > lowerBound && q != avoid){
lowerBound = COS;
found = *itp;
}
}
}
- if (found) return found;
+ if(found) return found;
else{
-
std::list<BDS_Face*> ts;
double x[2];
- const double L = 0.5*sqrt(3.0)*((*it2)->length() * (*it1)->length()) ;
+ const double L = 0.5*sqrt(3.0)*((*it2)->length() * (*it1)->length());
p->getTriangles(ts);
std::list<BDS_Face*>::iterator it = ts.begin();
std::list<BDS_Face*>::iterator ite = ts.end();
while(it != ite) {
BDS_Face *t = *it;
- if (!t->e4){
+ if(!t->e4){
BDS_Edge *e = t->oppositeEdge(p);
- if (e->numfaces() == 2){
+ if(e->numfaces() == 2){
BDS_Face *f = e->otherFace(t);
- if (!f->e4){
+ if(!f->e4){
BDS_Point *target = f->oppositeVertex(e);
// ONLY WORKS IN 2D for now !!!!!!!!!!!!!!!!!!!
- Intersect_Edges_2d ( e->p1->X,e->p1->Y,
- e->p2->X,e->p2->Y,
- p->X,p->Y,
- p->X+n.x(),p->Y + n.y(),x);
- if ( x[0] >= 0 && x[0] <= 1){
- SVector3 d (target->X-p->X,target->Y-p->Y,target->Z-p->Z);
+ Intersect_Edges_2d(e->p1->X,e->p1->Y,
+ e->p2->X,e->p2->Y,
+ p->X,p->Y,
+ p->X+n.x(),p->Y + n.y(),x);
+ if(x[0] >= 0 && x[0] <= 1){
+ SVector3 d(target->X-p->X,target->Y-p->Y,target->Z-p->Z);
d.normalize();
double COS = dot(n,d);
- double L2 = sqrt ((target->X - p->X) *(target->X - p->X) +
- (target->X - p->Y) *(target->X - p->Y) +
- (target->X - p->Z) *(target->X - p->Z) );
+ double L2 = sqrt((target->X - p->X) *(target->X - p->X) +
+ (target->X - p->Y) *(target->X - p->Y) +
+ (target->X - p->Z) *(target->X - p->Z) );
// swapping the edge alllow to find an edgge that has the right direction and
// right size
- if (COS > cos (M_PI/6.0) && L2 < L){
- m->swap_edge( e, BDS_SwapEdgeTestQuality(false,false));
+ if(COS > cos(M_PI/6.0) && L2 < L){
+ m->swap_edge(e, BDS_SwapEdgeTestQuality(false,false));
BDS_Edge *newE = m->find_edge(p,target);
// printf("swapping -> %p\n",newE);
return newE;
@@ -440,18 +444,17 @@ BDS_Edge *gmshEdgeFront::findOptimalEdge(BDS_Point *p, BDS_Point *avoid){
return 0;
}
-
-
-bool gmshEdgeFront::emptyFront (int tag){
+bool gmshEdgeFront::emptyFront(int tag)
+{
// front edges tagged "tag" is empty
- if (stat[tag].size() == 0)return true;
+ if(stat[tag].size() == 0) return true;
BDS_Edge *e = *(stat[tag].begin());
- BDS_Edge *left,*right=0;
- eiter it1,it2;
+ BDS_Edge *left, *right = 0;
+ eiter it1, it2;
- std::vector<eiter> fe1,fe2;
+ std::vector<eiter> fe1, fe2;
- printf("front edges %d %d tag %d\n",e->p1->iD,e->p2->iD,tag);
+ printf("front edges %d %d tag %d\n", e->p1->iD, e->p2->iD, tag);
switch(tag){
// both left and right neighboring edges are
@@ -459,11 +462,11 @@ bool gmshEdgeFront::emptyFront (int tag){
// form the quad
case 3:
{
- getFrontEdges (e->p1, it1, it2);
- if (*it1 == e)left=*it2;
+ getFrontEdges(e->p1, it1, it2);
+ if(*it1 == e) left=*it2;
else left = *it1;
- getFrontEdges (e->p2, it1, it2);
- if (*it1 == e)right=*it2;
+ getFrontEdges(e->p2, it1, it2);
+ if(*it1 == e) right = *it2;
else right = *it1;
// printf("case 3\n");
}
@@ -472,32 +475,32 @@ bool gmshEdgeFront::emptyFront (int tag){
// we therefore find a new front edge in the
// mesh that allows to move to tag 3
case 2:
- getFrontEdges (e->p1, it1, it2);
- if (*it1 == e)left=*it2;
+ getFrontEdges(e->p1, it1, it2);
+ if(*it1 == e) left = *it2;
else left = *it1;
// printf("case 2 left edge %p\n",left);
- right = findOptimalEdge(e->p2,left->othervertex(e->p1));
- if (right)getFrontEdges (right->othervertex(e->p2),fe2);
+ right = findOptimalEdge(e->p2, left->othervertex(e->p1));
+ if(right) getFrontEdges(right->othervertex(e->p2), fe2);
break;
// left edge is angled with current edge
// we therefore find a new front edge in the
// mesh that allows to move to tag 3
case 1:
- getFrontEdges (e->p2, it1, it2);
- if (*it2 == e)right=*it1;
+ getFrontEdges(e->p2, it1, it2);
+ if(*it2 == e) right = *it1;
else right = *it2;
// printf("case 1 right edge %p %p\n",e,right);
- left = findOptimalEdge( e->p1,right->othervertex(e->p2));
- if (left)getFrontEdges (left->othervertex(e->p1),fe1);
+ left = findOptimalEdge(e->p1, right->othervertex(e->p2));
+ if(left) getFrontEdges(left->othervertex(e->p1), fe1);
break;
// no neighboring edge is angled with current edge
// we therefore find a new front edge in the
// mesh that allows to move to tag 3
case 0:
- left = findOptimalEdge( e->p1,0);
- if (left)right= findOptimalEdge( e->p2,left->othervertex(e->p1));
- if(right)getFrontEdges (right->othervertex(e->p2),fe2);
- if(left)getFrontEdges (left->othervertex(e->p1),fe1);
+ left = findOptimalEdge(e->p1, 0);
+ if(left) right= findOptimalEdge(e->p2, left->othervertex(e->p1));
+ if(right) getFrontEdges(right->othervertex(e->p2), fe2);
+ if(left) getFrontEdges(left->othervertex(e->p1), fe1);
// printf("strange\n");
break;
default:
@@ -505,7 +508,7 @@ bool gmshEdgeFront::emptyFront (int tag){
return false;
}
- if ( fe1.size() || fe2.size() || !left || !right || !formQuad (e, left, right)){
+ if(fe1.size() || fe2.size() || !left || !right || !formQuad(e, left, right)){
// printf("front cloeses : algo has to be finished\n");
stat[tag].erase(stat[tag].begin());
stat[4].insert(e);
@@ -513,20 +516,21 @@ bool gmshEdgeFront::emptyFront (int tag){
return false;
}
-static int numQuads ( BDS_Mesh *m)
+static int numQuads(BDS_Mesh *m)
{
int N = 0;
- for (std::list<BDS_Face*> ::iterator it = m->triangles.begin() ; it != m->triangles.end() ; ++it){
- if ((*it)->e4) N++;
+ for(std::list<BDS_Face*>::iterator it = m->triangles.begin();
+ it != m->triangles.end(); ++it){
+ if((*it)->e4) N++;
}
return N;
}
-int gmshQMorph (GFace *gf)
+int gmshQMorph(GFace *gf)
{
// assert first that there exist a triangulation of
// the face
- if (!gf->triangles.size()){
+ if(!gf->triangles.size()){
Msg::Warning("Cannot Quadrilaterize a face that has not been triangulated");
return -1;
}
@@ -536,24 +540,24 @@ int gmshQMorph (GFace *gf)
BDS_Mesh *pm = gmsh2BDS(l);
// create the front
- gmshEdgeFront front (pm,gf);
+ gmshEdgeFront front(pm,gf);
front.initiate();
int ITER = 1;
// empty the front for front edges tagged 3, 2, 1 then 0
int _numQuads = 0;
- while (1){
- if (front.emptyFront ( 3 )){
- if (front.emptyFront ( 2 )){
- if (front.emptyFront ( 1 )){
- if (front.emptyFront ( 0 )){
+ while(1){
+ if(front.emptyFront(3)){
+ if(front.emptyFront(2)){
+ if(front.emptyFront(1)){
+ if(front.emptyFront(0)){
int ns;
smoothVertexPass(gf,*pm,ns,false);
printf("nex row iter %6d->>>\n",ITER);
front.initiate();
- int _numQuadsNew = numQuads (pm);
- if (front.edges.size() == 0 || _numQuads == _numQuadsNew)
+ int _numQuadsNew = numQuads(pm);
+ if(front.edges.size() == 0 || _numQuads == _numQuadsNew)
break;
_numQuads = _numQuadsNew;
}
@@ -561,15 +565,14 @@ int gmshQMorph (GFace *gf)
}
}
ITER++;
- if (1 || ITER%100 == 0){
+ if(1 || ITER%100 == 0){
char name[256];
sprintf(name,"qmorph-face%d-iter%d.pos",gf->tag(),ITER);
outputScalarField(pm->triangles, name, 0);
}
- // if (ITER == 1123)break;
+ // if(ITER == 1123)break;
}
// delete the BDS
delete pm;
return 1;
}
-
--
GitLab