Skip to content
Snippets Groups Projects
Commit d3a7255e authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

*** empty log message ***

parent 1b7bfcd8
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Mesh/BDS.cpp
+ 39
− 196
View file @ d3a7255e
// $Id: BDS.cpp,v 1.94 2008-01-24 09:35:41 remacle Exp $
// $Id: BDS.cpp,v 1.95 2008-01-26 17:47:58 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
......@@ -221,18 +221,23 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
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(WARNING," 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());
Msg(WARNING," 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));
......@@ -249,6 +254,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
// num1,num2,ix);
// ix = 0;
// }
// printf("%d %d\n",intersected.size(),ix);
if (!intersected.size() || ix > 1000)
{
......@@ -265,7 +271,9 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
}
int ichoice = ix++ % intersected.size();
swap_edge ( intersected[ichoice] , BDS_SwapEdgeTestQuality (false) );
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;
}
......@@ -375,24 +383,6 @@ BDS_Face *BDS_Mesh::add_quadrangle(BDS_Edge * e1, BDS_Edge * e2,
return t;
}
// BDS_Tet *BDS_Mesh::add_tet(int p1, int p2, int p3, int p4)
// {
// BDS_Edge *e1 = add_edge(p1, p2);
// BDS_Edge *e2 = add_edge(p2, p3);
// BDS_Edge *e3 = add_edge(p3, p1);
// BDS_Edge *e4 = add_edge(p1, p4);
// BDS_Edge *e5 = add_edge(p2, p4);
// BDS_Edge *e6 = add_edge(p3, p4);
// BDS_Face *t1 = add_triangle(e1, e2, e3);
// BDS_Face *t2 = add_triangle(e1, e4, e5);
// BDS_Face *t3 = add_triangle(e2, e6, e5);
// BDS_Face *t4 = add_triangle(e3, e4, e6);
// BDS_Tet *t = new BDS_Tet(t1, t2, t3, t4);
// tets.push_back(t);
// return t;
// }
void BDS_Mesh::del_face(BDS_Face * t)
{
t->e1->del(t);
......@@ -404,11 +394,9 @@ void BDS_Mesh::del_face(BDS_Face * t)
void BDS_Mesh::del_edge(BDS_Edge * e)
{
// edges.erase(e);
e->p1->del(e);
e->p2->del(e);
e->deleted = true;
// edges_to_delete.insert(e);
}
void BDS_Mesh::del_point(BDS_Point * p)
......@@ -477,12 +465,6 @@ void recur_tag(BDS_Face * t, BDS_GeomEntity * g)
double PointLessThanLexicographic::t = 0;
double BDS_Vector::t = 0;
bool BDS_Mesh::import_view(PView *view, const double tolerance)
{
Msg(GERROR, "View import must be reinterfaced");
return false;
}
template < class IT > void DESTROOOY(IT beg, IT end)
{
while(beg != end) {
......@@ -551,22 +533,6 @@ bool BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
e->oppositeof(op);
// double qt1 = qmTriangle(p1,op[0],mid,QMTRI_RHO);
// double qt2 = qmTriangle(p2,op[0],mid,QMTRI_RHO);
// double qt3 = qmTriangle(p1,op[1],mid,QMTRI_RHO);
// double qt4 = qmTriangle(p2,op[1],mid,QMTRI_RHO);
// if (qt1 < 1.e-5 || qt2 < 1.e-5 || qt3 < 1.e-5 || qt4 < 1.e-5)
// {
// return false;
// }
// double l1 = sqrt((op[0]->X-op[1]->X) *(op[0]->X-op[1]->X) +
// (op[0]->Y-op[1]->Y) *(op[0]->Y-op[1]->Y) +
// (op[0]->Z-op[1]->Z) *(op[0]->Z-op[1]->Z) );
// if (l1 < 0.1* mid->lc()) return false;
BDS_Point *pts1[4];
e->faces(0)->getNodes(pts1);
......@@ -613,14 +579,6 @@ bool BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
BDS_Edge *mid_op2 = new BDS_Edge(mid, op[1]);
edges.push_back(mid_op2);
// p1_mid->target_length = t_l;
// op1_mid->target_length = t_l;
// mid_p2->target_length = t_l;
// mid_op2->target_length = t_l;
// printf("split ends 1 %d (%d %d) %d %d \n",
// p1_op1->numfaces(), p1->iD, op[0]->iD,
// op1_mid->numfaces(),p1_mid->numfaces());
BDS_Face *t1, *t2, *t3, *t4;
if(orientation == 1) {
t1 = new BDS_Face(op1_mid, p1_op1, p1_mid);
......@@ -660,131 +618,6 @@ bool BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
return true;
}
void BDS_Mesh::saturate_edge(BDS_Edge * e, std::vector<BDS_Point *> &mids)
{
BDS_Point *op[2];
BDS_Point *p1 = e->p1;
BDS_Point *p2 = e->p2;
BDS_Point *pts1[4];
e->faces(0)->getNodes(pts1);
int orientation = 0;
for(int i = 0; i < 3; i++) {
if(pts1[i] == p1) {
if(pts1[(i + 1) % 3] == p2)
orientation = 1;
else
orientation = -1;
break;
}
}
// printf("sat %d\n",mids.size());
// we should project
e->oppositeof(op);
BDS_GeomEntity *g1 = 0, *g2 = 0, *ge = e->g;
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));
BDS_Edge *op2_p2 = find_edge(op[1], p2, e->faces(1));
if(e->faces(0)) {
g1 = e->faces(0)->g;
del_face(e->faces(0));
}
// not a bug !!!
if(e->faces(0)) {
g2 = e->faces(0)->g;
del_face(e->faces(0));
}
// double t_l = e->target_length;
del_edge(e);
// create all the sub-edges of e
std::vector<BDS_Edge*> subs;
for (unsigned int i = 0; i < mids.size() + 1; i++)
{
BDS_Point *a,*b;
if (i == 0)a = p1;
else a = mids[i-1];
if (i == mids.size())b = p2;
else b = mids[i];
BDS_Edge *sub = new BDS_Edge(a, b);
// printf("%g %g -> %g %g\n",a->X,a->Y,b->X,b->Y);
edges.push_back(sub);
subs.push_back(sub);
sub->g = ge;
}
// create edges that connect op1 and op2
std::vector<BDS_Edge*> conn1,conn2;
for (unsigned int i = 0; i < mids.size(); i++)
{
BDS_Edge *c1 = new BDS_Edge(mids[i], op[0]);
edges.push_back(c1);
conn1.push_back(c1);
BDS_Edge *c2 = new BDS_Edge(mids[i], op[1]);
edges.push_back(c2);
conn2.push_back(c2);
c1->g = g1;
c2->g = g2;
mids[i]->g = ge;
}
// create the triangles
for (unsigned int i = 0; i < mids.size() + 1; i++)
{
BDS_Edge *e1,*e2,*e3=subs[i];
// printf("--> %d %g %g ->%d %g %g\n",e3->p1->iD,e3->p1->X,e3->p1->Y,e3->p2->iD,e3->p2->X,e3->p2->Y);
if (i==0)e1 = p1_op1;
else e1 = conn1[i-1];
if (i==mids.size())e2 = op1_p2;
else e2 = conn1[i];
// printf("--> %d %g %g ->%d %g %g\n",e1->p1->iD,e1->p1->X,e1->p1->Y,e1->p2->iD,e1->p2->X,e1->p2->Y);
// printf("--> %d %g %g ->%d %g %g\n",e2->p1->iD,e2->p1->X,e2->p1->Y,e2->p2->iD,e2->p2->X,e2->p2->Y);
BDS_Face *t1;
if(orientation == 1)
t1 = new BDS_Face(e3, e2, e1);
else
t1 = new BDS_Face(e3, e1, e2);
if (i==0)e1 = p1_op2;
else e1 = conn2[i-1];
if (i==mids.size())e2 = op2_p2;
else e2 = conn2[i];
BDS_Face *t2;
if(orientation != 1)
t2 = new BDS_Face(e3, e2, e1);
else
t2 = new BDS_Face(e3, e1, e2);
t1->g = g1;
t2->g = g2;
triangles.push_back(t1);
triangles.push_back(t2);
}
// config has changed
p1->config_modified = true;
p2->config_modified = true;
op[0]->config_modified = true;
op[1]->config_modified = true;
}
// This function does actually the swap without taking into account
// the feasability of the operation. Those conditions have to be
// taken into account before doing the edge swap
......@@ -792,25 +625,32 @@ void BDS_Mesh::saturate_edge(BDS_Edge * e, std::vector<BDS_Point *> &mids)
bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1,BDS_Point *_p2,
BDS_Point *_q1,BDS_Point *_q2) const
{
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-10 * (s1+s2))return false;
if (s3 < .02 * (s1+s2) || s4 < .02 * (s1+s2))return false;
return true;
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 !
}
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,
......@@ -848,8 +688,6 @@ bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1,BDS_Point *_p2,BDS_Poi
}
void swap_config(BDS_Edge * e,
BDS_Point **p11,BDS_Point **p12,BDS_Point **p13,
BDS_Point **p21,BDS_Point **p22,BDS_Point **p23,
......@@ -1410,15 +1248,20 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point * p, GFace *gf, bool test_quality
// printf("%22.15E %22.15E\n",snew,sold);
if ( snew < .1 * sold) return false;
if (test_quality){
if (1 || test_quality){
p->X = gp.x();
p->Y = gp.y();
p->Z = gp.z();
newWorst = std::min(newWorst,qmTriangle(*it,QMTRI_RHO));
double norm1[3],norm2[3];
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;
}
++it;
}
......
Mesh/BDS.h
+ 3
− 5
View file @ d3a7255e
......@@ -379,9 +379,9 @@ class BDS_SwapEdgeTest
class BDS_SwapEdgeTestQuality : public BDS_SwapEdgeTest
{
bool testQuality;
bool testQuality, testSmallTriangles;
public:
BDS_SwapEdgeTestQuality (bool a) : testQuality(a){}
BDS_SwapEdgeTestQuality (bool a, bool b=true) : testQuality(a),testSmallTriangles(b){}
virtual bool operator() (BDS_Point *p1,BDS_Point *p2,
BDS_Point *q1,BDS_Point *q2) const ;
virtual bool operator() (BDS_Point *p1,BDS_Point *p2, BDS_Point *p3,
......@@ -466,14 +466,12 @@ public:
bool smooth_point_centroid(BDS_Point * p, GFace *gf, bool test_quality = false);
bool move_point(BDS_Point *p , double X, double Y, double Z);
bool split_edge(BDS_Edge *, BDS_Point *);
void saturate_edge(BDS_Edge *, std::vector<BDS_Point *>&);
bool split_face(BDS_Face *, BDS_Point *);
bool edge_constraint ( BDS_Point *p1, BDS_Point *p2 );
bool recombine_edge ( BDS_Edge *e );
// Global operators
void cleanup();
void recombineIntoQuads (const double angle, GFace *gf);
// io's
bool import_view(PView *view, const double tolerance);
};
void outputScalarField(std::list < BDS_Face * >t, const char *fn, int param);
......
Mesh/meshGFace.cpp
+ 4
− 4
View file @ d3a7255e
// $Id: meshGFace.cpp,v 1.112 2008-01-24 09:35:41 remacle Exp $
// $Id: meshGFace.cpp,v 1.113 2008-01-26 17:47:58 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
......@@ -1280,9 +1280,9 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
if (!AlgoDelaunay2D ( gf ))
{
gmshRefineMeshBDS (gf,*m,CTX.mesh.refine_steps,true);
gmshOptimizeMeshBDS(gf, *m, 2,&recover_map);
gmshOptimizeMeshBDS(gf, *m, 2);
gmshRefineMeshBDS (gf,*m,-CTX.mesh.refine_steps,false);
gmshOptimizeMeshBDS(gf, *m, -2,&recover_map);
gmshOptimizeMeshBDS(gf, *m, 2,&recover_map);
if (gf->meshAttributes.recombine)
{
......
Mesh/meshGFaceBDS.cpp
+ 159
− 84
View file @ d3a7255e
// $Id: meshGFaceBDS.cpp,v 1.1 2008-01-24 09:35:41 remacle Exp $
// $Id: meshGFaceBDS.cpp,v 1.2 2008-01-26 17:47:58 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
......@@ -65,6 +65,30 @@ inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f,
return l1+l2;
}
inline double computeEdgeMiddleCoord(BDS_Point *p1, BDS_Point *p2, GFace *f,
double SCALINGU, double SCALINGV)
{
if (f->geomType() == GEntity::Plane)
return 0.5;
GPoint GP = f->point (SPoint2(0.5*(p1->u+p2->u)*SCALINGU,0.5*(p1->v+p2->v)*SCALINGV));
const double dx1 = p1->X-GP.x();
const double dy1 = p1->Y-GP.y();
const double dz1 = p1->Z-GP.z();
const double l1 = sqrt(dx1*dx1+dy1*dy1+dz1*dz1);
const double dx2 = p2->X-GP.x();
const double dy2 = p2->Y-GP.y();
const double dz2 = p2->Z-GP.z();
const double l2 = sqrt(dx2*dx2+dy2*dy2+dz2*dz2);
if (l1 > l2)
return 0.25 * (l1+l2) / l1;
else
return 0.25 * (3*l2-l1)/l2;
}
inline double computeEdgeLinearLength(BDS_Edge *e, GFace *f, double SCALINGU, double SCALINGV)
{
if (f->geomType() == GEntity::Plane)
......@@ -386,6 +410,7 @@ void splitEdgePass ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
if ((*it)->numfaces() == 2 && (lone > MAXE_)){
const double coord = 0.5;
//const double coord = computeEdgeMiddleCoord((*it)->p1,(*it)->p2,gf,m.scalingU,m.scalingV);
BDS_Point *mid ;
mid = m.add_point(++m.MAXPOINTNUMBER,
coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
......@@ -438,7 +463,6 @@ void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q)
}
}
void gmshRefineMeshBDS (GFace *gf,
BDS_Mesh &m,
const int NIT,
......@@ -554,99 +578,90 @@ void gmshRefineMeshBDS (GFace *gf,
Msg(DEBUG1," ---------------------------------------");
}
void gmshOptimizeMeshBDS(GFace *gf,
BDS_Mesh &m,
const int NIT,
std::map<BDS_Point*,MVertex*> *recover_map = 0)
{
int nb_swap;
gmshDelaunayizeBDS ( gf, m, nb_swap );
void allowAppearanceofEdge (BDS_Point *p1, BDS_Point *p2){
}
for (int ITER = 0;ITER < 3;ITER++){
int nb_smooth;
smoothVertexPass ( gf,m,nb_smooth,true);
double LIMIT = .1;
for (int KK=0;KK<4;KK++){
// swap edges that provide a better configuration
int NN1 = m.edges.size();
int NN2 = 0;
std::list<BDS_Edge*>::iterator it = m.edges.begin();
while (1)
void invalidEdgesPeriodic(BDS_Mesh &m,
std::map<BDS_Point*,MVertex*> *recover_map,
std::set<BDS_Edge*> &toSplit)
{
if (NN2++ >= NN1)break;
if (evalSwapForOptimize(*it,gf,m))
m.swap_edge ( *it , BDS_SwapEdgeTestQuality(false));
++it;
}
m.cleanup();
}
// first look for degenerated vertices
// then collapse small edges (take care not to create overlapping triangles)
// in case of periodic surfaces, split all edges that are problematic
for (int KK=0;KK<1;KK++){
int NN1 = m.edges.size();
int NN2 = 0;
std::list<BDS_Edge*>::iterator it = m.edges.begin();
std::vector<BDS_Edge *> toSplit;
while (1){
if (NN2++ >= NN1)break;
if((*it)->numfaces() == 2){
if (recover_map){
std::map<BDS_Point*,MVertex*>::iterator itp1 = recover_map->find((*it)->p1);
std::map<BDS_Point*,MVertex*>::iterator itp2 = recover_map->find((*it)->p2);
BDS_Point *op[2];
(*it)->oppositeof (op);
std::map<BDS_Point*,MVertex*>::iterator itp3 = recover_map->find(op[0]);
std::map<BDS_Point*,MVertex*>::iterator itp4 = recover_map->find(op[1]);
// this edge goes from one side to the other of the periodic parametric space !
if (itp1 != recover_map->end() && itp2 != recover_map->end() && itp1->second == itp2->second)
toSplit.push_back(*it);
// this edge is internal but the 2 adjacent triangles are the same in the real space
if (itp3 != recover_map->end() && itp4 != recover_map->end() && itp3->second == itp4->second)
{
// the first point is internal, split both edges that go from this one to the two opposites (that are the same)
BDS_Edge *e1 , *e2 ;
// printf ("edge prob %d %d\n",(*it)->p1->iD,(*it)->p2->iD);
if (itp1 == recover_map->end()){
e1 = m.find_edge ((*it)->p1,itp3->first);
e2 = m.find_edge ((*it)->p1,itp4->first);
if (e1 && e1->numfaces() == 2)
toSplit.push_back(e1);
if (e2 && e2->numfaces() == 2)
toSplit.push_back(e2);
}
else if (itp2 == recover_map->end()){
e1 = m.find_edge ((*it)->p2,itp3->first);
e2 = m.find_edge ((*it)->p2,itp4->first);
if (e1 && e1->numfaces() == 2)
toSplit.push_back(e1);
if (e2 && e2->numfaces() == 2)
toSplit.push_back(e2);
std::set<MVertex *> degenerated;
while (it != m.edges.end()){
BDS_Edge *e = *it;
if (!e->deleted && e->numfaces() == 1){
std::map<BDS_Point*,MVertex*>::iterator itp1 = recover_map->find(e->p1);
std::map<BDS_Point*,MVertex*>::iterator itp2 = recover_map->find(e->p2);
if (itp1 != recover_map->end() &&
itp2 != recover_map->end() &&
itp1->second == itp2->second){
degenerated.insert(itp1->second);
}
else{
printf("zarbi\n");
}
++it;
}
// toSplit.push_back(*it);
// printf("%d degenerated\n",degenerated.size());
toSplit.clear();
it = m.edges.begin();
std::map< std::pair < MVertex*, BDS_Point*> , BDS_Edge *> touchPeriodic;
while (it != m.edges.end()){
BDS_Edge *e = *it;
if (!e->deleted && e->numfaces() == 2){
std::map<BDS_Point*,MVertex*>::iterator itp1 = recover_map->find(e->p1);
std::map<BDS_Point*,MVertex*>::iterator itp2 = recover_map->find(e->p2);
if (itp1 != recover_map->end() &&
itp2 != recover_map->end() &&
itp1->second == itp2->second) toSplit.insert(e);
else if (itp1 != recover_map->end() && itp2 == recover_map->end()){
std::pair<MVertex*,BDS_Point*> a ( itp1->second, e->p2 );
std::map< std::pair < MVertex*, BDS_Point*> , BDS_Edge *> :: iterator itf =
touchPeriodic.find (a);
if (itf == touchPeriodic.end()) touchPeriodic[a] = e;
else if (degenerated.find(itp1->second)==degenerated.end()){toSplit.insert(e);toSplit.insert(itf->second);}
}
else if (itp1 == recover_map->end() && itp2 != recover_map->end()){
std::pair<MVertex*,BDS_Point*> a ( itp2->second, e->p1 );
std::map< std::pair < MVertex*, BDS_Point*> , BDS_Edge *> :: iterator itf =
touchPeriodic.find (a);
if (itf == touchPeriodic.end()) touchPeriodic[a] = e;
else if (degenerated.find(itp2->second)==degenerated.end()){toSplit.insert(e);toSplit.insert(itf->second);}
}
}
++it;
}
// printf("%d edges to splitounette\n",toSplit.size());
}
// consider p1 and p2, 2 vertices that are different in the parametric
// plane BUT are the same in the real plane
// consider a vertex v that is internal
// if p1 is the start and the end of a degenerated edge, then allow edges p1 v and p2 v
// if not do not allow p1 v and p2 v, split them
// if p1 p2 exists and it is internal, split it
int gmshSolveInvalidPeriodic(GFace *gf,
BDS_Mesh &m,
std::map<BDS_Point*,MVertex*> *recover_map){
// return 0;
std::set<BDS_Edge*> toSplit;
invalidEdgesPeriodic(m,recover_map,toSplit);
std::set<BDS_Edge*>::iterator ite = toSplit.begin();
// printf("%d edges to split\n",toSplit.size());
for (int i=0;i<toSplit.size();i++){
BDS_Edge *e = toSplit[i];
if (!e->deleted){
for (;ite !=toSplit.end();++ite){
BDS_Edge *e = *ite;
if (!e->deleted && e->numfaces() == 2){
const double coord = 0.5;
BDS_Point *mid ;
// printf("%d %d\n",e->p1->iD,e->p2->iD);
mid = m.add_point(++m.MAXPOINTNUMBER,
coord * e->p1->u + (1 - coord) * e->p2->u,
coord * e->p1->v + (1 - coord) * e->p2->v,gf);
// printf("%d %d %d\n",e->p1->iD,e->p2->iD,mid->iD);
mid->lcBGM() = BGM_MeshSize(gf,
(coord * e->p1->u + (1 - coord) * e->p2->u)*m.scalingU,
(coord * e->p1->v + (1 - coord) * e->p2->v)*m.scalingV,
......@@ -655,8 +670,68 @@ void gmshOptimizeMeshBDS(GFace *gf,
if(!m.split_edge ( e, mid )) m.del_point(mid);
}
}
int nb_smooth;
if (toSplit.size())smoothVertexPass ( gf,m,nb_smooth,true);
m.cleanup();
return toSplit.size();
}
void gmshOptimizeMeshBDS(GFace *gf,
BDS_Mesh &m,
const int NIT,
std::map<BDS_Point*,MVertex*> *recover_map = 0)
{
int nb_swap;
gmshDelaunayizeBDS ( gf, m, nb_swap );
for (int ITER = 0;ITER < 3;ITER++){
double LIMIT = .1;
for (int KK=0;KK<4;KK++){
// swap edges that provide a better configuration
int NN1 = m.edges.size();
int NN2 = 0;
std::list<BDS_Edge*>::iterator it = m.edges.begin();
while (1){
if (NN2++ >= NN1)break;
if (evalSwapForOptimize(*it,gf,m))
m.swap_edge ( *it , BDS_SwapEdgeTestQuality(false));
++it;
}
m.cleanup();
int nb_smooth;
smoothVertexPass ( gf,m,nb_smooth,true);
}
}
int nbSplit = 0;
if (recover_map){
while(gmshSolveInvalidPeriodic(gf,m,recover_map)){}
}
}
// DELAUNAY BDS
// void delaunayPointInsertionBDS ( GFace *gf, BDS_Mesh &m, BDS_Vertex *v, BDS_Face *f){
// const double p[2] = {v->u,v->v};
// BDS_Edge *e1 = f->e1;
// BDS_Edge *e2 = f->e2;
// BDS_Edge *e3 = f->e3;
// m.split_face ( f , v );
// recur_delaunay_swap ( e1 );
// recur_delaunay_swap ( e2 );
// recur_delaunay_swap ( e3 );
// return;
// BDS_Point *n2[4];
// f1->getNodes(n1);
// const double a[2] = {v->u,v->v};
// const double b[2] = {v->u,v->v};
// const double c[2] = {v->u,v->v};
// }
benchmarks/testsuite/testsuite.sh
+ 1
− 0
View file @ d3a7255e
......@@ -9,5 +9,6 @@ date >> statreport.dat
../../bin/gmsh sphere_Rhino3D_default.igs -clscale 1 -v 0 -2 -append_statreport statreport.dat
../../bin/gmsh Cone_1.brep -clscale 1 -v 0 -2 -append_statreport statreport.dat
../../bin/gmsh Cylinder_1.brep -clscale 1 -v 0 -2 -append_statreport statreport.dat
../../bin/gmsh Torus_1.brep -clscale 300 -v 0 -2 -append_statreport statreport.dat
../../bin/gmsh A319.geo -clscale 3 -v 0 -2 -append_statreport statreport.dat
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment