*** empty log message ***

Mesh/BDS.cpp

+#include "BDS.h"
+#include <math.h>
+
+void BDS_Mesh :: add_point (int num , double x, double y,double z)
+{
+    points.insert ( BDS_Point ( num, x, y, z) );
+}
+
+BDS_Point *BDS_Mesh :: find_point  (int p)
+{
+    BDS_Point P ( p ) ;
+    std::set<BDS_Point>::iterator it = points.find(P);
+    if (it != points.end()) return (BDS_Point*) & (*it);
+    else return 0;
+}
+
+BDS_Edge *BDS_Mesh :: find_edge  (int num1, int num2)
+{
+    BDS_Point P1 ( num1 ) ;
+    BDS_Point P2 ( num2 ) ;
+    BDS_Edge  E (&P1, &P2);
+    std::set<BDS_Edge>::iterator it = edges.find(E);
+    if (it != edges.end()) return (BDS_Edge*) & (*it);
+    else return 0;
+}
+
+void BDS_Mesh :: add_edge  (int p1, int p2)
+{
+    BDS_Edge *efound = find_edge (p1,p2);
+    if (efound)return;
+
+    BDS_Point *pp1=find_point(p1);
+    BDS_Point *pp2=find_point(p2);
+    if (!pp1 || !pp2)throw;
+    edges.insert ( BDS_Edge ( pp1, pp2 ) );    
+}
+
+void BDS_Mesh :: add_triangle  (int p1, int p2, int p3, double nx, double ny, double nz)
+{
+    add_edge (p1,p2);
+    add_edge (p2,p3);
+    add_edge (p3,p1);
+    BDS_Edge *e1 = find_edge (p1,p2); 
+    BDS_Edge *e2 = find_edge (p2,p3); 
+    BDS_Edge *e3 = find_edge (p3,p1); 
+    triangles.insert ( BDS_Triangle ( e1, e2, e3 ,nx,ny,nz) );        
+}
+
+void BDS_Mesh :: add_geom  (int p1, int p2)
+{
+    geom.insert ( BDS_GeomEntity ( p1,p2 ) ) ;        
+}
+BDS_GeomEntity * BDS_Mesh :: get_geom  (int p1, int p2)
+{
+  std::set<BDS_GeomEntity>::iterator it = geom.find(BDS_GeomEntity ( p1,p2 ));
+  if (it == geom.end())return 0;
+  return (BDS_GeomEntity*) & (*it);
+}
+
+void BDS_Mesh :: save_gmsh_format ( const char *filename )
+{
+}
+
+void recur_tag ( BDS_Triangle *t , BDS_GeomEntity *g )
+{
+  if (!t->g)
+    {
+      t->g = g;
+      if ( ! t->e1->g && t->e1->numfaces() == 2)
+	{
+	  recur_tag (t->e1->otherFace (t) , g);
+	}
+      if ( ! t->e2->g && t->e2->numfaces() == 2)
+	{
+	  recur_tag (t->e2->otherFace (t) , g);
+	}
+      if ( ! t->e3->g && t->e3->numfaces() == 2)
+	{
+	  recur_tag (t->e3->otherFace (t) , g);
+	}
+    }
+}
+
+void BDS_Mesh :: classify ( double angle )
+{
+    BDS_GeomEntity EDGE_CLAS (1,0);
+    std::set<BDS_Edge>::iterator it  = edges.begin();
+    std::set<BDS_Edge>::iterator ite = edges.end();
+    while (it!=ite)
+    {
+	BDS_Edge &e = (BDS_Edge &) *it;
+	if ( e.numfaces() == 1) e.g = &EDGE_CLAS;
+	else if (e.numfaces() == 2)
+	  {
+	    double a[3] = { e.faces[0]->NX ,  e.faces[0]->NY ,  e.faces[0]->NZ };
+	    double b[3] = { e.faces[1]->NX ,  e.faces[1]->NY ,  e.faces[1]->NZ };
+	    double c[3];
+	    c[2] = a[0] * b[1] - a[1] * b[0];
+	    c[1] = -a[0] * b[2] + a[2] * b[0];
+	    c[0] = a[1] * b[2] - a[2] * b[1];
+	    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);
+	    if (fabs(ag) > angle)e.g = &EDGE_CLAS;
+	  }
+	++it;
+    }
+    {
+      int tag = 1;
+      while (1)
+	{
+	  std::set<BDS_Triangle>::iterator it  = triangles.begin();
+	  std::set<BDS_Triangle>::iterator ite = triangles.end();
+	  BDS_Triangle *start = 0;
+	  while (it!=ite)
+	    {
+	      if (!it->g)
+		{
+		  start = (BDS_Triangle*)& (*it);
+		}
+	      ++it;
+	    }
+	  if (!start)break;	  
+	  add_geom (tag, 2);
+	  BDS_GeomEntity *g = get_geom (tag,2);
+	  recur_tag ( start , g );
+	  tag++;
+	}
+    }
+
+}

Mesh/BDS.h

+// This is a 2D version of the Bidirectional Data Structure (BDS)
+// of shephard and beall
+// points may know the normals to the surface they are classified on
+// default values are 0,0,1
+
+#include <set>
+#include <list>
+
+class BDS_GeomEntity
+{
+public:
+    int classif_degree;
+    int classif_tag;
+    inline bool operator <  ( const BDS_GeomEntity & other ) const
+	{
+	    if (classif_degree < other.classif_degree)return true;
+	    if (classif_degree > other.classif_degree)return false;
+	    if (classif_tag < other.classif_tag)return true;
+	    return false;
+	}
+    BDS_GeomEntity (int a, int b)
+	: classif_degree (a),classif_tag(b)
+	{
+	}
+};
+
+
+class BDS_Edge;
+class BDS_Triangle;
+
+class BDS_Point
+{
+public:
+  int iD;
+  BDS_GeomEntity *g;
+  double X,Y,Z;
+  std::list<BDS_Edge*> edges;
+  inline bool operator <  ( const BDS_Point & other ) const
+  {
+    return iD < other.iD;
+  }
+  BDS_Point ( int id, double x=0, double y=0, double z=0 )
+    : iD(id),X(x),Y(y),Z(z),g(0)
+  {	    
+  }
+};
+
+class BDS_Edge
+{
+public:
+    BDS_Point *p1,*p2;
+    BDS_GeomEntity *g;
+    BDS_Triangle*faces[2];
+    inline int numfaces ( ) const 
+	{
+	    if (faces[1]) return 2;
+	    if (faces[0]) return 1;
+	    return 0;
+	}
+    inline void addface ( BDS_Triangle *f)
+	{
+	    if(faces [0])faces[1] = f;
+	    else faces[0] = f;
+	}
+    inline bool operator <  ( const BDS_Edge & other ) const
+	{
+	    if (*other.p1 < *p1) return true;
+	    if (*p1 < *other.p1) return false;
+	    if (*other.p2 < *p2) return true;
+	    return false;
+	}
+    inline BDS_Triangle * otherFace ( BDS_Triangle *f)
+      {
+	if (f == faces[0]) return faces[1];
+	if (f == faces[1]) return faces[0];
+	throw;
+      }
+    BDS_Edge ( BDS_Point *A, BDS_Point *B )
+	: g(0)
+	{	    
+	    faces[0] = faces[1] = 0;
+	    if (*A < *B) 
+	    {
+		p1=A;
+		p2=B;
+	    }
+	    else
+	    {
+		p1=B;
+		p2=A;
+	    }
+	    p1->edges.push_back(this);
+	    p2->edges.push_back(this);
+	}
+};
+
+class BDS_Triangle
+{
+public:
+    BDS_Edge *e1,*e2,*e3;
+    BDS_GeomEntity *g;
+    double NX,NY,NZ;
+    inline bool operator <  ( const BDS_Triangle & other ) const
+	{
+	    if (*other.e1 < *e1) return true;
+	    if (*e1 < *other.e1) return false;
+	    if (*other.e2 < *e2) return true;
+	    if (*e2 < *other.e2) return false;
+	    if (*other.e3 < *e3) return true;
+	    return false;
+	}
+    BDS_Triangle ( BDS_Edge *A, BDS_Edge *B, BDS_Edge *C , double nx=0, double ny=0, double nz=0)
+	: NX(nx),NY(ny),NZ(nz),g(0)
+	{	    
+	    if (*A < *B && *A < *C) 
+	    {
+		if (*B < *C) 
+		{
+		    e1=A;
+		    e2=B;
+		    e3=C;
+		}
+		else
+		{
+		    e1=A;
+		    e2=C;
+		    e3=B;
+		}
+	    }
+	    else if (*B < *A && *B < *C) 
+	    {
+		if (*A < *C) 
+		{
+		    e1=B;
+		    e2=A;
+		    e3=C;
+		}
+		else
+		{
+		    e1=B;
+		    e2=C;
+		    e3=A;
+		}
+	    }
+	    else
+	    {
+		if (*A < *B) 
+		{
+		    e1=C;
+		    e2=A;
+		    e3=B;
+		}
+		else
+		{
+		    e1=C;
+		    e2=B;
+		    e3=A;
+		}
+	    }
+	    e1->addface(this);
+	    e2->addface(this);
+	    e3->addface(this);
+	}
+};
+
+class BDS_Mesh 
+{
+ public:
+    std::set<BDS_GeomEntity> geom; 
+    std::set<BDS_Point>      points; 
+    std::set<BDS_Edge>       edges; 
+    std::set<BDS_Triangle>   triangles; 
+    void add_point (int num , double x, double y,double z);
+    void add_edge  (int p1, int p2);
+    void add_triangle  (int p1, int p2, int p3);
+    void add_triangle  (int p1, int p2, int p3, double nx, double ny, double nz);
+    void add_geom  (int degree, int tag);
+    BDS_Point *find_point (int num);
+    BDS_Edge  *find_edge (int p1, int p2);
+  BDS_GeomEntity *get_geom  (int p1, int p2);
+//    bool swap_edge ( BDS_Edge *);
+//    bool collapse_edge ( BDS_Edge *);
+//    bool split_edge ( BDS_Edge *, double coord);
+    void save_gmsh_format (const char *filename);
+    void classify ( double angle);
+};

Mesh/DiscreteSurface.cpp

-// $Id: DiscreteSurface.cpp,v 1.6 2005-03-14 18:12:29 geuzaine Exp $
+// $Id: DiscreteSurface.cpp,v 1.7 2005-04-11 08:53:15 remacle Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -32,6 +32,7 @@
 #include "Create.h"
 #include "Interpolation.h"
 #include "Context.h"
+#include "BDS.h"
 
 extern Mesh *THEM;
 extern Context_T CTX;
@@ -137,6 +138,35 @@ double SetLC(Vertex *v1, Vertex *v2, Vertex *v3, double factor)
   return lc;
 }
 
+void Mesh_To_BDS(Surface *s, BDS_Mesh *m)
+{
+  List_T *vertices = Tree2List ( s->Vertices ) ;
+  for (int i=0;i<List_Nbr ( vertices ) ;++i)
+    {
+      Vertex *v;
+      List_Read ( vertices, i, &v);
+      m->add_point (v->Num,v->Pos.X,v->Pos.Y,v->Pos.Z); 
+    }
+  List_Delete (vertices);
+  List_T *triangles = Tree2List ( s->Simplexes) ;
+  for (int i=0;i<List_Nbr ( triangles ) ;++i)
+    {
+      Simplex *simp;
+      List_Read ( triangles, i, &simp);
+      Vertex *v1 = simp->V[0];
+      Vertex *v2 = simp->V[1];
+      Vertex *v3 = simp->V[2];
+      double n[3];
+      normal3points ( v1->Pos.X , v1->Pos.Y , v1->Pos.Z,
+		      v2->Pos.X , v2->Pos.Y , v2->Pos.Z,
+		      v3->Pos.X , v3->Pos.Y , v3->Pos.Z,
+		      n);
+      m->add_triangle (v1->Num,v2->Num,v3->Num,n[0],n[1],n[2]); 
+    }
+  List_Delete (triangles);
+
+}
+
 void POLY_rep_To_Mesh(POLY_rep *prep, Surface *s)
 {  
   VertexBound = Tree_Create(sizeof(Vertex *), comparePosition);
@@ -579,6 +609,9 @@ int MeshDiscreteSurface(Surface *s)
     // routines to do that at the moment--so let's just use it and
     // hope for the best.
     POLY_rep_To_Mesh(s->thePolyRep, s);
+    BDS_Mesh bds;
+    Mesh_To_BDS(s,&bds);
+    bds.classify ( M_PI / 9 );
     return 1;
   }
   else if(s->Typ == MSH_SURF_DISCRETE){

Mesh/Makefile

-# $Id: Makefile,v 1.79 2005-02-20 06:36:54 geuzaine Exp $
+# $Id: Makefile,v 1.80 2005-04-11 08:53:15 remacle Exp $
 #
 # Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 #
@@ -53,6 +53,7 @@ SRC = 1D_Mesh.cpp \
         3D_Divide.cpp \
         3D_Bricks.cpp \
         3D_Mesh_Netgen.cpp \
+      BDS.cpp \
       MeshQuality.cpp \
       Create.cpp \
       Generator.cpp \