diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 85744518732ae067168057fad05ada7e3e7633e0..3847d4899c38fe307864b56773ee7317484de3ed 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.97 2008-02-05 14:40:30 remacle Exp $
+// $Id: BDS.cpp,v 1.98 2008-02-16 20:42:40 geuzaine Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -28,37 +28,39 @@
#include "GFace.h"
#include "qualityMeasures.h"
-bool test_move_point_parametric_triangle (BDS_Point * p, double u, double v, BDS_Face *t);
+bool test_move_point_parametric_triangle(BDS_Point *p, double u, double v, BDS_Face *t);
-void outputScalarField(std::list < BDS_Face * >t, const char *iii, int param, GFace *gf)
+void outputScalarField(std::list<BDS_Face*> t, const char *iii, int param, GFace *gf)
{
FILE *f = fopen(iii, "w");
if(!f) return;
fprintf(f, "View \"scalar\" {\n");
- std::list < BDS_Face * >::iterator tit = t.begin();
- std::list < BDS_Face * >::iterator tite = t.end();
+ std::list<BDS_Face*>::iterator tit = t.begin();
+ std::list<BDS_Face*>::iterator tite = t.end();
while(tit != tite) {
BDS_Point *pts[4];
(*tit)->getNodes(pts);
if (param)
- fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
pts[0]->u, pts[0]->v, 0.0,
pts[1]->u, pts[1]->v, 0.0,
- pts[2]->u, pts[2]->v, 0.0,(double)pts[0]->iD,(double)pts[1]->iD,(double)pts[2]->iD);
+ pts[2]->u, pts[2]->v, 0.0,
+ pts[0]->iD, pts[1]->iD, pts[2]->iD);
else{
if (!gf)
- fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
pts[0]->X, pts[0]->Y, pts[0]->Z,
pts[1]->X, pts[1]->Y, pts[1]->Z,
- pts[2]->X, pts[2]->Y, pts[2]->Z,(double)pts[0]->iD,(double)pts[1]->iD,(double)pts[2]->iD);
+ pts[2]->X, pts[2]->Y, pts[2]->Z,
+ pts[0]->iD, pts[1]->iD, pts[2]->iD);
else
fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
pts[0]->X, pts[0]->Y, pts[0]->Z,
pts[1]->X, pts[1]->Y, pts[1]->Z,
pts[2]->X, pts[2]->Y, pts[2]->Z,
- gf->curvature (SPoint2(pts[0]->u, pts[0]->v)),
- gf->curvature (SPoint2(pts[1]->u, pts[1]->v)),
- gf->curvature (SPoint2(pts[2]->u, pts[2]->v)));
+ gf->curvature(SPoint2(pts[0]->u, pts[0]->v)),
+ gf->curvature(SPoint2(pts[1]->u, pts[1]->v)),
+ gf->curvature(SPoint2(pts[2]->u, pts[2]->v)));
}
++tit;
}
@@ -110,18 +112,18 @@ double surface_triangle_param(BDS_Point * p1, BDS_Point * p2, BDS_Point * p3)
return (0.5 * c);
}
-void BDS_Point::getTriangles(std::list < BDS_Face * >&t) const
+void BDS_Point::getTriangles(std::list<BDS_Face*> &t) const
{
t.clear();
- std::list < BDS_Edge * >::const_iterator it = edges.begin();
- std::list < BDS_Edge * >::const_iterator ite = edges.end();
+ std::list<BDS_Edge*>::const_iterator it = edges.begin();
+ std::list<BDS_Edge*>::const_iterator ite = edges.end();
while(it != ite) {
int NF = (*it)->numfaces();
for(int i = 0; i < NF; ++i) {
BDS_Face *tt = (*it)->faces(i);
if(tt) {
- std::list < BDS_Face * >::iterator tit = t.begin();
- std::list < BDS_Face * >::iterator tite = t.end();
+ std::list<BDS_Face*>::iterator tit = t.begin();
+ std::list<BDS_Face*>::iterator tite = t.end();
int found = 0;
while(tit != tite) {
if(tt == *tit)
@@ -180,13 +182,13 @@ BDS_Edge *BDS_Mesh::find_edge(BDS_Point *p, int num2)
BDS_Edge *BDS_Mesh::find_edge(BDS_Point *p1, BDS_Point *p2)
{
- return find_edge (p1,p2->iD);
+ return find_edge(p1, p2->iD);
}
BDS_Edge *BDS_Mesh::find_edge(int num1, int num2)
{
BDS_Point *p = find_point(num1);
- return find_edge (p,num2);
+ return find_edge(p, num2);
}
int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
@@ -211,7 +213,6 @@ int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
// }
// return 0;
-
double mat[2][2];
double rhs[2], x[2];
mat[0][0] = (x2 - x1);
@@ -227,7 +228,8 @@ 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, std::set<EdgeToRecover> *not_recovered )
+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);
@@ -236,80 +238,70 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
BDS_Point *p1 = find_point(num1);
BDS_Point *p2 = find_point(num2);
- if (!p1 || !p2) throw;;
+ if (!p1 || !p2) throw;
- Msg(DEBUG," edge %d %d has to be recovered",num1,num2);
+ Msg(DEBUG, " edge %d %d has to be recovered", num1, num2);
int ix = 0;
int ixMax = 300;
- while(1)
- {
- std::vector<BDS_Edge *> intersected;
- std::list<BDS_Edge *>::iterator it = edges.begin();
- while (it!=edges.end())
- {
- e = (*it);
- // if (e->p1->iD >= 0 && e->p2->iD >= 0)Msg(INFO," %d %d %g %g - %g %g",e->p1->iD,e->p2->iD,e->p1->u,e->p1->v,e->p2->u,e->p2->v);
- if (!e->deleted && e->p1 != p1 && e->p1 != p2 && e->p2 != p1 && e->p2 != p2)
-// if (e->p1->iD == 99 || e->p2->iD == 99) {
-// printf("%d %d %d %d\n",e->p1->iD,e->p2->iD,p1->iD,p2->iD);
-// printf("%g %g,%g %g,%g %g,%g %g\n",e->p1->u, e->p1->v,e->p2->u, e->p2->v,p1->u,p1->v,p2->u,p2->v);
-// }
- if(Intersect_Edges_2d(e->p1->u, e->p1->v,
- e->p2->u, e->p2->v,
- p1->u, p1->v,
- p2->u, p2->v))
- {
- // printf("intersect\n");
- if (e2r && e2r->find(EdgeToRecover(e->p1->iD,e->p2->iD,0)) != e2r->end())
- {
- 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(DEBUG2," 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());
-
- // 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;
+ while(1){
+ std::vector<BDS_Edge*> intersected;
+ std::list<BDS_Edge*>::iterator it = edges.begin();
+ while(it != edges.end()){
+ e = (*it);
+ if(!e->deleted && e->p1 != p1 && e->p1 != p2 && e->p2 != p1 && e->p2 != p2)
+ if(Intersect_Edges_2d(e->p1->u, e->p1->v,
+ e->p2->u, e->p2->v,
+ p1->u, p1->v,
+ p2->u, p2->v)){
+ // intersect
+ if(e2r && e2r->find(EdgeToRecover(e->p1->iD, e->p2->iD, 0)) != e2r->end()){
+ 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(DEBUG2, " 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());
+ // 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 (ix > 300){
-// Msg(WARNING," edge %d %d cannot be recovered after %d iterations, trying again",
-// num1,num2,ix);
-// ix = 0;
-// }
-// printf("%d %d\n",intersected.size(),ix);
-
- if (!intersected.size() || ix > 1000)
- {
- BDS_Edge *eee = find_edge (num1, num2);
- if (!eee)
- {
- outputScalarField(triangles, "debugp.pos",1);
- outputScalarField(triangles, "debugr.pos",0);
- Msg(GERROR," edge %d %d cannot be recovered at all, look at debugp.pos and debugr.pos",
- num1,num2);
- return 0;
- }
- return eee;
- }
+// if (ix > 300){
+// Msg(WARNING," edge %d %d cannot be recovered after %d iterations, trying again",
+// num1, num2, ix);
+// ix = 0;
+// }
+// printf("%d %d\n",intersected.size(),ix);
- int ichoice = ix++ % intersected.size();
- bool success = swap_edge ( intersected[ichoice] , BDS_SwapEdgeTestQuality (false,false) );
-// printf("trying to swop %d %d = %d (%d %d)\n",intersected[ichoice]->p1->iD,intersected[ichoice]->p2->iD,success,
-// intersected[ichoice]->deleted,intersected[ichoice]->numfaces());
+ if(!intersected.size() || ix > 1000){
+ BDS_Edge *eee = find_edge (num1, num2);
+ if(!eee){
+ outputScalarField(triangles, "debugp.pos", 1);
+ outputScalarField(triangles, "debugr.pos", 0);
+ Msg(GERROR," edge %d %d cannot be recovered at all, look at debugp.pos "
+ "and debugr.pos", num1, num2);
+ return 0;
+ }
+ return eee;
}
+
+ int ichoice = ix++ % intersected.size();
+ bool success = swap_edge(intersected[ichoice], BDS_SwapEdgeTestQuality(false, false));
+ // printf("trying to swop %d %d = %d (%d %d)\n", intersected[ichoice]->p1->iD,
+ // intersected[ichoice]->p2->iD, success, intersected[ichoice]->deleted,
+ // intersected[ichoice]->numfaces());
+ }
return 0;
}
-BDS_Edge *BDS_Mesh::find_edge(BDS_Point * p1, BDS_Point * p2,
- BDS_Face * t) const
+BDS_Edge *BDS_Mesh::find_edge(BDS_Point *p1, BDS_Point *p2, BDS_Face *t) const
{
BDS_Point P1(p1->iD);
BDS_Point P2(p2->iD);
@@ -323,8 +315,7 @@ BDS_Edge *BDS_Mesh::find_edge(BDS_Point * p1, BDS_Point * p2,
return 0;
}
-BDS_Face *BDS_Mesh::find_triangle(BDS_Edge * e1, BDS_Edge * e2,
- BDS_Edge * e3)
+BDS_Face *BDS_Mesh::find_triangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3)
{
int i;
for(i = 0; i < e1->numfaces(); i++) {
@@ -372,8 +363,7 @@ BDS_Face *BDS_Mesh::find_triangle(BDS_Edge * e1, BDS_Edge * e2,
BDS_Edge *BDS_Mesh::add_edge(int p1, int p2)
{
BDS_Edge *efound = find_edge(p1, p2);
- if(efound)
- return efound;
+ if(efound) return efound;
BDS_Point *pp1 = find_point(p1);
BDS_Point *pp2 = find_point(p2);
@@ -392,21 +382,18 @@ BDS_Face *BDS_Mesh::add_triangle(int p1, int p2, int p3)
return add_triangle(e1, e2, e3);
}
-BDS_Face *BDS_Mesh::add_triangle(BDS_Edge * e1, BDS_Edge * e2,
- BDS_Edge * e3)
+BDS_Face *BDS_Mesh::add_triangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3)
{
-
- //BDS_Face *tfound = find_triangle(e1, e2, e3);
- // if(tfound)
- // return tfound;
+ // BDS_Face *tfound = find_triangle(e1, e2, e3);
+ // if(tfound) return tfound;
BDS_Face *t = new BDS_Face(e1, e2, e3);
triangles.push_back(t);
return t;
}
-BDS_Face *BDS_Mesh::add_quadrangle(BDS_Edge * e1, BDS_Edge * e2,
- BDS_Edge * e3,BDS_Edge * e4 )
+BDS_Face *BDS_Mesh::add_quadrangle(BDS_Edge *e1, BDS_Edge *e2,
+ BDS_Edge *e3, BDS_Edge *e4)
{
BDS_Face *t = new BDS_Face(e1, e2, e3, e4);
triangles.push_back(t);
@@ -418,7 +405,7 @@ void BDS_Mesh::del_face(BDS_Face * t)
t->e1->del(t);
t->e2->del(t);
t->e3->del(t);
- if(t->e4)t->e4->del(t);
+ if(t->e4) t->e4->del(t);
t->deleted = true;
}
@@ -469,8 +456,7 @@ BDS_GeomEntity *BDS_Mesh::get_geom(int p1, int p2)
{
BDS_GeomEntity ge(p1, p2);
std::set < BDS_GeomEntity *, GeomLessThan >::iterator it = geom.find(&ge);
- if(it == geom.end())
- return 0;
+ if(it == geom.end()) return 0;
return (BDS_GeomEntity *) (*it);
}
@@ -478,7 +464,7 @@ void recur_tag(BDS_Face * t, BDS_GeomEntity * g)
{
if(!t->g) {
t->g = g;
- // g->t.push_back(t);
+ // g->t.push_back(t);
if(!t->e1->g && t->e1->numfaces() == 2) {
recur_tag(t->e1->otherFace(t), g);
}
@@ -491,7 +477,6 @@ void recur_tag(BDS_Face * t, BDS_GeomEntity * g)
}
}
-
double PointLessThanLexicographic::t = 0;
double BDS_Vector::t = 0;
@@ -506,7 +491,7 @@ template < class IT > void DESTROOOY(IT beg, IT end)
void BDS_Mesh::cleanup()
{
{
- std::list<BDS_Face*> :: iterator it = triangles.begin();
+ std::list<BDS_Face*>::iterator it = triangles.begin();
while(it != triangles.end()){
if((*it)->deleted){
delete *it;
@@ -517,7 +502,7 @@ void BDS_Mesh::cleanup()
}
}
{
- std::list<BDS_Edge*> :: iterator it = edges.begin();
+ std::list<BDS_Edge*>::iterator it = edges.begin();
while(it != edges.end()){
if((*it)->deleted){
delete *it;
@@ -529,21 +514,20 @@ void BDS_Mesh::cleanup()
}
}
-BDS_Mesh ::~ BDS_Mesh ()
+BDS_Mesh::~BDS_Mesh()
{
- DESTROOOY ( geom.begin(),geom.end());
- DESTROOOY ( points.begin(),points.end());
+ DESTROOOY(geom.begin(), geom.end());
+ DESTROOOY(points.begin(), points.end());
cleanup();
- DESTROOOY ( edges.begin(),edges.end());
- DESTROOOY ( triangles.begin(),triangles.end());
+ DESTROOOY(edges.begin(), edges.end());
+ DESTROOOY(triangles.begin(), triangles.end());
}
-
-bool BDS_Mesh::split_face(BDS_Face * f, BDS_Point *mid)
+bool BDS_Mesh::split_face(BDS_Face *f, BDS_Point *mid)
{
- BDS_Point *p1 = f->e1->commonvertex (f->e2);
- BDS_Point *p2 = f->e3->commonvertex (f->e2);
- BDS_Point *p3 = f->e1->commonvertex (f->e3);
+ BDS_Point *p1 = f->e1->commonvertex(f->e2);
+ BDS_Point *p2 = f->e3->commonvertex(f->e2);
+ BDS_Point *p3 = f->e1->commonvertex(f->e3);
BDS_Edge *p1_mid = new BDS_Edge(p1, mid);
edges.push_back(p1_mid);
BDS_Edge *p2_mid = new BDS_Edge(p2, mid);
@@ -551,9 +535,9 @@ bool BDS_Mesh::split_face(BDS_Face * f, BDS_Point *mid)
BDS_Edge *p3_mid = new BDS_Edge(p3, mid);
edges.push_back(p2_mid);
BDS_Face *t1, *t2, *t3;
- t1 = new BDS_Face(f->e1, p1_mid,p3_mid);
- t2 = new BDS_Face(f->e2, p2_mid,p1_mid);
- t3 = new BDS_Face(f->e3, p3_mid,p2_mid);
+ t1 = new BDS_Face(f->e1, p1_mid, p3_mid);
+ t2 = new BDS_Face(f->e2, p2_mid, p1_mid);
+ t3 = new BDS_Face(f->e3, p3_mid, p2_mid);
t1->g = f->g;
t2->g = f->g;
@@ -597,7 +581,6 @@ bool BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
// p1,op2,mid +
*/
-
BDS_Point *op[2];
BDS_Point *p1 = e->p1;
BDS_Point *p2 = e->p2;
@@ -668,7 +651,6 @@ bool BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
t3->g = g1;
t4->g = g2;
-
p1_mid->g = ge;
mid_p2->g = ge;
op1_mid->g = g1;
@@ -693,51 +675,46 @@ bool BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
// the feasability of the operation. Those conditions have to be
// taken into account before doing the edge swap
-bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1,BDS_Point *_p2,
- BDS_Point *_q1,BDS_Point *_q2) const
+bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1, BDS_Point *_p2,
+ BDS_Point *_q1, BDS_Point *_q2) const
{
if (!testSmallTriangles){
-
double p1 [2] = {_p1->u,_p1->v};
double p2 [2] = {_p2->u,_p2->v};
double op1[2] = {_q1->u,_q1->v};
double op2[2] = {_q2->u,_q2->v};
-
double ori_t1 = gmsh::orient2d(op1, p1, op2);
double ori_t2 = gmsh::orient2d(op1, op2, p2);
- // printf("%d %d %d %d %g %g\n",_p1->iD,_p2->iD,_q1->iD,_q2->iD,ori_t1,ori_t2);
- return ( ori_t1 * ori_t2 > 0 ); // the quadrangle was strictly convex !
+ // printf("%d %d %d %d %g %g\n",_p1->iD,_p2->iD,_q1->iD,_q2->iD,ori_t1,ori_t2);
+ return (ori_t1 * ori_t2 > 0); // the quadrangle was strictly convex !
}
-
-
- double s1 = fabs(surface_triangle_param(_p1,_p2,_q1));
- double s2 = fabs(surface_triangle_param(_p1,_p2,_q2));
- double s3 = fabs(surface_triangle_param(_p1,_q1,_q2));
- double s4 = fabs(surface_triangle_param(_p2,_q1,_q2));
- if (fabs(s1+s2-s3-s4) > 1.e-12 * (s1+s2))return false;
- if (s3 < .02 * (s1+s2) || s4 < .02 * (s1+s2))
- return false;
- return true;
-
-
+
+ double s1 = fabs(surface_triangle_param(_p1, _p2, _q1));
+ double s2 = fabs(surface_triangle_param(_p1, _p2, _q2));
+ double s3 = fabs(surface_triangle_param(_p1, _q1, _q2));
+ double s4 = fabs(surface_triangle_param(_p2, _q1, _q2));
+ if (fabs(s1 + s2 - s3 - s4) > 1.e-12 * (s1 + s2)) return false;
+ if (s3 < .02 * (s1 + s2) || s4 < .02 * (s1 + s2))
+ return false;
+ return true;
}
-bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1,BDS_Point *_p2,BDS_Point *_p3,
- BDS_Point *_q1,BDS_Point *_q2,BDS_Point *_q3,
- BDS_Point *_op1,BDS_Point *_op2,BDS_Point *_op3,
- BDS_Point *_oq1,BDS_Point *_oq2,BDS_Point *_oq3 )const
+bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1, BDS_Point *_p2, BDS_Point *_p3,
+ BDS_Point *_q1, BDS_Point *_q2, BDS_Point *_q3,
+ BDS_Point *_op1, BDS_Point *_op2, BDS_Point *_op3,
+ BDS_Point *_oq1, BDS_Point *_oq2, BDS_Point *_oq3) const
{
if (!testQuality) return true;
- double n[3],q[3],on[3],oq[3];
- normal_triangle(_p1,_p2,_p3,n);
- normal_triangle(_q1,_q2,_q3,q);
- normal_triangle(_op1,_op2,_op3,on);
- normal_triangle(_oq1,_oq2,_oq3,oq);
+ double n[3], q[3], on[3], oq[3];
+ normal_triangle(_p1, _p2, _p3, n);
+ normal_triangle(_q1, _q2, _q3, q);
+ normal_triangle(_op1, _op2, _op3, on);
+ normal_triangle(_oq1, _oq2, _oq3, oq);
- double cosnq; prosca(n,q,&cosnq);
- double cosonq; prosca(on,oq,&cosonq);
+ double cosnq; prosca(n, q, &cosnq);
+ double cosonq; prosca(on, oq, &cosonq);
double qa1 = qmTriangle(_p1, _p2, _p3,QMTRI_RHO);
double qa2 = qmTriangle(_q1, _q2, _q3,QMTRI_RHO);
@@ -748,22 +725,21 @@ bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1,BDS_Point *_p2,BDS_Poi
double mina = std::min(qa1,qa2);
double minb = std::min(qb1,qb2);
-// if (cosnq < .3 && cosonq > .5 && minb > .1)
-// printf("mina = %g minb = %g cos %g %g\n",mina,minb,cosnq,cosonq);
+ // if (cosnq < .3 && cosonq > .5 && minb > .1)
+ // printf("mina = %g minb = %g cos %g %g\n",mina,minb,cosnq,cosonq);
- if (cosnq < .3 && cosonq > .5 && minb > .1) return true;
+ if (cosnq < .3 && cosonq > .5 && minb > .1) return true;
- if (minb > mina) return true;
+ if (minb > mina) return true;
// if (mina > minb && cosnq <= cosonq)return true;
return false;
-
}
void swap_config(BDS_Edge * e,
- BDS_Point **p11,BDS_Point **p12,BDS_Point **p13,
- BDS_Point **p21,BDS_Point **p22,BDS_Point **p23,
- BDS_Point **p31,BDS_Point **p32,BDS_Point **p33,
- BDS_Point **p41,BDS_Point **p42,BDS_Point **p43)
+ BDS_Point **p11, BDS_Point **p12, BDS_Point **p13,
+ BDS_Point **p21, BDS_Point **p22, BDS_Point **p23,
+ BDS_Point **p31, BDS_Point **p32, BDS_Point **p33,
+ BDS_Point **p41, BDS_Point **p42, BDS_Point **p43)
{
BDS_Point *op[2];
BDS_Point *p1 = e->p1;
@@ -819,7 +795,7 @@ void swap_config(BDS_Edge * e,
}
}
-bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
+bool BDS_Mesh::swap_edge(BDS_Edge *e, const BDS_SwapEdgeTest &theTest)
{
/*
@@ -848,7 +824,6 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
if(e->g && e->g->classif_degree == 1)
return false;
-
BDS_Point *op[2];
BDS_Point *p1 = e->p1;
BDS_Point *p2 = e->p2;
@@ -859,7 +834,6 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
BDS_Point *pts1[4];
e->faces(0)->getNodes(pts1);
-
// compute the orientation of the face
// with respect to the edge
int orientation = 0;
@@ -874,24 +848,23 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
}
if(orientation == 1) {
- if (!theTest ( p1, p2, op[0],
- p2,p1,op[1],
- p1,op[1],op[0],
- op[1],p2,op[0]))
- return false;
+ if (!theTest(p1, p2, op[0],
+ p2, p1, op[1],
+ p1, op[1], op[0],
+ op[1], p2, op[0]))
+ return false;
}
else{
- if (!theTest ( p2, p1, op[0],
- p1,p2,op[1],
- p1,op[0],op[1],
- op[1],op[0],p2))
+ if (!theTest(p2, p1, op[0],
+ p1, p2, op[1],
+ p1, op[0], op[1],
+ op[1], op[0], p2))
return false;
}
-
- if (!theTest ( p1, p2, op[0],op[1]))
+
+ if (!theTest(p1, p2, op[0], op[1]))
return false;
-
BDS_Edge *p1_op1 = find_edge(p1, op[0], e->faces(0));
BDS_Edge *op1_p2 = find_edge(op[0], p2, e->faces(0));
BDS_Edge *p1_op2 = find_edge(p1, op[1], e->faces(1));
@@ -908,7 +881,6 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
}
del_edge(e);
-
BDS_Edge *op1_op2 = new BDS_Edge(op[0], op[1]);
edges.push_back(op1_op2);
@@ -938,10 +910,10 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
return true;
}
-int BDS_Edge :: numTriangles() const
+int BDS_Edge::numTriangles() const
{
int NT = 0;
- for (unsigned int i=0;i<_faces.size();i++)
+ for (unsigned int i = 0; i < _faces.size(); i++)
if (faces(i)->numEdges() == 3) NT++ ;
return NT;
}
@@ -951,7 +923,6 @@ int BDS_Edge :: numTriangles() const
// taken into account before doing the edge swap
bool BDS_Mesh::recombine_edge(BDS_Edge * e)
{
-
/*
p1
/ | \
@@ -1012,10 +983,8 @@ bool BDS_Mesh::recombine_edge(BDS_Edge * e)
return true;
}
-
bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
{
-
if(e->numfaces() != 2)
return false;
if(p->g && p->g->classif_degree == 0)
@@ -1049,7 +1018,7 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
while(it != ite) {
BDS_Face *t = *it;
if(t->e1 != e && t->e2 != e && t->e3 != e) {
- if (!test_move_point_parametric_triangle ( p, o->u, o->v, t))
+ if (!test_move_point_parametric_triangle(p, o->u, o->v, t))
return false;
gs[nt] = t->g;
BDS_Point *pts[4];
@@ -1058,13 +1027,13 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
pt[1][nt] = (pts[1] == p) ? o : pts[1];
pt[2][nt] = (pts[2] == p) ? o : pts[2];
-// double qnew = qmTriangle(pt[0][nt],pt[1][nt],pt[2][nt],QMTRI_RHO);
-// double qold = qmTriangle(pts[0],pts[1],pts[2],QMTRI_RHO);
-// if ( qold > 1.e-4 && qnew < 1.e-4)return false;
+// double qnew = qmTriangle(pt[0][nt], pt[1][nt], pt[2][nt], QMTRI_RHO);
+// double qold = qmTriangle(pts[0], pts[1], pts[2], QMTRI_RHO);
+// if(qold > 1.e-4 && qnew < 1.e-4) return false;
nt++;
-// pt[0][nt] = (pts[0] == p) ? o->iD : pts[0]->iD;
-// pt[1][nt] = (pts[1] == p) ? o->iD : pts[1]->iD;
-// pt[2][nt++] = (pts[2] == p) ? o->iD : pts[2]->iD;
+// pt[0][nt] = (pts[0] == p) ? o->iD : pts[0]->iD;
+// pt[1][nt] = (pts[1] == p) ? o->iD : pts[1]->iD;
+// pt[2][nt++] = (pts[2] == p) ? o->iD : pts[2]->iD;
}
++it;
}
@@ -1085,11 +1054,11 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
int kk = 0;
{
- std::list < BDS_Edge * >edges(p->edges);
- std::list < BDS_Edge * >::iterator eit = edges.begin();
- std::list < BDS_Edge * >::iterator eite = edges.end();
+ std::list<BDS_Edge*> edges(p->edges);
+ std::list<BDS_Edge*>::iterator eit = edges.begin();
+ std::list<BDS_Edge*>::iterator eite = edges.end();
while(eit != eite) {
- (*eit)->p1->config_modified = (*eit)->p2->config_modified = true;
+ (*eit)->p1->config_modified = (*eit)->p2->config_modified = true;
ept[0][kk] = ((*eit)->p1 == p) ? o->iD : (*eit)->p1->iD;
ept[1][kk] = ((*eit)->p2 == p) ? o->iD : (*eit)->p2->iD;
egs[kk++] = (*eit)->g;
@@ -1113,24 +1082,22 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
e->g = egs[i];
}
-
return true;
}
-
-// use robust predicates for not allowing to revert a triangle
-// by moving one of its vertices
-bool test_move_point_parametric_triangle (BDS_Point * p, double u, double v, BDS_Face *t)
+// use robust predicates for not allowing to revert a triangle by
+// moving one of its vertices
+bool test_move_point_parametric_triangle(BDS_Point *p, double u, double v, BDS_Face *t)
{
BDS_Point *pts[4];
t->getNodes(pts);
- double pa[2] = { pts[0]->u,pts[0]->v};
- double pb[2] = { pts[1]->u,pts[1]->v};
- double pc[2] = { pts[2]->u,pts[2]->v};
+ double pa[2] = {pts[0]->u, pts[0]->v};
+ double pb[2] = {pts[1]->u, pts[1]->v};
+ double pc[2] = {pts[2]->u, pts[2]->v};
- double a[2] = {pb[0]-pa[0],pb[1]-pa[1]};
- double b[2] = {pc[0]-pa[0],pc[1]-pa[1]};
+ double a[2] = {pb[0] - pa[0], pb[1] - pa[1]};
+ double b[2] = {pc[0] - pa[0], pc[1] - pa[1]};
double area_init = fabs(a[0] * b[1] - a[1] * b[0]);
@@ -1138,42 +1105,46 @@ bool test_move_point_parametric_triangle (BDS_Point * p, double u, double v, BDS
double ori_init = gmsh::orient2d(pa, pb, pc);
- if (p == pts[0])
- {pa[0]=u;pa[1]=v;}
- else if (p == pts[1])
- {pb[0]=u;pb[1]=v;}
- else if (p == pts[2])
- {pc[0]=u;pc[1]=v;}
+ if (p == pts[0]){
+ pa[0] = u;
+ pa[1] = v;
+ }
+ else if (p == pts[1]){
+ pb[0] = u;
+ pb[1] = v;
+ }
+ else if (p == pts[2]){
+ pc[0] = u;
+ pc[1] = v;
+ }
else
return false;
-
+
double area_final = fabs(a[0] * b[1] - a[1] * b[0]);
- if (area_final < 0.1 * area_init)return false;
+ if (area_final < 0.1 * area_init) return false;
double ori_final = gmsh::orient2d(pa, pb, pc);
// allow to move a point when a triangle was flat
return ori_init*ori_final > 0;
}
-
/**
d^2_i = (x^2_i - x)^T M (x_i - x)
= M11 (x_i - x)^2 + 2 M21 (x_i-x)(y_i-y) + M22 (y_i-y)^2
-
~:-)
-
-
*/
struct smoothVertexData{
BDS_Point *p;
GFace *gf;
- double scalu,scalv;
- std::list < BDS_Face * >ts;
+ double scalu, scalv;
+ std::list<BDS_Face*> ts;
};
-double smoothing_objective_function(double U, double V, BDS_Point *v, std::list < BDS_Face * >&ts,double su, double sv,GFace *gf){
-
- GPoint gp = gf->point(U*su,V*sv);
+double smoothing_objective_function(double U, double V, BDS_Point *v,
+ std::list<BDS_Face*> &ts, double su, double sv,
+ GFace *gf)
+{
+ GPoint gp = gf->point(U * su, V * sv);
const double oldX = v->X;
const double oldY = v->Y;
@@ -1182,12 +1153,12 @@ double smoothing_objective_function(double U, double V, BDS_Point *v, std::list
v->Y = gp.y();
v->Z = gp.z();
- std::list < BDS_Face * >::iterator it = ts.begin();
- std::list < BDS_Face * >::iterator ite = ts.end();
- double qMin=1;
+ std::list<BDS_Face*>::iterator it = ts.begin();
+ std::list<BDS_Face*>::iterator ite = ts.end();
+ double qMin = 1.;
while(it != ite) {
BDS_Face *t = *it;
- qMin = std::min(qmTriangle(*it,QMTRI_RHO),qMin);
+ qMin = std::min(qmTriangle(*it, QMTRI_RHO), qMin);
++it;
}
v->X = oldX;
@@ -1197,66 +1168,69 @@ double smoothing_objective_function(double U, double V, BDS_Point *v, std::list
}
void deriv_smoothing_objective_function(double U, double V,
- double &F, double &dFdU, double &dFdV,void *data){
+ double &F, double &dFdU, double &dFdV,
+ void *data)
+{
smoothVertexData *svd = (smoothVertexData*)data;
BDS_Point *v = svd->p;
const double LARGE = 1.e5;
const double SMALL = 1./LARGE;
- F = smoothing_objective_function(U,V,v,svd->ts,svd->scalu,svd->scalv,svd->gf);
- double F_U = smoothing_objective_function(U+SMALL,V,v,svd->ts,svd->scalu,svd->scalv,svd->gf);
- double F_V = smoothing_objective_function(U,V+SMALL,v,svd->ts,svd->scalu,svd->scalv,svd->gf);
- dFdU = (F_U-F)*LARGE;
- dFdV = (F_V-F)*LARGE;
+ F = smoothing_objective_function(U, V, v, svd->ts,
+ svd->scalu, svd->scalv, svd->gf);
+ double F_U = smoothing_objective_function(U + SMALL, V, v, svd->ts,
+ svd->scalu, svd->scalv, svd->gf);
+ double F_V = smoothing_objective_function(U, V + SMALL, v, svd->ts,
+ svd->scalu, svd->scalv, svd->gf);
+ dFdU = (F_U - F) * LARGE;
+ dFdV = (F_V - F) * LARGE;
}
-double smooth_obj(double U, double V, void *data){
+double smooth_obj(double U, double V, void *data)
+{
smoothVertexData *svd = (smoothVertexData*)data;
- return smoothing_objective_function(U,V,svd->p,svd->ts,svd->scalu,svd->scalv,svd->gf);
+ return smoothing_objective_function(U, V, svd->p, svd->ts,
+ svd->scalu, svd->scalv, svd->gf);
}
-
-void optimize_vertex_position (GFace *GF, BDS_Point *data, double su, double sv){
+void optimize_vertex_position (GFace *GF, BDS_Point *data, double su, double sv)
+{
#ifdef HAVE_GSL
- if(data->g && data->g->classif_degree <= 1)
- return;
+ if(data->g && data->g->classif_degree <= 1) return;
smoothVertexData vd;
vd.p = data;
vd.scalu = su;
vd.scalv = sv;
vd.gf = GF;
data->getTriangles(vd.ts);
- double U=data->u,V=data->v,val;
+ double U = data->u, V = data->v, val;
val = smooth_obj(U, V, &vd);
- if (val < -.90)return;
+ if (val < -.90) return;
- minimize_2 ( smooth_obj, deriv_smoothing_objective_function, &vd, 5, U,V,val);
- std::list < BDS_Face * >::iterator it = vd.ts.begin();
- std::list < BDS_Face * >::iterator ite = vd.ts.end();
+ minimize_2(smooth_obj, deriv_smoothing_objective_function, &vd, 5, U,V,val);
+ std::list<BDS_Face*>::iterator it = vd.ts.begin();
+ std::list<BDS_Face*>::iterator ite = vd.ts.end();
while(it != ite) {
BDS_Face *t = *it;
- if (!test_move_point_parametric_triangle ( data, U, V, t)){
+ if (!test_move_point_parametric_triangle(data, U, V, t)){
return;
}
++it;
}
-
+
data->u = U;
data->v = V;
- GPoint gp = GF->point(U*su,V*sv);
+ GPoint gp = GF->point(U * su, V * sv);
data->X = gp.x();
data->Y = gp.y();
data->Z = gp.z();
#endif
}
-
-bool BDS_Mesh::smooth_point_centroid(BDS_Point * p, GFace *gf, bool test_quality)
+bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
{
-
- if (!p->config_modified)return false;
- if(p->g && p->g->classif_degree <= 1)
- return false;
+ if (!p->config_modified) return false;
+ if(p->g && p->g->classif_degree <= 1) return false;
std::list < BDS_Edge * >::iterator eit = p->edges.begin();
while(eit != p->edges.end()) {
@@ -1268,17 +1242,17 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point * p, GFace *gf, bool test_quality
double V = 0;
double LC = 0;
- std::list < BDS_Face * >ts;
+ std::list<BDS_Face*> ts;
p->getTriangles(ts);
- std::list < BDS_Face * >::iterator it = ts.begin();
- std::list < BDS_Face * >::iterator ite = ts.end();
+ std::list<BDS_Face*>::iterator it = ts.begin();
+ std::list<BDS_Face*>::iterator ite = ts.end();
double sTot = 0;
while(it != ite) {
BDS_Face *t = *it;
BDS_Point *n[4];
t->getNodes(n);
- // double S = fabs(surface_triangle(n[0],n[1],n[2]));
+ // double S = fabs(surface_triangle(n[0],n[1],n[2]));
double S = 1;
sTot += S;
U += (n[0]->u + n[1]->u + n[2]->u) *S;
@@ -1290,17 +1264,17 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point * p, GFace *gf, bool test_quality
V /= (3.*sTot);
LC /= (3.*sTot);
- GPoint gp = gf->point(U*scalingU,V*scalingV);
+ GPoint gp = gf->point(U * scalingU, V * scalingV);
const double oldX = p->X;
const double oldY = p->Y;
const double oldZ = p->Z;
- double oldU=p->u;
- double oldV=p->v;
+ double oldU = p->u;
+ double oldV = p->v;
it = ts.begin();
- double s1=0,s2=0;
+ double s1 = 0, s2 = 0;
double newWorst = 1.0;
double oldWorst = 1.0;
@@ -1310,40 +1284,39 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point * p, GFace *gf, bool test_quality
t->getNodes(n);
p->u = U;
p->v = V;
- double snew = fabs(surface_triangle_param(n[0],n[1],n[2]));
+ double snew = fabs(surface_triangle_param(n[0], n[1], n[2]));
s1 += snew;
p->u = oldU;
p->v = oldV;
- double sold = fabs(surface_triangle_param(n[0],n[1],n[2]));
+ double sold = fabs(surface_triangle_param(n[0], n[1], n[2]));
s2 += sold;
// printf("%22.15E %22.15E\n",snew,sold);
- if ( snew < .1 * sold) return false;
+ if (snew < .1 * sold) return false;
if (1 || test_quality){
p->X = gp.x();
p->Y = gp.y();
p->Z = gp.z();
- newWorst = std::min(newWorst,qmTriangle(*it,QMTRI_RHO));
+ newWorst = std::min(newWorst, qmTriangle(*it, QMTRI_RHO));
double norm1[3],norm2[3];
- normal_triangle(n[0],n[1],n[2],norm1);
+ normal_triangle(n[0], n[1], n[2], norm1);
p->X = oldX;
p->Y = oldY;
p->Z = oldZ;
- normal_triangle(n[0],n[1],n[2],norm2);
- oldWorst = std::min(oldWorst,qmTriangle(*it,QMTRI_RHO));
- double ps; prosca (norm1,norm2,&ps);
- if (ps < .5)return false;
- }
+ normal_triangle(n[0], n[1], n[2], norm2);
+ oldWorst = std::min(oldWorst, qmTriangle(*it, QMTRI_RHO));
+ double ps; prosca (norm1, norm2, &ps);
+ if (ps < .5) return false;
+ }
++it;
}
// printf("%22.15E %22.15E %22.15E\n",s1,s2,fabs(s2-s1));
- if (fabs(s2-s1) > 1.e-14 * (s2+s1))return false;
+ if (fabs(s2-s1) > 1.e-14 * (s2 + s1)) return false;
-
- if (test_quality && newWorst < oldWorst){
- return false;
- }
+ if (test_quality && newWorst < oldWorst){
+ return false;
+ }
p->u = U;
p->v = V;
@@ -1363,7 +1336,7 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
{
if (!p->config_modified)return false;
- if(p->g && p->g->classif_degree <= 1)
+ if (p->g && p->g->classif_degree <= 1)
return false;
double U = 0;
@@ -1376,8 +1349,8 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
if ((*eit)->numfaces() == 1) return false;
BDS_Point *op = ((*eit)->p1 == p) ? (*eit)->p2 : (*eit)->p1;
const double l_e = (*eit)->length();
- U += op->u*l_e;
- V += op->v*l_e;
+ U += op->u * l_e;
+ V += op->v * l_e;
tot_length += l_e;
LC += op->lc();
++eit;
@@ -1387,19 +1360,18 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
V /= tot_length;
LC /= p->edges.size();
- std::list < BDS_Face * >ts;
+ std::list<BDS_Face*> ts;
p->getTriangles(ts);
- std::list < BDS_Face * >::iterator it = ts.begin();
- std::list < BDS_Face * >::iterator ite = ts.end();
+ std::list<BDS_Face*>::iterator it = ts.begin();
+ std::list<BDS_Face*>::iterator ite = ts.end();
while(it != ite) {
BDS_Face *t = *it;
- if (!test_move_point_parametric_triangle ( p, U, V, t))
+ if (!test_move_point_parametric_triangle(p, U, V, t))
return false;
++it;
}
-
- GPoint gp = gf->point(U*scalingU,V*scalingV);
+ GPoint gp = gf->point(U * scalingU, V * scalingV);
p->u = U;
p->v = V;
p->lc() = LC;
@@ -1412,72 +1384,65 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
++eit;
}
-
return true;
}
-
struct recombine_T
{
- const BDS_Edge *e ;
+ const BDS_Edge *e;
double angle ;
- recombine_T (const BDS_Edge *_e );
- bool operator < ( const recombine_T & other ) const
+ recombine_T(const BDS_Edge *_e);
+ bool operator < (const recombine_T &other) const
{
return angle < other.angle;
}
};
-recombine_T::recombine_T (const BDS_Edge *_e )
+recombine_T::recombine_T(const BDS_Edge *_e)
: e(_e)
{
BDS_Point *op[2];
BDS_Point *p1 = e->p1;
BDS_Point *p2 = e->p2;
e->oppositeof(op);
- BDS_Vector v1 (*p1,*op[0]);
- BDS_Vector v2 (*op[0],*p2);
- BDS_Vector v3 (*p2,*op[1]);
- BDS_Vector v4 (*op[1],*p1);
+ BDS_Vector v1(*p1, *op[0]);
+ BDS_Vector v2(*op[0], *p2);
+ BDS_Vector v3(*p2, *op[1]);
+ BDS_Vector v4(*op[1], *p1);
angle = fabs(90.-v1.angle_deg(v2));
- angle = std::max(fabs(90.-v2.angle_deg(v3)),angle);
- angle = std::max(fabs(90.-v3.angle_deg(v4)),angle);
- angle = std::max(fabs(90.-v4.angle_deg(v1)),angle);
+ angle = std::max(fabs(90.-v2.angle_deg(v3)), angle);
+ angle = std::max(fabs(90.-v3.angle_deg(v4)), angle);
+ angle = std::max(fabs(90.-v4.angle_deg(v1)), angle);
}
void BDS_Mesh::recombineIntoQuads (const double angle_limit, GFace *gf)
{
Msg(INFO,"Recombining triangles for surface %d",gf->tag());
- for (int i=0;i<5;i++)
- {
+ for (int i = 0; i < 5 ; i++){
bool rec = false;
std::set<recombine_T> _pairs;
for(std::list < BDS_Edge * >::const_iterator it = edges.begin();
- it != edges.end(); ++it)
- {
- if (!(*it)->deleted && (*it)->numfaces () == 2)
- _pairs.insert ( recombine_T ( *it ) );
- }
+ it != edges.end(); ++it){
+ if (!(*it)->deleted && (*it)->numfaces () == 2)
+ _pairs.insert(recombine_T(*it));
+ }
+
+ std::set<recombine_T>::iterator itp = _pairs.begin();
+
+ while (itp != _pairs.end()){
+ rec |= recombine_edge ((BDS_Edge*)itp->e);
+ ++itp;
+ }
- std::set<recombine_T>::iterator itp = _pairs.begin();
-
- while (itp != _pairs.end())
- {
- rec |= recombine_edge ((BDS_Edge*)itp->e);
- ++itp;
- }
-
if (!rec) break;
std::set<BDS_Point*,PointLessThan>::iterator itpt = points.begin();
- while (itpt != points.end())
- {
-
- smooth_point_parametric(*itpt,gf);
- ++itpt;
- }
+ while (itpt != points.end()){
+ smooth_point_parametric(*itpt,gf);
+ ++itpt;
+ }
}
}