diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 3b40069c012dab3bdf257c122e4c4e12acdce7ed..acc06e8bbec1c98bf66c2d173b5e2049a17898a9 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -3,6 +3,28 @@
#include <math.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,
BDS_Point *p1, BDS_Point *p2, BDS_Point *p3 ,
double &x, double &y, double &z)
@@ -91,7 +113,7 @@ void BDS_Point::getTriangles (std::list<BDS_Triangle *> &t)
int NF = (*it)->numfaces();
for (int i=0;i<NF;++i)
{
- BDS_Triangle *tt = (*it)->faces[i];
+ BDS_Triangle *tt = (*it)->faces(i);
if (tt)
{
std::list<BDS_Triangle*>::iterator tit = t.begin();
@@ -198,22 +220,23 @@ void BDS_Mesh :: add_geom (int p1, int p2)
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];
- faces[0]->getNodes (pts);
- if (pts[0] != p1 && pts[0] != p2)oface[0] = pts[0];
- else if (pts[1] != p1 && pts[1] != p2)oface[0] = pts[1];
- else oface[0] = pts[2];
+ BDS_Point *pts[3];
+ faces(0)->getNodes (pts);
+ if (pts[0] != p1 && pts[0] != p2)oface[0] = pts[0];
+ else if (pts[1] != p1 && pts[1] != p2)oface[0] = pts[1];
+ else oface[0] = pts[2];
}
- if (faces[1])
+ if (faces(1))
{
- BDS_Point *pts[3];
- faces[1]->getNodes (pts);
- if (pts[0] != p1 && pts[0] != p2)oface[1] = pts[0];
- else if (pts[1] != p1 && pts[1] != p2)oface[1] = pts[1];
- else oface[1] = pts[2];
+ BDS_Point *pts[3];
+ faces(1)->getNodes (pts);
+ if (pts[0] != p1 && pts[0] != p2)oface[1] = pts[0];
+ else if (pts[1] != p1 && pts[1] != p2)oface[1] = pts[1];
+ else oface[1] = pts[2];
}
}
@@ -270,16 +293,21 @@ void BDS_Mesh :: classify ( double angle )
while (it!=ite)
{
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 &&
- e.faces[0]->g != e.faces[1]->g )
+ e.faces(0)->g != e.faces(1)->g )
{
e.g = &EDGE_CLAS;
}
+
else if (e.numfaces() == 2)
{
- BDS_Vector N0 = e.faces[0]->N();
- BDS_Vector N1 = e.faces[1]->N();
+ BDS_Vector N0 = e.faces(0)->N();
+ BDS_Vector N1 = e.faces(1)->N();
double a[3] = { N0.x , N0.y , N0.z };
double b[3] = { N1.x , N1.y , N1.z };
double c[3];
@@ -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 sina = sqrt (c[0]*c[0] + c[1]*c[1] + c[2]*c[2]);
double ag = atan2(sina,cosa);
-// 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)
+// 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;
}
+ else
+ {
+ e.g = &EDGE_CLAS;
+ }
++it;
}
}
@@ -347,14 +380,14 @@ void BDS_Mesh :: classify ( double angle )
BDS_GeomEntity *g;
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)
- found = edgetags.find (std::make_pair ( e.faces[1]->g->classif_tag,e.faces[0]->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));
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)
{
@@ -364,15 +397,15 @@ void BDS_Mesh :: classify ( double angle )
g = get_geom (edgetag,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)
{
- 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)]
+ 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)]
= edgetag;
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++;
@@ -486,6 +519,8 @@ bool BDS_Mesh :: read_stl ( const char *filename , const double tolerance)
(Min[1]-Max[1])*(Min[1]-Max[1])+
(Min[2]-Max[2])*(Min[2]-Max[2]));
+ printf("LC = %g\n",LC);
+
PointLessThanLexicographic::t = LC * tolerance;
rewind(f);
@@ -884,10 +919,10 @@ bool BDS_Mesh ::split_edge ( BDS_Edge *e, double coord)
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] );
+ 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) );
// BDS_Edge *p1_op1 = find_edge ( p1->iD, op[0]->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)
print_face(e->faces[0]);
print_face(e->faces[1]);
*/
- if (e->faces[0])
+ if (e->faces(0))
{
- g1 = e->faces[0]->g;
- del_triangle (e->faces[0]);
+ g1 = e->faces(0)->g;
+ del_triangle (e->faces(0));
}
// not a bug !!!
- if (e->faces[0])
+ if (e->faces(0))
{
- g2 = e->faces[0]->g;
- del_triangle (e->faces[0]);
+ g2 = e->faces(0)->g;
+ del_triangle (e->faces(0));
}
@@ -959,9 +994,9 @@ bool BDS_Mesh ::swap_edge ( BDS_Edge *e)
if (nbFaces != 2)return false;
BDS_Point *pts1[3];
- e->faces[0]->getNodes (pts1);
+ e->faces(0)->getNodes (pts1);
BDS_Point *pts2[3];
- e->faces[1]->getNodes (pts2);
+ e->faces(1)->getNodes (pts2);
if (e->g && e->g->classif_degree != 2)return false;
@@ -971,21 +1006,21 @@ bool BDS_Mesh ::swap_edge ( BDS_Edge *e)
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] );
+ 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])
+ if (e->faces(0))
{
- g1 = e->faces[0]->g;
- del_triangle (e->faces[0]);
+ g1 = e->faces(0)->g;
+ del_triangle (e->faces(0));
}
// not a bug !!!
- if (e->faces[0])
+ if (e->faces(0))
{
- g2 = e->faces[0]->g;
- del_triangle (e->faces[0]);
+ g2 = e->faces(0)->g;
+ del_triangle (e->faces(0));
}
del_edge (e);
@@ -1149,10 +1184,17 @@ int BDS_Mesh :: adapt_mesh ( double l)
double qb1 = quality_triangle ( (*it)->p1 , 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 qb = (qb1<qb2)?qb1:qb2;
// printf("qa %g qb %g ..\n",qa,qb);
- if (qb > qa)
+ if (qb > qa && d < 0.05 * dd)
{
// printf("swop ..\n");
nb_modif++;
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 4ef115473a06dd1609d436d3c9cf7476ffec0272..aa43cdf01d3d873a3611ff091d4449e308a97a6b 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -5,8 +5,10 @@
#include <set>
#include <map>
+#include <vector>
#include <list>
#include <math.h>
+#include <algorithm>
class BDS_GeomEntity
{
@@ -32,38 +34,6 @@ class BDS_Triangle;
class BDS_Mesh;
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
{
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
{
+ std::vector <BDS_Triangle *> _faces;
public:
bool deleted;
BDS_Point *p1,*p2;
BDS_GeomEntity *g;
- BDS_Triangle*faces[2];
+ inline BDS_Triangle* faces(int i) const
+ {
+ return _faces [i];
+ }
+
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));
}
inline int numfaces ( ) const
{
- if (faces[1]) return 2;
- if (faces[0]) return 1;
- return 0;
+ return _faces.size();
}
inline BDS_Point * commonvertex ( const BDS_Edge *other ) const
{
@@ -123,18 +144,7 @@ public:
inline void addface ( BDS_Triangle *f)
{
- if (faces[1])
- {
- 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;
+ _faces.push_back(f);
}
inline bool operator < ( const BDS_Edge & other ) const
{
@@ -145,22 +155,15 @@ public:
}
inline BDS_Triangle * otherFace ( const BDS_Triangle *f) const
{
- if (f == faces[0]) return faces[1];
- if (f == faces[1]) return faces[0];
- throw;
+ if (numfaces()!=2) throw;
+ if (f == _faces[0]) return _faces[1];
+ if (f == _faces[1]) return _faces[0];
+ throw;
}
inline void del (BDS_Triangle *t)
{
- if (faces[0] == t)
- {
- 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;
+ _faces.erase ( std::remove_if (_faces.begin(),_faces.end() , std::bind2nd(std::equal_to<BDS_Triangle*> (), t)) ,
+ _faces.end () );
}
inline void oppositeof (BDS_Point * oface[2]) const;
@@ -168,7 +171,6 @@ public:
BDS_Edge ( BDS_Point *A, BDS_Point *B )
: deleted(false), g(0)
{
- faces[0] = faces[1] = 0;
if (*A < *B)
{
p1=A;
diff --git a/Parser/OpenFile.cpp b/Parser/OpenFile.cpp
index c7d1083e09b7834f63da7a6e710571c4b60a7d63..48f073601981c75f7e8a986c5394e1f7db13ac04 100644
--- a/Parser/OpenFile.cpp
+++ b/Parser/OpenFile.cpp
@@ -1,4 +1,4 @@
-// $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
//
@@ -289,10 +289,10 @@ int MergeProblem(char *name, int warn_if_missing)
else
{
// 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->classify ( M_PI / 6 );
+ THEM->bds->classify ( M_PI / 4 );
THEM->bds->save_gmsh_format ( "2.msh" );
BDS_To_Mesh (THEM);
SetBoundingBox();