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

*** empty log message ***

parent 8acb8f52
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Mesh/BDS.cpp
+ 94
− 52
View file @ 34164fbb
...@@ -3,6 +3,28 @@ ...@@ -3,6 +3,28 @@
#include <math.h> #include <math.h>
#include "Numeric.h" #include "Numeric.h"
double dist_droites_gauches(BDS_Point *p1, BDS_Point *p2,
BDS_Point *p3 ,BDS_Point *p4)
{
BDS_Vector u1 ( *p2, *p1 );
BDS_Vector u2 ( *p4, *p3 );
BDS_Vector a ( *p2, *p4 );
BDS_Vector u1xu2 = u1%u2;
double x = sqrt (u1xu2 * u1xu2 );
// les droites sont paralleles
if (x == 0)
{
throw;
}
// les droites sont gauches
else
{
double y = fabs((a % u1) * u2);
return y/x;
}
}
void proj_point_plane ( double xa, double ya, double za, void proj_point_plane ( double xa, double ya, double za,
BDS_Point *p1, BDS_Point *p2, BDS_Point *p3 , BDS_Point *p1, BDS_Point *p2, BDS_Point *p3 ,
double &x, double &y, double &z) double &x, double &y, double &z)
...@@ -91,7 +113,7 @@ void BDS_Point::getTriangles (std::list<BDS_Triangle *> &t) ...@@ -91,7 +113,7 @@ void BDS_Point::getTriangles (std::list<BDS_Triangle *> &t)
int NF = (*it)->numfaces(); int NF = (*it)->numfaces();
for (int i=0;i<NF;++i) for (int i=0;i<NF;++i)
{ {
BDS_Triangle *tt = (*it)->faces[i]; BDS_Triangle *tt = (*it)->faces(i);
if (tt) if (tt)
{ {
std::list<BDS_Triangle*>::iterator tit = t.begin(); std::list<BDS_Triangle*>::iterator tit = t.begin();
...@@ -198,22 +220,23 @@ void BDS_Mesh :: add_geom (int p1, int p2) ...@@ -198,22 +220,23 @@ void BDS_Mesh :: add_geom (int p1, int p2)
void BDS_Edge::oppositeof (BDS_Point * oface[2]) const void BDS_Edge::oppositeof (BDS_Point * oface[2]) const
{ {
oface[0] = oface[1] = 0;
if (faces[0]) oface[0] = oface[1] = 0;
if (faces(0))
{ {
BDS_Point *pts[3]; BDS_Point *pts[3];
faces[0]->getNodes (pts); faces(0)->getNodes (pts);
if (pts[0] != p1 && pts[0] != p2)oface[0] = pts[0]; if (pts[0] != p1 && pts[0] != p2)oface[0] = pts[0];
else if (pts[1] != p1 && pts[1] != p2)oface[0] = pts[1]; else if (pts[1] != p1 && pts[1] != p2)oface[0] = pts[1];
else oface[0] = pts[2]; else oface[0] = pts[2];
} }
if (faces[1]) if (faces(1))
{ {
BDS_Point *pts[3]; BDS_Point *pts[3];
faces[1]->getNodes (pts); faces(1)->getNodes (pts);
if (pts[0] != p1 && pts[0] != p2)oface[1] = pts[0]; if (pts[0] != p1 && pts[0] != p2)oface[1] = pts[0];
else if (pts[1] != p1 && pts[1] != p2)oface[1] = pts[1]; else if (pts[1] != p1 && pts[1] != p2)oface[1] = pts[1];
else oface[1] = pts[2]; else oface[1] = pts[2];
} }
} }
...@@ -270,16 +293,21 @@ void BDS_Mesh :: classify ( double angle ) ...@@ -270,16 +293,21 @@ void BDS_Mesh :: classify ( double angle )
while (it!=ite) while (it!=ite)
{ {
BDS_Edge &e = *((BDS_Edge *) *it); BDS_Edge &e = *((BDS_Edge *) *it);
if ( e.numfaces() == 1) e.g = &EDGE_CLAS; if ( e.numfaces() == 1)
{
e.g = &EDGE_CLAS;
}
else if (e.numfaces() == 2 && else if (e.numfaces() == 2 &&
e.faces[0]->g != e.faces[1]->g ) e.faces(0)->g != e.faces(1)->g )
{ {
e.g = &EDGE_CLAS; e.g = &EDGE_CLAS;
} }
else if (e.numfaces() == 2) else if (e.numfaces() == 2)
{ {
BDS_Vector N0 = e.faces[0]->N(); BDS_Vector N0 = e.faces(0)->N();
BDS_Vector N1 = e.faces[1]->N(); BDS_Vector N1 = e.faces(1)->N();
double a[3] = { N0.x , N0.y , N0.z }; double a[3] = { N0.x , N0.y , N0.z };
double b[3] = { N1.x , N1.y , N1.z }; double b[3] = { N1.x , N1.y , N1.z };
double c[3]; double c[3];
...@@ -289,10 +317,15 @@ void BDS_Mesh :: classify ( double angle ) ...@@ -289,10 +317,15 @@ void BDS_Mesh :: classify ( double angle )
double cosa = a[0]*b[0] +a[1]*b[1] +a[2]*b[2]; double cosa = a[0]*b[0] +a[1]*b[1] +a[2]*b[2];
double sina = sqrt (c[0]*c[0] + c[1]*c[1] + c[2]*c[2]); double sina = sqrt (c[0]*c[0] + c[1]*c[1] + c[2]*c[2]);
double ag = atan2(sina,cosa); double ag = atan2(sina,cosa);
// printf("edge with angle = %g sin %g cos %g n1 %g %g %g n2 %g %g %g\n", // if (fabs(ag) > angle)
// ag,sina,cosa,a[0],a[1],a[2],b[0],b[1],b[2]); // printf("edge with angle = %g sin %g cos %g n1 %g %g %g n2 %g %g %g\n",
// ag,sina,cosa,a[0],a[1],a[2],b[0],b[1],b[2]);
if (fabs(ag) > angle)e.g = &EDGE_CLAS; if (fabs(ag) > angle)e.g = &EDGE_CLAS;
} }
else
{
e.g = &EDGE_CLAS;
}
++it; ++it;
} }
} }
...@@ -347,14 +380,14 @@ void BDS_Mesh :: classify ( double angle ) ...@@ -347,14 +380,14 @@ void BDS_Mesh :: classify ( double angle )
BDS_GeomEntity *g; BDS_GeomEntity *g;
if ( e.numfaces() == 1) if ( e.numfaces() == 1)
{ {
found = edgetags.find (std::make_pair ( e.faces[0]->g->classif_tag,-1)); found = edgetags.find (std::make_pair ( e.faces(0)->g->classif_tag,-1));
} }
else if (e.numfaces() == 2 && e.faces[0]->g->classif_tag != e.faces[1]->g->classif_tag) else if (e.numfaces() == 2 && e.faces(0)->g->classif_tag != e.faces(1)->g->classif_tag)
{ {
if (e.faces[0]->g->classif_tag > e.faces[1]->g->classif_tag) if (e.faces(0)->g->classif_tag > e.faces(1)->g->classif_tag)
found = edgetags.find (std::make_pair ( e.faces[1]->g->classif_tag,e.faces[0]->g->classif_tag)); found = edgetags.find (std::make_pair ( e.faces(1)->g->classif_tag,e.faces(0)->g->classif_tag));
else else
found = edgetags.find (std::make_pair ( e.faces[0]->g->classif_tag,e.faces[1]->g->classif_tag)); found = edgetags.find (std::make_pair ( e.faces(0)->g->classif_tag,e.faces(1)->g->classif_tag));
} }
if (e.g) if (e.g)
{ {
...@@ -364,15 +397,15 @@ void BDS_Mesh :: classify ( double angle ) ...@@ -364,15 +397,15 @@ void BDS_Mesh :: classify ( double angle )
g = get_geom (edgetag,1); g = get_geom (edgetag,1);
if ( e.numfaces() == 1) if ( e.numfaces() == 1)
{ {
edgetags[std::make_pair ( e.faces[0]->g->classif_tag,-1)] = edgetag; edgetags[std::make_pair ( e.faces(0)->g->classif_tag,-1)] = edgetag;
} }
else if (e.numfaces() == 2) else if (e.numfaces() == 2)
{ {
if (e.faces[0]->g->classif_tag > e.faces[1]->g->classif_tag) if (e.faces(0)->g->classif_tag > e.faces(1)->g->classif_tag)
edgetags[std::make_pair ( e.faces[1]->g->classif_tag,e.faces[0]->g->classif_tag)] edgetags[std::make_pair ( e.faces(1)->g->classif_tag,e.faces(0)->g->classif_tag)]
= edgetag; = edgetag;
else else
edgetags[std::make_pair ( e.faces[0]->g->classif_tag,e.faces[1]->g->classif_tag)] edgetags[std::make_pair ( e.faces(0)->g->classif_tag,e.faces(1)->g->classif_tag)]
= edgetag; = edgetag;
} }
edgetag++; edgetag++;
...@@ -486,6 +519,8 @@ bool BDS_Mesh :: read_stl ( const char *filename , const double tolerance) ...@@ -486,6 +519,8 @@ bool BDS_Mesh :: read_stl ( const char *filename , const double tolerance)
(Min[1]-Max[1])*(Min[1]-Max[1])+ (Min[1]-Max[1])*(Min[1]-Max[1])+
(Min[2]-Max[2])*(Min[2]-Max[2])); (Min[2]-Max[2])*(Min[2]-Max[2]));
printf("LC = %g\n",LC);
PointLessThanLexicographic::t = LC * tolerance; PointLessThanLexicographic::t = LC * tolerance;
rewind(f); rewind(f);
...@@ -884,10 +919,10 @@ bool BDS_Mesh ::split_edge ( BDS_Edge *e, double coord) ...@@ -884,10 +919,10 @@ bool BDS_Mesh ::split_edge ( BDS_Edge *e, double coord)
BDS_GeomEntity *g1=0,*g2=0,*ge=e->g; BDS_GeomEntity *g1=0,*g2=0,*ge=e->g;
BDS_Edge *p1_op1 = find_edge ( p1, op[0], e->faces[0] ); 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 *op1_p2 = find_edge ( op[0],p2,e->faces(0) );
BDS_Edge *p1_op2 = find_edge ( p1, op[1],e->faces[1] ); BDS_Edge *p1_op2 = find_edge ( p1, op[1],e->faces(1) );
BDS_Edge *op2_p2 = find_edge ( op[1],p2,e->faces[1] ); BDS_Edge *op2_p2 = find_edge ( op[1],p2,e->faces(1) );
// BDS_Edge *p1_op1 = find_edge ( p1->iD, op[0]->iD ); // BDS_Edge *p1_op1 = find_edge ( p1->iD, op[0]->iD );
// BDS_Edge *op1_p2 = find_edge ( op[0]->iD,p2->iD ); // BDS_Edge *op1_p2 = find_edge ( op[0]->iD,p2->iD );
...@@ -903,16 +938,16 @@ bool BDS_Mesh ::split_edge ( BDS_Edge *e, double coord) ...@@ -903,16 +938,16 @@ bool BDS_Mesh ::split_edge ( BDS_Edge *e, double coord)
print_face(e->faces[0]); print_face(e->faces[0]);
print_face(e->faces[1]); print_face(e->faces[1]);
*/ */
if (e->faces[0]) if (e->faces(0))
{ {
g1 = e->faces[0]->g; g1 = e->faces(0)->g;
del_triangle (e->faces[0]); del_triangle (e->faces(0));
} }
// not a bug !!! // not a bug !!!
if (e->faces[0]) if (e->faces(0))
{ {
g2 = e->faces[0]->g; g2 = e->faces(0)->g;
del_triangle (e->faces[0]); del_triangle (e->faces(0));
} }
...@@ -959,9 +994,9 @@ bool BDS_Mesh ::swap_edge ( BDS_Edge *e) ...@@ -959,9 +994,9 @@ bool BDS_Mesh ::swap_edge ( BDS_Edge *e)
if (nbFaces != 2)return false; if (nbFaces != 2)return false;
BDS_Point *pts1[3]; BDS_Point *pts1[3];
e->faces[0]->getNodes (pts1); e->faces(0)->getNodes (pts1);
BDS_Point *pts2[3]; BDS_Point *pts2[3];
e->faces[1]->getNodes (pts2); e->faces(1)->getNodes (pts2);
if (e->g && e->g->classif_degree != 2)return false; if (e->g && e->g->classif_degree != 2)return false;
...@@ -971,21 +1006,21 @@ bool BDS_Mesh ::swap_edge ( BDS_Edge *e) ...@@ -971,21 +1006,21 @@ bool BDS_Mesh ::swap_edge ( BDS_Edge *e)
e->oppositeof (op); e->oppositeof (op);
BDS_GeomEntity *g1=0,*g2=0,*ge=e->g; BDS_GeomEntity *g1=0,*g2=0,*ge=e->g;
BDS_Edge *p1_op1 = find_edge ( p1, op[0], e->faces[0] ); 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 *op1_p2 = find_edge ( op[0],p2,e->faces(0) );
BDS_Edge *p1_op2 = find_edge ( p1, op[1],e->faces[1] ); BDS_Edge *p1_op2 = find_edge ( p1, op[1],e->faces(1) );
BDS_Edge *op2_p2 = find_edge ( op[1],p2,e->faces[1] ); BDS_Edge *op2_p2 = find_edge ( op[1],p2,e->faces(1) );
if (e->faces[0]) if (e->faces(0))
{ {
g1 = e->faces[0]->g; g1 = e->faces(0)->g;
del_triangle (e->faces[0]); del_triangle (e->faces(0));
} }
// not a bug !!! // not a bug !!!
if (e->faces[0]) if (e->faces(0))
{ {
g2 = e->faces[0]->g; g2 = e->faces(0)->g;
del_triangle (e->faces[0]); del_triangle (e->faces(0));
} }
del_edge (e); del_edge (e);
...@@ -1149,10 +1184,17 @@ int BDS_Mesh :: adapt_mesh ( double l) ...@@ -1149,10 +1184,17 @@ int BDS_Mesh :: adapt_mesh ( double l)
double qb1 = quality_triangle ( (*it)->p1 , op[0] , op[1] ); double qb1 = quality_triangle ( (*it)->p1 , op[0] , op[1] );
double qb2 = quality_triangle ( (*it)->p2 , op[0] , op[1] ); double qb2 = quality_triangle ( (*it)->p2 , op[0] , op[1] );
double d = dist_droites_gauches((*it)->p1, (*it)->p2,
op[0],op[1]);
BDS_Vector v1 (*((*it)->p1), *((*it)->p2));
BDS_Vector v2 (*(op[0]),*(op[1]));
double dd = sqrt (v1*v1 + v2*v2);
double qa = (qa1<qa2)?qa1:qa2; double qa = (qa1<qa2)?qa1:qa2;
double qb = (qb1<qb2)?qb1:qb2; double qb = (qb1<qb2)?qb1:qb2;
// printf("qa %g qb %g ..\n",qa,qb); // printf("qa %g qb %g ..\n",qa,qb);
if (qb > qa) if (qb > qa && d < 0.05 * dd)
{ {
// printf("swop ..\n"); // printf("swop ..\n");
nb_modif++; nb_modif++;
......
Mesh/BDS.h
+ 64
− 62
View file @ 34164fbb
...@@ -5,8 +5,10 @@ ...@@ -5,8 +5,10 @@
#include <set> #include <set>
#include <map> #include <map>
#include <vector>
#include <list> #include <list>
#include <math.h> #include <math.h>
#include <algorithm>
class BDS_GeomEntity class BDS_GeomEntity
{ {
...@@ -32,38 +34,6 @@ class BDS_Triangle; ...@@ -32,38 +34,6 @@ class BDS_Triangle;
class BDS_Mesh; class BDS_Mesh;
void print_face (BDS_Triangle *t); void print_face (BDS_Triangle *t);
class BDS_Vector
{
public:
double x,y,z;
BDS_Vector operator + (const BDS_Vector &v)
{
return BDS_Vector (x+v.x,y+v.y,z+v.z);
}
BDS_Vector& operator += (const BDS_Vector &v)
{
x+=v.x;
y+=v.y;
z+=v.z;
return *this;
}
BDS_Vector operator / (const double &v)
{
return BDS_Vector (x/v,y/v,z/v);
}
BDS_Vector operator * (const double &v)
{
return BDS_Vector (x*v,y*v,z*v);
}
BDS_Vector (double ix=0, double iy=0, double iz=0)
// : x(ix),(iy),z(iz)
{
x=ix;
y=iy;
z=iz;
}
};
class BDS_Point class BDS_Point
{ {
public: public:
...@@ -97,22 +67,73 @@ public: ...@@ -97,22 +67,73 @@ public:
} }
}; };
class BDS_Vector
{
public:
double x,y,z;
BDS_Vector operator + (const BDS_Vector &v)
{
return BDS_Vector (x+v.x,y+v.y,z+v.z);
}
inline BDS_Vector operator % (const BDS_Vector &other) const
{
return BDS_Vector(y*other.z-z*other.y,
z*other.x-x*other.z,
x*other.y-y*other.x);
}
BDS_Vector& operator += (const BDS_Vector &v)
{
x+=v.x;
y+=v.y;
z+=v.z;
return *this;
}
BDS_Vector operator / (const double &v)
{
return BDS_Vector (x/v,y/v,z/v);
}
BDS_Vector operator * (const double &v)
{
return BDS_Vector (x*v,y*v,z*v);
}
double operator * (const BDS_Vector &v) const
{
return (x*v.x+y*v.y+z*v.z);
}
BDS_Vector (const BDS_Point &p2,const BDS_Point &p1)
: x(p2.X-p1.X),y(p2.Y-p1.Y),z(p2.Z-p1.Z)
{
}
BDS_Vector (double ix=0, double iy=0, double iz=0)
// : x(ix),(iy),z(iz)
{
x=ix;
y=iy;
z=iz;
}
};
class BDS_Edge class BDS_Edge
{ {
std::vector <BDS_Triangle *> _faces;
public: public:
bool deleted; bool deleted;
BDS_Point *p1,*p2; BDS_Point *p1,*p2;
BDS_GeomEntity *g; BDS_GeomEntity *g;
BDS_Triangle*faces[2]; inline BDS_Triangle* faces(int i) const
{
return _faces [i];
}
inline double length () const inline double length () const
{ {
return sqrt ((p1->X-p2->X)*(p1->X-p2->X)+(p1->Y-p2->Y)*(p1->Y-p2->Y)+(p1->Z-p2->Z)*(p1->Z-p2->Z)); return sqrt ((p1->X-p2->X)*(p1->X-p2->X)+(p1->Y-p2->Y)*(p1->Y-p2->Y)+(p1->Z-p2->Z)*(p1->Z-p2->Z));
} }
inline int numfaces ( ) const inline int numfaces ( ) const
{ {
if (faces[1]) return 2; return _faces.size();
if (faces[0]) return 1;
return 0;
} }
inline BDS_Point * commonvertex ( const BDS_Edge *other ) const inline BDS_Point * commonvertex ( const BDS_Edge *other ) const
{ {
...@@ -123,18 +144,7 @@ public: ...@@ -123,18 +144,7 @@ public:
inline void addface ( BDS_Triangle *f) inline void addface ( BDS_Triangle *f)
{ {
if (faces[1]) _faces.push_back(f);
{
printf("Non Manifold model not done yet\n");
print_face (f);
print_face (faces[0]);
print_face (faces[1]);
throw 1;
}
if(faces [0])faces[1] = f;
else faces[0] = f;
} }
inline bool operator < ( const BDS_Edge & other ) const inline bool operator < ( const BDS_Edge & other ) const
{ {
...@@ -145,22 +155,15 @@ public: ...@@ -145,22 +155,15 @@ public:
} }
inline BDS_Triangle * otherFace ( const BDS_Triangle *f) const inline BDS_Triangle * otherFace ( const BDS_Triangle *f) const
{ {
if (f == faces[0]) return faces[1]; if (numfaces()!=2) throw;
if (f == faces[1]) return faces[0]; if (f == _faces[0]) return _faces[1];
throw; if (f == _faces[1]) return _faces[0];
throw;
} }
inline void del (BDS_Triangle *t) inline void del (BDS_Triangle *t)
{ {
if (faces[0] == t) _faces.erase ( std::remove_if (_faces.begin(),_faces.end() , std::bind2nd(std::equal_to<BDS_Triangle*> (), t)) ,
{ _faces.end () );
faces[0] = faces[1];
}
else if (faces[1]!=t)
{
printf("edge with faces %p %p : cannot delete %p\n",faces[0],faces[1],t);
throw;
}
faces [1] = 0;
} }
inline void oppositeof (BDS_Point * oface[2]) const; inline void oppositeof (BDS_Point * oface[2]) const;
...@@ -168,7 +171,6 @@ public: ...@@ -168,7 +171,6 @@ public:
BDS_Edge ( BDS_Point *A, BDS_Point *B ) BDS_Edge ( BDS_Point *A, BDS_Point *B )
: deleted(false), g(0) : deleted(false), g(0)
{ {
faces[0] = faces[1] = 0;
if (*A < *B) if (*A < *B)
{ {
p1=A; p1=A;
......
Parser/OpenFile.cpp
+ 3
− 3
View file @ 34164fbb
// $Id: OpenFile.cpp,v 1.76 2005-04-21 09:26:11 remacle Exp $ // $Id: OpenFile.cpp,v 1.77 2005-04-22 14:36:43 remacle Exp $
// //
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
// //
...@@ -289,10 +289,10 @@ int MergeProblem(char *name, int warn_if_missing) ...@@ -289,10 +289,10 @@ int MergeProblem(char *name, int warn_if_missing)
else else
{ {
// STEREO LITOGRAPHY (STL) FILE // STEREO LITOGRAPHY (STL) FILE
THEM->bds->read_stl ( name , 1.e-6 ); THEM->bds->read_stl ( name , 5.e-7 );
} }
THEM->bds->save_gmsh_format ( "1.msh" ); THEM->bds->save_gmsh_format ( "1.msh" );
THEM->bds->classify ( M_PI / 6 ); THEM->bds->classify ( M_PI / 4 );
THEM->bds->save_gmsh_format ( "2.msh" ); THEM->bds->save_gmsh_format ( "2.msh" );
BDS_To_Mesh (THEM); BDS_To_Mesh (THEM);
SetBoundingBox(); SetBoundingBox();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment