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

*** empty log message ***

parent 271d2e03
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 276 additions and 311 deletions
Mesh/BDS.cpp
+ 276
− 311
View file @ a869aef8
// $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
//
......@@ -41,16 +41,18 @@ void outputScalarField(std::list < BDS_Face * >t, const char *iii, int param, GF
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,
......@@ -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,39 +238,31 @@ 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);
int ix = 0;
int ixMax = 300;
while(1)
{
while(1){
std::vector<BDS_Edge*> intersected;
std::list<BDS_Edge*>::iterator it = edges.begin();
while (it!=edges.end())
{
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(),
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));
......@@ -286,15 +280,13 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
// }
// printf("%d %d\n",intersected.size(),ix);
if (!intersected.size() || ix > 1000)
{
if(!intersected.size() || ix > 1000){
BDS_Edge *eee = find_edge (num1, num2);
if (!eee)
{
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);
Msg(GERROR," edge %d %d cannot be recovered at all, look at debugp.pos "
"and debugr.pos", num1, num2);
return 0;
}
return eee;
......@@ -302,14 +294,14 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
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());
// 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,13 +382,10 @@ 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;
// if(tfound) return tfound;
BDS_Face *t = new BDS_Face(e1, e2, e3);
triangles.push_back(t);
......@@ -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);
}
......@@ -491,7 +477,6 @@ void recur_tag(BDS_Face * t, BDS_GeomEntity * g)
}
}
double PointLessThanLexicographic::t = 0;
double BDS_Vector::t = 0;
......@@ -538,7 +523,6 @@ BDS_Mesh ::~ BDS_Mesh ()
DESTROOOY(triangles.begin(), triangles.end());
}
bool BDS_Mesh::split_face(BDS_Face *f, BDS_Point *mid)
{
BDS_Point *p1 = f->e1->commonvertex(f->e2);
......@@ -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;
......@@ -697,13 +679,11 @@ 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);
......@@ -711,7 +691,6 @@ bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1,BDS_Point *_p2,
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));
......@@ -720,8 +699,6 @@ bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1,BDS_Point *_p2,
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,
......@@ -756,7 +733,6 @@ bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1,BDS_Point *_p2,BDS_Poi
if (minb > mina) return true;
// if (mina > minb && cosnq <= cosonq)return true;
return false;
}
void swap_config(BDS_Edge * e,
......@@ -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;
......@@ -891,7 +865,6 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
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);
......@@ -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)
......@@ -1113,13 +1082,11 @@ 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
// 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];
......@@ -1138,12 +1105,18 @@ 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;
......@@ -1154,14 +1127,10 @@ bool test_move_point_parametric_triangle (BDS_Point * p, double u, double v, BDS
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{
......@@ -1171,8 +1140,10 @@ struct smoothVertexData{
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){
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;
......@@ -1184,7 +1155,7 @@ double smoothing_objective_function(double U, double V, BDS_Point *v, std::list
std::list<BDS_Face*>::iterator it = ts.begin();
std::list<BDS_Face*>::iterator ite = ts.end();
double qMin=1;
double qMin = 1.;
while(it != ite) {
BDS_Face *t = *it;
qMin = std::min(qmTriangle(*it, QMTRI_RHO), qMin);
......@@ -1197,28 +1168,34 @@ 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);
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;
......@@ -1250,13 +1227,10 @@ void optimize_vertex_position (GFace *GF, BDS_Point *data, double su, double sv)
#endif
}
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->g && p->g->classif_degree <= 1) return false;
std::list < BDS_Edge * >::iterator eit = p->edges.begin();
while(eit != p->edges.end()) {
......@@ -1340,7 +1314,6 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point * p, GFace *gf, bool test_quality
// 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 (test_quality && newWorst < oldWorst){
return false;
}
......@@ -1398,7 +1371,6 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
++it;
}
GPoint gp = gf->point(U * scalingU, V * scalingV);
p->u = U;
p->v = V;
......@@ -1412,11 +1384,9 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
++eit;
}
return true;
}
struct recombine_T
{
const BDS_Edge *e;
......@@ -1449,22 +1419,19 @@ 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)
{
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())
{
while (itp != _pairs.end()){
rec |= recombine_edge ((BDS_Edge*)itp->e);
++itp;
}
......@@ -1472,9 +1439,7 @@ void BDS_Mesh::recombineIntoQuads (const double angle_limit, GFace *gf)
if (!rec) break;
std::set<BDS_Point*,PointLessThan>::iterator itpt = points.begin();
while (itpt != points.end())
{
while (itpt != points.end()){
smooth_point_parametric(*itpt,gf);
++itpt;
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment