From 14e0fc6be6baf9a7e524b4e79b0d4e5523d5a32a Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Mon, 9 May 2005 06:56:13 +0000
Subject: [PATCH] *** empty log message ***
---
Mesh/BDS.cpp | 211 +++++++++++++++++++++++++++++++--------
Mesh/BDS.h | 27 +++++
Mesh/DiscreteSurface.cpp | 4 +-
3 files changed, 199 insertions(+), 43 deletions(-)
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index c92d99023f..9fea688a09 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -3,6 +3,13 @@
#include <math.h>
#include "Numeric.h"
+/*
+ (X-Xc)^2 = R^2
+
+
+
+*/
+
BDS_Sphere :: BDS_Sphere( BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, BDS_Point *p4)
@@ -431,7 +438,7 @@ void BDS_Mesh :: reverseEngineerCAD ( )
if (fabs(dist) > distmax) distmax = fabs(dist);
++pit;
}
- if (distmax < 1.e-2 * LC )
+ if (distmax < 1.e-116 * LC )
{
(*it)->surf = sphere;
printf("spherical surface %d found %g (R=%g,LC=%g)\n",(*it)->classif_tag,distmax,sphere->R,LC);
@@ -474,18 +481,7 @@ void BDS_Mesh :: classify ( double angle )
{
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);
-// 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]);
+ double ag = N0.angle (N1);
if (fabs(ag) > angle)e.g = &EDGE_CLAS;
}
else
@@ -622,6 +618,42 @@ void BDS_Mesh :: classify ( double angle )
++it;
}
}
+ {
+ std::set<BDS_Point*,PointLessThan>::iterator it = points.begin();
+ std::set<BDS_Point*,PointLessThan>::iterator ite = points.end();
+ while (it != ite)
+ {
+ if ((*it)->g->classif_degree > 1)
+ {
+ std::list<BDS_Triangle *> t;
+ (*it)->getTriangles (t);
+ std::list<BDS_Triangle *>::iterator tit = t.begin();
+ std::list<BDS_Triangle *>::iterator tite = t.end();
+ BDS_Vector SumN;
+ while (tit!=tite)
+ {
+ SumN += (*tit)->N();
+ ++tit;
+ }
+ double max_angle = 0;
+ double NORM = sqrt(SumN*SumN);
+ SumN /= NORM;
+ tit = t.begin();
+ while (tit!=tite)
+ {
+ double ag = SumN.angle ((*tit)->N());
+ if (fabs(ag) > max_angle) max_angle = fabs(ag);
+ ++tit;
+ }
+ if (fabs(max_angle) > angle)
+ {
+ add_geom (vertextag, 0);
+ (*it)->g = get_geom(vertextag++,0);
+ }
+ }
+ ++it;
+ }
+ }
reverseEngineerCAD ( ) ;
printf(" end classifying %d edgetags %d vertextags\n",edgetag-1,vertextag-1);
@@ -977,9 +1009,33 @@ bool BDS_Mesh :: read_mesh ( const char *filename )
}
}
}
+ else if (!strcmp (name,"Quadrilaterals"))
+ {
+
+ int nbt,cl,n1,n2,n3,n4;
+ fgets (buffer,255,f);
+ sscanf (buffer,"%d",&nbt);
+ printf("%d Quadrilaterals\n",nbt);
+ for (int i=0;i<nbt;i++)
+ {
+ fgets (buffer,255,f);
+ sscanf (buffer,"%d %d %d %d %d",&n1,&n2,&n3,&n4,&cl);
+ BDS_Triangle *t1 = add_triangle(n1,n2,n3);
+ // pas 12
+ // pas 13
+ BDS_Triangle *t2 = add_triangle(n1,n3,n4);
+ t1->g = get_geom (cl,2);
+ if (!t1->g)
+ {
+ add_geom (cl,2);
+ t1->g = get_geom (cl,2);
+ }
+ t2->g = t1->g;
+ }
+ }
}
}
-
+
LC = sqrt ((Min[0]-Max[0])*(Min[0]-Max[0])+
(Min[1]-Max[1])*(Min[1]-Max[1])+
(Min[2]-Max[2])*(Min[2]-Max[2]));
@@ -1122,6 +1178,23 @@ bool BDS_Mesh ::split_edge ( BDS_Edge *e, double coord)
{
+ /*
+
+ p1
+ / | \
+ / | \
+ op1/ 0mid1 \op2
+ \ | /
+ \ | /
+ \ p2/
+
+// p1,op1,mid -
+// p2,op2,mid -
+// p2,op1,mid +
+// p1,op2,mid +
+
+ */
+
int nbFaces = e->numfaces();
if (nbFaces != 2)return false;
@@ -1129,6 +1202,26 @@ bool BDS_Mesh ::split_edge ( BDS_Edge *e, double coord)
BDS_Point *op[2];
BDS_Point *p1 = e->p1;
BDS_Point *p2 = e->p2;
+
+ BDS_Point *pts1[3];
+ e->faces(0)->getNodes (pts1);
+
+
+ int orientation = 0;
+ for (int i=0;i<3;i++)
+ {
+ if (pts1[i] == p1)
+ {
+ if (pts1[(i+1)%3] == p2)
+ orientation = 1;
+ else
+ orientation = -1;
+
+ break;
+ }
+ }
+
+
MAXPOINTNUMBER++;
BDS_Point *mid = add_point (MAXPOINTNUMBER ,
coord * p1->X + (1-coord) * p2->X,
@@ -1145,20 +1238,6 @@ bool BDS_Mesh ::split_edge ( BDS_Edge *e, double coord)
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;
@@ -1185,11 +1264,21 @@ bool BDS_Mesh ::split_edge ( BDS_Edge *e, double coord)
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 );
- BDS_Triangle*t2 = new BDS_Triangle ( op2_p2, mid_op2, mid_p2 );
- BDS_Triangle*t3 = new BDS_Triangle ( op1_p2, op1_mid, mid_p2 );
- BDS_Triangle*t4 = new BDS_Triangle ( p1_op2, mid_op2, p1_mid );
-
+ BDS_Triangle *t1,*t2,*t3,*t4;
+ if (orientation == 1)
+ {
+ t1 = new BDS_Triangle ( op1_mid, p1_op1,p1_mid );
+ t2 = new BDS_Triangle ( mid_op2, op2_p2, mid_p2 );
+ t3 = new BDS_Triangle ( op1_p2, op1_mid, mid_p2 );
+ t4 = new BDS_Triangle ( p1_op2, mid_op2, p1_mid );
+ }
+ else
+ {
+ t1 = new BDS_Triangle ( p1_op1,op1_mid,p1_mid );
+ t2 = new BDS_Triangle ( op2_p2, mid_op2,mid_p2 );
+ t3 = new BDS_Triangle ( op1_mid,op1_p2, mid_p2 );
+ t4 = new BDS_Triangle ( mid_op2, p1_op2,p1_mid );
+ }
t1->g = g1;
t2->g = g2;
t3->g = g1;
@@ -1212,16 +1301,29 @@ bool BDS_Mesh ::split_edge ( BDS_Edge *e, double coord)
bool BDS_Mesh ::swap_edge ( BDS_Edge *e)
{
+
+ /*
+
+ p1
+ / | \
+ / | \
+ op1/ 0 | 1 \op2
+ \ | /
+ \ | /
+ \ p2/
+
+ // op1 p1 op2
+ // op1 op2 p2
+
+ */
+
+
if (e->deleted)return false;
int nbFaces = e->numfaces();
if (nbFaces != 2)return false;
- BDS_Point *pts1[3];
- e->faces(0)->getNodes (pts1);
- BDS_Point *pts2[3];
- e->faces(1)->getNodes (pts2);
if (e->g && e->g->classif_degree != 2)return false;
@@ -1231,6 +1333,24 @@ bool BDS_Mesh ::swap_edge ( BDS_Edge *e)
e->oppositeof (op);
BDS_GeomEntity *g1=0,*g2=0,*ge=e->g;
+ BDS_Point *pts1[3];
+ e->faces(0)->getNodes (pts1);
+
+
+ int orientation = 0;
+ for (int i=0;i<3;i++)
+ {
+ if (pts1[i] == p1)
+ {
+ if (pts1[(i+1)%3] == p2)
+ orientation = 1;
+ else
+ orientation = -1;
+
+ break;
+ }
+ }
+
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) );
@@ -1252,9 +1372,18 @@ bool BDS_Mesh ::swap_edge ( BDS_Edge *e)
BDS_Edge *op1_op2 = new BDS_Edge ( op[0], op[1] );
edges.insert ( op1_op2 );
- BDS_Triangle*t1 = new BDS_Triangle ( p1_op1, p1_op2, op1_op2 );
- BDS_Triangle*t2 = new BDS_Triangle ( op2_p2, op1_op2, op1_p2 );
-
+ BDS_Triangle*t1,*t2;
+ if (orientation == 1)
+ {
+ t1 = new BDS_Triangle ( p1_op1, p1_op2, op1_op2 );
+ t2 = new BDS_Triangle ( op1_op2, op2_p2, op1_p2 );
+ }
+ else
+ {
+ t1 = new BDS_Triangle ( p1_op2,p1_op1, op1_op2 );
+ t2 = new BDS_Triangle ( op2_p2,op1_op2, op1_p2 );
+ }
+
t1->g = g1;
t2->g = g2;
@@ -1492,7 +1621,7 @@ int BDS_Mesh :: adapt_mesh ( double l, bool smooth, BDS_Mesh *geom_mesh)
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/0.7 && ll > 0.5 * l){
+ if (!(*it)->deleted && (*it)->length() > l/0.7 && ll > 0.25 * l){
split_edge (*it, 0.5 );
nb_modif++;
}
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 0e321c155c..0ea17a2b6d 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -102,6 +102,20 @@ public:
z+=v.z;
return *this;
}
+ BDS_Vector& operator *= (const double &v)
+ {
+ x*=v;
+ y*=v;
+ z*=v;
+ return *this;
+ }
+ BDS_Vector& operator /= (const double &v)
+ {
+ x/=v;
+ y/=v;
+ z/=v;
+ return *this;
+ }
BDS_Vector operator / (const double &v)
{
return BDS_Vector (x/v,y/v,z/v);
@@ -110,6 +124,19 @@ public:
{
return BDS_Vector (x*v,y*v,z*v);
}
+ double angle (const BDS_Vector &v) const
+ {
+ double a[3] = { x , y , z };
+ double b[3] = { v.x , v.y , v.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);
+ return ag;
+ }
double operator * (const BDS_Vector &v) const
{
return (x*v.x+y*v.y+z*v.z);
diff --git a/Mesh/DiscreteSurface.cpp b/Mesh/DiscreteSurface.cpp
index ba20185f56..eba9abed8e 100644
--- a/Mesh/DiscreteSurface.cpp
+++ b/Mesh/DiscreteSurface.cpp
@@ -1,4 +1,4 @@
-// $Id: DiscreteSurface.cpp,v 1.12 2005-05-06 16:06:42 remacle Exp $
+// $Id: DiscreteSurface.cpp,v 1.13 2005-05-09 06:56:13 remacle Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -501,7 +501,7 @@ 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 / 100, true))
+ while (iter < 20 && THEM->bds_mesh -> adapt_mesh ( CTX.mesh.lc_factor * THEM->bds->LC, true,THEM->bds))
{
printf("iter %d done\n",iter);
iter ++;
--
GitLab