From e5c3f30ae1384bba3e75cb5f217c871c6c5954c9 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Mon, 11 Apr 2005 08:53:15 +0000
Subject: [PATCH] *** empty log message ***
---
Mesh/BDS.cpp | 131 +++++++++++++++++++++++++++
Mesh/BDS.h | 186 +++++++++++++++++++++++++++++++++++++++
Mesh/DiscreteSurface.cpp | 35 +++++++-
Mesh/Makefile | 3 +-
4 files changed, 353 insertions(+), 2 deletions(-)
create mode 100644 Mesh/BDS.cpp
create mode 100644 Mesh/BDS.h
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
new file mode 100644
index 0000000000..47e0bbb039
--- /dev/null
+++ b/Mesh/BDS.cpp
@@ -0,0 +1,131 @@
+#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++;
+ }
+ }
+
+}
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
new file mode 100644
index 0000000000..8404642032
--- /dev/null
+++ b/Mesh/BDS.h
@@ -0,0 +1,186 @@
+// 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);
+};
diff --git a/Mesh/DiscreteSurface.cpp b/Mesh/DiscreteSurface.cpp
index 3f510bc311..346c74af64 100644
--- a/Mesh/DiscreteSurface.cpp
+++ b/Mesh/DiscreteSurface.cpp
@@ -1,4 +1,4 @@
-// $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){
diff --git a/Mesh/Makefile b/Mesh/Makefile
index e5d7ba738f..063830fc3d 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $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 \
--
GitLab