From c6376f6c46a171978971a7e68c2f0187b32ff12b Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 15 Apr 2005 14:32:40 +0000
Subject: [PATCH] *** empty log message ***
---
Mesh/2D_Mesh.cpp | 11 +-
Mesh/2D_Recombine.cpp | 108 +++++-
Mesh/BDS.cpp | 784 ++++++++++++++++++++++++++++++++++++---
Mesh/BDS.h | 259 +++++++++----
Mesh/DiscreteSurface.cpp | 23 +-
Mesh/Makefile | 64 ++--
6 files changed, 1084 insertions(+), 165 deletions(-)
diff --git a/Mesh/2D_Mesh.cpp b/Mesh/2D_Mesh.cpp
index adb0163777..125ea29dea 100644
--- a/Mesh/2D_Mesh.cpp
+++ b/Mesh/2D_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: 2D_Mesh.cpp,v 1.73 2005-01-20 19:05:09 geuzaine Exp $
+// $Id: 2D_Mesh.cpp,v 1.74 2005-04-15 14:32:40 remacle Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -914,7 +914,7 @@ void Maillage_Surface(void *data, void *dum)
delete_NXE(tnxe);
}
- if(s->Recombine){
+ if(s->Recombine && CTX.mesh.algo_recombine ==1 ){
Msg(STATUS3, "Recombining surface %d", s->Num);
Recombine(s->Vertices, s->Simplexes, s->Quadrangles, s->RecombineAngle);
}
@@ -947,4 +947,11 @@ void Maillage_Surface(void *data, void *dum)
ReOrientSurfaceMesh(s);
}
+
+ if(CTX.mesh.algo_recombine ==2 )
+ {
+ printf("coucou\n");
+ Recombine_All (THEM);
+ }
+
}
diff --git a/Mesh/2D_Recombine.cpp b/Mesh/2D_Recombine.cpp
index 313499018e..753833d63c 100644
--- a/Mesh/2D_Recombine.cpp
+++ b/Mesh/2D_Recombine.cpp
@@ -1,4 +1,4 @@
-// $Id: 2D_Recombine.cpp,v 1.22 2005-01-01 19:35:30 geuzaine Exp $
+// $Id: 2D_Recombine.cpp,v 1.23 2005-04-15 14:32:40 remacle Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -166,3 +166,109 @@ int Recombine(Tree_T * Vertices, Tree_T * Simplexes, Tree_T * Quadrangles,
return ntot;
}
+
+/*
+ This function recombines EVERYTHING.
+
+ First, all characteristic lengths are divided by 2.
+
+ This is done in 3 steps.
+ -) the previous technique is applied to every surface
+ -) remainder triangles are split into 3 quads and all recombined quads
+ are split into 4 quads
+ -) Enhancements are performned on the quad mesh.
+ */
+
+int Recombine_All (Mesh *THEM)
+{
+ if(Tree_Nbr(THEM->Surfaces))
+ {
+ // recombine all surfaces
+ List_T *surfaces = Tree2List(THEM->Surfaces);
+ for(int i = 0; i < List_Nbr(surfaces); i++){
+ Surface *s;
+ List_Read(surfaces, i, &s);
+ if(s->Recombine){
+ Msg(STATUS3, "Recombining surface %d", s->Num);
+ Recombine(s->Vertices, s->Simplexes, s->Quadrangles, s->RecombineAngle);
+ }
+ }
+
+ // add 2nd order nodes to all elements
+
+ Degre2(2);
+
+ // then split everybody
+
+ for(int i = 0; i < List_Nbr(surfaces); i++)
+ {
+ Surface *s;
+ List_Read(surfaces, i, &s);
+ List_T *Quadrangles = Tree2List ( s->Quadrangles );
+
+ for (int j=0 ; j<List_Nbr (Quadrangles); j++)
+ {
+ Quadrangle *q;
+ List_Read (Quadrangles, j, &q);
+ Quadrangle *q1 = Create_Quadrangle(q->V[0], q->VSUP[0], q->VSUP[4], q->VSUP[3]);
+ Quadrangle *q2 = Create_Quadrangle(q->V[1], q->VSUP[1], q->VSUP[4], q->VSUP[0]);
+ Quadrangle *q3 = Create_Quadrangle(q->V[2], q->VSUP[2], q->VSUP[4], q->VSUP[1]);
+ Quadrangle *q4 = Create_Quadrangle(q->V[3], q->VSUP[3], q->VSUP[4], q->VSUP[2]);
+ Tree_Add ( s->Quadrangles, &q1);
+ Tree_Add ( s->Quadrangles, &q2);
+ Tree_Add ( s->Quadrangles, &q3);
+ Tree_Add ( s->Quadrangles, &q4);
+ Tree_Suppress (s->Quadrangles, &q);
+ }
+
+ List_Delete(Quadrangles);
+
+ List_T *Triangles = Tree2List ( s->Simplexes );
+
+ if(s->Recombine)
+ {
+ for (int j=0 ; j<List_Nbr (Triangles); j++)
+ {
+ Simplex *t;
+ List_Read (Triangles, j, &t);
+ Vertex *c = Create_Vertex(++THEM->MaxPointNum,
+ (t->V[0]->Pos.X+t->V[1]->Pos.X+t->V[2]->Pos.X)/3.,
+ (t->V[0]->Pos.Y+t->V[1]->Pos.Y+t->V[2]->Pos.Y)/3.,
+ (t->V[0]->Pos.Z+t->V[1]->Pos.Z+t->V[2]->Pos.Z)/3.,
+ (t->V[0]->lc+t->V[1]->lc+t->V[2]->lc)/3.,
+ (t->V[0]->u+t->V[1]->u+t->V[2]->u)/3.);
+ Quadrangle *q1 = Create_Quadrangle(t->V[0], t->VSUP[0], c, t->VSUP[2]);
+ Quadrangle *q2 = Create_Quadrangle(t->V[1], t->VSUP[1], c, t->VSUP[0]);
+ Quadrangle *q3 = Create_Quadrangle(t->V[2], t->VSUP[2], c, t->VSUP[1]);
+ Tree_Add ( s->Quadrangles, &q1);
+ Tree_Add ( s->Quadrangles, &q2);
+ Tree_Add ( s->Quadrangles, &q3);
+ Tree_Suppress (s->Simplexes, &t);
+ Free_Simplex(&t,0);
+ }
+ }
+ else
+ {
+ for (int j=0 ; j<List_Nbr (Triangles); j++)
+ {
+ Simplex *t;
+ List_Read (Triangles, j, &t);
+ Simplex *t1 = Create_Simplex(t->V[0], t->VSUP[0],t->VSUP[2],0);
+ Simplex *t2 = Create_Simplex(t->V[1], t->VSUP[1],t->VSUP[0],0);
+ Simplex *t3 = Create_Simplex(t->V[2], t->VSUP[2],t->VSUP[1],0);
+ Simplex *t4 = Create_Simplex(t->VSUP[0], t->VSUP[1],t->VSUP[2],0);
+ Tree_Add ( s->Simplexes, &t1);
+ Tree_Add ( s->Simplexes, &t2);
+ Tree_Add ( s->Simplexes, &t3);
+ Tree_Add ( s->Simplexes, &t4);
+ Tree_Suppress (s->Simplexes, &t);
+ Free_Simplex(&t,0);
+ }
+ }
+ List_Delete(Triangles);
+ }
+ List_Delete(surfaces);
+ Degre1();
+ }
+ return 1;
+}
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 47e0bbb039..b2a5b0b4d8 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,16 +1,126 @@
#include "BDS.h"
+#include <map>
#include <math.h>
+#include "Numeric.h"
-void BDS_Mesh :: add_point (int num , double x, double y,double z)
+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)
{
- points.insert ( BDS_Point ( num, x, y, z) );
+ /// plane ax+by+cz+1 = 0;
+
+ 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 b[3] = {-1,-1,-1};
+ double pl[3] , det ;
+ sys3x3 (mat,b,pl,&det);
+ double n[3] = {pl[0],pl[1],pl[2]};
+// norme(n);
+ // pp(x,y,z) = p + k n
+ // pp belongs to the plane
+ // a x + b y + c z + 1 = 0
+ // x = xa + k nx
+ // y = ya + k ny
+ // z = za + k nz
+ // a (xa + k nx) + b ( ya + k ny ) + c ( za + k nz ) + 1 = 0
+ // k = - ( a xa + b ya + c za + 1) / (a nx + b ny + c nz )
+ double k = - ( pl[0] * xa + pl[1] * ya + pl[2] * za + 1 ) / ( pl[0] * n[0] + pl[1] * n[1] + pl[2] * n[2] );
+ x = xa + k * n[0];
+ y = ya + k * n[1];
+ z = za + k * n[2];
+}
+
+void print_face (BDS_Triangle *t)
+{
+ BDS_Point *pts[3];
+ t->getNodes (pts);
+ printf("face %p with nodes %d %d %d and edges %p (%d %d) %p (%d %d) %p (%d %d)\n",t,pts[0]->iD,pts[1]->iD,pts[2]->iD,
+ t->e1,t->e1->p1->iD,t->e1->p2->iD,t->e2,t->e2->p1->iD,t->e2->p2->iD,t->e3,t->e3->p1->iD,t->e3->p2->iD);
+}
+
+void vector_triangle (BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3])
+{
+ double a[3] = { p1->X - p2->X,p1->Y - p2->Y,p1->Z - p2->Z};
+ double b[3] = { p1->X - p3->X,p1->Y - p3->Y,p1->Z - p3->Z};
+ 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];
+}
+
+void normal_triangle (BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3])
+{
+ vector_triangle (p1,p2,p3,c);
+ double l = sqrt (c[0]*c[0]+c[1]*c[1]+c[2]*c[2]);
+ if (l==0)return;
+ c[0]/=l;
+ c[1]/=l;
+ c[2]/=l;
+}
+
+double surface_triangle (BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
+{
+ double c[3];
+ vector_triangle (p1,p2,p3,c);
+ return 0.5*sqrt (c[0]*c[0]+c[1]*c[1]+c[2]*c[2]);
+}
+double quality_triangle (BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
+{
+ double e12 = ((p1->X - p2->X)*(p1->X - p2->X)+
+ (p1->Y - p2->Y)*(p1->Y - p2->Y)+
+ (p1->Z - p2->Z)*(p1->Z - p2->Z));
+ double e22 = ((p1->X - p3->X)*(p1->X - p3->X)+
+ (p1->Y - p3->Y)*(p1->Y - p3->Y)+
+ (p1->Z - p3->Z)*(p1->Z - p3->Z));
+ double e32 = ((p3->X - p2->X)*(p3->X - p2->X)+
+ (p3->Y - p2->Y)*(p3->Y - p2->Y)+
+ (p3->Z - p2->Z)*(p3->Z - p2->Z));
+ double a = surface_triangle ( p1,p2,p3);
+ return a / (e12+e22+e32);
+}
+
+
+void BDS_Point::getTriangles (std::list<BDS_Triangle *> &t)
+{
+ t.clear ();
+ std::list<BDS_Edge*>::iterator it = edges.begin();
+ std::list<BDS_Edge*>::iterator ite = edges.end();
+ while(it!=ite)
+ {
+ for (int i=0;i<2;++i)
+ {
+ BDS_Triangle *tt = (*it)->faces[i];
+ if (tt)
+ {
+ std::list<BDS_Triangle*>::iterator tit = t.begin();
+ std::list<BDS_Triangle*>::iterator tite = t.end();
+ int found = 0;
+ while (tit!=tite)
+ {
+ if (tt == *tit)found = 1;
+ ++tit;
+ }
+ if (!found) t.push_back (tt);
+ }
+ }
+ ++it;
+ }
+}
+
+BDS_Point* BDS_Mesh :: add_point (int num , double x, double y,double z)
+{
+ BDS_Point *pp = new BDS_Point ( num, x, y, z);
+ points.insert ( pp );
+ MAXPOINTNUMBER = ( MAXPOINTNUMBER < num )? num : MAXPOINTNUMBER;
+ return pp;
}
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);
+ std::set<BDS_Point*,PointLessThan>::iterator it = points.find(&P);
+ if (it != points.end()) return (BDS_Point*) (*it);
else return 0;
}
@@ -19,20 +129,33 @@ 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);
+ std::set<BDS_Edge*, EdgeLessThan>::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 *BDS_Mesh :: find_edge (BDS_Point *p1, BDS_Point *p2, BDS_Triangle *t) const
+{
+ BDS_Point P1 ( p1->iD ) ;
+ BDS_Point P2 ( p2->iD ) ;
+ BDS_Edge E (&P1, &P2);
+ if (t->e1->p1->iD == E.p1->iD && t->e1->p2->iD == E.p2->iD)return t->e1;
+ if (t->e2->p1->iD == E.p1->iD && t->e2->p2->iD == E.p2->iD)return t->e2;
+ if (t->e3->p1->iD == E.p1->iD && t->e3->p2->iD == E.p2->iD)return t->e3;
+ return 0;
+}
+
+BDS_Edge *BDS_Mesh :: add_edge (int p1, int p2)
{
BDS_Edge *efound = find_edge (p1,p2);
- if (efound)return;
+ if (efound)return efound;
BDS_Point *pp1=find_point(p1);
BDS_Point *pp2=find_point(p2);
if (!pp1 || !pp2)throw;
- edges.insert ( BDS_Edge ( pp1, pp2 ) );
+ BDS_Edge *e = new BDS_Edge ( pp1, pp2 );
+ edges.insert ( e );
+ return e;
}
void BDS_Mesh :: add_triangle (int p1, int p2, int p3, double nx, double ny, double nz)
@@ -43,24 +166,61 @@ void BDS_Mesh :: add_triangle (int p1, int p2, int p3, double nx, double ny, do
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) );
+ BDS_Triangle *t = new BDS_Triangle ( e1, e2, e3 ,find_point(p1));
+ triangles.push_back ( t );
+}
+void BDS_Mesh :: del_triangle (BDS_Triangle *t)
+{
+ t->e1->del ( t );
+ t->e2->del ( t );
+ t->e3->del ( t );
+ t->deleted = true;
+}
+void BDS_Mesh :: del_edge (BDS_Edge *e)
+{
+ edges.erase ( e );
+ e->p1->del ( e );
+ e->p2->del ( e );
+ e->deleted = true;
+ edges_to_delete.insert (e);
}
-
void BDS_Mesh :: add_geom (int p1, int p2)
{
- geom.insert ( BDS_GeomEntity ( p1,p2 ) ) ;
+ geom.insert ( new BDS_GeomEntity ( p1,p2 ) ) ;
}
-BDS_GeomEntity * BDS_Mesh :: get_geom (int p1, int p2)
+
+void BDS_Edge::oppositeof (BDS_Point * oface[2]) const
{
- std::set<BDS_GeomEntity>::iterator it = geom.find(BDS_GeomEntity ( p1,p2 ));
- if (it == geom.end())return 0;
- return (BDS_GeomEntity*) & (*it);
+ 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];
+ }
+ 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];
+ }
}
-void BDS_Mesh :: save_gmsh_format ( const char *filename )
+
+
+BDS_GeomEntity * BDS_Mesh :: get_geom (int p1, int p2)
{
+ BDS_GeomEntity ge ( p1,p2 );
+ std::set<BDS_GeomEntity*,GeomLessThan>::iterator it = geom.find(&ge);
+ if (it == geom.end())return 0;
+ return (BDS_GeomEntity*) (*it);
}
+
void recur_tag ( BDS_Triangle *t , BDS_GeomEntity *g )
{
if (!t->g)
@@ -81,51 +241,565 @@ void recur_tag ( BDS_Triangle *t , BDS_GeomEntity *g )
}
}
+BDS_Vector BDS_Triangle :: N () const
+{
+ BDS_Point *pts[3];
+ getNodes (pts);
+ double c[3];
+ normal_triangle (pts[0],pts[1],pts[2],c);
+ return BDS_Vector ( c[0], c[1], c[2] );
+}
+
+
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)
+
+ static BDS_GeomEntity EDGE_CLAS (0,1);
+ {
+ std::set<BDS_Edge*, EdgeLessThan>::iterator it = edges.begin();
+ std::set<BDS_Edge*, EdgeLessThan>::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)
{
- if (!it->g)
+ 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];
+ 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);
+// 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;
+ }
+ ++it;
+ }
+ }
+ {
+ int tag = 1;
+ while (1)
+ {
+ printf("tag %d\n",tag);
+ std::list<BDS_Triangle*>::iterator it = triangles.begin();
+ std::list<BDS_Triangle*>::iterator ite = triangles.end();
+ BDS_Triangle *start = 0;
+ while (it!=ite)
+ {
+ if (!(*it)->deleted && !(*it)->g)
+ {
+ start = (BDS_Triangle*) (*it);
+ }
+ if (start)break;
+ ++it;
+ }
+ if (!start)break;
+ add_geom (tag, 2);
+ BDS_GeomEntity *g = get_geom (tag,2);
+ recur_tag ( start , g );
+ tag++;
+ }
+ }
+// return;
+ {
+ std::map< std::pair<int,int> , int > edgetags;
+ std::set<BDS_Edge*, EdgeLessThan>::iterator it = edges.begin();
+ std::set<BDS_Edge*, EdgeLessThan>::iterator ite = edges.end();
+ int edgetag = 1;
+ while (it!=ite)
+ {
+ BDS_Edge &e = *((BDS_Edge *) *it);
+
+ if (e.g)
+ {
+ std::map< std::pair<int,int> , int >::iterator found;
+ BDS_GeomEntity *g;
+ if ( e.numfaces() == 1)
{
- start = (BDS_Triangle*)& (*it);
+ found = edgetags.find (std::make_pair ( e.faces[0]->g->classif_tag,-1));
}
- ++it;
+ 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));
+ else
+ found = edgetags.find (std::make_pair ( e.faces[0]->g->classif_tag,e.faces[1]->g->classif_tag));
+ }
+ if (e.g)
+ {
+ if (found == edgetags.end())
+ {
+ add_geom (edgetag, 1);
+ g = get_geom (edgetag,1);
+ if ( e.numfaces() == 1)
+ {
+ 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)]
+ = edgetag;
+ else
+ edgetags[std::make_pair ( e.faces[0]->g->classif_tag,e.faces[1]->g->classif_tag)]
+ = edgetag;
+ }
+ edgetag++;
+ }
+ else
+ {
+ g = get_geom(found->second,1);
+ }
+ e.g = g;
+ }
+ }
+ ++it;
+ }
+ }
+}
+
+void BDS_Mesh :: save_gmsh_format ( const char *filename )
+{
+ cleanup();
+ FILE *f = fopen ( filename, "w");
+ {
+ std::set<BDS_Point*,PointLessThan>::iterator it = points.begin();
+ std::set<BDS_Point*,PointLessThan>::iterator ite = points.end();
+
+ fprintf(f,"$NOD\n");
+ fprintf(f,"%d\n",points.size());
+ while(it!=ite)
+ {
+ fprintf(f,"%d %g %g %g\n",(*it)->iD,(*it)->X,(*it)->Y,(*it)->Z);
+ ++it;
+ }
+ fprintf(f,"$ENDNOD\n");
+ }
+ {
+ fprintf(f,"$ELM\n");
+
+ int nbClasEdges = 0;
+ {
+ std::set<BDS_Edge*, EdgeLessThan>::iterator it = edges.begin();
+ std::set<BDS_Edge*, EdgeLessThan>::iterator ite = edges.end();
+ while(it!=ite)
+ {
+ if ((*it)->g)nbClasEdges++;
+ ++it;
+ }
+ }
+
+ fprintf(f,"%d\n",nbClasEdges+triangles.size());
+
+ int k=1;
+ {
+ std::set<BDS_Edge*, EdgeLessThan>::iterator it = edges.begin();
+ std::set<BDS_Edge*, EdgeLessThan>::iterator ite = edges.end();
+ while(it!=ite)
+ {
+ if ((*it)->g)
+ fprintf(f,"%d %d %d %d %d %d %d\n",
+ k++,1,(*it)->g->classif_tag,(*it)->g->classif_tag,2,
+ (*it)->p1->iD, (*it)->p2->iD);
+ ++it;
+
+ }
+ }
+ {
+ std::list<BDS_Triangle*>::iterator it = triangles.begin();
+ std::list<BDS_Triangle*>::iterator ite = triangles.end();
+ while(it!=ite)
+ {
+ if ((*it)->deleted)
+ {
+ printf("aahg\n");
+ throw;
+ }
+ BDS_Point *nod[3];
+ (*it)->getNodes (nod);
+ if ((*it)->g)
+ fprintf(f,"%d %d %d %d %d %d %d %d\n",
+ k++,2,(*it)->g->classif_tag,(*it)->g->classif_tag,3,
+ nod[0]->iD,nod[1]->iD,nod[2]->iD);
+ else
+ fprintf(f,"%d %d %d %d %d %d %d %d\n",
+ k++,2,1,1,3,
+ nod[0]->iD,nod[1]->iD,nod[2]->iD);
+ ++it;
+
+ }
+ }
+ fprintf(f,"$ENDELM\n");
+ }
+ fclose (f);
+}
+
+template <class IT>
+void DESTROOOY ( IT beg , IT end)
+{
+ while (beg != end)
+ {
+ delete *beg;
+ beg ++;
+ }
+}
+void BDS_Mesh :: cleanup()
+{
+ for (std::list<BDS_Triangle*> :: iterator it = triangles.begin();
+ it != triangles.end();
+ it++)
+ {
+ while (it != triangles.end() && (*it)->deleted)
+ {
+ delete *it;
+ it = triangles.erase (it);
+ }
+ }
+ DESTROOOY ( edges_to_delete.begin(),edges_to_delete.end());
+ edges_to_delete.clear();
+}
+
+BDS_Mesh ::~ BDS_Mesh ()
+{
+ DESTROOOY ( geom.begin(),geom.end());
+ DESTROOOY ( points.begin(),points.end());
+ DESTROOOY ( edges.begin(),edges.end());
+ cleanup();
+ DESTROOOY ( triangles.begin(),triangles.end());
+}
+
+bool BDS_Mesh ::split_edge ( BDS_Edge *e, double coord, const BDS_Projector &)
+{
+
+// printf("split start %d\n",e->deleted);
+
+ int nbFaces = e->numfaces();
+
+ if (nbFaces == 1)return false;
+
+ BDS_Point *op[2];
+ BDS_Point *p1 = e->p1;
+ BDS_Point *p2 = e->p2;
+ MAXPOINTNUMBER++;
+ BDS_Point *mid = add_point (MAXPOINTNUMBER ,
+ coord * p1->X + (1-coord) * p2->X,
+ coord * p1->Y + (1-coord) * p2->Y,
+ coord * p1->Z + (1-coord) * p2->Z);
+
+ // we should project
+ 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->iD, op[0]->iD );
+// BDS_Edge *op1_p2 = find_edge ( op[0]->iD,p2->iD );
+// BDS_Edge *p1_op2 = find_edge ( p1->iD, op[1]->iD );
+// BDS_Edge *op2_p2 = find_edge ( op[1]->iD,p2->iD );
+
+// printf("%p %p %p %p\n",p1_op1,op1_p2,p1_op2,op2_p2);
+
+/* printf("edge %p (%d %d) with %d neighbors\n",p1_op1,p1->iD, op[0]->iD,p1_op1->numfaces());
+ print_face(p1_op1->faces[0]);
+ print_face(p1_op1->faces[1]);
+ printf("edge %p (%d %d) with %d neighbors\n",e,p1->iD, p2->iD,e->numfaces());
+ print_face(e->faces[0]);
+ print_face(e->faces[1]);
+*/
+ if (e->faces[0])
+ {
+ g1 = e->faces[0]->g;
+ del_triangle (e->faces[0]);
+ }
+// not a bug !!!
+ if (e->faces[0])
+ {
+ g2 = e->faces[0]->g;
+ del_triangle (e->faces[0]);
+ }
+
+
+ del_edge (e);
+
+
+ BDS_Edge *p1_mid = new BDS_Edge ( p1, mid );
+ edges.insert ( p1_mid );
+ BDS_Edge *mid_p2 = new BDS_Edge ( mid, p2 );
+ edges.insert ( mid_p2 );
+ BDS_Edge *op1_mid = new BDS_Edge ( op[0], mid );
+ edges.insert ( op1_mid );
+ BDS_Edge *mid_op2 = new BDS_Edge ( mid, op[1] );
+ edges.insert ( mid_op2 );
+
+// printf("split ends 1 %d (%d %d) %d %d \n",p1_op1->numfaces(), p1->iD, op[0]->iD, op1_mid->numfaces(),p1_mid->numfaces());
+ BDS_Triangle*t1 = new BDS_Triangle ( p1_op1, op1_mid, p1_mid ,op[0] );
+ BDS_Triangle*t2 = new BDS_Triangle ( op2_p2, mid_op2, mid_p2 ,op[1] );
+ BDS_Triangle*t3 = new BDS_Triangle ( op1_p2, op1_mid, mid_p2 ,p2);
+ BDS_Triangle*t4 = new BDS_Triangle ( p1_op2, mid_op2, p1_mid ,p1);
+
+ t1->g = g1;
+ t2->g = g2;
+ t3->g = g1;
+ t4->g = g2;
+
+ p1_mid->g = ge;
+ mid_p2->g = ge;
+
+ triangles.push_back(t1);
+ triangles.push_back (t2);
+ triangles.push_back (t3);
+ triangles.push_back (t4);
+
+ return true;
+}
+
+bool BDS_Mesh ::swap_edge ( BDS_Edge *e)
+{
+ if (e->deleted)return false;
+
+ int nbFaces = e->numfaces();
+
+
+ BDS_Point *pts1[3];
+ e->faces[0]->getNodes (pts1);
+ BDS_Point *pts2[3];
+ e->faces[1]->getNodes (pts2);
+
+
+ if (nbFaces != 2)return false;
+ if (e->g && e->g->classif_degree != 2)return false;
+
+ BDS_Point *op[2];
+ BDS_Point *p1 = e->p1;
+ BDS_Point *p2 = e->p2;
+ 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] );
+
+ if (e->faces[0])
+ {
+ g1 = e->faces[0]->g;
+ del_triangle (e->faces[0]);
+ }
+// not a bug !!!
+ if (e->faces[0])
+ {
+ g2 = e->faces[0]->g;
+ del_triangle (e->faces[0]);
+ }
+ del_edge (e);
+
+ BDS_Edge *op1_op2 = new BDS_Edge ( op[0], op[1] );
+ edges.insert ( op1_op2 );
+
+ BDS_Point *first_in_0 = 0;
+ for (int k=0;k<3;k++)
+ if (pts1[k] == p1)
+ {
+ if (pts1[(k+1)%3] == op[0])first_in_0 = p1;
+ else first_in_0 = op[0];
+ }
+ BDS_Point *first_in_1 = 0;
+ for (int k=0;k<3;k++)
+ if (pts2[k] == p2)
+ {
+ if (pts2[(k+1)%3] == op[1])first_in_1 = p2;
+ else first_in_1 = op[1];
+ }
+
+ BDS_Triangle*t1 = new BDS_Triangle ( p1_op1, p1_op2, op1_op2 , first_in_0 );
+ BDS_Triangle*t2 = new BDS_Triangle ( op2_p2, op1_op2, op1_p2 , first_in_1 );
+
+ t1->g = g1;
+ t2->g = g2;
+
+ op1_op2->g = ge;
+
+ triangles.push_back(t1);
+ triangles.push_back(t2);
+ return true;
+}
+
+
+bool BDS_Mesh ::collapse_edge ( BDS_Edge *e, BDS_Point *p, const BDS_Projector &)
+{
+
+}
+
+bool BDS_Mesh ::smooth_point ( BDS_Point *p, const BDS_Projector &)
+{
+ std::list<BDS_Edge *>::iterator eit = p->edges.begin();
+ std::list<BDS_Edge *>::iterator eite = p->edges.end();
+ while (eit!=eite)
+ {
+ if ((*eit)->g)
+ if((*eit)->g->classif_degree == 1)return false;
+ ++eit;
+ }
+ double X = 0;
+ double Y = 0;
+ double Z = 0;
+
+ eit = p->edges.begin();
+ while (eit!=eite)
+ {
+ BDS_Point *op = ((*eit)->p1 == p)?(*eit)->p2:(*eit)->p1;
+ X+=op->X;
+ Y+=op->Y;
+ Z+=op->Z;
+ ++eit;
+ }
+
+ X /= p->edges.size();
+ Y /= p->edges.size();
+ Z /= p->edges.size();
+
+ std::list<BDS_Triangle *> t;
+ p->getTriangles (t);
+ {
+ std::list<BDS_Triangle *>::iterator it = t.begin();
+ std::list<BDS_Triangle *>::iterator ite = t.end();
+ double dist = 1.e22;
+ double XX = X,YY=Y,ZZ=Z;
+ while (it != ite)
+ {
+ 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;
}
- if (!start)break;
- add_geom (tag, 2);
- BDS_GeomEntity *g = get_geom (tag,2);
- recur_tag ( start , g );
- tag++;
+ ++it;
}
+ X = XX;
+ Y = YY;
+ Z = ZZ;
}
+ p->X = X;
+ p->Y = Y;
+ p->Z = Z;
+
+ return true;
+}
+
+void BDS_Mesh :: adapt_mesh ( double l, const BDS_Projector &proj)
+{
+ std::set<BDS_Edge*, EdgeLessThan> small_to_long;
+ {
+ std::set<BDS_Edge*, EdgeLessThan>::iterator it = edges.begin();
+ std::set<BDS_Edge*, EdgeLessThan>::iterator ite = edges.end();
+ while (it != ite)
+ {
+ if ((*it)->numfaces()==1)printf("one face\n");
+ small_to_long.insert (*it);
+ ++it;
+ }
+ }
+ {
+ std::set<BDS_Edge*, EdgeLessThan>::iterator it = small_to_long.begin();
+ std::set<BDS_Edge*, EdgeLessThan>::iterator ite = small_to_long.end();
+
+ int nbspl = 0;
+
+ while (it != ite)
+ {
+ BDS_Point *op[2];
+ (*it)->oppositeof (op);
+ double ll = sqrt((op[0]->X-op[1]->X)*(op[0]->X-op[1]->X)+
+ (op[0]->Y-op[1]->Y)*(op[0]->Y-op[1]->Y)+
+ (op[0]->Z-op[1]->Z)*(op[0]->Z-op[1]->Z));
+ if (!(*it)->deleted && (*it)->length() > l && ll > 0.5 * l){
+ split_edge (*it, 0.5 ,proj);
+ nbspl++;
+ }
+ ++it;
+ }
+ printf("nbspl = %d\n",nbspl);
+ }
+// return;
+ cleanup();
+
+ {
+ small_to_long.clear();
+ std::set<BDS_Edge*, EdgeLessThan>::iterator it = edges.begin();
+ std::set<BDS_Edge*, EdgeLessThan>::iterator ite = edges.end();
+ while (it != ite)
+ {
+ small_to_long.insert (*it);
+ ++it;
+ }
+ }
+ {
+ std::set<BDS_Edge*, EdgeLessThan>::iterator it = small_to_long.begin();
+ std::set<BDS_Edge*, EdgeLessThan>::iterator ite = small_to_long.end();
+ while (it != ite)
+ {
+ if (!(*it)->deleted && (*it)->numfaces()==2)
+ {
+ BDS_Point *op[2];
+ (*it)->oppositeof (op);
+
+ double a1 = surface_triangle ( (*it)->p1 , (*it)->p2 , op[0] );
+ double a2 = surface_triangle ( (*it)->p1 , (*it)->p2 , op[1] );
+ double b1 = surface_triangle ( (*it)->p1 , op[0] , op[1] );
+ double b2 = surface_triangle ( (*it)->p2 , op[0] , op[1] );
+
+ double cb1[3],cb2[3];
+ normal_triangle ( (*it)->p1 , op[0] , op[1],cb1 );
+ normal_triangle ( (*it)->p2 , op[0] , op[1],cb2 );
+
+ double prosc = cb1[0]*cb2[0]+cb1[1]*cb2[1]+cb1[2]*cb2[2];
+
+ if (fabs(a1+a2-b1-b2) < 0.1 * (a1+a2) && prosc < 0 && fabs (prosc) > 0.7)
+ {
+ double qa1 = quality_triangle ( (*it)->p1 , (*it)->p2 , op[0] );
+ double qa2 = quality_triangle ( (*it)->p1 , (*it)->p2 , op[1] );
+ double qb1 = quality_triangle ( (*it)->p1 , op[0] , op[1] );
+ double qb2 = quality_triangle ( (*it)->p2 , op[0] , op[1] );
+
+ double qa = (qa1<qa2)?qa1:qa2;
+ double qb = (qb1<qb2)?qb1:qb2;
+// printf("qa %g qb %g ..\n",qa,qb);
+ if (qb > qa)
+ {
+// printf("swop ..\n");
+ swap_edge ( *it );
+ }
+ }
+ }
+ ++it;
+ }
+ }
+ cleanup();
+ {
+ for (int i=0;i<4;++i)
+ {
+ std::set<BDS_Point*, PointLessThan>::iterator it = points.begin();
+ std::set<BDS_Point*, PointLessThan>::iterator ite = points.end();
+ while (it != ite)
+ {
+ smooth_point(*it,proj);
+ ++it;
+ }
+ }
+ }
}
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 10268ec373..251d206607 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -4,13 +4,15 @@
// default values are 0,0,1
#include <set>
+#include <map>
#include <list>
+#include <math.h>
class BDS_GeomEntity
{
public:
- int classif_degree;
int classif_tag;
+ int classif_degree;
inline bool operator < ( const BDS_GeomEntity & other ) const
{
if (classif_degree < other.classif_degree)return true;
@@ -19,7 +21,7 @@ public:
return false;
}
BDS_GeomEntity (int a, int b)
- : classif_degree (a),classif_tag(b)
+ : classif_tag (a),classif_degree(b)
{
}
};
@@ -27,30 +29,84 @@ public:
class BDS_Edge;
class BDS_Triangle;
+class BDS_Mesh;
-class BDS_Point
+class BDS_Vector
{
public:
- int iD;
- BDS_GeomEntity *g;
- double X,Y,Z;
- std::list<BDS_Edge*> edges;
- inline bool operator < ( const BDS_Point & other ) const
+ double x,y,z;
+ BDS_Vector operator + (const BDS_Vector &v)
{
- return iD < other.iD;
+ return BDS_Vector (x+v.x,y+v.y,z+v.z);
}
- BDS_Point ( int id, double x=0, double y=0, double z=0 )
- : iD(id),g(0),X(x),Y(y),Z(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:
+ int iD;
+ double X,Y,Z;
+ BDS_GeomEntity *g;
+ std::list<BDS_Edge*> edges;
+ inline bool operator < ( const BDS_Point & other ) const
+ {
+ return iD < other.iD;
+ }
+ inline void del (BDS_Edge *e)
+ {
+ std::list<BDS_Edge*>::iterator it = edges.begin();
+ std::list<BDS_Edge*>::iterator ite = edges.end();
+ while(it!=ite)
+ {
+ if (*it == e)
+ {
+ edges.erase(it);
+ break;
+ }
+ ++it;
+ }
+ }
+ void getTriangles (std::list<BDS_Triangle *> &t);
+
+ 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:
+ bool deleted;
BDS_Point *p1,*p2;
BDS_GeomEntity *g;
BDS_Triangle*faces[2];
+ 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;
@@ -59,6 +115,11 @@ public:
}
inline void addface ( BDS_Triangle *f)
{
+ if (faces[1])
+ {
+ printf("Non Manifold model not done yet\n");
+ throw 1;
+ }
if(faces [0])faces[1] = f;
else faces[0] = f;
}
@@ -75,8 +136,24 @@ public:
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;
+ }
+
+ inline void oppositeof (BDS_Point * oface[2]) const;
+
BDS_Edge ( BDS_Point *A, BDS_Point *B )
- : g(0)
+ : deleted(false), g(0)
{
faces[0] = faces[1] = 0;
if (*A < *B)
@@ -97,90 +174,122 @@ public:
class BDS_Triangle
{
public:
+ bool deleted;
+ int ori_first_edge;
BDS_Edge *e1,*e2,*e3;
+ BDS_Vector N() const ;
BDS_GeomEntity *g;
- double NX,NY,NZ;
- inline bool operator < ( const BDS_Triangle & other ) const
+ inline void getNodes (BDS_Point *n[3]) 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)
- : g(0),NX(nx),NY(ny),NZ(nz)
- {
- if (*A < *B && *A < *C)
+ if (ori_first_edge == 1)
{
- 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;
- }
+ n[0] = e1->p1;
+ n[1] = e1->p2;
}
else
{
- if (*A < *B)
- {
- e1=C;
- e2=A;
- e3=B;
- }
- else
- {
- e1=C;
- e2=B;
- e3=A;
- }
+ n[0] = e1->p2;
+ n[1] = e1->p1;
}
+ if (e2->p1 != n[0] && e2->p1 != n[1])n[2] = e2->p1;
+ else n[2] = e2->p2;
+ }
+ BDS_Triangle ( BDS_Edge *A, BDS_Edge *B, BDS_Edge *C, BDS_Point *first_vertex)
+ : deleted (false) , e1(A),e2(B),e3(C),g(0)
+ {
+ if (first_vertex == A->p1)
+ ori_first_edge = 1;
+ else if (first_vertex == A->p2)
+ ori_first_edge = -1;
+ else
+ throw;
e1->addface(this);
e2->addface(this);
e3->addface(this);
}
};
+class GeomLessThan
+{
+ public:
+ bool operator()(const BDS_GeomEntity* ent1, const BDS_GeomEntity* ent2) const
+ {
+ return *ent1 < *ent2;
+ }
+};
+class PointLessThan
+{
+ public:
+ bool operator()(const BDS_Point* ent1, const BDS_Point* ent2) const
+ {
+ return *ent1 < *ent2;
+ }
+};
+class EdgeLessThan
+{
+ public:
+ bool operator()(const BDS_Edge* ent1, const BDS_Edge* ent2) const
+ {
+ return *ent1 < *ent2;
+ }
+};
+
+class BDS_Projector
+{
+ public:
+ virtual bool project ( const double &X, const double &Y, const double &Z, const BDS_Vector & n,
+ double &Xp, double &Yp, double &Zp ) const = 0;
+};
+
+class BDS_NoProjector : public BDS_Projector
+{
+ public:
+ bool project ( const double &X, const double &Y, const double &Z, const BDS_Vector & n,
+ double &Xp, double &Yp, double &Zp ) const
+ {
+ Xp = X;
+ Yp = Y;
+ Zp = Z;
+ return true;
+ }
+};
+
+class BDS_BDSProjector : public BDS_Projector
+{
+ BDS_Mesh &m;
+ public:
+ BDS_BDSProjector (BDS_Mesh &_m) : m(_m){}
+ bool project ( const double &X, const double &Y, const double &Z, const BDS_Vector & n,
+ double &Xp, double &Yp, double &Zp ) const;
+};
class BDS_Mesh
{
+ int MAXPOINTNUMBER;
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);
+ std::set<BDS_GeomEntity*,GeomLessThan> geom;
+ std::set<BDS_Point*,PointLessThan> points;
+ std::set<BDS_Edge*, EdgeLessThan> edges;
+ std::set<BDS_Edge*> edges_to_delete;
+ std::list<BDS_Triangle*> triangles;
+ BDS_Point * add_point (int num , double x, double y,double z);
+ BDS_Edge * add_edge (int p1, int p2);
+ void del_edge (BDS_Edge *e);
void add_triangle (int p1, int p2, int p3);
+ void del_triangle (BDS_Triangle *t);
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);
+ BDS_Edge *find_edge (BDS_Point *p1, BDS_Point *p2, BDS_Triangle *t)const;
+ BDS_GeomEntity *get_geom (int p1, int p2);
+ bool swap_edge ( BDS_Edge *);
+ bool collapse_edge ( BDS_Edge *, BDS_Point*, const BDS_Projector &);
+ bool smooth_point ( BDS_Point*, const BDS_Projector &);
+ bool split_edge ( BDS_Edge *, double coord, const BDS_Projector &);
void save_gmsh_format (const char *filename);
void classify ( double angle);
+ BDS_Mesh() : MAXPOINTNUMBER (0){}
+ virtual ~BDS_Mesh ();
+ void adapt_mesh(double, const BDS_Projector &);
+ void cleanup();
};
diff --git a/Mesh/DiscreteSurface.cpp b/Mesh/DiscreteSurface.cpp
index 346c74af64..a83f4d6d26 100644
--- a/Mesh/DiscreteSurface.cpp
+++ b/Mesh/DiscreteSurface.cpp
@@ -1,4 +1,4 @@
-// $Id: DiscreteSurface.cpp,v 1.7 2005-04-11 08:53:15 remacle Exp $
+// $Id: DiscreteSurface.cpp,v 1.8 2005-04-15 14:32:40 remacle Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -169,6 +169,7 @@ void Mesh_To_BDS(Surface *s, BDS_Mesh *m)
void POLY_rep_To_Mesh(POLY_rep *prep, Surface *s)
{
+ printf ("compareposition : %22.15e\n",CTX.lc);
VertexBound = Tree_Create(sizeof(Vertex *), comparePosition);
Tree_Action(THEM->Vertices, InsertInVertexBound);
@@ -379,6 +380,7 @@ void STLEndSolid()
Tree_Add(THEM->Surfaces, &STL_Surf);
Tree_Action(VertexBound, Free_Vertex);
}
+
Tree_Delete(VertexBound);
}
@@ -545,6 +547,7 @@ double SetLC(Vertex *v1, Vertex *v2, double factor)
void SEGM_rep_To_Mesh(SEGM_rep *srep, Curve *c)
{
+
VertexBound = Tree_Create(sizeof(Vertex *), comparePosition);
Tree_Action(THEM->Vertices, InsertInVertexBound);
@@ -610,8 +613,24 @@ int MeshDiscreteSurface(Surface *s)
// hope for the best.
POLY_rep_To_Mesh(s->thePolyRep, s);
BDS_Mesh bds;
+ BDS_NoProjector nop;
Mesh_To_BDS(s,&bds);
- bds.classify ( M_PI / 9 );
+ bds.save_gmsh_format ( "1.msh" );
+ bds.classify ( M_PI / 6 );
+ bds.save_gmsh_format ( "2.msh" );
+ try{
+ for (int i=0;i<15;i++)
+ {
+ bds.adapt_mesh (CTX.lc/30 , nop );
+ printf("%d \n",i);
+ }
+ }
+ catch(...)
+ {
+ printf("coucou\n");
+ }
+ bds.save_gmsh_format ( "3.msh" );
+
return 1;
}
else if(s->Typ == MSH_SURF_DISCRETE){
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 417dc5f10e..63b9a83b9e 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.81 2005-04-12 17:05:39 geuzaine Exp $
+# $Id: Makefile,v 1.82 2005-04-15 14:32:40 remacle Exp $
#
# Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
#
@@ -83,7 +83,7 @@ ${LIB}: ${OBJ}
${RANLIB} ${LIB}
.cpp.o:
- ${CXX} ${CFLAGS} -c $<
+ ${CXX} ${CFLAGS} -c -pg $<
# Don't optimize 3D_Mesh: it sometimes mysteriously crashes on Linux
3D_Mesh.o:
@@ -114,8 +114,8 @@ depend:
../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
../Mesh/DiscreteSurface.h ../Common/VertexArray.h \
- ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h Utils.h \
- Create.h 2D_Mesh.h ../Common/Context.h Interpolation.h
+ ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h Mesh.h \
+ Utils.h Vertex.h Create.h 2D_Mesh.h ../Common/Context.h Interpolation.h
2D_Transfinite.o: 2D_Transfinite.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../Geo/Geo.h Mesh.h Vertex.h \
@@ -128,7 +128,7 @@ depend:
../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
../Mesh/DiscreteSurface.h ../Common/VertexArray.h \
- ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h
+ ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h Mesh.h
2D_BGMesh.o: 2D_BGMesh.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
@@ -196,21 +196,22 @@ depend:
../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
../Mesh/DiscreteSurface.h ../Common/VertexArray.h \
../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h \
- Interpolation.h 2D_Mesh.h Create.h ../Common/Context.h
+ Interpolation.h Vertex.h Mesh.h 2D_Mesh.h Create.h ../Common/Context.h
2D_Mesh_Aniso.o: 2D_Mesh_Aniso.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
../Mesh/DiscreteSurface.h ../Common/VertexArray.h \
- ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h \
- Interpolation.h Create.h ../Common/Context.h
+ ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h Mesh.h \
+ Interpolation.h Vertex.h Create.h ../Common/Context.h
2D_Mesh_Triangle.o: 2D_Mesh_Triangle.cpp ../Common/Gmsh.h \
../Common/Message.h ../DataStr/Malloc.h ../DataStr/List.h \
../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h Mesh.h Vertex.h \
Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h \
DiscreteSurface.h ../Common/VertexArray.h ../Common/SmoothNormals.h \
- Metric.h Matrix.h ../Numeric/Numeric.h ../Common/Context.h
+ Metric.h Matrix.h ../Numeric/Numeric.h ../Common/Context.h \
+ ../Triangle/triangle.h
3D_Mesh.o: 3D_Mesh.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
@@ -237,16 +238,16 @@ depend:
../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
../Mesh/DiscreteSurface.h ../Common/VertexArray.h \
- ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h \
- ../Common/Context.h Create.h
+ ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h Mesh.h \
+ ../Common/Context.h Create.h Vertex.h
3D_Extrude_Old.o: 3D_Extrude_Old.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
../Mesh/DiscreteSurface.h ../Common/VertexArray.h \
- ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h \
- ../Common/Context.h Create.h
+ ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h Mesh.h \
+ ../Common/Context.h Create.h Vertex.h
3D_Coherence.o: 3D_Coherence.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
@@ -270,8 +271,9 @@ depend:
../DataStr/avl.h ../DataStr/Tools.h Mesh.h Vertex.h Element.h Simplex.h \
Face.h Edge.h ../Geo/ExtrudeParams.h DiscreteSurface.h \
../Common/VertexArray.h ../Common/SmoothNormals.h Metric.h Matrix.h \
- Create.h ../Numeric/Numeric.h ../Common/Context.h
-BDS.o: BDS.cpp BDS.h
+ Create.h ../Numeric/Numeric.h ../Common/Context.h \
+ ../Netgen/libsrc/interface/nglib.h ../Netgen/nglib_addon.h
+BDS.o: BDS.cpp BDS.h ../Numeric/Numeric.h
MeshQuality.o: MeshQuality.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
@@ -284,8 +286,8 @@ Create.o: Create.cpp ../Common/Gmsh.h ../Common/Message.h \
../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
../Mesh/DiscreteSurface.h ../Common/VertexArray.h \
- ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h Utils.h \
- ../Common/Context.h Create.h
+ ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h Mesh.h \
+ Utils.h Vertex.h ../Common/Context.h Create.h
Generator.o: Generator.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
@@ -300,23 +302,24 @@ Print_Mesh.o: Print_Mesh.cpp ../Common/Gmsh.h ../Common/Message.h \
../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
../Mesh/DiscreteSurface.h ../Common/VertexArray.h \
- ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h Create.h \
- ../Common/Context.h
+ ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h Mesh.h \
+ Create.h Vertex.h ../Common/Context.h
Read_Mesh.o: Read_Mesh.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../Geo/Geo.h ../Geo/CAD.h \
../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h \
../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
../Mesh/DiscreteSurface.h ../Common/VertexArray.h \
- ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h 3D_Mesh.h \
- Create.h ../Geo/MinMax.h ../Common/Context.h
+ ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h Mesh.h \
+ 3D_Mesh.h Create.h Vertex.h ../Geo/MinMax.h ../Common/Context.h
DiscreteSurface.o: DiscreteSurface.cpp ../Common/Gmsh.h \
../Common/Message.h ../DataStr/Malloc.h ../DataStr/List.h \
../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
../Numeric/Numeric.h Mesh.h Vertex.h Element.h Simplex.h Face.h Edge.h \
../Geo/ExtrudeParams.h DiscreteSurface.h ../Common/VertexArray.h \
- ../Common/SmoothNormals.h Metric.h Matrix.h ../Geo/CAD.h ../Geo/Geo.h \
- Create.h Interpolation.h ../Common/Context.h BDS.h
+ ../Common/SmoothNormals.h Metric.h Matrix.h ../Geo/CAD.h ../Mesh/Mesh.h \
+ ../Mesh/Vertex.h ../Geo/Geo.h Create.h Interpolation.h \
+ ../Common/Context.h BDS.h
SwapEdge.o: SwapEdge.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
@@ -329,29 +332,30 @@ Utils.o: Utils.cpp ../Common/Gmsh.h ../Common/Message.h \
../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
../Mesh/DiscreteSurface.h ../Common/VertexArray.h \
- ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h \
- Interpolation.h ../Common/Context.h
+ ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h Mesh.h \
+ Interpolation.h Vertex.h ../Common/Context.h
Metric.o: Metric.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
../Mesh/DiscreteSurface.h ../Common/VertexArray.h \
- ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h \
- Interpolation.h
+ ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Matrix.h Mesh.h \
+ Matrix.h Interpolation.h Vertex.h
Nurbs.o: Nurbs.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h Nurbs.h Vertex.h Mesh.h Element.h \
Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h DiscreteSurface.h \
../Common/VertexArray.h ../Common/SmoothNormals.h Metric.h Matrix.h \
- ../Geo/Geo.h ../Geo/GeoUtils.h Create.h ../Geo/CAD.h
+ ../Geo/Geo.h ../Geo/GeoUtils.h ../Mesh/Mesh.h Create.h ../Geo/CAD.h \
+ ../Mesh/Vertex.h
Interpolation.o: Interpolation.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
Nurbs.h Vertex.h Mesh.h Element.h Simplex.h Face.h Edge.h \
../Geo/ExtrudeParams.h DiscreteSurface.h ../Common/VertexArray.h \
- ../Common/SmoothNormals.h Metric.h Matrix.h ../Geo/CAD.h Utils.h \
- Interpolation.h
+ ../Common/SmoothNormals.h Metric.h Matrix.h ../Geo/CAD.h ../Mesh/Mesh.h \
+ ../Mesh/Vertex.h Utils.h Interpolation.h
SecondOrder.o: SecondOrder.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../Geo/Geo.h Mesh.h Vertex.h \
--
GitLab