diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index bc41ce924eb93486824ef358ba8e4517815f94c2..c92d99023f400955d5e85362b59b7647d8d9b094 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -3,6 +3,55 @@
#include <math.h>
#include "Numeric.h"
+
+
+BDS_Sphere :: BDS_Sphere( BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, BDS_Point *p4)
+{
+ double mat[3][3], b[3], dum,res[3];
+ int i;
+ double X[4] = {p1->X,p2->X,p3->X,p4->X};
+ double Y[4] = {p1->Y,p2->Y,p3->Y,p4->Y};
+ double Z[4] = {p1->Z,p2->Z,p3->Z,p4->Z};
+
+ b[0] = X[1] * X[1] - X[0] * X[0] +
+ Y[1] * Y[1] - Y[0] * Y[0] + Z[1] * Z[1] - Z[0] * Z[0];
+ b[1] = X[2] * X[2] - X[1] * X[1] +
+ Y[2] * Y[2] - Y[1] * Y[1] + Z[2] * Z[2] - Z[1] * Z[1];
+ b[2] = X[3] * X[3] - X[2] * X[2] +
+ Y[3] * Y[3] - Y[2] * Y[2] + Z[3] * Z[3] - Z[2] * Z[2];
+
+ for(i = 0; i < 3; i++)
+ b[i] *= 0.5;
+
+ mat[0][0] = X[1] - X[0];
+ mat[0][1] = Y[1] - Y[0];
+ mat[0][2] = Z[1] - Z[0];
+ mat[1][0] = X[2] - X[1];
+ mat[1][1] = Y[2] - Y[1];
+ mat[1][2] = Z[2] - Z[1];
+ mat[2][0] = X[3] - X[2];
+ mat[2][1] = Y[3] - Y[2];
+ mat[2][2] = Z[3] - Z[2];
+
+
+// printf("X = %g %g %g %g\n",X[0],X[1],X[2],X[3]);
+// printf("Y = %g %g %g %g\n",Y[0],Y[1],Y[2],Y[3]);
+// printf("Z = %g %g %g %g\n",Z[0],Z[1],Z[2],Z[3]);
+
+ if(!sys3x3(mat, b, res, &dum)) {
+ R = 1.e22;
+ xc = yc = zc = 1.e33;
+// printf("ahrgh\n");
+ }
+ else
+ {
+ xc = res[0];
+ yc = res[1];
+ zc = res[2];
+ R = sqrt ((p1->X-xc)*(p1->X-xc)+(p1->Y-yc)*(p1->Y-yc)+(p1->Z-zc)*(p1->Z-zc));
+ }
+}
+
double dist_droites_gauches(BDS_Point *p1, BDS_Point *p2,
BDS_Point *p3 ,BDS_Point *p4)
{
@@ -111,6 +160,20 @@ double quality_triangle (BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
return a / (e12+e22+e32);
}
+BDS_Plane::BDS_Plane (BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
+{
+ double mat[3][3];
+ mat[0][0] = p1->X; mat[0][1] = p1->Y; mat[0][2] = p1->Z;
+ mat[1][0] = p2->X; mat[1][1] = p2->Y; mat[1][2] = p2->Z;
+ mat[2][0] = p3->X; mat[2][1] = p3->Y; mat[2][2] = p3->Z;
+ double rhs[3] = {-1,-1,-1};
+ double pl[3] , det ;
+ sys3x3 (mat,rhs,pl,&det);
+ a = pl[0];
+ b = pl[1];
+ c = pl[2];
+}
+
void BDS_Point::getTriangles (std::list<BDS_Triangle *> &t)
{
@@ -295,6 +358,95 @@ BDS_Vector BDS_Triangle :: N () const
}
+void BDS_Mesh :: reverseEngineerCAD ( )
+{
+ for (std::set<BDS_GeomEntity*,GeomLessThan>::iterator it = geom.begin();
+ it != geom.end();
+ ++it)
+ {
+ if ((*it)->classif_degree == 2)
+ {
+ std::set<BDS_Point*,PointLessThan> pts;
+
+ std::list<BDS_Triangle*>::iterator tit = triangles.begin();
+ std::list<BDS_Triangle*>::iterator tite = triangles.end();
+ while (tit!=tite)
+ {
+ if ((*tit)->g == (*it))
+ {
+ BDS_Point *nod[3];
+ (*tit)->getNodes (nod);
+ pts.insert(nod[0]);
+ pts.insert(nod[1]);
+ pts.insert(nod[2]);
+ }
+ tit++;
+ }
+ if(pts.size() > 3)
+ {
+ // TEST PLANE
+ {
+ std::set<BDS_Point*>::iterator pit = pts.begin();
+ std::set<BDS_Point*>::iterator pite = pts.end();
+
+ BDS_Point *p1 = *pit;++pit;
+ BDS_Point *p2 = *pit;++pit;
+ BDS_Point *p3 = *pit;++pit;
+ BDS_Plane *plane = new BDS_Plane ( p1 , p2 , p3 );
+ double distmax = 0;
+ while (pit!=pite)
+ {
+ double dist = plane->signedDistanceTo (*pit);
+ if (fabs(dist) > distmax) distmax = fabs(dist);
+ ++pit;
+ }
+ if (distmax < 1.e-6 * LC )
+ {
+ (*it)->surf = plane;
+ printf("plane surface found %g\n",distmax);
+ }
+ else
+ {
+// printf("non plane surface found %g\n",distmax);
+ delete plane;
+ }
+ }
+ }
+ if(pts.size() > 4 && !(*it)->surf)
+ {
+ // TEST SPHERE
+ {
+ std::set<BDS_Point*,PointLessThan>::iterator pit = pts.begin();
+ std::set<BDS_Point*,PointLessThan>::iterator pite = pts.end();
+
+ BDS_Point *p1 = *pit;++pit;
+ BDS_Point *p2 = *pit;++pit;
+ BDS_Point *p3 = *pit;++pit;
+ BDS_Point *p4 = *pit;++pit;
+ BDS_Sphere *sphere = new BDS_Sphere ( p1 , p2 , p3 , p4);
+ double distmax = 0;
+ while (pit!=pite)
+ {
+ double dist = sphere->signedDistanceTo (*pit);
+ if (fabs(dist) > distmax) distmax = fabs(dist);
+ ++pit;
+ }
+ if (distmax < 1.e-2 * LC )
+ {
+ (*it)->surf = sphere;
+ printf("spherical surface %d found %g (R=%g,LC=%g)\n",(*it)->classif_tag,distmax,sphere->R,LC);
+ }
+ else
+ {
+// printf("non spherical surface found %g\n",distmax);
+ delete sphere;
+ }
+ }
+ }
+ }
+ }
+}
+
void BDS_Mesh :: classify ( double angle )
{
@@ -470,6 +622,7 @@ void BDS_Mesh :: classify ( double angle )
++it;
}
}
+ reverseEngineerCAD ( ) ;
printf(" end classifying %d edgetags %d vertextags\n",edgetag-1,vertextag-1);
}
@@ -1253,50 +1406,58 @@ bool BDS_Mesh ::smooth_point ( BDS_Point *p , BDS_Mesh *geom_mesh )
Y /= p->edges.size();
Z /= p->edges.size();
- std::list<BDS_Triangle *> t;
- std::list<BDS_Triangle *>::iterator it ;
- std::list<BDS_Triangle *>::iterator ite;
- if (geom_mesh)
+ if (p->g->surf)
{
- it = geom_mesh->triangles.begin();
- ite = geom_mesh->triangles.end();
+ p->g->surf->projection ( X,Y,Z,p->X,p->Y,p->Z);
+// printf("coucou\n");
}
else
{
- p->getTriangles (t);
- it = t.begin();
- ite = t.end();
- }
-
- double dist = 1.e22;
- double XX = X,YY=Y,ZZ=Z;
-
- int p_classif = p->g->classif_tag;
-
- while (it != ite)
- {
- if ((*it)->g->classif_tag == p_classif)
+ std::list<BDS_Triangle *> t;
+ std::list<BDS_Triangle *>::iterator it ;
+ std::list<BDS_Triangle *>::iterator ite;
+ if (geom_mesh)
{
- BDS_Point *pts[3];
- (*it)->getNodes (pts);
- double xp,yp,zp;
- proj_point_plane ( X,Y,Z,pts[0],pts[1],pts[2],xp,yp,zp);
- double d = sqrt ((xp-X)*(xp-X)+(yp-Y)*(yp-Y)+(zp-Z)*(zp-Z));
- if (d < dist)
+ it = geom_mesh->triangles.begin();
+ ite = geom_mesh->triangles.end();
+ }
+ else
+ {
+ p->getTriangles (t);
+ it = t.begin();
+ ite = t.end();
+ }
+
+ double dist = 1.e22;
+ double XX = X,YY=Y,ZZ=Z;
+
+ int p_classif = p->g->classif_tag;
+
+ while (it != ite)
+ {
+ if ((*it)->g->classif_tag == p_classif)
{
- XX = xp; YY = yp; ZZ = zp;
- dist = d;
+ BDS_Point *pts[3];
+ (*it)->getNodes (pts);
+ double xp,yp,zp;
+ proj_point_plane ( X,Y,Z,pts[0],pts[1],pts[2],xp,yp,zp);
+ double d = sqrt ((xp-X)*(xp-X)+(yp-Y)*(yp-Y)+(zp-Z)*(zp-Z));
+ if (d < dist)
+ {
+ XX = xp; YY = yp; ZZ = zp;
+ dist = d;
+ }
}
+ ++it;
}
- ++it;
+ X = XX;
+ Y = YY;
+ Z = ZZ;
+
+ p->X = X;
+ p->Y = Y;
+ p->Z = Z;
}
- X = XX;
- Y = YY;
- Z = ZZ;
-
- p->X = X;
- p->Y = Y;
- p->Z = Z;
return true;
}
@@ -1448,6 +1609,8 @@ BDS_Mesh::BDS_Mesh (const BDS_Mesh &other)
++it)
{
add_geom((*it)->classif_tag,(*it)->classif_degree);
+ BDS_GeomEntity* g = get_geom((*it)->classif_tag,(*it)->classif_degree);
+ g->surf = (*it)->surf;
}
for (std::set<BDS_Point*,PointLessThan>::iterator it = other.points.begin();
it != other.points.end();
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 0848601a90d569d62272e32ec68597405d577f32..0e321c155c4297eb6185d6f92e68a43c97aca671 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -10,12 +10,28 @@
#include <math.h>
#include <algorithm>
+class BDS_Edge;
+class BDS_Triangle;
+class BDS_Mesh;
+class BDS_Point;
+
+class BDS_Surface
+{
+public :
+ virtual double signedDistanceTo ( BDS_Point * ) const = 0;
+ virtual void projection ( double xa, double ya, double za,
+ double &x, double &y, double &z) const =0;
+};
+
+
class BDS_GeomEntity
{
public:
int classif_tag;
int classif_degree;
- bool is_plane_surface;
+
+ BDS_Surface *surf;
+
inline bool operator < ( const BDS_GeomEntity & other ) const
{
if (classif_degree < other.classif_degree)return true;
@@ -24,15 +40,12 @@ public:
return false;
}
BDS_GeomEntity (int a, int b)
- : classif_tag (a),classif_degree(b)
+ : classif_tag (a),classif_degree(b),surf(0)
{
}
};
-class BDS_Edge;
-class BDS_Triangle;
-class BDS_Mesh;
void print_face (BDS_Triangle *t);
class BDS_Point
@@ -216,6 +229,45 @@ public:
e3->addface(this);
}
};
+
+
+class BDS_Plane : public BDS_Surface
+{
+ double a,b,c;
+public :
+ BDS_Plane (BDS_Point *p1, BDS_Point *p2, BDS_Point *p3);
+ virtual double signedDistanceTo ( BDS_Point * p) const {return a*p->X + b*p->Y + c*p->Z + 1;}
+ virtual void projection ( double xa, double ya, double za,
+ double &x, double &y, double &z) const
+ {
+ double k = - ( a * xa + b * ya + c * za + 1 ) / ( a * a + b * b + c * c );
+ x = xa + k * a;
+ y = ya + k * b;
+ z = za + k * c;
+ }
+};
+
+class BDS_Sphere : public BDS_Surface
+{
+public :
+ double xc,yc,zc,R;
+ BDS_Sphere (BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, BDS_Point *p4);
+ virtual double signedDistanceTo ( BDS_Point * p) const {
+ return sqrt((p->X-xc)*(p->X-xc)+
+ (p->Y-yc)*(p->Y-yc)+
+ (p->Z-zc)*(p->Z-zc)) - R;
+ }
+ virtual void projection ( double xa, double ya, double za,
+ double &x, double &y, double &z) const
+ {
+ double k = ( R * R )/((xa-xc)*(xa-xc)+(ya-yc)*(ya-yc)+(za-zc)*(za-zc)) ;
+ x = xc + k * (xa-xc);
+ y = yc + k * (ya-yc);
+ z = zc + k * (za-zc);
+ }
+};
+
+
class GeomLessThan
{
public:
@@ -287,6 +339,7 @@ class BDS_Mesh
bool smooth_point ( BDS_Point* , BDS_Mesh *geom = 0);
bool split_edge ( BDS_Edge *, double coord);
void classify ( double angle);
+ void reverseEngineerCAD ( ) ;
int adapt_mesh(double,bool smooth = false,BDS_Mesh *geom = 0);
void cleanup();
// io's
diff --git a/Mesh/DiscreteSurface.cpp b/Mesh/DiscreteSurface.cpp
index 30c8afed806ee6c39fb25ea14b48821e1bce7a30..ba20185f5691978e52bf26090a7bd4a946d28aee 100644
--- a/Mesh/DiscreteSurface.cpp
+++ b/Mesh/DiscreteSurface.cpp
@@ -1,4 +1,4 @@
-// $Id: DiscreteSurface.cpp,v 1.11 2005-05-04 14:42:22 remacle Exp $
+// $Id: DiscreteSurface.cpp,v 1.12 2005-05-06 16:06:42 remacle Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -501,20 +501,11 @@ int MeshDiscreteSurface(Surface *s)
{
THEM->bds_mesh = new BDS_Mesh (*(THEM->bds));
int iter = 0;
- while (iter < 20 && THEM->bds_mesh -> adapt_mesh ( THEM->bds->LC / 170, true))
+ while (iter < 20 && THEM->bds_mesh -> adapt_mesh ( THEM->bds->LC / 100, true))
{
printf("iter %d done\n",iter);
iter ++;
}
- printf("smoothing 1/4\n");
- THEM->bds_mesh -> adapt_mesh ( THEM->bds->LC / 170, true,THEM->bds);
- printf("smoothing 2/4 \n");
- THEM->bds_mesh -> adapt_mesh ( THEM->bds->LC / 170, true,THEM->bds);
- printf("smoothing 3/4 \n");
- THEM->bds_mesh -> adapt_mesh ( THEM->bds->LC / 170, true,THEM->bds);
- printf("smoothing 4/4\n");
- THEM->bds_mesh -> adapt_mesh ( THEM->bds->LC / 170, true,THEM->bds);
- printf("smoothing done \n");
THEM->bds_mesh->save_gmsh_format ( "3.msh" );
}
return 1;