diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index cdc4db20fc889035031f9597703e7d5bfa495fca..e3aed1d68704cf6ab22dcd40f3bb06e958a32af4 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.455 2006-08-27 23:10:35 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.456 2006-08-29 10:39:48 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -3521,9 +3521,7 @@ void mesh_optimize_cb(CALLBACK_ARGS)
void mesh_remesh_cb(CALLBACK_ARGS)
{
- ReMesh();
- Draw();
- Msg(STATUS2N, " ");
+ Msg(GERROR, "Surface ReMeshing must be reinterfaced");
}
void mesh_update_edges_cb(CALLBACK_ARGS)
diff --git a/Geo/GFace.h b/Geo/GFace.h
index def562358c674cfeb8d923fde60df42092db2e39..2daf19571c65f1b128305a1c5c8c65fafcc466bf 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -48,7 +48,10 @@ class GFace : public GEntity
std::list<GVertex *> embedded_vertices;
public:
- GFace(GModel *model, int tag) : GEntity(model, tag), r1(0), r2(0) {}
+ GFace(GModel *model, int tag) : GEntity(model, tag), r1(0), r2(0)
+ {
+ meshAttributes.recombine = 0;
+ }
virtual ~GFace();
void addRegion(GRegion *r){ r1 ? r2 = r : r1 = r; }
@@ -126,6 +129,12 @@ class GFace : public GEntity
std::vector<MTriangle*> triangles;
std::vector<MQuadrangle*> quadrangles;
+
+ struct {
+ int recombine;
+ double recombineAngle;
+ } meshAttributes ;
+
};
#endif
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 3aa0ff290fb72259430f87256f22b7d926b40d1b..332f886bd60228ef1bb70a7a3e0022bc324f5d11 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshFace.cpp,v 1.16 2006-08-23 19:53:39 geuzaine Exp $
+// $Id: gmshFace.cpp,v 1.17 2006-08-29 10:39:48 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -68,6 +68,9 @@ gmshFace::gmshFace(GModel *m, Surface *face)
embedded_vertices.push_back(gv);
}
}
+
+ meshAttributes.recombine = s->Recombine;
+ meshAttributes.recombineAngle = s->RecombineAngle;
}
gmshFace::gmshFace(GModel *m, int num)
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 46a2dbeda897c10c9077acffe854b98e9186693d..6896150dcab317ee479c72ae2bf2c62040051cf1 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.59 2006-08-26 15:13:22 remacle Exp $
+// $Id: BDS.cpp,v 1.60 2006-08-29 10:39:48 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -27,7 +27,7 @@
#include "Message.h"
#include "GFace.h"
-bool test_move_point_parametric_triangle (BDS_Point * p, double u, double v, BDS_Triangle *t);
+bool test_move_point_parametric_triangle (BDS_Point * p, double u, double v, BDS_Face *t);
double BDS_Point::min_edge_length()
{
@@ -42,14 +42,14 @@ double BDS_Point::min_edge_length()
return L;
}
-void outputScalarField(std::list < BDS_Triangle * >t, const char *iii)
+void outputScalarField(std::list < BDS_Face * >t, const char *iii)
{
FILE *f = fopen(iii, "w");
fprintf(f, "View \"scalar\" {\n");
- std::list < BDS_Triangle * >::iterator tit = t.begin();
- std::list < BDS_Triangle * >::iterator tite = t.end();
+ std::list < BDS_Face * >::iterator tit = t.begin();
+ std::list < BDS_Face * >::iterator tite = t.end();
while(tit != tite) {
- BDS_Point *pts[3];
+ BDS_Point *pts[4];
(*tit)->getNodes(pts);
fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
pts[0]->X, pts[0]->Y, pts[0]->Z,
@@ -63,148 +63,11 @@ void outputScalarField(std::list < BDS_Triangle * >t, const char *iii)
fclose(f);
}
-double BDS_Quadric::normalCurv(double x, double y, double z) const
-{
- // K = div n = div ( Grad(f) / ||Grad(f)||)
- // = Laplacian (f) / ||Grad(f)|| - Grad^t(f) Hessian Grad(f) / ||Grad(f)||^3
- BDS_Vector g = Gradient(x, y, z);
- const double normG2 = 1. / (g.x * g.x + g.y * g.y + g.z * g.z);
- double normG = sqrt(normG2);
- double c1 = 2 * (a + b + c);
- double c2 =
- 2 * a * g.x * g.x +
- 2 * b * g.y * g.y +
- 2 * c * g.z * g.z +
- 4 * d * g.x * g.y + 4 * e * g.x * g.z + 4 * f * g.y * g.z;
-
- // printf("%g %g %g %g %g %g\n",a,b,c,normG,c1, c2 / (normG*normG*normG));
- return fabs(c1 - c2 * (normG2)) * normG;
-}
-
-void BDS_Quadric::projection(double xa, double ya, double za,
- double &x, double &y, double &z) const
-{
- x = xa;
- y = ya;
- z = za;
- const double residual = signedDistanceTo(xa, ya, za);
-
- int ITER = 0;
- while(1) {
- // the equation to solve is
- // X = XA + z * N (XA)
- // N = Gradient F (XA)
- // F(X) is the quadric :: X^t A X + X^t B + C
- // z^2 ( N^t A N ) + z ( XA^t A N + N^t A XA + N^t B) + F ( XA )
- BDS_Vector grad = Gradient(x, y, z);
- double coef1 =
- a * grad.x * grad.x +
- b * grad.y * grad.y +
- c * grad.z * grad.z +
- 2 * d * grad.x * grad.y +
- 2 * e * grad.x * grad.z + 2 * f * grad.y * grad.z;
- double coef2 =
- 2 * a * xa * grad.x +
- 2 * b * ya * grad.y +
- 2 * c * za * grad.z +
- 2 * d * (xa * grad.y + ya * grad.x) +
- 2 * e * (xa * grad.z + za * grad.x) +
- 2 * f * (za * grad.y + ya * grad.z) +
- grad.x * g + grad.y * h + grad.z * i;
- double delta = coef2 * coef2 - 4. * coef1 * residual;
- if(delta < 0) {
- // printf("aaargh error projection delta = %g\n",delta);
- }
- else {
- double alpha;
- double alpha1 = (-coef2 + sqrt(delta)) / (2 * coef1);
- double alpha2 = (-coef2 - sqrt(delta)) / (2 * coef1);
- if(fabs(alpha1) > fabs(alpha2))
- alpha = alpha2;
- else
- alpha = alpha1;
- x = xa + alpha * grad.x;
- y = ya + alpha * grad.y;
- z = za + alpha * grad.z;
- }
- // printf("point (%g,%g,%g) ==> (%g,%g,%g) dist %g\n",
- // xa,ya,za,x,y,z,signedDistanceTo(x,y,z));;
- if(ITER++ == 0)
- break;
- }
-}
-
-void BDS_GeomEntity::getClosestTriangles(double x, double y, double z,
- std::list < BDS_Triangle * >&l,
- double &radius ,
- double &X, double &Y, double &Z)
-{
-#ifdef HAVE_ANN_
- l.clear();
-
- if(kdTree == 0) {
- l = t;
- return;
- }
-
- queryPt[0] = x;
- queryPt[1] = y;
- queryPt[2] = z;
-
- double eps = 0.;
- kdTree->annkSearch(queryPt, // query point
- nbK, // number of near neighbors
- nnIdx, // nearest neighbors (returned)
- dists, // distance (returned)
- eps); // error bound
-
- std::list < BDS_Triangle * >l1, l2, l3;
-
- X = kdTree->thePoints()[nnIdx[0]][0];
- Y = kdTree->thePoints()[nnIdx[0]][1];
- Z = kdTree->thePoints()[nnIdx[0]][2];
-
- radius = sR[nnIdx[0]];
-
- for(int K = 0; K < nbK; K++) {
- BDS_Edge *e = sE[nnIdx[K]];
- e->p1->getTriangles(l1);
- e->p2->getTriangles(l2);
- l.insert(l.begin(), l1.begin(), l1.end());
- l.insert(l.begin(), l2.begin(), l2.end());
- l.sort();
- l.unique();
- }
-
-#else
- {
- l = t;
- }
-#endif
-}
-
-
BDS_Vector::BDS_Vector(const BDS_Point & p2, const BDS_Point & p1)
:x(p2.X - p1.X), y(p2.Y - p1.Y), z(p2.Z - p1.Z)
{
}
-BDS_Vector BDS_Point::N() const
-{
- std::list < BDS_Triangle * >t;
- 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 NORM = sqrt(SumN * SumN);
- SumN /= NORM;
- return SumN;
-}
-
double dist_droites_gauches(BDS_Point * p1, BDS_Point * p2,
BDS_Point * p3, BDS_Point * p4)
{
@@ -226,7 +89,7 @@ double dist_droites_gauches(BDS_Point * p1, BDS_Point * p2,
bool proj_point_triangle(double xa, double ya, double za,
const BDS_Vector & n,
- BDS_Triangle * t, double &x, double &y, double &z)
+ BDS_Face * t, double &x, double &y, double &z)
{
const double eps_prec = 1.e-10;
// BDS_Vector n = t->N();
@@ -234,7 +97,7 @@ bool proj_point_triangle(double xa, double ya, double za,
double b[3];
double res[2];
double det;
- BDS_Point *pts[3];
+ BDS_Point *pts[4];
t->getNodes(pts);
mat[0][0] = pts[1]->X - pts[0]->X;
@@ -307,9 +170,9 @@ void proj_point_plane(double xa, double ya, double za,
z = za + k * n[2];
}
-void print_face(BDS_Triangle * t)
+void print_face(BDS_Face * t)
{
- BDS_Point *pts[3];
+ BDS_Point *pts[4];
t->getNodes(pts);
printf("face %p with nodes %d %d %d and edges %p (%d %d) %p (%d %d) %p (%d %d)\n",
(void *)t, pts[0]->iD, pts[1]->iD, pts[2]->iD, (void *)t->e1,
@@ -387,7 +250,7 @@ double quality_triangle(BDS_Point * p1, BDS_Point * p2, BDS_Point * p3)
return a / (e12 + e22 + e32);
}
-void BDS_Point::getTriangles(std::list < BDS_Triangle * >&t) const
+void BDS_Point::getTriangles(std::list < BDS_Face * >&t) const
{
t.clear();
std::list < BDS_Edge * >::const_iterator it = edges.begin();
@@ -395,10 +258,10 @@ void BDS_Point::getTriangles(std::list < BDS_Triangle * >&t) const
while(it != ite) {
int NF = (*it)->numfaces();
for(int i = 0; i < NF; ++i) {
- BDS_Triangle *tt = (*it)->faces(i);
+ BDS_Face *tt = (*it)->faces(i);
if(tt) {
- std::list < BDS_Triangle * >::iterator tit = t.begin();
- std::list < BDS_Triangle * >::iterator tite = t.end();
+ std::list < BDS_Face * >::iterator tit = t.begin();
+ std::list < BDS_Face * >::iterator tite = t.end();
int found = 0;
while(tit != tite) {
if(tt == *tit)
@@ -429,9 +292,9 @@ BDS_Point *BDS_Mesh::add_point(int num, double u, double v, GFace *gf)
pp->u = u;
pp->v = v;
points.insert(pp);
- double curvature = gf->curvature (SPoint2(u,v));
+ // double curvature = gf->curvature (SPoint2(u,v));
// pp->radius() = (curvature ==0.0) ? 1.e22:1./curvature;
- pp->radius() = curvature;
+ // pp->radius() = curvature;
MAXPOINTNUMBER = (MAXPOINTNUMBER < num) ? num : MAXPOINTNUMBER;
return pp;
}
@@ -515,13 +378,13 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2)
}
int ichoice = ix++ % intersected.size();
- swap_edge ( intersected[ichoice] , BDS_SwapEdgeTestPlanar () );
+ swap_edge ( intersected[ichoice] , BDS_SwapEdgeTestParametric () );
}
return 0;
}
BDS_Edge *BDS_Mesh::find_edge(BDS_Point * p1, BDS_Point * p2,
- BDS_Triangle * t) const
+ BDS_Face * t) const
{
BDS_Point P1(p1->iD);
BDS_Point P2(p2->iD);
@@ -535,12 +398,12 @@ BDS_Edge *BDS_Mesh::find_edge(BDS_Point * p1, BDS_Point * p2,
return 0;
}
-BDS_Triangle *BDS_Mesh::find_triangle(BDS_Edge * e1, BDS_Edge * e2,
+BDS_Face *BDS_Mesh::find_triangle(BDS_Edge * e1, BDS_Edge * e2,
BDS_Edge * e3)
{
int i;
for(i = 0; i < e1->numfaces(); i++) {
- BDS_Triangle *t = e1->faces(i);
+ BDS_Face *t = e1->faces(i);
BDS_Edge *o1 = t->e1;
BDS_Edge *o2 = t->e2;
BDS_Edge *o3 = t->e3;
@@ -553,7 +416,7 @@ BDS_Triangle *BDS_Mesh::find_triangle(BDS_Edge * e1, BDS_Edge * e2,
return t;
}
for(i = 0; i < e2->numfaces(); i++) {
- BDS_Triangle *t = e2->faces(i);
+ BDS_Face *t = e2->faces(i);
BDS_Edge *o1 = t->e1;
BDS_Edge *o2 = t->e2;
BDS_Edge *o3 = t->e3;
@@ -566,7 +429,7 @@ BDS_Triangle *BDS_Mesh::find_triangle(BDS_Edge * e1, BDS_Edge * e2,
return t;
}
for(i = 0; i < e3->numfaces(); i++) {
- BDS_Triangle *t = e3->faces(i);
+ BDS_Face *t = e3->faces(i);
BDS_Edge *o1 = t->e1;
BDS_Edge *o2 = t->e2;
BDS_Edge *o3 = t->e3;
@@ -596,7 +459,7 @@ BDS_Edge *BDS_Mesh::add_edge(int p1, int p2)
return e;
}
-BDS_Triangle *BDS_Mesh::add_triangle(int p1, int p2, int p3)
+BDS_Face *BDS_Mesh::add_triangle(int p1, int p2, int p3)
{
BDS_Edge *e1 = add_edge(p1, p2);
BDS_Edge *e2 = add_edge(p2, p3);
@@ -604,42 +467,46 @@ BDS_Triangle *BDS_Mesh::add_triangle(int p1, int p2, int p3)
return add_triangle(e1, e2, e3);
}
-BDS_Triangle *BDS_Mesh::add_triangle(BDS_Edge * e1, BDS_Edge * e2,
+BDS_Face *BDS_Mesh::add_triangle(BDS_Edge * e1, BDS_Edge * e2,
BDS_Edge * e3)
{
- BDS_Triangle *tfound = find_triangle(e1, e2, e3);
- if(tfound)
- return tfound;
-
- try {
- BDS_Triangle *t = new BDS_Triangle(e1, e2, e3);
- triangles.push_back(t);
- return t;
- }
- catch(...) {
- return 0;
- }
+ //BDS_Face *tfound = find_triangle(e1, e2, e3);
+ // if(tfound)
+ // return tfound;
+
+ BDS_Face *t = new BDS_Face(e1, e2, e3);
+ triangles.push_back(t);
+ return t;
}
-BDS_Tet *BDS_Mesh::add_tet(int p1, int p2, int p3, int p4)
+BDS_Face *BDS_Mesh::add_quadrangle(BDS_Edge * e1, BDS_Edge * e2,
+ BDS_Edge * e3,BDS_Edge * e4 )
{
- BDS_Edge *e1 = add_edge(p1, p2);
- BDS_Edge *e2 = add_edge(p2, p3);
- BDS_Edge *e3 = add_edge(p3, p1);
- BDS_Edge *e4 = add_edge(p1, p4);
- BDS_Edge *e5 = add_edge(p2, p4);
- BDS_Edge *e6 = add_edge(p3, p4);
- BDS_Triangle *t1 = add_triangle(e1, e2, e3);
- BDS_Triangle *t2 = add_triangle(e1, e4, e5);
- BDS_Triangle *t3 = add_triangle(e2, e6, e5);
- BDS_Triangle *t4 = add_triangle(e3, e4, e6);
- BDS_Tet *t = new BDS_Tet(t1, t2, t3, t4);
- tets.push_back(t);
+ BDS_Face *t = new BDS_Face(e1, e2, e3, e4);
+ triangles.push_back(t);
return t;
}
-void BDS_Mesh::del_triangle(BDS_Triangle * t)
+
+// BDS_Tet *BDS_Mesh::add_tet(int p1, int p2, int p3, int p4)
+// {
+// BDS_Edge *e1 = add_edge(p1, p2);
+// BDS_Edge *e2 = add_edge(p2, p3);
+// BDS_Edge *e3 = add_edge(p3, p1);
+// BDS_Edge *e4 = add_edge(p1, p4);
+// BDS_Edge *e5 = add_edge(p2, p4);
+// BDS_Edge *e6 = add_edge(p3, p4);
+// BDS_Face *t1 = add_triangle(e1, e2, e3);
+// BDS_Face *t2 = add_triangle(e1, e4, e5);
+// BDS_Face *t3 = add_triangle(e2, e6, e5);
+// BDS_Face *t4 = add_triangle(e3, e4, e6);
+// BDS_Tet *t = new BDS_Tet(t1, t2, t3, t4);
+// tets.push_back(t);
+// return t;
+// }
+
+void BDS_Mesh::del_triangle(BDS_Face * t)
{
t->e1->del(t);
t->e2->del(t);
@@ -671,7 +538,7 @@ void BDS_Edge::oppositeof(BDS_Point * oface[2]) const
{
oface[0] = oface[1] = 0;
if(faces(0)) {
- BDS_Point *pts[3];
+ BDS_Point *pts[4];
faces(0)->getNodes(pts);
if(pts[0] != p1 && pts[0] != p2)
oface[0] = pts[0];
@@ -681,7 +548,7 @@ void BDS_Edge::oppositeof(BDS_Point * oface[2]) const
oface[0] = pts[2];
}
if(faces(1)) {
- BDS_Point *pts[3];
+ BDS_Point *pts[4];
faces(1)->getNodes(pts);
if(pts[0] != p1 && pts[0] != p2)
oface[1] = pts[0];
@@ -701,11 +568,11 @@ BDS_GeomEntity *BDS_Mesh::get_geom(int p1, int p2)
return (BDS_GeomEntity *) (*it);
}
-void recur_tag(BDS_Triangle * t, BDS_GeomEntity * g)
+void recur_tag(BDS_Face * t, BDS_GeomEntity * g)
{
if(!t->g) {
t->g = g;
- g->t.push_back(t);
+ // g->t.push_back(t);
if(!t->e1->g && t->e1->numfaces() == 2) {
recur_tag(t->e1->otherFace(t), g);
}
@@ -718,664 +585,6 @@ void recur_tag(BDS_Triangle * t, BDS_GeomEntity * g)
}
}
-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_Vector > pts;
- std::list < BDS_Triangle * >::iterator tit = (*it)->t.begin();
- std::list < BDS_Triangle * >::iterator tite = (*it)->t.end();
- BDS_Vector::t = LC / 1.e6;
- while(tit != tite) {
- BDS_Point *nod[3];
- (*tit)->getNodes(nod);
- BDS_Vector p1 = BDS_Vector(nod[0]->X, nod[0]->Y, nod[0]->Z);
- BDS_Vector p2 = BDS_Vector(nod[1]->X, nod[1]->Y, nod[1]->Z);
- BDS_Vector p3 = BDS_Vector(nod[2]->X, nod[2]->Y, nod[2]->Z);
- pts.insert(p1);
- pts.insert(p2);
- pts.insert(p3);
- pts.insert((p1 + p2) * 0.5);
- pts.insert((p1 + p3) * 0.5);
- pts.insert((p2 + p3) * 0.5);
- tit++;
- }
-
- // We consider the following implicit quadrics surface
- // f(x) = x^T A x + b^T x + c = 0
- // A symmetric ( 6 coeffs ), b a vector ( 3 coeffs ) and c a scalar
- // 10 coefficients
- // we try to fit that using least squares
- // we use both points and edges mid points for the fitting
-
- if(pts.size() > 3) {
- // TEST PLANE (easier than quadrics...)
- {
- Double_Matrix PLANE(pts.size(), 3);
- Double_Vector ONES(pts.size());
- Double_Vector RSLT(3);
- std::set < BDS_Vector >::iterator pit = pts.begin();
- std::set < BDS_Vector >::iterator pite = pts.end();
- int k = 0;
- while(pit != pite) {
- PLANE(k, 0) = (*pit).x + (rand() % 1000) * LC / 1.e12;
- PLANE(k, 1) = (*pit).y + (rand() % 1000) * LC / 1.e12;
- PLANE(k, 2) = (*pit).z + (rand() % 1000) * LC / 1.e12;
- ONES(k) = -1;
- k++;
- ++pit;
- }
-
- PLANE.least_squares(ONES, RSLT);
- // printf("%d points plane surface %d ?? : %g %g %g\n",
- // pts.size(),(*it)->classif_tag,RSLT(0),RSLT(1),RSLT(2));
-
- bool plane = true;
- for(unsigned int i = 0; i < pts.size(); i++) {
- const double dist = PLANE(i, 0) * RSLT(0) + PLANE(i, 1) * RSLT(1) +
- PLANE(i, 2) * RSLT(2) + 1;
- if(fabs(dist) >
- 1.e-5 * LC * sqrt(RSLT(0) * RSLT(0) + RSLT(1) * RSLT(1) +
- RSLT(2) * RSLT(2)))
- plane = false;
- }
-
- if(plane) {
- // printf("plane surface %d found : %g %g %g\n",
- // (*it)->classif_tag,RSLT(0),RSLT(1),RSLT(2));
- (*it)->surf = new BDS_Plane(RSLT(0), RSLT(1), RSLT(2));
- }
- }
- if(!(*it)->surf && pts.size() > 20) {
- Double_Matrix QUADRIC(pts.size(), 9);
- Double_Vector ONES(pts.size());
- Double_Vector RSLT(9);
- std::set < BDS_Vector >::iterator pit = pts.begin();
- std::set < BDS_Vector >::iterator pite = pts.end();
- int k = 0;
- while(pit != pite) {
- QUADRIC(k, 0) = (*pit).x * (*pit).x;
- QUADRIC(k, 1) = (*pit).y * (*pit).y;
- QUADRIC(k, 2) = (*pit).z * (*pit).z;
- QUADRIC(k, 3) = 2 * (*pit).x * (*pit).y;
- QUADRIC(k, 4) = 2 * (*pit).x * (*pit).z;
- QUADRIC(k, 5) = 2 * (*pit).y * (*pit).z;
- QUADRIC(k, 6) = (*pit).x;
- QUADRIC(k, 7) = (*pit).y;
- QUADRIC(k, 8) = (*pit).z;
- ONES(k) = 1;
- k++;
- ++pit;
- }
- QUADRIC.least_squares(ONES, RSLT);
- bool quad = true;
- BDS_Quadric temp(RSLT(0), RSLT(1), RSLT(2), RSLT(3), RSLT(4),
- RSLT(5), RSLT(6), RSLT(7), RSLT(8));
- pit = pts.begin();
- while(pit != pite) {
- double x, y, z;
- temp.projection((*pit).x, (*pit).y, (*pit).z, x, y, z);
- const double dist = sqrt(((*pit).x - x) * ((*pit).x - x) +
- ((*pit).y - y) * ((*pit).y - y) +
- ((*pit).z - z) * ((*pit).z - z));
- // printf("%g %g %g ... %g %g %g\n",(*pit).x , (*pit).y , (*pit).z ,x,y,z);
-
- if(dist > 1.e-3 * LC)
- quad = false;
- ++pit;
- }
- if(quad) {
- // printf("quadric surface %d found : %g %g %g %g %g %g %g %g %g\n",
- // (*it)->classif_tag, RSLT(0),RSLT(1),RSLT(2),RSLT(3),RSLT(4),
- // RSLT(5),RSLT(6),RSLT(7),RSLT(8));
- (*it)->surf =
- new BDS_Quadric(RSLT(0), RSLT(1), RSLT(2), RSLT(3), RSLT(4),
- RSLT(5), RSLT(6), RSLT(7), RSLT(8));
- //test
- /*
- FILE *f = fopen ("QUADRIC.pos","w");
- fprintf(f,"View \"quadric\" {\n");
- const int NNN = 20;
- for (int I=0;I<NNN;I++){
- for (int J=0;J<NNN;J++){
- for (int K=0;K<NNN;K++){
- double u1 = (double)I/NNN;
- double v1 = (double)J/NNN;
- double w1 = (double)K/NNN;
- double X1 = Min[0] + u1 * (Max[0]-Min[0]);
- double Y1 = Min[1] + v1 * (Max[1]-Min[1]);
- double Z1 = Min[2] + w1 * (Max[2]-Min[2]);
- double u2 = (double)(I+1)/NNN;
- double v2 = (double)(J+1)/NNN;
- double w2 = (double)(K+1)/NNN;
- double X2 = Min[0] + u2 * (Max[0]-Min[0]);
- double Y2 = Min[1] + v2 * (Max[1]-Min[1]);
- double Z2 = Min[2] + w2 * (Max[2]-Min[2]);
- fprintf(f,"SH(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,"
- "%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g};\n",
- X1,Y1,Z1,X2,Y1,Z1,X2,Y2,Z1,X1,Y2,Z1,X1,Y1,Z2,X2,
- Y1,Z2,X2,Y2,Z2,X1,Y2,Z2,
- (*it)->surf->signedDistanceTo (X1,Y1,Z1),
- (*it)->surf->signedDistanceTo (X2,Y1,Z1),
- (*it)->surf->signedDistanceTo (X2,Y2,Z1),
- (*it)->surf->signedDistanceTo (X1,Y2,Z1),
- (*it)->surf->signedDistanceTo (X1,Y1,Z2),
- (*it)->surf->signedDistanceTo (X2,Y1,Z2),
- (*it)->surf->signedDistanceTo (X2,Y2,Z2),
- (*it)->surf->signedDistanceTo (X1,Y2,Z2));
- }
- }
- }
- fprintf(f,"};\n");
- fclose(f);
- */
- }
- }
- }
- }
- }
-}
-
-void BDS_Mesh::createSearchStructures()
-{
-#ifdef HAVE_ANN_
-
- Msg(INFO, "Creating the ANN search structure");
-
- const double LC_SEARCH = LC * 1e-3;
-
- for(std::set < BDS_GeomEntity *, GeomLessThan >::iterator it = geom.begin();
- it != geom.end(); ++it) {
- if((*it)->classif_degree == 2 && !(*it)->surf) {
- if((*it)->t.size() > 5) {
- int maxPts = 0;
-
- std::set < BDS_Edge * >edg;
-
- std::list < BDS_Triangle * >::iterator tit = (*it)->t.begin();
- std::list < BDS_Triangle * >::iterator tite = (*it)->t.end();
-
- while(tit != tite) {
- edg.insert((*tit)->e1);
- edg.insert((*tit)->e2);
- edg.insert((*tit)->e3);
- tit++;
- }
-
- std::set < BDS_Edge * >::iterator eit = edg.begin();
- std::set < BDS_Edge * >::iterator eite = edg.end();
- while(eit != eite) {
- double l = (*eit)->length();
- maxPts += (int)(l / (LC_SEARCH) + 2);
- eit++;
- }
-
- // printf("%d pts found\n", maxPts);
-
- const int dim = 3;
- (*it)->queryPt = annAllocPt(dim); // allocate query point
- (*it)->dataPts = annAllocPts(maxPts, dim); // allocate data points
- (*it)->nnIdx = new ANNidx[(*it)->nbK]; // allocate near neigh indices
- (*it)->dists = new ANNdist[(*it)->nbK]; // allocate near neighbor dists
- (*it)->sE.resize(maxPts);
- (*it)->sR.resize(maxPts);
- int I = 0;
-
- eit = edg.begin();
- eite = edg.end();
- while(eit != eite) {
- double l = (*eit)->length();
- int N = (int)(l / (LC_SEARCH) + 2);
- BDS_Point *p1 = (*eit)->p1;
- BDS_Point *p2 = (*eit)->p2;
-
- for(int i = 0; i < N; i++) {
- double u = (double)i / (N - 1);
- (*it)->dataPts[I][0] = p1->X + (p2->X - p1->X) * (u);
- (*it)->dataPts[I][1] = p1->Y + (p2->Y - p1->Y) * (u);
- (*it)->dataPts[I][2] = p1->Z + (p2->Z - p1->Z) * (u);
- (*it)->sE[I] = (*eit);
- (*it)->sR[I] =
- p1->radius() + (p2->radius() -
- p1->radius()) * (u);;
- I++;
- }
- eit++;
- }
- // printf("building kd Tree for surface %d -- %d points\n",(*it)->classif_tag,maxPts);
- (*it)->kdTree = new ANNkd_tree( // build search structure
- (*it)->dataPts, // the data points
- maxPts, // number of points
- dim);
- }
- }
- }
-#endif
-}
-
-void BDS_Point::compute_curvature()
-{
- std::list < BDS_Triangle * >t;
- getTriangles(t);
-
- // printf("curvature using %d triangles\n",t.size());
-
- if(g && g->classif_degree != 2) {
- radius() = 1.e22;
- return;
- }
-
- if(t.size() > 4) {
- Double_Matrix M(t.size(), 4);
- Double_Vector Nx(t.size());
- Double_Vector Ny(t.size());
- Double_Vector Nz(t.size());
- Double_Vector Cx(4);
- Double_Vector Cy(4);
- Double_Vector Cz(4);
-
- std::list < BDS_Triangle * >::iterator tit = t.begin();
- std::list < BDS_Triangle * >::iterator tite = t.end();
- int k = 0;
- while(tit != tite) {
- const BDS_Vector cog = (*tit)->cog();
- M(k, 0) = X - cog.x;
- M(k, 1) = Y - cog.y;
- M(k, 2) = Z - cog.z;
- M(k, 3) = 1.0;
- BDS_Vector N = (*tit)->N();
- Nx(k) = N.x;
- Ny(k) = N.y;
- Nz(k) = N.z;
- k++;
- ++tit;
- }
- M.least_squares(Nx, Cx);
- M.least_squares(Ny, Cy);
- M.least_squares(Nz, Cz);
-
- // curvature = divergence of n
-
- double curvature = fabs(Cx(0)) + fabs(Cy(1)) + fabs(Cz(2));
-
- if(curvature != 0.0)
- radius() = fabs(1. / curvature);
- else
- radius() = 1.e22;
- // printf(" R = %g\n",radius());
- }
-}
-
-void compute_curvatures(std::list < BDS_Edge * >&edges)
-{
- {
- std::list < BDS_Edge * >::iterator it = edges.begin();
- std::list < BDS_Edge * >::iterator ite = edges.end();
- while(it != ite) {
- (*it)->target_length = 1.e22;
- (*it)->p1->radius() = 1.e22;
- (*it)->p2->radius() = 1.e22;
- ++it;
- }
-
- }
- {
- std::list < BDS_Edge * >::iterator it = edges.begin();
- std::list < BDS_Edge * >::iterator ite = edges.end();
- while(it != ite) {
- if((*it)->numfaces() == 2) {
- if((*it)->faces(0)->g == (*it)->faces(1)->g) {
- double l1 = 2 * (*it)->faces(0)->inscribed_radius();
- double l2 = 2 * (*it)->faces(1)->inscribed_radius();
- BDS_Vector N1 = (*it)->faces(0)->N();
- BDS_Vector N2 = (*it)->faces(1)->N();
- BDS_Vector DIFFN = N2 - N1;
- BDS_Vector DIST = l1 + l2;
- double crv = 1. / sqrt((DIFFN * DIFFN) / (DIST * DIST));
-
- if((*it)->p1->radius() > crv)
- (*it)->p1->radius() = crv;
- if((*it)->p2->radius() > crv)
- (*it)->p2->radius() = crv;
- }
- }
- ++it;
- }
- }
-}
-
-
-void recur_color_plane_surf(const double eps,
- BDS_Triangle * t,
- const BDS_Vector & n,
- std::set < BDS_Triangle * >&all,
- std::list < BDS_Triangle * >&plane)
-{
- if(all.find(t) == all.end()) {
- BDS_Vector n2 = t->N();
- double angle = fabs(n2.angle(n));
- // printf("triangle %p angle %g\n",t,angle);
- if(angle < eps) {
- all.insert(t);
- plane.push_back(t);
- if(t->e1->numfaces() == 2)
- recur_color_plane_surf(eps, t->e1->otherFace(t), n, all, plane);
- if(t->e2->numfaces() == 2)
- recur_color_plane_surf(eps, t->e2->otherFace(t), n, all, plane);
- if(t->e3->numfaces() == 2)
- recur_color_plane_surf(eps, t->e3->otherFace(t), n, all, plane);
- }
- }
-}
-
-
-
-
-void BDS_Mesh::color_plane_surf(double eps, int NB_T)
-{
- int current_status = 100000;
- std::set < BDS_Triangle * >all;
- std::list < BDS_Triangle * >::iterator it = triangles.begin();
- std::list < BDS_Triangle * >::iterator ite = triangles.end();
- while(it != ite)
- {
- if(all.find(*it) == all.end())
- {
- std::list < BDS_Triangle * >plane;
- BDS_Triangle *start = (BDS_Triangle *) (*it);
- recur_color_plane_surf(eps, start, start->N(), all, plane);
- if((int)plane.size() > NB_T)
- {
- // printf("plane surface found %d triangles\n",plane.size());
- std::list < BDS_Triangle * >::iterator xit = plane.begin();
- std::list < BDS_Triangle * >::iterator xite = plane.end();
- while(xit != xite)
- {
- (*xit)->status = current_status;
- ++xit;
- }
- current_status++;
- }
- }
- ++it;
- }
-}
-
-bool BDS_Mesh::extractVolumes()
-{
- bool closed = true;
- {
- std::list < BDS_Edge * >::iterator it = edges.begin();
- std::list < BDS_Edge * >::iterator ite = edges.end();
- while(it != ite)
- {
- if ( (*it)->numfaces() !=2 ) closed = false;
- ++it;
- }
- }
-
- // printf("the domain is closed ? : %d\n",closed);
- return closed;
-}
-
-
-void BDS_Mesh::classify(double angle, int NB_T)
-{
- // printf(" classifying \n");
-
- static BDS_GeomEntity EDGE_CLAS(0, 1);
- // clear everything
-
- {
- {
- std::set < BDS_Point *, PointLessThan >::iterator it = points.begin();
- std::set < BDS_Point *, PointLessThan >::iterator ite = points.end();
- while(it != ite) {
- (*it)->g = 0;
- ++it;
- }
- }
- {
- std::list < BDS_Edge * >::iterator it = edges.begin();
- std::list < BDS_Edge * >::iterator ite = edges.end();
- while(it != ite) {
- (*it)->g = 0;
- ++it;
- }
- }
- {
- std::list < BDS_Triangle * >::iterator it = triangles.begin();
- std::list < BDS_Triangle * >::iterator ite = triangles.end();
- while(it != ite) {
- (*it)->g = 0;
- if((*it)->status > 10000)
- (*it)->status = 0;
- ++it;
- }
- }
- geom.clear();
- }
-
- color_plane_surf(3.1415/200, 40);
-
- {
- std::list < BDS_Edge * >::iterator it = edges.begin();
- std::list < BDS_Edge * >::iterator ite = edges.end();
- while(it != ite) {
- BDS_Edge & e = *((BDS_Edge *) * it);
- if(e.status == 1) {
- e.g = &EDGE_CLAS;
- }
- else if(e.status == -1) {
- }
- else if(e.numfaces() == 1) {
- e.g = &EDGE_CLAS;
- }
- else if(e.numfaces() == 2 && e.faces(0)->status != e.faces(1)->status) {
- e.g = &EDGE_CLAS;
- }
- else if(e.numfaces() == 2) {
- BDS_Vector N0 = e.faces(0)->N();
- BDS_Vector N1 = e.faces(1)->N();
- double ag = N0.angle(N1);
- if(fabs(ag) > angle) {
- e.g = &EDGE_CLAS;
- }
- }
- else {
- e.g = &EDGE_CLAS;
- }
- ++it;
- }
- }
-
- /*
- {
- std::list<BDS_Triangle*>::iterator it = triangles.begin();
- std::list<BDS_Triangle*>::iterator ite = triangles.end();
- while(it!=ite){
- (*it)->g = 0;
- ++it;
- }
- geom.clear();
- }
- */
-
- // printf(" classifying faces\n");
-
- {
- int tag = 1;
- while(1) {
- 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++;
- }
- // printf(" classifying edges (%d face tags found)\n",tag-1);
- }
-
- int edgetag = 1;
- {
- std::map < std::pair < int, int >, int >edgetags;
- std::list < BDS_Edge * >::iterator it = edges.begin();
- std::list < BDS_Edge * >::iterator ite = edges.end();
- while(it != ite) {
- BDS_Edge & e = *((BDS_Edge *) * it);
-
- if(e.g) {
- int MIN_FAC = 100000;
- int MAX_FAC = -100000;
- std::map < std::pair < int, int >, int >::iterator found;
- BDS_GeomEntity *g;
- if(e.numfaces() == 1) {
- found =
- edgetags.find(std::make_pair(e.faces(0)->g->classif_tag, -1));
- }
- 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));
- }
- else {
- for(int K = 0; K < e.numfaces(); K++) {
- if(MIN_FAC > e.faces(K)->g->classif_tag)
- MIN_FAC = e.faces(K)->g->classif_tag;
- if(MAX_FAC < e.faces(K)->g->classif_tag)
- MAX_FAC = e.faces(K)->g->classif_tag;
- found = edgetags.find(std::make_pair(MIN_FAC, MAX_FAC));
- }
- }
-
- 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;
- }
- else {
- edgetags[std::make_pair(MIN_FAC, MAX_FAC)] = edgetag;
- }
- edgetag++;
- }
- else {
- g = get_geom(found->second, 1);
- }
- e.g = g;
- g->e.push_back(&e);
- }
- }
- else {
- e.g = e.faces(0)->g;
- e.g->e.push_back(&e);
- }
- ++it;
- }
- }
-
- // printf(" classifying points\n");
-
- int vertextag = 1;
- {
- std::set < BDS_Point *, PointLessThan >::iterator it = points.begin();
- std::set < BDS_Point *, PointLessThan >::iterator ite = points.end();
- while(it != ite) {
- std::list < BDS_Edge * >::iterator eit = (*it)->edges.begin();
- std::list < BDS_Edge * >::iterator eite = (*it)->edges.end();
- std::set < BDS_GeomEntity * >geoms;
- BDS_GeomEntity *vg = 0;
- while(eit != eite) {
- if((*eit)->g && (*eit)->g->classif_degree == 1)
- geoms.insert((*eit)->g);
- else
- vg = (*eit)->g;
- ++eit;
- }
- if(geoms.size() == 0) {
- (*it)->g = vg;
- }
- else if(geoms.size() == 1) {
- (*it)->g = (*(geoms.begin()));
- }
- else {
- add_geom(vertextag, 0);
- (*it)->g = get_geom(vertextag++, 0);
- (*it)->g->p = (*it);
- }
- if(!(*it)->g)
- Msg(GERROR, "Error in BDS");
- ++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);
- BDS_Vector SumN = (*it)->N();
- std::list < BDS_Triangle * >::iterator tit = t.begin();
- std::list < BDS_Triangle * >::iterator tite = t.end();
- double max_angle = 0;
- 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)->g->p = (*it);
- }
- }
- ++it;
- }
- }
-
- Msg(INFO, "Computing curvatures");
- compute_curvatures(edges);
- Msg(INFO, "Reverse engineering surfaces");
-
- reverseEngineerCAD();
- Msg(INFO, "Creating search structures");
- createSearchStructures();
- Msg(INFO, "End classifying %d edgetags %d vertextags", edgetag - 1, vertextag - 1);
- // outputScalarField (triangles,"R_curvature.pos");
- extractVolumes();
-}
double PointLessThanLexicographic::t = 0;
double BDS_Vector::t = 0;
@@ -1448,7 +657,7 @@ template < class IT > void DESTROOOY(IT beg, IT end)
void BDS_Mesh::cleanup()
{
{
- for (std::list<BDS_Triangle*> :: iterator it = triangles.begin();
+ for (std::list<BDS_Face*> :: iterator it = triangles.begin();
it != triangles.end();
it++)
{
@@ -1480,7 +689,6 @@ BDS_Mesh ::~ BDS_Mesh ()
cleanup();
DESTROOOY ( edges.begin(),edges.end());
DESTROOOY ( triangles.begin(),triangles.end());
- DESTROOOY ( tets.begin(),tets.end());
}
@@ -1506,7 +714,7 @@ void BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
BDS_Point *p1 = e->p1;
BDS_Point *p2 = e->p2;
- BDS_Point *pts1[3];
+ BDS_Point *pts1[4];
e->faces(0)->getNodes(pts1);
int orientation = 0;
@@ -1540,7 +748,7 @@ void BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
del_triangle(e->faces(0));
}
- double t_l = e->target_length;
+ // double t_l = e->target_length;
del_edge(e);
@@ -1553,26 +761,26 @@ void BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
BDS_Edge *mid_op2 = new BDS_Edge(mid, op[1]);
edges.push_back(mid_op2);
- p1_mid->target_length = t_l;
- op1_mid->target_length = t_l;
- mid_p2->target_length = t_l;
- mid_op2->target_length = t_l;
+ // p1_mid->target_length = t_l;
+ // op1_mid->target_length = t_l;
+ // mid_p2->target_length = t_l;
+ // mid_op2->target_length = t_l;
// 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, *t2, *t3, *t4;
+ BDS_Face *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);
+ t1 = new BDS_Face(op1_mid, p1_op1, p1_mid);
+ t2 = new BDS_Face(mid_op2, op2_p2, mid_p2);
+ t3 = new BDS_Face(op1_p2, op1_mid, mid_p2);
+ t4 = new BDS_Face(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 = new BDS_Face(p1_op1, op1_mid, p1_mid);
+ t2 = new BDS_Face(op2_p2, mid_op2, mid_p2);
+ t3 = new BDS_Face(op1_mid, op1_p2, mid_p2);
+ t4 = new BDS_Face(mid_op2, p1_op2, p1_mid);
}
t1->g = g1;
t2->g = g2;
@@ -1598,45 +806,6 @@ void BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
op[1]->config_modified = true;
}
-bool BDS_SwapEdgeTestNonPlanar::operator () (BDS_Edge *e,
- BDS_Point *p1,BDS_Point *p2,BDS_Point *p3,
- BDS_Point *q1,BDS_Point *q2,BDS_Point *q3) const
-{
- double cb1[3],cb2[3];
- normal_triangle ( p1,p2,p3,cb1);
- normal_triangle ( q1,q2,q3,cb2);
- BDS_Vector n1 = e->faces(0)->N();
- BDS_Vector n2 = e->faces(1)->N();
-
- double psc = n1.x * cb1[0] +n1.y * cb1[1] +n1.z * cb1[2];
- if (psc < 0) return false;
- psc = n2.x * cb1[0] +n2.y * cb1[1] +n2.z * cb1[2];
- if (psc < 0) return false;
- psc = n2.x * cb2[0] +n2.y * cb2[1] +n2.z * cb2[2];
- if (psc < 0) return false;
- psc = n1.x * cb2[0] +n1.y * cb2[1] +n1.z * cb2[2];
- if (psc < 0) return false;
- return true;
-}
-
-bool BDS_SwapEdgeTestPlanar::operator () (BDS_Edge *e,
- BDS_Point *p1,BDS_Point *p2,BDS_Point *p3,
- BDS_Point *q1,BDS_Point *q2,BDS_Point *q3) const
-{
-
-
- double s1 = surface_triangle ( p1,p2,p3);
- double s2 = surface_triangle ( q1,q2,q3);
- double n1 = e->faces(0)->S();
- double n2 = e->faces(1)->S();
-
- if (s1 < 1.e-4 * (n1+n2))return false;
- if (s2 < 1.e-4 * (n1+n2))return false;
- if (fabs(s1+s2 - n1-n2) > 1.e-10 * (s1+s2))return false;
-
- return true;
-}
-
bool BDS_SwapEdgeTestParametric::operator () (BDS_Edge *e,
BDS_Point *p1,BDS_Point *p2,BDS_Point *p3,
BDS_Point *q1,BDS_Point *q2,BDS_Point *q3) const
@@ -1644,8 +813,8 @@ bool BDS_SwapEdgeTestParametric::operator () (BDS_Edge *e,
// return false;
double s1 = fabs(surface_triangle_param ( p1,p2,p3));
double s2 = fabs(surface_triangle_param ( q1,q2,q3));
- BDS_Point * pf1[3];
- BDS_Point * pf2[3];
+ BDS_Point * pf1[4];
+ BDS_Point * pf2[4];
e->faces(0)->getNodes(pf1);
e->faces(1)->getNodes(pf2);
double n1 = fabs(surface_triangle_param ( pf1[0],pf1[1],pf1[2]));
@@ -1697,7 +866,7 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
BDS_GeomEntity *g1 = 0, *g2 = 0, *ge = e->g;
- BDS_Point *pts1[3];
+ BDS_Point *pts1[4];
e->faces(0)->getNodes(pts1);
@@ -1743,14 +912,14 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
BDS_Edge *op1_op2 = new BDS_Edge(op[0], op[1]);
edges.push_back(op1_op2);
- BDS_Triangle *t1, *t2;
+ BDS_Face *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);
+ t1 = new BDS_Face(p1_op1, p1_op2, op1_op2);
+ t2 = new BDS_Face(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 = new BDS_Face(p1_op2, p1_op1, op1_op2);
+ t2 = new BDS_Face(op2_p2, op1_op2, op1_p2);
}
t1->g = g1;
@@ -1769,6 +938,77 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
return true;
}
+int BDS_Edge :: numTriangles() const
+{
+ int NT = 0;
+ for (int i=0;i<_faces.size();i++)
+ if (faces (i)->numEdges() == 3) NT++ ;
+ return NT;
+}
+
+// This function does actually the swap without taking into account
+// the feasability of the operation. Those conditions have to be
+// taken into account before doing the edge swap
+bool BDS_Mesh::recombine_edge(BDS_Edge * e)
+{
+
+ /*
+ p1
+ / | \
+ / | \
+ op1/ 0 | 1 \op2
+ \ | /
+ \ | /
+ \ p2/
+
+ // op1 p1 op2
+ // op1 op2 p2
+ */
+
+ // we test if the edge is deleted
+ // return false;
+ if(e->deleted)
+ return false;
+ if(e->numfaces() != 2 || e->numTriangles() != 2)
+ return false;
+ if(e->g && e->g->classif_degree == 1)
+ return false;
+
+ BDS_Point *op[2];
+ BDS_Point *p1 = e->p1;
+ BDS_Point *p2 = e->p2;
+ e->oppositeof(op);
+
+ BDS_Point *pts1[4];
+ e->faces(0)->getNodes(pts1);
+
+ 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_GeomEntity *g=0;
+ if(e->faces(0)) {
+ g = e->faces(0)->g;
+ del_triangle(e->faces(0));
+ }
+ // not a bug !!!
+ if(e->faces(0)) {
+ del_triangle(e->faces(0));
+ }
+ del_edge(e);
+
+ BDS_Face *f = add_quadrangle (p1_op1, op1_p2, op2_p2, p1_op2);
+ f->g = g;
+
+ p1->config_modified = true;
+ p2->config_modified = true;
+ op[0]->config_modified = true;
+ op[1]->config_modified = true;
+
+ return true;
+}
+
bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
{
@@ -1785,7 +1025,7 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
return false;
}
- std::list < BDS_Triangle * >t;
+ std::list < BDS_Face * >t;
BDS_Point *o = e->othervertex(p);
if(o->g != p->g)
@@ -1804,16 +1044,16 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
int nt = 0;
{
p->getTriangles(t);
- std::list < BDS_Triangle * >::iterator it = t.begin();
- std::list < BDS_Triangle * >::iterator ite = t.end();
+ std::list < BDS_Face * >::iterator it = t.begin();
+ std::list < BDS_Face * >::iterator ite = t.end();
while(it != ite) {
- BDS_Triangle *t = *it;
+ BDS_Face *t = *it;
// print_face(t);
if(t->e1 != e && t->e2 != e && t->e3 != e) {
if (!test_move_point_parametric_triangle ( p, o->u, o-> v, t))
return false;
gs[nt] = t->g;
- BDS_Point *pts[3];
+ BDS_Point *pts[4];
t->getNodes(pts);
pt[0][nt] = (pts[0] == p) ? o->iD : pts[0]->iD;
pt[1][nt] = (pts[1] == p) ? o->iD : pts[1]->iD;
@@ -1825,11 +1065,11 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
{
- std::list < BDS_Triangle * >::iterator it = t.begin();
- std::list < BDS_Triangle * >::iterator ite = t.end();
+ std::list < BDS_Face * >::iterator it = t.begin();
+ std::list < BDS_Face * >::iterator ite = t.end();
while(it != ite) {
- BDS_Triangle *t = *it;
+ BDS_Face *t = *it;
del_triangle(t);
++it;
}
@@ -1856,7 +1096,7 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
{
for(int k = 0; k < nt; k++) {
- BDS_Triangle *t = add_triangle(pt[0][k], pt[1][k], pt[2][k]);
+ BDS_Face *t = add_triangle(pt[0][k], pt[1][k], pt[2][k]);
t->g = gs[k];
}
}
@@ -1871,185 +1111,11 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
return true;
}
-bool project_point_on_a_list_of_triangles(BDS_Point * p,
- const std::list <
- BDS_Triangle * >&t, double &X,
- double &Y, double &Z)
-{
- bool global_ok = false;
-
- int p_classif = p->g->classif_tag;
- std::list < BDS_Triangle * >::const_iterator it = t.begin();
- std::list < BDS_Triangle * >::const_iterator ite = t.end();
- double dist2 = 1.e22;
- double XX = p->X;
- double YY = p->Y;
- double ZZ = p->Z;
- while(it != ite) {
- if((*it)->g->classif_tag == p_classif) {
- {
- double xp, yp, zp;
- bool ok =
- proj_point_triangle(p->X, p->Y, p->Z, (*it)->N(), *it, xp, yp, zp);
- if(ok) {
- global_ok = true;
- double d2 =
- ((xp - X) * (xp - X) + (yp - Y) * (yp - Y) + (zp - Z) * (zp - Z));
- if(d2 < dist2) {
- // printf("one found among %d\n",t.size());
- XX = xp;
- YY = yp;
- ZZ = zp;
- dist2 = d2;
- }
- }
- // ok = proj_point_triangle ( p->X,p->Y,p->Z,p->N(),*it,xp,yp,zp);
- // if(ok){
- // global_ok = true;
- // double d2 = ((xp-X)*(xp-X)+(yp-Y)*(yp-Y)+(zp-Z)*(zp-Z));
- // if(d2 < dist2 ){
- // // printf("one found among %d\n",t.size());
- // XX = xp; YY = yp; ZZ = zp;
- // dist2 = d2;
- // }
- // }
-
- }
- }
- ++it;
- }
-
- X = XX;
- Y = YY;
- Z = ZZ;
- return global_ok;
-}
-
-bool BDS_Mesh::move_point(BDS_Point * p, double X, double Y, double Z)
-{
- std::list < BDS_Triangle * >t;
- p->getTriangles(t);
- std::list < BDS_Triangle * >::iterator tit = t.begin();
- std::list < BDS_Triangle * >::iterator tite = t.end();
- const double oldx = p->X;
- const double oldy = p->Y;
- const double oldz = p->Z;
- while(tit != tite) {
- BDS_Vector n1 = (*tit)->N();
- p->X = X;
- p->Y = Y;
- p->Z = Z;
- BDS_Vector n2 = (*tit)->N_on_the_fly();
- p->X = oldx;
- p->Y = oldy;
- p->Z = oldz;
-
- if (n1*n2 < 0) return false;
-
- double angle = n1.angle(n2);
- if(fabs(angle) > M_PI / 3)
- return false;
- ++tit;
- }
- p->X = X;
- p->Y = Y;
- p->Z = Z;
-
- tit = t.begin();
- while(tit != tite) {
- (*tit)->_update();
- ++tit;
- }
- return true;
-}
-
-void BDS_Mesh::snap_point(BDS_Point * p, BDS_Mesh * geom_mesh)
-{
- double X, Y, Z;
- if(p->g && p->g->surf) {
- p->g->surf->projection(p->X, p->Y, p->Z, X, Y, Z);
- if(move_point(p, X, Y, Z)) {
- SNAP_SUCCESS++;
- }
- else {
- SNAP_FAILURE++;
- }
- }
- else if(p->g && p->g->classif_degree == 2 && geom_mesh) {
- double xx, yy, zz;
- std::list < BDS_Triangle * >l;
- BDS_GeomEntity *gg =
- geom_mesh->get_geom(p->g->classif_tag, p->g->classif_degree);
- gg->getClosestTriangles(p->X, p->Y, p->Z, l, p->radius(), xx,
- yy, zz);
-
- bool ok = project_point_on_a_list_of_triangles(p, l, X, Y, Z);
-
- if(!ok) {
- if(move_point(p, xx, yy, zz)) {
- SNAP_SUCCESS++;
- }
- else {
- SNAP_FAILURE++;
- }
- }
- else {
- if(move_point(p, X, Y, Z)) {
- SNAP_SUCCESS++;
- }
- else {
- SNAP_FAILURE++;
- }
- }
- }
- else {
- return;
- }
-}
-bool BDS_Mesh::smooth_point(BDS_Point * p, BDS_Mesh * geom_mesh)
-{
-
-
-
- if(p->g && p->g->classif_degree <= 1)
- return false;
- // if(!p->g->surf)
- // return false;
-
-
- double X = 0;
- double Y = 0;
- double Z = 0;
-
- std::list < BDS_Edge * >::iterator eit = p->edges.begin();
-
- while(eit != p->edges.end()) {
- if ((*eit)->numfaces() == 1) return false;
- BDS_Point *op = ((*eit)->p1 == p) ? (*eit)->p2 : (*eit)->p1;
- X += op->X;
- Y += op->Y;
- Z += op->Z;
- // lengths are wrong;
- (*eit)->target_length = -1;
- ++eit;
- }
-
- bool success =
- move_point(p, X / p->edges.size(), Y / p->edges.size(),
- Z / p->edges.size());
-
- if(success)
- snap_point(p, geom_mesh);
-
- return true;
-}
-
-
-bool test_move_point_parametric_triangle (BDS_Point * p, double u, double v, BDS_Triangle *t)
+bool test_move_point_parametric_triangle (BDS_Point * p, double u, double v, BDS_Face *t)
{
- BDS_Point *pts[3];
+ BDS_Point *pts[4];
t->getNodes(pts);
double U = p->u;
double V = p->v;
@@ -2089,12 +1155,12 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
V /= p->edges.size();
LC /= p->edges.size();
- std::list < BDS_Triangle * >ts;
+ std::list < BDS_Face * >ts;
p->getTriangles(ts);
- std::list < BDS_Triangle * >::iterator it = ts.begin();
- std::list < BDS_Triangle * >::iterator ite = ts.end();
+ std::list < BDS_Face * >::iterator it = ts.begin();
+ std::list < BDS_Face * >::iterator ite = ts.end();
while(it != ite) {
- BDS_Triangle *t = *it;
+ BDS_Face *t = *it;
if (!test_move_point_parametric_triangle ( p, U, V, t))
return false;
++it;
@@ -2108,7 +1174,7 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
p->X = gp.x();
p->Y = gp.y();
p->Z = gp.z();
- p->radius() = gf->curvature (SPoint2(U,V));
+ // p->radius() = gf->curvature (SPoint2(U,V));
eit = p->edges.begin();
while(eit != p->edges.end()) {
@@ -2120,304 +1186,70 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
return true;
}
-void BDS_Mesh::compute_metric_edge_lengths(const BDS_Metric & metric)
-{
-
-
- // printf("computation of metric lengths\n");
- {
- std::list < BDS_Edge * >::iterator it = edges.begin();
- std::list < BDS_Edge * >::iterator ite = edges.end();
- while(it != ite) {
- if(!(*it)->deleted) {
- (*it)->target_length =
- metric.target_length(0.5 * ((*it)->p2->X + (*it)->p1->X),
- 0.5 * ((*it)->p2->Y + (*it)->p1->Y),
- 0.5 * ((*it)->p2->Z + (*it)->p1->Z));
- BDS_Edge *e = (*it);
- BDS_GeomEntity *g = e->g;
- if (g && g->classif_degree == 1)
- {
- double l = (*it)->length() * 1000;
- if (l < (*it)->target_length)
- (*it)->target_length =l;
- }
- }
- ++it;
- }
- }
-
-
+struct recombine_T
+{
+ const BDS_Edge *e ;
+ double angle ;
+ recombine_T (const BDS_Edge *_e );
+ bool operator < ( const recombine_T & other ) const
{
- std::list < BDS_Edge * >::iterator it = edges.begin();
- std::list < BDS_Edge * >::iterator ite = edges.end();
- while(it != ite) {
- BDS_Edge *e = (*it);
- BDS_GeomEntity *g = e->g;
-
- // do not unrefined curves
- // minimal length is the actual length
-
- if(g && g->surf) {
- double curvature = g->surf->normalCurv(0.5 * (e->p1->X + e->p2->X),
- 0.5 * (e->p1->Y + e->p2->Y),
- 0.5 * (e->p1->Z + e->p2->Z));
- double radius = 1. / curvature;
- double target = radius / metric.nb_elements_per_radius_of_curvature;
- e->target_length =
- metric.update_target_length(target, e->target_length);
- // printf("e1 radius %g target %g length %g mlp %g\n",
- // radius, target,e->length(),e->length()/target);
- }
- else {
- double radius =
- 0.5 * (e->p1->radius() + e->p2->radius());
- double target =
- 3.14159 * radius / metric.nb_elements_per_radius_of_curvature;
- e->target_length =
- metric.update_target_length(target, e->target_length);
- }
- ++it;
- }
- }
- // printf("smoothing\n");
- const int NITER = 3;
- const double BETA = metric.beta;
-
- std::vector<BDS_Point *> temp(points.size());
- std::copy( points.begin(),points.end(),temp.begin());
-
- for(int I = 0; I < NITER; ++I) {
-
- std::vector < BDS_Point * >::const_iterator it = temp.begin();
- std::vector < BDS_Point * >::const_iterator ite = temp.end();
-
- while(it != ite) {
-
- std::list < BDS_Edge * >::iterator eit = (*it)->edges.begin();
- std::list < BDS_Edge * >::iterator eite = (*it)->edges.end();
-
- double l_min = metric._max;
-
- while(eit != eite) {
- BDS_Edge *ee = (*eit);
- if(l_min > ee->target_length)
- l_min = ee->target_length;
- ++eit;
- }
-
- l_min /= BETA;
-
- eit = (*it)->edges.begin();
- while(eit != eite) {
- BDS_Edge *ee = (*eit);
- if(ee->target_length > l_min)
- ee->target_length = l_min;
- ++eit;
- }
- ++it;
- }
+ return angle < other.angle;
}
+};
- // printf("end computation of metric lengths\n");
+recombine_T::recombine_T (const BDS_Edge *_e )
+ : e(_e)
+{
+ BDS_Point *op[2];
+ BDS_Point *p1 = e->p1;
+ BDS_Point *p2 = e->p2;
+ e->oppositeof(op);
+ BDS_Vector v1 (*p1,*op[0]);
+ BDS_Vector v2 (*op[0],*p2);
+ BDS_Vector v3 (*p2,*op[1]);
+ BDS_Vector v4 (*op[1],*p1);
+ angle = fabs(90.-v1.angle_deg(v2));
+ angle = std::max(fabs(90.-v2.angle_deg(v3)),angle);
+ angle = std::max(fabs(90.-v3.angle_deg(v4)),angle);
+ angle = std::max(fabs(90.-v4.angle_deg(v1)),angle);
}
-int BDS_Mesh::adapt_mesh(const BDS_Metric & metric, bool smooth,
- BDS_Mesh * geom_mesh)
+void BDS_Mesh::recombineIntoQuads (const double angle_limit, GFace *gf)
{
- int nb_modif = 0;
- SNAP_SUCCESS = 0;
- SNAP_FAILURE = 0;
-
- std::list < BDS_Edge * >small_to_long(edges);
-
- // split edges
- for(int ii=0;ii<2;ii++){
- std::list < BDS_Edge * >::iterator it = small_to_long.begin();
- std::list < BDS_Edge * >::iterator ite = small_to_long.end();
- compute_metric_edge_lengths(metric);
-
- while(it != ite) {
- if((*it)->numfaces() == 2) {
- double length = (*it)->length();
- if(!(*it)->deleted && length > (*it)->target_length / 0.7) {
- MAXPOINTNUMBER++;
- BDS_Point *mid ;
- double coord = 0.5;
- mid = add_point(MAXPOINTNUMBER,
- coord * (*it)->p1->X + (1 - coord) * (*it)->p2->X,
- coord * (*it)->p1->Y + (1 - coord) * (*it)->p2->Y,
- coord * (*it)->p1->Z + (1 - coord) * (*it)->p2->Z);
- mid->radius() =
- coord * (*it)->p1->radius() + (1 - coord) * (*it)->p2->radius();
- split_edge(*it,mid);
- if(geom_mesh || (mid->g && mid->g->surf))
- snap_point(mid, geom_mesh);
- nb_modif++;
- }
- }
- ++it;
- }
- cleanup();
- small_to_long = edges;
- }
- // collapse
- for(int ii=0;ii<1;ii++){
- std::list<BDS_Edge*>::iterator it = small_to_long.begin();
- std::list<BDS_Edge*>::iterator ite = small_to_long.end();
- // compute_metric_edge_lengths (metric);
+ Msg(INFO,"Recombining triangles for surface %d",gf->tag());
+ for (int i=0;i<5;i++)
+ {
+ bool rec = false;
+ std::set<recombine_T> _pairs;
- while (it != ite)
+ for(std::list < BDS_Edge * >::const_iterator it = edges.begin();
+ it != edges.end(); ++it)
{
- double length = (*it)->length();
- if (!(*it)->deleted && length < (*it)->target_length * 0.7 ){
-// if (nb_modif %2 )
-// {
-// if (collapse_edge (*it, (*it)->p1, 0.1 ))
-// nb_modif++;
-// }
-// else
-// {
-// if (collapse_edge (*it, (*it)->p2, 0.1 ))
-// nb_modif++;
-// }
- }
- ++it;
- }
- cleanup();
- small_to_long = edges;
- }
-
- {
- std::list<BDS_Edge*>::iterator it = small_to_long.begin();
- std::list<BDS_Edge*>::iterator ite = small_to_long.end();
- while (it != ite)
+ if (!(*it)->deleted && (*it)->numfaces () == 2)
+ _pairs.insert ( recombine_T ( *it ) );
+ }
+
+ std::set<recombine_T>::iterator itp = _pairs.begin();
+
+ while (itp != _pairs.end())
{
- if (!(*it)->deleted && (*it)->numfaces()==2)
- {
- BDS_Point *op[2];
- (*it)->oppositeof (op);
-
- 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;
-
- double dd = dist_droites_gauches((*it)->p1 , (*it)->p2,
- op[0],op[1]);
-
- double ll = (*it)->length();
-
- if ((qb > qa && dd < 0.1 * ll) || (qb > 5 * qa))
- {
- nb_modif++;
- swap_edge ( *it , BDS_SwapEdgeTestNonPlanar()) ;
- }
- }
- ++it;
+ rec |= recombine_edge ((BDS_Edge*)itp->e);
+ ++itp;
}
- }
- cleanup();
- if (smooth ){
- Msg(INFO,"smoothing %d points\n",points.size());
+ if (!rec) break;
- std::vector <BDS_Point *> temp_l(points.size());
- std::copy( points.begin(),points.end(),temp_l.begin());
- std::vector < BDS_Point * >::iterator itx = temp_l.begin();
- std::vector < BDS_Point * >::iterator itxe = temp_l.end();
- while (itx != itxe)
+ std::set<BDS_Point*,PointLessThan>::iterator itpt = points.begin();
+ while (itpt != points.end())
{
- smooth_point(*itx,geom_mesh);
- ++itx;
+
+ smooth_point_parametric(*itpt,gf);
+ ++itpt;
}
- printf("coucouc1\n");
- }
-
- Msg(INFO,"%d snaps have succeeded , %d have failed\n",SNAP_SUCCESS,SNAP_FAILURE);
- // outputScalarField (triangles,"b.pos");
- cleanup();
- return nb_modif;
-}
-
-BDS_Mesh::BDS_Mesh(const BDS_Mesh & other)
-{
- MAXPOINTNUMBER = other.MAXPOINTNUMBER;
- LC = other.LC;
- for(int i = 0; i < 3; i++) {
- Min[i] = other.Min[i];
- Max[i] = other.Max[i];
- }
- for(std::set < BDS_GeomEntity *, GeomLessThan >::const_iterator it =
- other.geom.begin(); it != other.geom.end(); ++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 >::const_iterator it =
- other.points.begin(); it != other.points.end(); ++it) {
- BDS_Point *p = add_point((*it)->iD, (*it)->X, (*it)->Y, (*it)->Z);
- p->g =
- ((*it)->g) ? get_geom((*it)->g->classif_tag,
- (*it)->g->classif_degree) : 0;
- p->radius() = (*it)->radius();
- if(p->g->classif_degree == 0)
- p->g->p = p;
- }
- for(std::list < BDS_Edge * >::const_iterator it = other.edges.begin();
- it != other.edges.end(); ++it) {
- BDS_Edge *e = add_edge((*it)->p1->iD, (*it)->p2->iD);
- e->g =
- ((*it)->g) ? get_geom((*it)->g->classif_tag,
- (*it)->g->classif_degree) : 0;
- if(e->g->classif_degree == 1)
- e->g->e.push_back(e);
- e->partition = (*it)->partition;
}
- for(std::list < BDS_Triangle * >::const_iterator it =
- other.triangles.begin(); it != other.triangles.end(); ++it) {
- BDS_Point *n[3];
- (*it)->getNodes(n);
- BDS_Triangle *t = add_triangle(n[0]->iD, n[1]->iD, n[2]->iD);
- t->g = get_geom((*it)->g->classif_tag, (*it)->g->classif_degree);
- t->g->t.push_back(t);
- t->partition = (*it)->partition;
- }
-
- for(std::list < BDS_Tet * >::const_iterator it = other.tets.begin();
- it != other.tets.end(); ++it) {
- BDS_Point *n[4];
- (*it)->getNodes(n);
- BDS_Tet *t = add_tet(n[0]->iD, n[1]->iD, n[2]->iD, n[3]->iD);
- t->g = get_geom((*it)->g->classif_tag, (*it)->g->classif_degree);
- // t->g->t.push_back(t);
- t->partition = (*it)->partition;
- }
-
}
-/*
-void BDS_Mesh::recombineIntoQuads (const double angle)
-{
- for(std::list < BDS_Edge * >::const_iterator it = edges.begin();
- it != edges.end(); ++it)
- {
- BDS_Edge *e = *it;
- if (e->numfaces() == 2
- && !e->faces(0)->q
- && !e->faces(1)->q)
- {
- BDS_Point *op[2];
- e->oppositeof(op);
-
- }
- }
-
-}
-*/
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 83da04f8ff1a37e5b1cf85646c5a5674d4e622db..3ba289e297616828ea0b6f94efd5c70f082b8ca1 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -22,9 +22,6 @@
// points may know the normals to the surface they are classified on
// default values are 0,0,1
-#ifdef HAVE_ANN_
-#include "ANN/ANN.h"
-#endif
#include <string>
#include <set>
#include <map>
@@ -35,10 +32,8 @@
#include "GFace.h"
#include "Views.h"
-class BDS_Tet;
class BDS_Edge;
-class BDS_Triangle;
-class BDS_Quad;
+class BDS_Face;
class BDS_Mesh;
class BDS_Point;
class BDS_Vector;
@@ -50,65 +45,13 @@ double surface_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3);
double surface_triangle_param(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3);
double quality_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3);
-class BDS_Metric
-{
-public:
- const double target,_min,_max,treshold, beta;
- const double nb_elements_per_radius_of_curvature;
- BDS_Metric(double _target , double _mmin, double _mmax, double _b, double cc, double _tres = 0.7)
- : target(_target),_min(_mmin),_max(_mmax),treshold(_tres),beta(_b),nb_elements_per_radius_of_curvature(cc)
- {}
- inline double update_target_length(double _target, double old_target_length) const
- {
- if(_target <= _min) return _min;
- if(_target >= _max && old_target_length > _max) return _max;
- if(old_target_length > _target)return _target ;
- return old_target_length;
- }
- inline double target_length(double x, double y, double z) const
- {
- return target;
- }
-};
-
-class BDS_Surface
-{
-public :
- virtual ~BDS_Surface(){}
- virtual double signedDistanceTo(double x, double y, double z) const = 0;
- virtual void projection(double xa, double ya, double za,
- double &x, double &y, double &z) const =0;
- virtual std::string nameOf() const = 0;
- virtual BDS_Vector Gradient(double x, double y, double z) const = 0;
- virtual double normalCurv(double x, double y, double z) const = 0;
-};
-
-
class BDS_GeomEntity
{
public:
- int nbK;
+
int classif_tag;
int classif_degree;
-
-#ifdef HAVE_ANN_
- ANNpointArray dataPts; // data points
- ANNpoint queryPt; // query point
- ANNidxArray nnIdx; // near neighbor indices
- ANNdistArray dists; // near neighbor distances
- ANNkd_tree* kdTree; // search structure
- std::vector<BDS_Edge *> sE;
- std::vector<double> sR;
-#endif
- std::list<BDS_Triangle *> t;
- std::list<BDS_Edge *> e;
- BDS_Point *p;
- BDS_Surface *surf;
- void getClosestTriangles(double x, double y, double z,
- std::list<BDS_Triangle*> &l ,
- double &radius,
- double &X, double &Y, double &Z);
- inline bool operator < (const BDS_GeomEntity & other) const
+ inline bool operator < (const BDS_GeomEntity & other) const
{
if(classif_degree < other.classif_degree)return true;
if(classif_degree > other.classif_degree)return false;
@@ -116,28 +59,15 @@ public:
return false;
}
BDS_GeomEntity(int a, int b)
- : classif_tag(a),classif_degree(b),p(0),surf(0)
+ : classif_tag(a),classif_degree(b)
{
- nbK=1;
-#ifdef HAVE_ANN_
- kdTree = 0;
-#endif
}
~BDS_GeomEntity()
{
-#ifdef HAVE_ANN_
- if (kdTree)
- {
- annDeallocPts(dataPts);
- delete [] nnIdx; // clean things up
- delete [] dists;
- delete kdTree;
- }
-#endif
}
};
-void print_face(BDS_Triangle *t);
+void print_face(BDS_Face *t);
class BDS_Vector
{
@@ -208,6 +138,10 @@ public:
double ag = atan2(sina,cosa);
return ag;
}
+ double angle_deg(const BDS_Vector &v) const
+ {
+ return angle(v) * 180 / M_PI;
+ }
double operator * (const BDS_Vector &v) const
{
return (x*v.x+y*v.y+z*v.z);
@@ -223,7 +157,7 @@ public:
class BDS_Point
{
- double radius_of_curvature,_lc;
+ double _lc;
public:
double X,Y,Z; // Real COORDINATES
double u,v; // Parametric COORDINATES
@@ -233,11 +167,9 @@ public:
std::list<BDS_Edge*> edges;
// just a transition
- double & radius () {return radius_of_curvature;}
+ double & radius () {return _lc;}
double & lc () {return _lc;}
- BDS_Vector N() const;
-
inline bool operator < (const BDS_Point & other) const
{
return iD < other.iD;
@@ -255,10 +187,10 @@ public:
}
}
double min_edge_length();
- void getTriangles(std::list<BDS_Triangle *> &t) const;
+ void getTriangles(std::list<BDS_Face *> &t) const;
void compute_curvature();
BDS_Point(int id, double x=0, double y=0, double z=0)
- : radius_of_curvature(1.e22),_lc(1.e22),X(x),Y(y),Z(z),u(0),v(0),config_modified(true),iD(id),g(0)
+ : _lc(1.e22),X(x),Y(y),Z(z),u(0),v(0),config_modified(true),iD(id),g(0)
{
}
};
@@ -266,15 +198,12 @@ public:
class BDS_Edge
{
double _length;
- std::vector <BDS_Triangle *> _faces;
+ std::vector <BDS_Face *> _faces;
public:
bool deleted;
- int status;
- int partition;
- double target_length;
BDS_Point *p1,*p2;
BDS_GeomEntity *g;
- inline BDS_Triangle* faces(int i) const
+ inline BDS_Face* faces(int i) const
{
return _faces [i];
}
@@ -286,6 +215,7 @@ public:
{
return _faces.size();
}
+ int numTriangles() const ;
inline BDS_Point * commonvertex(const BDS_Edge *other) const
{
if(p1 == other->p1 || p1 == other->p2) return p1;
@@ -298,7 +228,7 @@ public:
if(p2 == p) return p1;
return 0;
}
- inline void addface(BDS_Triangle *f)
+ inline void addface(BDS_Face *f)
{
_faces.push_back(f);
}
@@ -309,16 +239,16 @@ public:
if(*other.p2 < *p2) return true;
return false;
}
- inline BDS_Triangle * otherFace(const BDS_Triangle *f) const
+ inline BDS_Face * otherFace(const BDS_Face *f) const
{
if(numfaces()!=2) throw;
if(f == _faces[0]) return _faces[1];
if(f == _faces[1]) return _faces[0];
throw;
}
- inline void del(BDS_Triangle *t)
+ inline void del(BDS_Face *t)
{
- _faces.erase(std::remove_if(_faces.begin(),_faces.end() , std::bind2nd(std::equal_to<BDS_Triangle*>(), t)) ,
+ _faces.erase(std::remove_if(_faces.begin(),_faces.end() , std::bind2nd(std::equal_to<BDS_Face*>(), t)) ,
_faces.end());
}
@@ -330,7 +260,7 @@ public:
}
BDS_Edge(BDS_Point *A, BDS_Point *B)
- : deleted(false), status(0),partition(0),target_length(-1.0),g(0)
+ : deleted(false),g(0)
{
if(*A < *B){
p1=A;
@@ -347,248 +277,42 @@ public:
};
-class BDS_Triangle
+class BDS_Face
{
public:
bool deleted;
- int status;
- int partition;
- BDS_Tet *t1,*t2;
- BDS_Edge *e1,*e2,*e3;
- BDS_Quad *q;
- BDS_Vector NORMAL;
- double surface;
- inline BDS_Vector N() const {return NORMAL;}
- inline double S() const {return surface;}
+ BDS_Edge *e1,*e2,*e3,*e4;
BDS_GeomEntity *g;
-
- inline double inscribed_radius() const
- {
- double l1 = e1->length();
- double l2 = e2->length();
- double l3 = e3->length();
- return (2 * S() / (l1+l2+l3));
- }
-
- inline BDS_Tet * opposite_tet(BDS_Tet *t)
- {
- if(t == t1)return t2;
- if(t == t2)return t1;
- throw;
- }
-
- inline int numtets() const
- {
- return ((t1!=0) + (t2 !=0));
- }
-
- inline BDS_Vector cog() const
- {
- BDS_Point *n[3];
- getNodes(n);
- return BDS_Vector((n[0]->X+n[1]->X+n[2]->X)/3.,
- (n[0]->Y+n[1]->Y+n[2]->Y)/3.,
- (n[0]->Z+n[1]->Z+n[2]->Z)/3.);
- }
-
- inline double longest_edge_length() const
- {
- double l1 = e1->length();
- double l2 = e2->length();
- double l3 = e3->length();
-
- if (l1 > l2 && l1 > l3) return l1;
- if (l2 > l1 && l2 > l3) return l2;
- return l3;
- }
-
- inline void _update()
- {
- BDS_Point *pts[3];
- getNodes(pts);
- double c[3];
- vector_triangle(pts[0],pts[1],pts[2],c);
- surface = 0.5 * sqrt(c[0]*c[0]+c[1]*c[1]+c[2]*c[2]);
- NORMAL.x = 2*c[0]/surface;
- NORMAL.y = 2*c[1]/surface;
- NORMAL.z = 2*c[2]/surface;
- }
-
- inline void getNodes(BDS_Point *n[3]) const
- {
- n[0] = e1->commonvertex(e3);
- n[1] = e1->commonvertex(e2);
- n[2] = e2->commonvertex(e3);
- }
-
- inline BDS_Vector N_on_the_fly() const
- {
- double nn[3];
- BDS_Point *pp[3];
- getNodes(pp);
- normal_triangle(pp[0], pp[1], pp[2],nn);
- return BDS_Vector(nn[0],nn[1],nn[2]);
- }
-
- inline void addtet(BDS_Tet *t)
- {
- if(!t1) t1 = t;
- else if(!t2) t2 = t;
- else throw;
- }
-
- inline void del(BDS_Tet *t)
+ inline int numEdges () const {return e4?4:3;}
+ inline void getNodes(BDS_Point *n[4]) const
{
- if(t == t1){
- t1 = t2;
- t2 = 0;
- }
- else if(t == t2){
- t2 = 0;
- }
+ if (!e4)
+ {
+ n[0] = e1->commonvertex(e3);
+ n[1] = e1->commonvertex(e2);
+ n[2] = e2->commonvertex(e3);
+ n[3] = 0;
+ }
else
- throw;
+ {
+ n[0] = e1->commonvertex(e4);
+ n[1] = e1->commonvertex(e2);
+ n[2] = e2->commonvertex(e3);
+ n[3] = e3->commonvertex(e4);
+ }
}
- BDS_Triangle(BDS_Edge *A, BDS_Edge *B, BDS_Edge *C)
- : deleted(false) , status(0), partition(0),t1(0),t2(0),e1(A),e2(B),e3(C),q(0),g(0)
+ BDS_Face(BDS_Edge *A, BDS_Edge *B, BDS_Edge *C,BDS_Edge *D = 0)
+ : deleted(false) ,e1(A),e2(B),e3(C),e4(D),g(0)
{
e1->addface(this);
e2->addface(this);
e3->addface(this);
- _update();
+ if(e4)e4->addface(this);
}
};
-class BDS_Quad
-{
-public:
- BDS_Triangle *t1,*t2;
- inline void getNodes(BDS_Point *n[4]) const
- {
- }
-
- BDS_Quad(BDS_Triangle *A, BDS_Triangle *B)
- : t1(A), t2(B)
- {
- t1->q = this;
- t2->q = this;
- }
-};
-
-
-class BDS_Tet
-{
-public:
- bool deleted;
- int status;
- int partition;
- BDS_Triangle *f1,*f2,*f3,*f4;
- double volume;
- inline double V() const {return volume;}
- BDS_GeomEntity *g;
-
- inline BDS_Vector cog() const
- {
- BDS_Point *n[4];
- getNodes(n);
- return BDS_Vector((n[0]->X+n[1]->X+n[2]->X+n[3]->X)/4.,
- (n[0]->Y+n[1]->Y+n[2]->Y+n[3]->Y)/4.,
- (n[0]->Z+n[1]->Z+n[2]->Z+n[3]->Z)/4.);
- }
-
- inline void _update()
- {
- }
-
- inline void getNodes(BDS_Point *n[4]) const
- {
- BDS_Point *o[3];
- f1->getNodes(n);
- f2->getNodes(o);
- n[3] = 0; //for stupid gcc warning
- if(o[0] != n[0] && o[0] != n[1] &&o[0] != n[2])n[3] = o[0];
- if(o[1] != n[0] && o[1] != n[1] &&o[1] != n[2])n[3] = o[1];
- if(o[2] != n[0] && o[2] != n[1] &&o[2] != n[2])n[3] = o[2];
- }
-
- BDS_Tet(BDS_Triangle *A, BDS_Triangle *B, BDS_Triangle *C, BDS_Triangle *D)
- : deleted(false) , status(0), partition(0),f1(A),f2(B),f3(C),f4(D),g(0)
- {
- f1->addtet(this);
- f2->addtet(this);
- f3->addtet(this);
- f4->addtet(this);
- _update();
- }
-};
-
-
-class BDS_Plane : public BDS_Surface
-{
- double a,b,c,d;
-public :
- BDS_Plane(const double &A, const double &B, const double &C)
- : a(A),b(B),c(C)
- {
- }
- virtual double signedDistanceTo(double x, double y, double z) const {return a*x + b*y + c*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;
- }
- virtual std::string nameOf() const {return std::string("Plane");}
- virtual BDS_Vector Gradient(double x, double y, double z) const
- {
- return BDS_Vector(a , b , c);
- }
- virtual double normalCurv(double x, double y, double z) const
- {
- return 1.e-22;
- }
-};
-
-class BDS_Quadric : public BDS_Surface
-{
-public :
- double a,b,c,d,e,f,g,h,i;
- BDS_Quadric(double A,double B,double C, double D, double E, double F, double G, double H, double I)
- : a(A),b(B),c(C),d(D),e(E),f(F),g(G),h(H),i(I)
- {
- }
-
- virtual BDS_Vector Gradient(double x, double y, double z) const
- {
- return BDS_Vector(2 * (a * x + d * y + e * z) + g ,
- 2 * (d * x + b * y + f * z) + h ,
- 2 * (e * x + f * y + c * z) + i);
- }
-
- virtual double normalCurv(double x, double y, double z) const;
-
- virtual double signedDistanceTo(double x, double y, double z) const {
- const double q =
- a * x * x +
- b * y * y +
- c * z * z +
- 2 * d * x * y +
- 2 * e * x * z +
- 2 * f * y * z +
- g * x +
- h * y +
- i * z - 1.0;
- return q;
- }
- virtual void projection(double xa, double ya, double za,
- double &x, double &y, double &z) const ;
- virtual std::string nameOf() const {return std::string("Quadric");}
-};
-
class GeomLessThan
{
public:
@@ -656,30 +380,18 @@ class BDS_SwapEdgeTestParametric : public BDS_SwapEdgeTest
BDS_Point *q1,BDS_Point *q2,BDS_Point *q3) const;
};
-class BDS_SwapEdgeTestNonPlanar : public BDS_SwapEdgeTest
-{
- public:
- virtual bool operator() (BDS_Edge *e,
- BDS_Point *p1,BDS_Point *p2,BDS_Point *p3,
- BDS_Point *q1,BDS_Point *q2,BDS_Point *q3) const;
-};
-
class BDS_Mesh
{
public:
int MAXPOINTNUMBER,SNAP_SUCCESS,SNAP_FAILURE;
double Min[3],Max[3],LC;
-
- void projection(double &x, double &y, double &z);
-
BDS_Mesh(int _MAXX = 0) : MAXPOINTNUMBER(_MAXX){}
virtual ~BDS_Mesh();
BDS_Mesh(const BDS_Mesh &other);
std::set<BDS_GeomEntity*,GeomLessThan> geom;
std::set<BDS_Point*,PointLessThan> points;
std::list<BDS_Edge*> edges;
- std::list<BDS_Triangle*> triangles;
- std::list<BDS_Tet*> tets;
+ std::list<BDS_Face*> triangles;
// Points
BDS_Point * add_point(int num , double x, double y,double z);
BDS_Point * add_point(int num , double u, double v , GFace *gf);
@@ -689,14 +401,15 @@ public:
BDS_Edge * add_edge(int p1, int p2);
void del_edge(BDS_Edge *e);
BDS_Edge *find_edge(int p1, int p2);
- BDS_Edge *find_edge(BDS_Point *p1, BDS_Point *p2, BDS_Triangle *t)const;
- // Triangles
- BDS_Triangle *add_triangle(int p1, int p2, int p3);
- BDS_Triangle *add_triangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3);
- void del_triangle(BDS_Triangle *t);
- BDS_Triangle *find_triangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3);
- // Tets
- BDS_Tet *add_tet(int p1, int p2, int p3,int p4);
+ BDS_Edge *find_edge(BDS_Point *p1, BDS_Point *p2, BDS_Face *t)const;
+ // Triangles & Quadrangles
+ BDS_Face *add_triangle(int p1, int p2, int p3);
+ BDS_Face *add_quadrangle(int p1, int p2, int p3,int p4);
+ BDS_Face *add_triangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3);
+ BDS_Face *add_quadrangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3,BDS_Edge *e4);
+ void del_triangle(BDS_Face *t);
+ BDS_Face *find_triangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3);
+ BDS_Face *find_quadrangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3,BDS_Edge *e4);
// Geom entities
void add_geom(int degree, int tag);
BDS_GeomEntity *get_geom(int p1, int p2);
@@ -710,21 +423,13 @@ public:
bool move_point(BDS_Point *p , double X, double Y, double Z);
void split_edge(BDS_Edge *, BDS_Point *);
bool edge_constraint ( BDS_Point *p1, BDS_Point *p2 );
+ bool recombine_edge ( BDS_Edge *e );
// Global operators
- void classify(double angle, int nb = -1);
- void color_plane_surf(double eps , int nb);
- void reverseEngineerCAD() ;
- void createSearchStructures() ;
- int adapt_mesh(const BDS_Metric & ,bool smooth = false,BDS_Mesh *geom = 0);
- void compute_metric_edge_lengths(const BDS_Metric & metric);
void cleanup();
- bool extractVolumes();
+ void recombineIntoQuads (const double angle, GFace *gf);
// io's
bool import_view(Post_View *view, const double tolerance);
- void applyOptimizationPatterns();
};
-void outputScalarField(std::list < BDS_Triangle * >t, const char *fn);
-void recur_tag(BDS_Triangle * t, BDS_GeomEntity * g);
-bool project_point_on_a_list_of_triangles(BDS_Point *p , const std::list<BDS_Triangle*> &t,
- double &X, double &Y, double &Z);
+void outputScalarField(std::list < BDS_Face * >t, const char *fn);
+void recur_tag(BDS_Face * t, BDS_GeomEntity * g);
diff --git a/Mesh/DiscreteSurface.cpp b/Mesh/DiscreteSurface.cpp
index 3652282f77610d02861f24588e39ac769cc940bd..81205b446365e23dc5a8661e742639d5051104a0 100644
--- a/Mesh/DiscreteSurface.cpp
+++ b/Mesh/DiscreteSurface.cpp
@@ -1,4 +1,4 @@
-// $Id: DiscreteSurface.cpp,v 1.45 2006-08-22 01:58:34 geuzaine Exp $
+// $Id: DiscreteSurface.cpp,v 1.46 2006-08-29 10:39:48 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -19,449 +19,16 @@
//
// Please report all bugs and problems to "gmsh@geuz.org".
-#include "Gmsh.h"
-#include "Numeric.h"
-#include "Mesh.h"
-#include "CAD.h"
-#include "Geo.h"
-#include "Create.h"
-#include "Interpolation.h"
-#include "Context.h"
-#include "BDS.h"
-#include "PartitionMesh.h"
-#include "OS.h"
-
-extern Mesh *THEM;
-extern Context_T CTX;
-
-void Mesh_To_BDS(Mesh *M)
-{
- Msg(STATUS2, "Moving mesh into new structure");
- Move_SimplexBaseToSimplex(3);
- // create the structure
- if(M->bds)
- delete M->bds;
- M->bds = new BDS_Mesh;
- PhysicalGroup *p;
-
- M->bds->Min[0] = M->bds->Min[1] = M->bds->Min[2] = 1.e12;
- M->bds->Max[0] = M->bds->Max[1] = M->bds->Max[2] = -1.e12;
-
- Msg(STATUS2, "Moving nodes");
- // copy the nodes
- List_T *vertices = Tree2List(M->Vertices);
- for(int i = 0; i < List_Nbr(vertices); ++i) {
- Vertex *v;
- List_Read(vertices, i, &v);
- if(v->Pos.X < M->bds->Min[0])
- M->bds->Min[0] = v->Pos.X;
- if(v->Pos.Y < M->bds->Min[1])
- M->bds->Min[1] = v->Pos.Y;
- if(v->Pos.Z < M->bds->Min[2])
- M->bds->Min[2] = v->Pos.Z;
- if(v->Pos.X > M->bds->Max[0])
- M->bds->Max[0] = v->Pos.X;
- if(v->Pos.Y > M->bds->Max[1])
- M->bds->Max[1] = v->Pos.Y;
- if(v->Pos.Z > M->bds->Max[2])
- M->bds->Max[2] = v->Pos.Z;
- M->bds->add_point(v->Num, v->Pos.X, v->Pos.Y, v->Pos.Z);
- }
- M->bds->LC =
- sqrt((M->bds->Min[0] - M->bds->Max[0]) * (M->bds->Min[0] -
- M->bds->Max[0]) +
- (M->bds->Min[1] - M->bds->Max[1]) * (M->bds->Min[1] -
- M->bds->Max[1]) +
- (M->bds->Min[2] - M->bds->Max[2]) * (M->bds->Min[2] -
- M->bds->Max[2]));
- List_Delete(vertices);
-
- for(int i = 0; i < List_Nbr(M->PhysicalGroups); i++) {
- List_Read(M->PhysicalGroups, i, &p);
- if(p->Typ == MSH_PHYSICAL_POINT) {
- M->bds->add_geom(p->Num, 0);
- BDS_GeomEntity *g = M->bds->get_geom(p->Num, 0);
- for(int j = 0; j < List_Nbr(p->Entities); j++) {
- int Num;
- List_Read(p->Entities, j, &Num);
- BDS_Point *ppp = M->bds->find_point(Num);
- ppp->g = g;
- g->p = ppp;
- }
- }
- }
-
- Msg(STATUS2, "Moving curves");
- List_T *curves = Tree2List(M->Curves);
- for(int i = 0; i < List_Nbr(curves); ++i) {
- Curve *c;
- List_Read(curves, i, &c);
- M->bds->add_geom(c->Num, 1);
- BDS_GeomEntity *g = M->bds->get_geom(c->Num, 1);
- List_T *simplices = Tree2List(c->Simplexes);
- Simplex *simp;
- for(int j = 0; j < List_Nbr(simplices); ++j) {
- List_Read(simplices, j, &simp);
- BDS_Edge *edge = M->bds->add_edge(simp->V[0]->Num, simp->V[1]->Num);
- edge->g = g;
- if(!edge->p1->g)
- edge->p1->g = g;
- if(!edge->p2->g)
- edge->p2->g = g;
- }
- List_Delete(simplices);
- }
- List_Delete(curves);
-
- Msg(STATUS2, "Moving surfaces");
- List_T *surfaces = Tree2List(M->Surfaces);
- for(int i = 0; i < List_Nbr(surfaces); ++i) {
- Surface *s;
- List_Read(surfaces, i, &s);
- M->bds->add_geom(s->Num, 2);
- BDS_GeomEntity *g = M->bds->get_geom(s->Num, 2);
-
- Msg(INFO, "Created new surface %d %d", g->classif_tag, g->classif_degree);
-
- List_T *simplices = Tree2List(s->Simplexes);
- Simplex *simp;
- for(int j = 0; j < List_Nbr(simplices); ++j) {
- List_Read(simplices, j, &simp);
- BDS_Triangle *t =
- M->bds->add_triangle(simp->V[0]->Num, simp->V[1]->Num,
- simp->V[2]->Num);
- t->g = g;
- BDS_Point *n[3];
- t->getNodes(n);
- for(int K = 0; K < 3; K++)
- if(!n[K]->g)
- n[K]->g = g;
- if(!t->e1->g)
- t->e1->g = g;
- if(!t->e2->g)
- t->e2->g = g;
- if(!t->e3->g)
- t->e3->g = g;
- }
- List_Delete(simplices);
- }
- List_Delete(surfaces);
-
- Msg(STATUS2, "Moving %d volumes", Tree_Nbr(M->Volumes));
- List_T *volumes = Tree2List(M->Volumes);
- for(int i = 0; i < List_Nbr(volumes); ++i) {
- Volume *v;
- List_Read(volumes, i, &v);
- M->bds->add_geom(v->Num, 3);
- BDS_GeomEntity *g = M->bds->get_geom(v->Num, 3);
- List_T *simplices = Tree2List(v->Simplexes);
- Simplex *simp;
-
- for(int j = 0; j < List_Nbr(simplices); ++j) {
- List_Read(simplices, j, &simp);
- BDS_Tet *t =
- M->bds->add_tet(simp->V[0]->Num, simp->V[1]->Num, simp->V[2]->Num,
- simp->V[3]->Num);
- t->g = g;
- BDS_Point *n[4];
- t->getNodes(n);
- for(int K = 0; K < 4; K++)
- if(!n[K]->g)
- n[K]->g = g;
- if(!t->f1->g)
- t->f1->g = g;
- if(!t->f2->g)
- t->f2->g = g;
- if(!t->f3->g)
- t->f3->g = g;
- if(!t->f4->g)
- t->f4->g = g;
- if(!t->f1->e1->g)
- t->f1->e1->g = g;
- if(!t->f2->e1->g)
- t->f2->e1->g = g;
- if(!t->f3->e1->g)
- t->f3->e1->g = g;
- if(!t->f4->e1->g)
- t->f4->e1->g = g;
- if(!t->f1->e2->g)
- t->f1->e2->g = g;
- if(!t->f2->e2->g)
- t->f2->e2->g = g;
- if(!t->f3->e2->g)
- t->f3->e2->g = g;
- if(!t->f4->e2->g)
- t->f4->e2->g = g;
- if(!t->f1->e3->g)
- t->f1->e3->g = g;
- if(!t->f2->e3->g)
- t->f2->e3->g = g;
- if(!t->f3->e3->g)
- t->f3->e3->g = g;
- if(!t->f4->e3->g)
- t->f4->e3->g = g;
-
- }
- List_Delete(simplices);
- }
- List_Delete(volumes);
-
-}
-
-void BDS_To_Mesh_2(Mesh *M)
-{
- Msg(STATUS2, "Moving surface mesh into old structure");
-
- Tree_Action(M->Vertices, Free_Vertex);
-
- Tree_Delete(M->Vertices);
- M->Vertices = Tree_Create(sizeof(Vertex *), compareVertex);
-
- {
- std::set < BDS_Point *, PointLessThan >::iterator it =
- M->bds_mesh->points.begin();
- std::set < BDS_Point *, PointLessThan >::iterator ite =
- M->bds_mesh->points.end();
-
-
- while(it != ite) {
- // double dx = 1.e-3 * M->bds_mesh->LC * (double)rand() / (double)RAND_MAX;
- // double dy = 1.e-3 * M->bds_mesh->LC * (double)rand() / (double)RAND_MAX;
- // double dz = 1.e-3 * M->bds_mesh->LC * (double)rand() / (double)RAND_MAX;
- double dx = 0;
- double dy = 0;
- double dz = 0;
- Vertex *vert =
- Create_Vertex((*it)->iD, (*it)->X+dx, (*it)->Y+dy, (*it)->Z+dz, (*it)->min_edge_length(), 0.0);
- Tree_Add(M->Vertices, &vert);
- ++it;
- }
- }
- {
- std::list < BDS_Edge * >::iterator it = M->bds_mesh->edges.begin();
- std::list < BDS_Edge * >::iterator ite = M->bds_mesh->edges.end();
- while(it != ite) {
- BDS_GeomEntity *g = (*it)->g;
- if(g && g->classif_degree == 1) {
- Vertex *v1 = FindVertex((*it)->p1->iD, M);
- Vertex *v2 = FindVertex((*it)->p2->iD, M);
- Simplex *simp = Create_Simplex(v1, v2, NULL, NULL);
- Curve *c = FindCurve(g->classif_tag, M);
- if(c)
- simp->iEnt = g->classif_tag;
- Tree_Insert(c->Simplexes, &simp);
- }
- ++it;
- }
- }
- {
- std::list < BDS_Triangle * >::iterator it =
- M->bds_mesh->triangles.begin();
- std::list < BDS_Triangle * >::iterator ite = M->bds_mesh->triangles.end();
- while(it != ite) {
- BDS_GeomEntity *g = (*it)->g;
- if(g && g->classif_degree == 2) {
- BDS_Point *nod[3];
- (*it)->getNodes(nod);
- Vertex *v1 = FindVertex(nod[0]->iD, M);
- Vertex *v2 = FindVertex(nod[1]->iD, M);
- Vertex *v3 = FindVertex(nod[2]->iD, M);
- Simplex *simp = Create_Simplex(v1, v2, v3, NULL);
- BDS_GeomEntity *g = (*it)->g;
- Surface *s = FindSurface(g->classif_tag, M);
- if(s) {
- simp->iEnt = g->classif_tag;
- }
- else
- Msg(GERROR, "Impossible to find surface %d", g->classif_tag);
- Tree_Add(s->Simplexes, &simp);
- Tree_Add(s->Vertices, &v1);
- Tree_Add(s->Vertices, &v2);
- Tree_Add(s->Vertices, &v3);
- }
- ++it;
- }
- }
- {
- std::list < BDS_Tet * >::iterator it = M->bds_mesh->tets.begin();
- std::list < BDS_Tet * >::iterator ite = M->bds_mesh->tets.end();
- while(it != ite) {
- BDS_Point *nod[4];
- (*it)->getNodes(nod);
- Vertex *v1 = FindVertex(nod[0]->iD, M);
- Vertex *v2 = FindVertex(nod[1]->iD, M);
- Vertex *v3 = FindVertex(nod[2]->iD, M);
- Vertex *v4 = FindVertex(nod[3]->iD, M);
- Simplex *simp = Create_Simplex(v1, v2, v3, v4);
- BDS_GeomEntity *g = (*it)->g;
- Volume *v = FindVolume(g->classif_tag, M);
- if(v) {
- simp->iEnt = g->classif_tag;
- }
- else
- Msg(GERROR, "Error in BDS");
- Tree_Add(v->Simplexes, &simp);
- ++it;
- }
- }
-
- Msg(STATUS2N, " ");
-}
-
-void BDS_To_Mesh(Mesh *M)
-{
- Tree_Action(M->Points, Free_Vertex);
- Tree_Delete(M->Points);
- Tree_Action(M->Volumes, Free_Volume);
- Tree_Action(M->Surfaces, Free_Surface);
- Tree_Action(M->Curves, Free_Curve);
- Tree_Delete(M->Surfaces);
- Tree_Delete(M->Curves);
- Tree_Delete(M->Volumes);
- M->Points = Tree_Create(sizeof(Vertex *), compareVertex);
- M->Curves = Tree_Create(sizeof(Curve *), compareCurve);
- M->Surfaces = Tree_Create(sizeof(Surface *), compareSurface);
- M->Volumes = Tree_Create(sizeof(Volume *), compareVolume);
-
- std::set < BDS_GeomEntity *, GeomLessThan >::iterator it =
- M->bds->geom.begin();
- std::set < BDS_GeomEntity *, GeomLessThan >::iterator ite =
- M->bds->geom.end();
-
- while(it != ite) {
- if((*it)->classif_degree == 3) {
- Volume *_Vol = 0;
- _Vol = FindVolume((*it)->classif_tag, M);
- if(!_Vol)
- _Vol = Create_Volume((*it)->classif_tag, MSH_VOLUME_DISCRETE);
- Tree_Add(M->Volumes, &_Vol);
- }
- else if((*it)->classif_degree == 2) {
- Surface *_Surf = 0;
- _Surf = FindSurface((*it)->classif_tag, M);
- if(!_Surf)
- _Surf = Create_Surface((*it)->classif_tag, MSH_SURF_DISCRETE);
- // _Surf->bds = M->bds;
- End_Surface(_Surf);
- Tree_Add(M->Surfaces, &_Surf);
- }
- else if((*it)->classif_degree == 1) {
- Curve *_c = 0;
- _c = FindCurve((*it)->classif_tag, M);
- if(!_c)
- _c =
- Create_Curve((*it)->classif_tag, MSH_SEGM_DISCRETE, 1, NULL, NULL,
- -1, -1, 0., 1.);
- // _c->bds = M->bds;
- End_Curve(_c);
- Tree_Add(M->Curves, &_c);
- }
- else if((*it)->classif_degree == 0) {
- BDS_Point *p = (*it)->p;
- if(p) {
- Vertex *_v = Create_Vertex(p->iD, p->X, p->Y, p->Z, 1, 0);
- Tree_Add(M->Points, &_v);
- }
- }
- ++it;
- }
-
- CTX.mesh.changed = ENT_SURFACE;
-}
-
-
-void CreateVolumeWithAllSurfaces(Mesh *M)
-{
- // Volume *vol = Create_Volume(1, MSH_VOLUME_DISCRETE);
- // Volume *vol = Create_Volume(1, 99999);
- Volume *vol2 = Create_Volume(2, MSH_VOLUME);
- List_T *surfaces = Tree2List(M->Surfaces);
- for(int i = 0; i < List_Nbr(surfaces); ++i) {
- Surface *s;
- List_Read(surfaces, i, &s);
- //List_Add (vol->Surfaces,&s);
- List_Add (vol2->Surfaces,&s);
- int ori = 1;
- Move_SimplexBaseToSimplex(&s->SimplexesBase, s->Simplexes);
- List_Add (vol2->SurfacesOrientations,&ori);
- }
- // Tree_Add(M->Volumes, &vol);
- Tree_Add(M->Volumes, &vol2);
-}
-
-int ReMesh()
-{
- if(THEM->status != 2)
- return 0;
-
- if(!THEM->bds) {
- Mesh_To_BDS(THEM);
- THEM->bds->classify(CTX.mesh.dihedral_angle_tol * M_PI / 180);
- BDS_To_Mesh(THEM);
- }
-
- DeleteMesh(THEM);
-
- if(THEM->bds_mesh) {
- delete THEM->bds_mesh;
- THEM->bds_mesh = 0;
- }
-
-
- MeshDiscreteSurface((Surface *) 0);
- CreateVolumeWithAllSurfaces(THEM);
- CTX.mesh.changed = ENT_SURFACE;
- return 1;
-}
-
// Public interface for discrete surface/curve mesh algo
+#include "Mesh.h"
+
int MeshDiscreteSurface(Surface * s)
{
- static int NITER = 10;
- if(THEM->bds) {
- BDS_Metric metric(THEM->bds->LC / CTX.mesh.target_elem_size_fact,
- THEM->bds->LC / CTX.mesh.min_elem_size_fact,
- THEM->bds->LC,
- CTX.mesh.beta_smooth_metric, CTX.mesh.nb_elem_per_rc);
- if(!THEM->bds_mesh) {
- Msg(STATUS1, "Remeshing 2D...");
- double t1 = Cpu();
-
- THEM->bds_mesh = new BDS_Mesh(*(THEM->bds));
- int iter = 0;
- while(iter < NITER && THEM->bds_mesh->adapt_mesh(metric, true, THEM->bds)) {
- Msg(STATUS2, "Iter=%d/%d Tri=%d", iter, NITER, THEM->bds_mesh->triangles.size());
- iter++;
- }
- BDS_To_Mesh_2(THEM);
- Msg(INFO, "Mesh has %d vertices (%d)", Tree_Nbr(THEM->Vertices),
- THEM->bds->points.size());
- // THEM->bds_mesh->save_gmsh_format("3.msh");
-
- double t2 = Cpu();
- Msg(INFO, "Remesh 2D complete (%g s)", t2 - t1);
- // NITER++;
- return 1;
- }
- return 2;
- }
- else if(s->Typ == MSH_SURF_DISCRETE) {
- // nothing else to do: we assume that the surface is represented
- // by a mesh that will not be modified
- return 1;
- }
- else
- return 0;
+ throw;
}
int MeshDiscreteCurve(Curve * c)
{
- if(c->Typ == MSH_SEGM_DISCRETE) {
- // nothing else to do: we assume that the curve is represented by
- // a mesh that will not be modified
- return 1;
- }
- else
- return 0;
+ throw;
}
diff --git a/Mesh/PartitionMesh.cpp b/Mesh/PartitionMesh.cpp
index 3cb9c976c709f1f5c16ec60c1f727ce083200a84..25e4c0ee6fd2ebc385cee16c266bdddcc16739f5 100644
--- a/Mesh/PartitionMesh.cpp
+++ b/Mesh/PartitionMesh.cpp
@@ -1,4 +1,4 @@
-// $Id: PartitionMesh.cpp,v 1.8 2006-01-16 17:55:44 geuzaine Exp $
+// $Id: PartitionMesh.cpp,v 1.9 2006-08-29 10:39:49 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -19,203 +19,15 @@
//
// Please report all bugs and problems to "gmsh@geuz.org".
-#include "Gmsh.h"
-#include "Numeric.h"
#include "Mesh.h"
-#include "CAD.h"
-#include "Geo.h"
-#include "Create.h"
-#include "Interpolation.h"
-#include "Context.h"
-#include "BDS.h"
#include "PartitionMesh.h"
-#include "OpenFile.h"
-
-#ifdef HAVE_METIS
-extern "C"
-{
-#include "metis.h"
-}
-#endif
-
-extern void Mesh_To_BDS(Mesh * m);
-extern void BDS_To_Mesh_2(Mesh * m);
-extern void BDS_To_Mesh(Mesh * m);
-
-void DeleteMesh(Mesh * M)
-{
- List_T *Curves = Tree2List(M->Curves);
- for(int i = 0; i < List_Nbr(Curves); i++) {
- Curve *c;
- List_Read(Curves, i, &c);
- Tree_Action(c->Simplexes, Free_Simplex);
- Tree_Delete(c->Simplexes);
- Tree_Action(c->SimplexesBase, Free_SimplexBase);
- Tree_Delete(c->SimplexesBase);
- c->Simplexes = Tree_Create(sizeof(Simplex *), compareSimplex);
- c->SimplexesBase = Tree_Create(sizeof(SimplexBase *), compareSimplexBase);
- }
- List_T *Surfaces = Tree2List(M->Surfaces);
- for(int i = 0; i < List_Nbr(Surfaces); i++) {
- Surface *s;
- List_Read(Surfaces, i, &s);
- Tree_Action(s->Simplexes, Free_Simplex);
- Tree_Delete(s->Simplexes);
- Tree_Delete(s->Vertices);
- Tree_Action(s->SimplexesBase, Free_SimplexBase);
- Tree_Delete(s->SimplexesBase);
- s->Simplexes = Tree_Create(sizeof(Simplex *), compareSimplex);
- s->Vertices = Tree_Create(sizeof(Simplex *), compareVertex);
- s->SimplexesBase = Tree_Create(sizeof(SimplexBase *), compareSimplexBase);
- }
- List_Delete(Surfaces);
-
- List_T *Volumes = Tree2List(M->Volumes);
- for(int i = 0; i < List_Nbr(Volumes); i++) {
- Volume *v;
- List_Read(Volumes, i, &v);
- Tree_Action(v->Simplexes, Free_Simplex);
- Tree_Delete(v->Simplexes);
- Tree_Delete(v->Vertices);
- Tree_Action(v->SimplexesBase, Free_SimplexBase);
- Tree_Delete(v->SimplexesBase);
- v->Simplexes = Tree_Create(sizeof(Simplex *), compareSimplex);
- v->Vertices = Tree_Create(sizeof(Simplex *), compareVertex);
- v->SimplexesBase = Tree_Create(sizeof(SimplexBase *), compareSimplexBase);
- }
- List_Delete(Volumes);
-}
void PartitionMesh(Mesh * M, int NP)
{
- Msg(INFO, "moving the mesh to BDS");
- Mesh_To_BDS(M);
- BDS_Mesh *m = M->bds;
- Msg(INFO, "Partitioning");
- PartitionMesh(m, NP);
- Msg(INFO, "Moving back to the old data str");
- DeleteMesh(M);
- //BDS_To_Mesh(M);
- M->bds_mesh = new BDS_Mesh(*M->bds);
- BDS_To_Mesh_2(M);
- delete M->bds;
- //delete M->bds_mesh;
- M->bds = 0;
- //M->bds_mesh = 0;
- SetBoundingBox();
+ throw;
}
void PartitionMesh(BDS_Mesh * m, int NP)
{
-#ifdef HAVE_METIS
- //NN = number of nodes of the graph
- int dim = (m->tets.size() == 0) ? 2 : 3;
- int NN = (dim == 2) ? m->triangles.size() : m->tets.size();
-
- Msg(INFO, "%d nodes in the graph", NN);
-
- int *partitionVector = new int[NN];
- int *xadj = new int[NN + 2];
-
- int totCount = 0;
-
- std::list< BDS_Triangle*>::iterator it2 = m->triangles.begin();
- std::list< BDS_Tet*>::iterator it3 = m->tets.begin();
-
- xadj[0] = 0;
- for(int i = 0; i < NN; i++) {
- int nbAdj = 0;
- if(dim == 2) {
- BDS_Triangle *t = *it2;
- t->partition = i;
- ++it2;
- nbAdj = (t->e1->numfaces() + t->e2->numfaces() + t->e3->numfaces() - 3);
- totCount += nbAdj;
- }
- else if(dim == 3) {
- BDS_Tet *t = *it3;
- t->partition = i;
- ++it3;
- nbAdj = (t->f1->numtets() + t->f2->numtets() +
- t->f3->numtets() + t->f4->numtets() - 4);
- totCount += nbAdj;
- }
- xadj[i + 1] = xadj[i] + nbAdj;
- }
-
- Msg(INFO, "Tot Count %d", totCount);
-
- it2 = m->triangles.begin();
- it3 = m->tets.begin();
-
- int *adjncy = new int[totCount + 1];
-
- int count = 0;
-
- for(int i = 0; i < NN; i++) {
- if(dim == 2) {
- BDS_Triangle *t = *it2;
- for(int j = 0; j < t->e1->numfaces(); j++) {
- BDS_Triangle *f = t->e1->faces(j);
- if(f != t)
- adjncy[count++] = f->partition;
- }
- for(int j = 0; j < t->e2->numfaces(); j++) {
- BDS_Triangle *f = t->e2->faces(j);
- if(f != t)
- adjncy[count++] = f->partition;
- }
- for(int j = 0; j < t->e3->numfaces(); j++) {
- BDS_Triangle *f = t->e3->faces(j);
- if(f != t)
- adjncy[count++] = f->partition;
- }
- ++it2;
- }
- else if(dim == 3) {
- BDS_Tet *t = *it3;
- BDS_Tet *o = t->f1->opposite_tet(t);
- if(o)
- adjncy[count++] = o->partition;
- o = t->f2->opposite_tet(t);
- if(o)
- adjncy[count++] = o->partition;
- o = t->f3->opposite_tet(t);
- if(o)
- adjncy[count++] = o->partition;
- o = t->f4->opposite_tet(t);
- if(o)
- adjncy[count++] = o->partition;
- ++it3;
- }
- }
-
- int wgtflag = 0;
- int numflag = 0;
- int options[4];
- options[0] = 0;
- int edgecut;
- METIS_PartGraphRecursive(&NN, xadj, adjncy, 0, 0, &wgtflag,
- &numflag, &NP, options, &edgecut, partitionVector);
- delete[]xadj;
- delete[]adjncy;
-
- it2 = m->triangles.begin();
- it3 = m->tets.begin();
-
- for(int i = 0; i < NN; i++) {
- if(dim == 2) {
- BDS_Triangle *t = *it2;
- t->partition = partitionVector[i];
- ++it2;
- }
- else if(dim == 3) {
- BDS_Tet *t = *it3;
- t->partition = partitionVector[i];
- ++it3;
- }
- }
-
- delete[]partitionVector;
-#endif
+ throw;
}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 9438f876fa2245db1744bace53cb9282c27b6326..2812775edba75bd66e68b9ac1c98ebfe653e4b2c 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.8 2006-08-26 15:13:22 remacle Exp $
+// $Id: meshGFace.cpp,v 1.9 2006-08-29 10:39:49 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -172,49 +172,6 @@ void computeEdgeLoops (const GFace *gf,
}
-
-
-void computeEdgeParameters ( double x1, double y1, double x2, double y2, GFace *gf , const int numberOfTestPoints, double &coordMiddle, double &edgeLength )
-{
- std::vector<GPoint> pts;
- for (int i=0;i<numberOfTestPoints;++i)
- {
- double xi = (double)i/(double)(numberOfTestPoints-1);
- double u = x1 * (1-xi) + x2 * xi;
- double v = y1 * (1-xi) + y2 * xi;
- pts.push_back(gf->point(u,v));
- }
- edgeLength = 0;
- for (int i=1;i<numberOfTestPoints;++i)
- {
- GPoint p1 = pts[i-1];
- GPoint p2 = pts[i];
- edgeLength += 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()) );
- }
- double ll = 0;
- for (int i=1;i<numberOfTestPoints;++i)
- {
- double xi = (double)(i-1)/(double)(numberOfTestPoints-1);
- GPoint p1 = pts[i-1];
- GPoint p2 = pts[i];
-
- double oldll = ll;
-
- ll += 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()) );
-
- if (oldll <= 0.5*edgeLength && ll >= 0.5*edgeLength)
- {
- double xi2 = (0.5*edgeLength - oldll)/(ll-oldll);
- coordMiddle = xi + xi2/ (double)(numberOfTestPoints-1);
- return;
- }
- }
-}
-
extern double F_LC_ANALY (double xx, double yy, double zz);
double NewGetLc ( BDS_Point *p )
@@ -303,7 +260,7 @@ void OptimizeMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
std::set<BDS_Point*,PointLessThan>::iterator itp = PTS.begin();
while (itp != PTS.end())
{
- std::list < BDS_Triangle * >t;
+ std::list < BDS_Face * >t;
(*itp)->getTriangles(t);
if ((t.size()==3 && (*itp)->edges.size() == 3)||
(t.size()==4 && (*itp)->edges.size() == 4))
@@ -368,12 +325,12 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
while (1)
{
// double stot = 0;
-// std::list<BDS_Triangle *>::iterator ittt = m.triangles.begin();
+// std::list<BDS_Face *>::iterator ittt = m.triangles.begin();
// while (ittt!= m.triangles.end())
// {
// if (!(*ittt)->deleted)
// {
-// BDS_Point *pts[3];
+// BDS_Point *pts[4];
// (*ittt)->getNodes(pts);
// stot += fabs( surface_triangle_param(pts[0], pts[1], pts[2]));
// }
@@ -522,7 +479,6 @@ void OldMeshGenerator ( GFace *gf,
extern PointRecord *gPointArray;
-
// FIX THAT !!
double LC2D = 1;
if (gf->geomType () == GEntity::Plane)
@@ -711,8 +667,6 @@ bool recover_medge ( BDS_Mesh *m, GEdge *ge)
void gmsh2DMeshGenerator ( GFace *gf )
{
- // if(gf->tag() != 136) return;
-
std::set<MVertex*>all_vertices;
std::map<int, MVertex*>numbered_vertices;
std::list<GEdge*> edges = gf->edges();
@@ -911,10 +865,10 @@ void gmsh2DMeshGenerator ( GFace *gf )
// delete useless stuff
{
- std::list<BDS_Triangle*>::iterator itt = m->triangles.begin();
+ std::list<BDS_Face*>::iterator itt = m->triangles.begin();
while (itt != m->triangles.end())
{
- BDS_Triangle *t = *itt;
+ BDS_Face *t = *itt;
if (!t->g)
m->del_triangle (t);
++itt;
@@ -948,7 +902,13 @@ void gmsh2DMeshGenerator ( GFace *gf )
// start mesh generation
RefineMesh (gf,*m,20);
- OptimizeMesh (gf,*m,2);
+ // OptimizeMesh (gf,*m,2);
+
+
+ if (gf->meshAttributes.recombine)
+ {
+ m->recombineIntoQuads (gf->meshAttributes.recombineAngle,gf);
+ }
// fill the small gmsh structures
@@ -959,32 +919,40 @@ void gmsh2DMeshGenerator ( GFace *gf )
BDS_Point *p = *itp;
if (numbered_vertices.find(p->iD) == numbered_vertices.end())
{
- //MVertex *v = new MFaceVertex (p->X,p->Y,p->Z,gf,p->u,p->v);
- MVertex *v = new MFaceVertex (p->u,p->v,0,gf,p->u,p->v);
+ MVertex *v = new MFaceVertex (p->X,p->Y,p->Z,gf,p->u,p->v);
+ //MVertex *v = new MFaceVertex (p->u,p->v,0,gf,p->u,p->v);
numbered_vertices[p->iD]=v;
gf->mesh_vertices.push_back(v);
}
++itp;
}
}
+
{
- std::list<BDS_Triangle*>::iterator itt = m->triangles.begin();
+ std::list<BDS_Face*>::iterator itt = m->triangles.begin();
while (itt != m->triangles.end())
{
- BDS_Triangle *t = *itt;
- BDS_Point *n[3];
- t->getNodes(n);
+ BDS_Face *t = *itt;
if (!t->deleted)
{
+ BDS_Point *n[4];
+ t->getNodes(n);
MVertex *v1 = numbered_vertices[n[0]->iD];
MVertex *v2 = numbered_vertices[n[1]->iD];
MVertex *v3 = numbered_vertices[n[2]->iD];
- gf->triangles.push_back(new MTriangle (v1,v2,v3) );
+ if (!n[3])
+ gf->triangles.push_back(new MTriangle (v1,v2,v3) );
+ else
+ {
+ MVertex *v4 = numbered_vertices[n[3]->iD];
+ gf->quadrangles.push_back(new MQuadrangle (v1,v2,v3,v4) );
+ }
}
++itt;
}
}
+
// delete the mesh
// char name[245];
@@ -1008,9 +976,9 @@ void gmsh2DMeshGenerator ( GFace *gf )
delete [] Y_;
delete [] Z_;
- fromParametricToCartesian p2c ( gf );
+ // fromParametricToCartesian p2c ( gf );
// std::for_each(all_vertices.begin(),all_vertices.end(),p2c);
- std::for_each(gf->mesh_vertices.begin(),gf->mesh_vertices.end(),p2c);
+ // std::for_each(gf->mesh_vertices.begin(),gf->mesh_vertices.end(),p2c);
}
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index fda11e2ec70b17cf704726e412ae6fa2a5bb73e7..bc0497f1a27c57b3d49ede24ec77c8885a57bd24 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -300,7 +300,7 @@
/* Copy the first part of user declarations. */
#line 1 "Gmsh.y"
-// $Id: Gmsh.tab.cpp,v 1.273 2006-08-19 08:26:47 remacle Exp $
+// $Id: Gmsh.tab.cpp,v 1.274 2006-08-29 10:39:49 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -837,27 +837,27 @@ static const unsigned short int yyrline[] =
1414, 1439, 1458, 1477, 1493, 1513, 1530, 1547, 1568, 1573,
1578, 1583, 1591, 1592, 1593, 1598, 1601, 1605, 1621, 1637,
1653, 1674, 1688, 1694, 1700, 1712, 1721, 1731, 1745, 1764,
- 1778, 1784, 1790, 1799, 1813, 1858, 1873, 1884, 1903, 1913,
- 1935, 1939, 1944, 1949, 1961, 1976, 1992, 2018, 2045, 2077,
- 2084, 2089, 2095, 2099, 2108, 2116, 2124, 2133, 2132, 2146,
- 2145, 2159, 2158, 2173, 2180, 2187, 2194, 2201, 2208, 2215,
- 2222, 2229, 2237, 2236, 2249, 2248, 2261, 2260, 2273, 2272,
- 2285, 2284, 2297, 2296, 2309, 2308, 2321, 2320, 2333, 2332,
- 2348, 2351, 2357, 2383, 2407, 2416, 2434, 2452, 2470, 2499,
- 2534, 2561, 2588, 2602, 2621, 2627, 2633, 2636, 2643, 2647,
- 2657, 2658, 2659, 2660, 2661, 2662, 2663, 2664, 2665, 2672,
- 2673, 2674, 2675, 2676, 2677, 2678, 2679, 2680, 2681, 2682,
- 2683, 2684, 2685, 2686, 2687, 2688, 2689, 2690, 2691, 2692,
- 2693, 2694, 2695, 2696, 2697, 2698, 2699, 2700, 2701, 2702,
- 2703, 2705, 2706, 2707, 2708, 2709, 2710, 2711, 2712, 2713,
- 2714, 2715, 2716, 2717, 2718, 2719, 2720, 2721, 2722, 2723,
- 2724, 2725, 2734, 2735, 2736, 2737, 2738, 2739, 2740, 2744,
- 2760, 2775, 2795, 2809, 2822, 2845, 2863, 2881, 2899, 2917,
- 2925, 2929, 2933, 2937, 2941, 2948, 2952, 2956, 2960, 2968,
- 2970, 2974, 2981, 2986, 2994, 2999, 3003, 3008, 3012, 3024,
- 3030, 3041, 3061, 3071, 3081, 3091, 3108, 3127, 3151, 3180,
- 3185, 3189, 3193, 3206, 3210, 3222, 3229, 3251, 3255, 3270,
- 3275, 3282, 3286, 3294, 3302, 3316, 3330, 3334, 3353, 3375
+ 1778, 1784, 1790, 1799, 1813, 1858, 1873, 1884, 1904, 1914,
+ 1936, 1940, 1945, 1950, 1962, 1977, 1993, 2019, 2046, 2078,
+ 2085, 2090, 2096, 2100, 2109, 2117, 2125, 2134, 2133, 2147,
+ 2146, 2160, 2159, 2174, 2181, 2188, 2195, 2202, 2209, 2216,
+ 2223, 2230, 2238, 2237, 2250, 2249, 2262, 2261, 2274, 2273,
+ 2286, 2285, 2298, 2297, 2310, 2309, 2322, 2321, 2334, 2333,
+ 2349, 2352, 2358, 2384, 2408, 2417, 2435, 2453, 2471, 2500,
+ 2535, 2562, 2589, 2603, 2622, 2628, 2634, 2637, 2644, 2648,
+ 2658, 2659, 2660, 2661, 2662, 2663, 2664, 2665, 2666, 2673,
+ 2674, 2675, 2676, 2677, 2678, 2679, 2680, 2681, 2682, 2683,
+ 2684, 2685, 2686, 2687, 2688, 2689, 2690, 2691, 2692, 2693,
+ 2694, 2695, 2696, 2697, 2698, 2699, 2700, 2701, 2702, 2703,
+ 2704, 2706, 2707, 2708, 2709, 2710, 2711, 2712, 2713, 2714,
+ 2715, 2716, 2717, 2718, 2719, 2720, 2721, 2722, 2723, 2724,
+ 2725, 2726, 2735, 2736, 2737, 2738, 2739, 2740, 2741, 2745,
+ 2761, 2776, 2796, 2810, 2823, 2846, 2864, 2882, 2900, 2918,
+ 2926, 2930, 2934, 2938, 2942, 2949, 2953, 2957, 2961, 2969,
+ 2971, 2975, 2982, 2987, 2995, 3000, 3004, 3009, 3013, 3025,
+ 3031, 3042, 3062, 3072, 3082, 3092, 3109, 3128, 3152, 3181,
+ 3186, 3190, 3194, 3207, 3211, 3223, 3230, 3252, 3256, 3271,
+ 3276, 3283, 3287, 3295, 3303, 3317, 3331, 3335, 3354, 3376
};
#endif
@@ -5380,7 +5380,8 @@ yyreduce:
SleepInSeconds((yyvsp[-1].d));
}
else if(!strcmp((yyvsp[-2].c), "Remesh")){
- ReMesh();
+ Msg(GERROR, "Surface ReMeshing must be reinterfaced");
+ // ReMesh();
}
else if(!strcmp((yyvsp[-2].c), "Mesh")){
yymsg(GERROR, "Mesh directives are not (yet) allowed in scripts");
@@ -5396,7 +5397,7 @@ yyreduce:
break;
case 138:
-#line 1904 "Gmsh.y"
+#line 1905 "Gmsh.y"
{
try {
GMSH_PluginManager::instance()->action((yyvsp[-4].c), (yyvsp[-1].c), 0);
@@ -5409,7 +5410,7 @@ yyreduce:
break;
case 139:
-#line 1914 "Gmsh.y"
+#line 1915 "Gmsh.y"
{
if(!strcmp((yyvsp[-1].c), "ElementsFromAllViews"))
CombineViews(0, 1, CTX.post.combine_remove_orig);
@@ -5434,14 +5435,14 @@ yyreduce:
break;
case 140:
-#line 1936 "Gmsh.y"
+#line 1937 "Gmsh.y"
{
exit(0);
;}
break;
case 141:
-#line 1940 "Gmsh.y"
+#line 1941 "Gmsh.y"
{
CTX.forced_bbox = 0;
SetBoundingBox();
@@ -5449,7 +5450,7 @@ yyreduce:
break;
case 142:
-#line 1945 "Gmsh.y"
+#line 1946 "Gmsh.y"
{
CTX.forced_bbox = 1;
SetBoundingBox((yyvsp[-12].d), (yyvsp[-10].d), (yyvsp[-8].d), (yyvsp[-6].d), (yyvsp[-4].d), (yyvsp[-2].d));
@@ -5457,7 +5458,7 @@ yyreduce:
break;
case 143:
-#line 1950 "Gmsh.y"
+#line 1951 "Gmsh.y"
{
#if defined(HAVE_FLTK)
Draw();
@@ -5466,7 +5467,7 @@ yyreduce:
break;
case 144:
-#line 1962 "Gmsh.y"
+#line 1963 "Gmsh.y"
{
LoopControlVariablesTab[ImbricatedLoop][0] = (yyvsp[-3].d);
LoopControlVariablesTab[ImbricatedLoop][1] = (yyvsp[-1].d);
@@ -5484,7 +5485,7 @@ yyreduce:
break;
case 145:
-#line 1977 "Gmsh.y"
+#line 1978 "Gmsh.y"
{
LoopControlVariablesTab[ImbricatedLoop][0] = (yyvsp[-5].d);
LoopControlVariablesTab[ImbricatedLoop][1] = (yyvsp[-3].d);
@@ -5503,7 +5504,7 @@ yyreduce:
break;
case 146:
-#line 1993 "Gmsh.y"
+#line 1994 "Gmsh.y"
{
LoopControlVariablesTab[ImbricatedLoop][0] = (yyvsp[-3].d);
LoopControlVariablesTab[ImbricatedLoop][1] = (yyvsp[-1].d);
@@ -5532,7 +5533,7 @@ yyreduce:
break;
case 147:
-#line 2019 "Gmsh.y"
+#line 2020 "Gmsh.y"
{
LoopControlVariablesTab[ImbricatedLoop][0] = (yyvsp[-5].d);
LoopControlVariablesTab[ImbricatedLoop][1] = (yyvsp[-3].d);
@@ -5562,7 +5563,7 @@ yyreduce:
break;
case 148:
-#line 2046 "Gmsh.y"
+#line 2047 "Gmsh.y"
{
if(ImbricatedLoop <= 0){
yymsg(GERROR, "Invalid For/EndFor loop");
@@ -5597,7 +5598,7 @@ yyreduce:
break;
case 149:
-#line 2078 "Gmsh.y"
+#line 2079 "Gmsh.y"
{
if(!FunctionManager::Instance()->createFunction((yyvsp[0].c), yyin, yyname, yylineno))
yymsg(GERROR, "Redefinition of function %s", (yyvsp[0].c));
@@ -5607,7 +5608,7 @@ yyreduce:
break;
case 150:
-#line 2085 "Gmsh.y"
+#line 2086 "Gmsh.y"
{
if(!FunctionManager::Instance()->leaveFunction(&yyin, yyname, yylineno))
yymsg(GERROR, "Error while exiting function");
@@ -5615,7 +5616,7 @@ yyreduce:
break;
case 151:
-#line 2090 "Gmsh.y"
+#line 2091 "Gmsh.y"
{
if(!FunctionManager::Instance()->enterFunction((yyvsp[-1].c), &yyin, yyname, yylineno))
yymsg(GERROR, "Unknown function %s", (yyvsp[-1].c));
@@ -5624,20 +5625,20 @@ yyreduce:
break;
case 152:
-#line 2096 "Gmsh.y"
+#line 2097 "Gmsh.y"
{
if(!(yyvsp[-1].d)) skip_until("If", "EndIf");
;}
break;
case 153:
-#line 2100 "Gmsh.y"
+#line 2101 "Gmsh.y"
{
;}
break;
case 154:
-#line 2109 "Gmsh.y"
+#line 2110 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShapes(TRANSLATE, (yyvsp[-1].l),
@@ -5648,7 +5649,7 @@ yyreduce:
break;
case 155:
-#line 2117 "Gmsh.y"
+#line 2118 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShapes(ROTATE, (yyvsp[-1].l),
@@ -5659,7 +5660,7 @@ yyreduce:
break;
case 156:
-#line 2125 "Gmsh.y"
+#line 2126 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShapes(TRANSLATE_ROTATE, (yyvsp[-1].l),
@@ -5670,7 +5671,7 @@ yyreduce:
break;
case 157:
-#line 2133 "Gmsh.y"
+#line 2134 "Gmsh.y"
{
extr.mesh.ExtrudeMesh = false;
extr.mesh.Recombine = false;
@@ -5678,7 +5679,7 @@ yyreduce:
break;
case 158:
-#line 2138 "Gmsh.y"
+#line 2139 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShapes(TRANSLATE, (yyvsp[-3].l),
@@ -5689,7 +5690,7 @@ yyreduce:
break;
case 159:
-#line 2146 "Gmsh.y"
+#line 2147 "Gmsh.y"
{
extr.mesh.ExtrudeMesh = false;
extr.mesh.Recombine = false;
@@ -5697,7 +5698,7 @@ yyreduce:
break;
case 160:
-#line 2151 "Gmsh.y"
+#line 2152 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShapes(ROTATE, (yyvsp[-3].l),
@@ -5708,7 +5709,7 @@ yyreduce:
break;
case 161:
-#line 2159 "Gmsh.y"
+#line 2160 "Gmsh.y"
{
extr.mesh.ExtrudeMesh = false;
extr.mesh.Recombine = false;
@@ -5716,7 +5717,7 @@ yyreduce:
break;
case 162:
-#line 2164 "Gmsh.y"
+#line 2165 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShapes(TRANSLATE_ROTATE, (yyvsp[-3].l),
@@ -5727,7 +5728,7 @@ yyreduce:
break;
case 163:
-#line 2174 "Gmsh.y"
+#line 2175 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(TRANSLATE, MSH_POINT, (int)(yyvsp[-4].d),
@@ -5737,7 +5738,7 @@ yyreduce:
break;
case 164:
-#line 2181 "Gmsh.y"
+#line 2182 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(TRANSLATE, MSH_SEGM_LINE, (int)(yyvsp[-4].d),
@@ -5747,7 +5748,7 @@ yyreduce:
break;
case 165:
-#line 2188 "Gmsh.y"
+#line 2189 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(TRANSLATE, MSH_SURF_PLAN, (int)(yyvsp[-4].d),
@@ -5757,7 +5758,7 @@ yyreduce:
break;
case 166:
-#line 2195 "Gmsh.y"
+#line 2196 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(ROTATE, MSH_POINT, (int)(yyvsp[-8].d),
@@ -5767,7 +5768,7 @@ yyreduce:
break;
case 167:
-#line 2202 "Gmsh.y"
+#line 2203 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(ROTATE, MSH_SEGM_LINE, (int)(yyvsp[-8].d),
@@ -5777,7 +5778,7 @@ yyreduce:
break;
case 168:
-#line 2209 "Gmsh.y"
+#line 2210 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(ROTATE, MSH_SURF_PLAN, (int)(yyvsp[-8].d),
@@ -5787,7 +5788,7 @@ yyreduce:
break;
case 169:
-#line 2216 "Gmsh.y"
+#line 2217 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(TRANSLATE_ROTATE, MSH_POINT, (int)(yyvsp[-10].d),
@@ -5797,7 +5798,7 @@ yyreduce:
break;
case 170:
-#line 2223 "Gmsh.y"
+#line 2224 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(TRANSLATE_ROTATE, MSH_SEGM_LINE, (int)(yyvsp[-10].d),
@@ -5807,7 +5808,7 @@ yyreduce:
break;
case 171:
-#line 2230 "Gmsh.y"
+#line 2231 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(TRANSLATE_ROTATE, MSH_SURF_PLAN, (int)(yyvsp[-10].d),
@@ -5817,7 +5818,7 @@ yyreduce:
break;
case 172:
-#line 2237 "Gmsh.y"
+#line 2238 "Gmsh.y"
{
extr.mesh.ExtrudeMesh = false;
extr.mesh.Recombine = false;
@@ -5825,7 +5826,7 @@ yyreduce:
break;
case 173:
-#line 2242 "Gmsh.y"
+#line 2243 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(TRANSLATE, MSH_POINT, (int)(yyvsp[-8].d),
@@ -5835,7 +5836,7 @@ yyreduce:
break;
case 174:
-#line 2249 "Gmsh.y"
+#line 2250 "Gmsh.y"
{
extr.mesh.ExtrudeMesh = false;
extr.mesh.Recombine = false;
@@ -5843,7 +5844,7 @@ yyreduce:
break;
case 175:
-#line 2254 "Gmsh.y"
+#line 2255 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(TRANSLATE, MSH_SEGM_LINE, (int)(yyvsp[-8].d),
@@ -5853,7 +5854,7 @@ yyreduce:
break;
case 176:
-#line 2261 "Gmsh.y"
+#line 2262 "Gmsh.y"
{
extr.mesh.ExtrudeMesh = false;
extr.mesh.Recombine = false;
@@ -5861,7 +5862,7 @@ yyreduce:
break;
case 177:
-#line 2266 "Gmsh.y"
+#line 2267 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(TRANSLATE, MSH_SURF_PLAN, (int)(yyvsp[-8].d),
@@ -5871,7 +5872,7 @@ yyreduce:
break;
case 178:
-#line 2273 "Gmsh.y"
+#line 2274 "Gmsh.y"
{
extr.mesh.ExtrudeMesh = false;
extr.mesh.Recombine = false;
@@ -5879,7 +5880,7 @@ yyreduce:
break;
case 179:
-#line 2278 "Gmsh.y"
+#line 2279 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(ROTATE, MSH_POINT, (int)(yyvsp[-12].d),
@@ -5889,7 +5890,7 @@ yyreduce:
break;
case 180:
-#line 2285 "Gmsh.y"
+#line 2286 "Gmsh.y"
{
extr.mesh.ExtrudeMesh = false;
extr.mesh.Recombine = false;
@@ -5897,7 +5898,7 @@ yyreduce:
break;
case 181:
-#line 2290 "Gmsh.y"
+#line 2291 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(ROTATE, MSH_SEGM_LINE, (int)(yyvsp[-12].d),
@@ -5907,7 +5908,7 @@ yyreduce:
break;
case 182:
-#line 2297 "Gmsh.y"
+#line 2298 "Gmsh.y"
{
extr.mesh.ExtrudeMesh = false;
extr.mesh.Recombine = false;
@@ -5915,7 +5916,7 @@ yyreduce:
break;
case 183:
-#line 2302 "Gmsh.y"
+#line 2303 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(ROTATE, MSH_SURF_PLAN, (int)(yyvsp[-12].d),
@@ -5925,7 +5926,7 @@ yyreduce:
break;
case 184:
-#line 2309 "Gmsh.y"
+#line 2310 "Gmsh.y"
{
extr.mesh.ExtrudeMesh = false;
extr.mesh.Recombine = false;
@@ -5933,7 +5934,7 @@ yyreduce:
break;
case 185:
-#line 2314 "Gmsh.y"
+#line 2315 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(TRANSLATE_ROTATE, MSH_POINT, (int)(yyvsp[-14].d),
@@ -5943,7 +5944,7 @@ yyreduce:
break;
case 186:
-#line 2321 "Gmsh.y"
+#line 2322 "Gmsh.y"
{
extr.mesh.ExtrudeMesh = false;
extr.mesh.Recombine = false;
@@ -5951,7 +5952,7 @@ yyreduce:
break;
case 187:
-#line 2326 "Gmsh.y"
+#line 2327 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(TRANSLATE_ROTATE, MSH_SEGM_LINE, (int)(yyvsp[-14].d),
@@ -5961,7 +5962,7 @@ yyreduce:
break;
case 188:
-#line 2333 "Gmsh.y"
+#line 2334 "Gmsh.y"
{
extr.mesh.ExtrudeMesh = false;
extr.mesh.Recombine = false;
@@ -5969,7 +5970,7 @@ yyreduce:
break;
case 189:
-#line 2338 "Gmsh.y"
+#line 2339 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(Shape));
ExtrudeShape(TRANSLATE_ROTATE, MSH_SURF_PLAN, (int)(yyvsp[-14].d),
@@ -5979,19 +5980,19 @@ yyreduce:
break;
case 190:
-#line 2349 "Gmsh.y"
+#line 2350 "Gmsh.y"
{
;}
break;
case 191:
-#line 2352 "Gmsh.y"
+#line 2353 "Gmsh.y"
{
;}
break;
case 192:
-#line 2358 "Gmsh.y"
+#line 2359 "Gmsh.y"
{
double d;
extr.mesh.ExtrudeMesh = true;
@@ -6020,7 +6021,7 @@ yyreduce:
break;
case 193:
-#line 2384 "Gmsh.y"
+#line 2385 "Gmsh.y"
{
double d;
extr.mesh.ExtrudeMesh = true;
@@ -6047,14 +6048,14 @@ yyreduce:
break;
case 194:
-#line 2408 "Gmsh.y"
+#line 2409 "Gmsh.y"
{
extr.mesh.Recombine = true;
;}
break;
case 195:
-#line 2417 "Gmsh.y"
+#line 2418 "Gmsh.y"
{
for(int i = 0; i < List_Nbr((yyvsp[-3].l)); i++){
double d;
@@ -6075,7 +6076,7 @@ yyreduce:
break;
case 196:
-#line 2435 "Gmsh.y"
+#line 2436 "Gmsh.y"
{
for(int i = 0; i < List_Nbr((yyvsp[-6].l)); i++){
double d;
@@ -6096,7 +6097,7 @@ yyreduce:
break;
case 197:
-#line 2453 "Gmsh.y"
+#line 2454 "Gmsh.y"
{
for(int i = 0; i < List_Nbr((yyvsp[-6].l)); i++){
double d;
@@ -6117,7 +6118,7 @@ yyreduce:
break;
case 198:
-#line 2471 "Gmsh.y"
+#line 2472 "Gmsh.y"
{
Surface *s = FindSurface((int)(yyvsp[-4].d), THEM);
if(!s)
@@ -6149,7 +6150,7 @@ yyreduce:
break;
case 199:
-#line 2500 "Gmsh.y"
+#line 2501 "Gmsh.y"
{
Surface *s = FindSurface((int)(yyvsp[-5].d), THEM);
if(!s)
@@ -6187,7 +6188,7 @@ yyreduce:
break;
case 200:
-#line 2535 "Gmsh.y"
+#line 2536 "Gmsh.y"
{
Surface *s = FindSurface((int)(yyvsp[-4].d), THEM);
if(!s)
@@ -6217,7 +6218,7 @@ yyreduce:
break;
case 201:
-#line 2562 "Gmsh.y"
+#line 2563 "Gmsh.y"
{
Volume *v = FindVolume((int)(yyvsp[-4].d), THEM);
if(!v)
@@ -6247,7 +6248,7 @@ yyreduce:
break;
case 202:
-#line 2589 "Gmsh.y"
+#line 2590 "Gmsh.y"
{
for(int i = 0; i < List_Nbr((yyvsp[-3].l)); i++){
double d;
@@ -6264,7 +6265,7 @@ yyreduce:
break;
case 203:
-#line 2603 "Gmsh.y"
+#line 2604 "Gmsh.y"
{
for(int i = 0; i < List_Nbr((yyvsp[-1].l)); i++){
double d;
@@ -6280,7 +6281,7 @@ yyreduce:
break;
case 204:
-#line 2622 "Gmsh.y"
+#line 2623 "Gmsh.y"
{
Surface *s = FindSurface((int)(yyvsp[-1].d), THEM);
if(s)
@@ -6289,7 +6290,7 @@ yyreduce:
break;
case 205:
-#line 2628 "Gmsh.y"
+#line 2629 "Gmsh.y"
{
Surface *s = FindSurface((int)(yyvsp[-1].d), THEM);
if(s)
@@ -6298,73 +6299,73 @@ yyreduce:
break;
case 206:
-#line 2634 "Gmsh.y"
+#line 2635 "Gmsh.y"
{
;}
break;
case 207:
-#line 2637 "Gmsh.y"
+#line 2638 "Gmsh.y"
{
;}
break;
case 208:
-#line 2644 "Gmsh.y"
+#line 2645 "Gmsh.y"
{
ReplaceAllDuplicates(THEM);
;}
break;
case 209:
-#line 2648 "Gmsh.y"
+#line 2649 "Gmsh.y"
{
IntersectAllSegmentsTogether();
;}
break;
case 210:
-#line 2657 "Gmsh.y"
+#line 2658 "Gmsh.y"
{ (yyval.d) = (yyvsp[0].d); ;}
break;
case 211:
-#line 2658 "Gmsh.y"
+#line 2659 "Gmsh.y"
{ (yyval.d) = (yyvsp[-1].d); ;}
break;
case 212:
-#line 2659 "Gmsh.y"
+#line 2660 "Gmsh.y"
{ (yyval.d) = -(yyvsp[0].d); ;}
break;
case 213:
-#line 2660 "Gmsh.y"
+#line 2661 "Gmsh.y"
{ (yyval.d) = (yyvsp[0].d); ;}
break;
case 214:
-#line 2661 "Gmsh.y"
+#line 2662 "Gmsh.y"
{ (yyval.d) = !(yyvsp[0].d); ;}
break;
case 215:
-#line 2662 "Gmsh.y"
+#line 2663 "Gmsh.y"
{ (yyval.d) = (yyvsp[-2].d) - (yyvsp[0].d); ;}
break;
case 216:
-#line 2663 "Gmsh.y"
+#line 2664 "Gmsh.y"
{ (yyval.d) = (yyvsp[-2].d) + (yyvsp[0].d); ;}
break;
case 217:
-#line 2664 "Gmsh.y"
+#line 2665 "Gmsh.y"
{ (yyval.d) = (yyvsp[-2].d) * (yyvsp[0].d); ;}
break;
case 218:
-#line 2666 "Gmsh.y"
+#line 2667 "Gmsh.y"
{
if(!(yyvsp[0].d))
yymsg(GERROR, "Division by zero in '%g / %g'", (yyvsp[-2].d), (yyvsp[0].d));
@@ -6374,307 +6375,307 @@ yyreduce:
break;
case 219:
-#line 2672 "Gmsh.y"
+#line 2673 "Gmsh.y"
{ (yyval.d) = (int)(yyvsp[-2].d) % (int)(yyvsp[0].d); ;}
break;
case 220:
-#line 2673 "Gmsh.y"
+#line 2674 "Gmsh.y"
{ (yyval.d) = pow((yyvsp[-2].d), (yyvsp[0].d)); ;}
break;
case 221:
-#line 2674 "Gmsh.y"
+#line 2675 "Gmsh.y"
{ (yyval.d) = (yyvsp[-2].d) < (yyvsp[0].d); ;}
break;
case 222:
-#line 2675 "Gmsh.y"
+#line 2676 "Gmsh.y"
{ (yyval.d) = (yyvsp[-2].d) > (yyvsp[0].d); ;}
break;
case 223:
-#line 2676 "Gmsh.y"
+#line 2677 "Gmsh.y"
{ (yyval.d) = (yyvsp[-2].d) <= (yyvsp[0].d); ;}
break;
case 224:
-#line 2677 "Gmsh.y"
+#line 2678 "Gmsh.y"
{ (yyval.d) = (yyvsp[-2].d) >= (yyvsp[0].d); ;}
break;
case 225:
-#line 2678 "Gmsh.y"
+#line 2679 "Gmsh.y"
{ (yyval.d) = (yyvsp[-2].d) == (yyvsp[0].d); ;}
break;
case 226:
-#line 2679 "Gmsh.y"
+#line 2680 "Gmsh.y"
{ (yyval.d) = (yyvsp[-2].d) != (yyvsp[0].d); ;}
break;
case 227:
-#line 2680 "Gmsh.y"
+#line 2681 "Gmsh.y"
{ (yyval.d) = (yyvsp[-2].d) && (yyvsp[0].d); ;}
break;
case 228:
-#line 2681 "Gmsh.y"
+#line 2682 "Gmsh.y"
{ (yyval.d) = (yyvsp[-2].d) || (yyvsp[0].d); ;}
break;
case 229:
-#line 2682 "Gmsh.y"
+#line 2683 "Gmsh.y"
{ (yyval.d) = (yyvsp[-4].d)? (yyvsp[-2].d) : (yyvsp[0].d); ;}
break;
case 230:
-#line 2683 "Gmsh.y"
+#line 2684 "Gmsh.y"
{ (yyval.d) = exp((yyvsp[-1].d)); ;}
break;
case 231:
-#line 2684 "Gmsh.y"
+#line 2685 "Gmsh.y"
{ (yyval.d) = log((yyvsp[-1].d)); ;}
break;
case 232:
-#line 2685 "Gmsh.y"
+#line 2686 "Gmsh.y"
{ (yyval.d) = log10((yyvsp[-1].d)); ;}
break;
case 233:
-#line 2686 "Gmsh.y"
+#line 2687 "Gmsh.y"
{ (yyval.d) = sqrt((yyvsp[-1].d)); ;}
break;
case 234:
-#line 2687 "Gmsh.y"
+#line 2688 "Gmsh.y"
{ (yyval.d) = sin((yyvsp[-1].d)); ;}
break;
case 235:
-#line 2688 "Gmsh.y"
+#line 2689 "Gmsh.y"
{ (yyval.d) = asin((yyvsp[-1].d)); ;}
break;
case 236:
-#line 2689 "Gmsh.y"
+#line 2690 "Gmsh.y"
{ (yyval.d) = cos((yyvsp[-1].d)); ;}
break;
case 237:
-#line 2690 "Gmsh.y"
+#line 2691 "Gmsh.y"
{ (yyval.d) = acos((yyvsp[-1].d)); ;}
break;
case 238:
-#line 2691 "Gmsh.y"
+#line 2692 "Gmsh.y"
{ (yyval.d) = tan((yyvsp[-1].d)); ;}
break;
case 239:
-#line 2692 "Gmsh.y"
+#line 2693 "Gmsh.y"
{ (yyval.d) = atan((yyvsp[-1].d)); ;}
break;
case 240:
-#line 2693 "Gmsh.y"
+#line 2694 "Gmsh.y"
{ (yyval.d) = atan2((yyvsp[-3].d), (yyvsp[-1].d));;}
break;
case 241:
-#line 2694 "Gmsh.y"
+#line 2695 "Gmsh.y"
{ (yyval.d) = sinh((yyvsp[-1].d)); ;}
break;
case 242:
-#line 2695 "Gmsh.y"
+#line 2696 "Gmsh.y"
{ (yyval.d) = cosh((yyvsp[-1].d)); ;}
break;
case 243:
-#line 2696 "Gmsh.y"
+#line 2697 "Gmsh.y"
{ (yyval.d) = tanh((yyvsp[-1].d)); ;}
break;
case 244:
-#line 2697 "Gmsh.y"
+#line 2698 "Gmsh.y"
{ (yyval.d) = fabs((yyvsp[-1].d)); ;}
break;
case 245:
-#line 2698 "Gmsh.y"
+#line 2699 "Gmsh.y"
{ (yyval.d) = floor((yyvsp[-1].d)); ;}
break;
case 246:
-#line 2699 "Gmsh.y"
+#line 2700 "Gmsh.y"
{ (yyval.d) = ceil((yyvsp[-1].d)); ;}
break;
case 247:
-#line 2700 "Gmsh.y"
+#line 2701 "Gmsh.y"
{ (yyval.d) = fmod((yyvsp[-3].d), (yyvsp[-1].d)); ;}
break;
case 248:
-#line 2701 "Gmsh.y"
+#line 2702 "Gmsh.y"
{ (yyval.d) = fmod((yyvsp[-3].d), (yyvsp[-1].d)); ;}
break;
case 249:
-#line 2702 "Gmsh.y"
+#line 2703 "Gmsh.y"
{ (yyval.d) = sqrt((yyvsp[-3].d)*(yyvsp[-3].d)+(yyvsp[-1].d)*(yyvsp[-1].d)); ;}
break;
case 250:
-#line 2703 "Gmsh.y"
+#line 2704 "Gmsh.y"
{ (yyval.d) = (yyvsp[-1].d)*(double)rand()/(double)RAND_MAX; ;}
break;
case 251:
-#line 2705 "Gmsh.y"
+#line 2706 "Gmsh.y"
{ (yyval.d) = exp((yyvsp[-1].d)); ;}
break;
case 252:
-#line 2706 "Gmsh.y"
+#line 2707 "Gmsh.y"
{ (yyval.d) = log((yyvsp[-1].d)); ;}
break;
case 253:
-#line 2707 "Gmsh.y"
+#line 2708 "Gmsh.y"
{ (yyval.d) = log10((yyvsp[-1].d)); ;}
break;
case 254:
-#line 2708 "Gmsh.y"
+#line 2709 "Gmsh.y"
{ (yyval.d) = sqrt((yyvsp[-1].d)); ;}
break;
case 255:
-#line 2709 "Gmsh.y"
+#line 2710 "Gmsh.y"
{ (yyval.d) = sin((yyvsp[-1].d)); ;}
break;
case 256:
-#line 2710 "Gmsh.y"
+#line 2711 "Gmsh.y"
{ (yyval.d) = asin((yyvsp[-1].d)); ;}
break;
case 257:
-#line 2711 "Gmsh.y"
+#line 2712 "Gmsh.y"
{ (yyval.d) = cos((yyvsp[-1].d)); ;}
break;
case 258:
-#line 2712 "Gmsh.y"
+#line 2713 "Gmsh.y"
{ (yyval.d) = acos((yyvsp[-1].d)); ;}
break;
case 259:
-#line 2713 "Gmsh.y"
+#line 2714 "Gmsh.y"
{ (yyval.d) = tan((yyvsp[-1].d)); ;}
break;
case 260:
-#line 2714 "Gmsh.y"
+#line 2715 "Gmsh.y"
{ (yyval.d) = atan((yyvsp[-1].d)); ;}
break;
case 261:
-#line 2715 "Gmsh.y"
+#line 2716 "Gmsh.y"
{ (yyval.d) = atan2((yyvsp[-3].d), (yyvsp[-1].d));;}
break;
case 262:
-#line 2716 "Gmsh.y"
+#line 2717 "Gmsh.y"
{ (yyval.d) = sinh((yyvsp[-1].d)); ;}
break;
case 263:
-#line 2717 "Gmsh.y"
+#line 2718 "Gmsh.y"
{ (yyval.d) = cosh((yyvsp[-1].d)); ;}
break;
case 264:
-#line 2718 "Gmsh.y"
+#line 2719 "Gmsh.y"
{ (yyval.d) = tanh((yyvsp[-1].d)); ;}
break;
case 265:
-#line 2719 "Gmsh.y"
+#line 2720 "Gmsh.y"
{ (yyval.d) = fabs((yyvsp[-1].d)); ;}
break;
case 266:
-#line 2720 "Gmsh.y"
+#line 2721 "Gmsh.y"
{ (yyval.d) = floor((yyvsp[-1].d)); ;}
break;
case 267:
-#line 2721 "Gmsh.y"
+#line 2722 "Gmsh.y"
{ (yyval.d) = ceil((yyvsp[-1].d)); ;}
break;
case 268:
-#line 2722 "Gmsh.y"
+#line 2723 "Gmsh.y"
{ (yyval.d) = fmod((yyvsp[-3].d), (yyvsp[-1].d)); ;}
break;
case 269:
-#line 2723 "Gmsh.y"
+#line 2724 "Gmsh.y"
{ (yyval.d) = fmod((yyvsp[-3].d), (yyvsp[-1].d)); ;}
break;
case 270:
-#line 2724 "Gmsh.y"
+#line 2725 "Gmsh.y"
{ (yyval.d) = sqrt((yyvsp[-3].d)*(yyvsp[-3].d)+(yyvsp[-1].d)*(yyvsp[-1].d)); ;}
break;
case 271:
-#line 2725 "Gmsh.y"
+#line 2726 "Gmsh.y"
{ (yyval.d) = (yyvsp[-1].d)*(double)rand()/(double)RAND_MAX; ;}
break;
case 272:
-#line 2734 "Gmsh.y"
+#line 2735 "Gmsh.y"
{ (yyval.d) = (yyvsp[0].d); ;}
break;
case 273:
-#line 2735 "Gmsh.y"
+#line 2736 "Gmsh.y"
{ (yyval.d) = 3.141592653589793; ;}
break;
case 274:
-#line 2736 "Gmsh.y"
+#line 2737 "Gmsh.y"
{ (yyval.d) = ParUtil::Instance()->rank(); ;}
break;
case 275:
-#line 2737 "Gmsh.y"
+#line 2738 "Gmsh.y"
{ (yyval.d) = ParUtil::Instance()->size(); ;}
break;
case 276:
-#line 2738 "Gmsh.y"
+#line 2739 "Gmsh.y"
{ (yyval.d) = Get_GmshMajorVersion(); ;}
break;
case 277:
-#line 2739 "Gmsh.y"
+#line 2740 "Gmsh.y"
{ (yyval.d) = Get_GmshMinorVersion(); ;}
break;
case 278:
-#line 2740 "Gmsh.y"
+#line 2741 "Gmsh.y"
{ (yyval.d) = Get_GmshPatchVersion(); ;}
break;
case 279:
-#line 2745 "Gmsh.y"
+#line 2746 "Gmsh.y"
{
Symbol TheSymbol;
TheSymbol.Name = (yyvsp[0].c);
@@ -6690,7 +6691,7 @@ yyreduce:
break;
case 280:
-#line 2761 "Gmsh.y"
+#line 2762 "Gmsh.y"
{
char tmpstring[1024];
sprintf(tmpstring, "%s_%d", (yyvsp[-4].c), (int)(yyvsp[-1].d)) ;
@@ -6708,7 +6709,7 @@ yyreduce:
break;
case 281:
-#line 2776 "Gmsh.y"
+#line 2777 "Gmsh.y"
{
Symbol TheSymbol;
TheSymbol.Name = (yyvsp[-3].c);
@@ -6731,7 +6732,7 @@ yyreduce:
break;
case 282:
-#line 2796 "Gmsh.y"
+#line 2797 "Gmsh.y"
{
Symbol TheSymbol;
TheSymbol.Name = (yyvsp[-2].c);
@@ -6748,7 +6749,7 @@ yyreduce:
break;
case 283:
-#line 2810 "Gmsh.y"
+#line 2811 "Gmsh.y"
{
Symbol TheSymbol;
TheSymbol.Name = (yyvsp[-1].c);
@@ -6764,7 +6765,7 @@ yyreduce:
break;
case 284:
-#line 2823 "Gmsh.y"
+#line 2824 "Gmsh.y"
{
Symbol TheSymbol;
TheSymbol.Name = (yyvsp[-4].c);
@@ -6787,7 +6788,7 @@ yyreduce:
break;
case 285:
-#line 2846 "Gmsh.y"
+#line 2847 "Gmsh.y"
{
double (*pNumOpt)(int num, int action, double value);
StringXNumber *pNumCat;
@@ -6808,7 +6809,7 @@ yyreduce:
break;
case 286:
-#line 2864 "Gmsh.y"
+#line 2865 "Gmsh.y"
{
double (*pNumOpt)(int num, int action, double value);
StringXNumber *pNumCat;
@@ -6829,7 +6830,7 @@ yyreduce:
break;
case 287:
-#line 2882 "Gmsh.y"
+#line 2883 "Gmsh.y"
{
double (*pNumOpt)(int num, int action, double value);
StringXNumber *pNumCat;
@@ -6850,7 +6851,7 @@ yyreduce:
break;
case 288:
-#line 2900 "Gmsh.y"
+#line 2901 "Gmsh.y"
{
double (*pNumOpt)(int num, int action, double value);
StringXNumber *pNumCat;
@@ -6871,7 +6872,7 @@ yyreduce:
break;
case 289:
-#line 2918 "Gmsh.y"
+#line 2919 "Gmsh.y"
{
(yyval.d) = GetValue((yyvsp[-3].c), (yyvsp[-1].d));
Free((yyvsp[-3].c));
@@ -6879,90 +6880,90 @@ yyreduce:
break;
case 290:
-#line 2926 "Gmsh.y"
+#line 2927 "Gmsh.y"
{
memcpy((yyval.v), (yyvsp[0].v), 5*sizeof(double));
;}
break;
case 291:
-#line 2930 "Gmsh.y"
+#line 2931 "Gmsh.y"
{
for(int i = 0; i < 5; i++) (yyval.v)[i] = -(yyvsp[0].v)[i];
;}
break;
case 292:
-#line 2934 "Gmsh.y"
+#line 2935 "Gmsh.y"
{
for(int i = 0; i < 5; i++) (yyval.v)[i] = (yyvsp[0].v)[i];
;}
break;
case 293:
-#line 2938 "Gmsh.y"
+#line 2939 "Gmsh.y"
{
for(int i = 0; i < 5; i++) (yyval.v)[i] = (yyvsp[-2].v)[i] - (yyvsp[0].v)[i];
;}
break;
case 294:
-#line 2942 "Gmsh.y"
+#line 2943 "Gmsh.y"
{
for(int i = 0; i < 5; i++) (yyval.v)[i] = (yyvsp[-2].v)[i] + (yyvsp[0].v)[i];
;}
break;
case 295:
-#line 2949 "Gmsh.y"
+#line 2950 "Gmsh.y"
{
(yyval.v)[0] = (yyvsp[-9].d); (yyval.v)[1] = (yyvsp[-7].d); (yyval.v)[2] = (yyvsp[-5].d); (yyval.v)[3] = (yyvsp[-3].d); (yyval.v)[4] = (yyvsp[-1].d);
;}
break;
case 296:
-#line 2953 "Gmsh.y"
+#line 2954 "Gmsh.y"
{
(yyval.v)[0] = (yyvsp[-7].d); (yyval.v)[1] = (yyvsp[-5].d); (yyval.v)[2] = (yyvsp[-3].d); (yyval.v)[3] = (yyvsp[-1].d); (yyval.v)[4] = 1.0;
;}
break;
case 297:
-#line 2957 "Gmsh.y"
+#line 2958 "Gmsh.y"
{
(yyval.v)[0] = (yyvsp[-5].d); (yyval.v)[1] = (yyvsp[-3].d); (yyval.v)[2] = (yyvsp[-1].d); (yyval.v)[3] = 0.0; (yyval.v)[4] = 1.0;
;}
break;
case 298:
-#line 2961 "Gmsh.y"
+#line 2962 "Gmsh.y"
{
(yyval.v)[0] = (yyvsp[-5].d); (yyval.v)[1] = (yyvsp[-3].d); (yyval.v)[2] = (yyvsp[-1].d); (yyval.v)[3] = 0.0; (yyval.v)[4] = 1.0;
;}
break;
case 299:
-#line 2968 "Gmsh.y"
+#line 2969 "Gmsh.y"
{
;}
break;
case 300:
-#line 2971 "Gmsh.y"
+#line 2972 "Gmsh.y"
{
(yyval.l) = (yyvsp[-1].l);
;}
break;
case 301:
-#line 2975 "Gmsh.y"
+#line 2976 "Gmsh.y"
{
(yyval.l) = (yyvsp[-1].l);
;}
break;
case 302:
-#line 2982 "Gmsh.y"
+#line 2983 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(List_T*));
List_Add((yyval.l), &((yyvsp[0].l)));
@@ -6970,14 +6971,14 @@ yyreduce:
break;
case 303:
-#line 2987 "Gmsh.y"
+#line 2988 "Gmsh.y"
{
List_Add((yyval.l), &((yyvsp[0].l)));
;}
break;
case 304:
-#line 2995 "Gmsh.y"
+#line 2996 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(double));
List_Add((yyval.l), &((yyvsp[0].d)));
@@ -6985,14 +6986,14 @@ yyreduce:
break;
case 305:
-#line 3000 "Gmsh.y"
+#line 3001 "Gmsh.y"
{
(yyval.l) = (yyvsp[0].l);
;}
break;
case 306:
-#line 3004 "Gmsh.y"
+#line 3005 "Gmsh.y"
{
// creates an empty list
(yyval.l) = List_Create(2, 1, sizeof(double));
@@ -7000,14 +7001,14 @@ yyreduce:
break;
case 307:
-#line 3009 "Gmsh.y"
+#line 3010 "Gmsh.y"
{
(yyval.l) = (yyvsp[-1].l);
;}
break;
case 308:
-#line 3013 "Gmsh.y"
+#line 3014 "Gmsh.y"
{
(yyval.l) = (yyvsp[-1].l);
double *pd;
@@ -7019,7 +7020,7 @@ yyreduce:
break;
case 309:
-#line 3025 "Gmsh.y"
+#line 3026 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(double));
for(double d = (yyvsp[-2].d); ((yyvsp[-2].d) < (yyvsp[0].d)) ? (d <= (yyvsp[0].d)) : (d >= (yyvsp[0].d)); ((yyvsp[-2].d) < (yyvsp[0].d)) ? (d += 1.) : (d -= 1.))
@@ -7028,7 +7029,7 @@ yyreduce:
break;
case 310:
-#line 3031 "Gmsh.y"
+#line 3032 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(double));
if(!(yyvsp[0].d) || ((yyvsp[-4].d) < (yyvsp[-2].d) && (yyvsp[0].d) < 0) || ((yyvsp[-4].d) > (yyvsp[-2].d) && (yyvsp[0].d) > 0)){
@@ -7042,7 +7043,7 @@ yyreduce:
break;
case 311:
-#line 3042 "Gmsh.y"
+#line 3043 "Gmsh.y"
{
// Returns the coordinates of a point and fills a list with it.
// This allows to ensure e.g. that relative point positions are
@@ -7065,7 +7066,7 @@ yyreduce:
break;
case 312:
-#line 3062 "Gmsh.y"
+#line 3063 "Gmsh.y"
{
(yyval.l) = List_Create(List_Nbr((yyvsp[0].l)), 1, sizeof(double));
for(int i = 0; i < List_Nbr((yyvsp[0].l)); i++){
@@ -7078,7 +7079,7 @@ yyreduce:
break;
case 313:
-#line 3072 "Gmsh.y"
+#line 3073 "Gmsh.y"
{
(yyval.l) = List_Create(List_Nbr((yyvsp[0].l)), 1, sizeof(double));
for(int i = 0; i < List_Nbr((yyvsp[0].l)); i++){
@@ -7091,7 +7092,7 @@ yyreduce:
break;
case 314:
-#line 3082 "Gmsh.y"
+#line 3083 "Gmsh.y"
{
(yyval.l) = List_Create(List_Nbr((yyvsp[0].l)), 1, sizeof(double));
for(int i = 0; i < List_Nbr((yyvsp[0].l)); i++){
@@ -7104,7 +7105,7 @@ yyreduce:
break;
case 315:
-#line 3092 "Gmsh.y"
+#line 3093 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(double));
Symbol TheSymbol;
@@ -7124,7 +7125,7 @@ yyreduce:
break;
case 316:
-#line 3109 "Gmsh.y"
+#line 3110 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(double));
Symbol TheSymbol;
@@ -7146,7 +7147,7 @@ yyreduce:
break;
case 317:
-#line 3128 "Gmsh.y"
+#line 3129 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(double));
Symbol TheSymbol;
@@ -7173,7 +7174,7 @@ yyreduce:
break;
case 318:
-#line 3152 "Gmsh.y"
+#line 3153 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(double));
Symbol TheSymbol;
@@ -7202,7 +7203,7 @@ yyreduce:
break;
case 319:
-#line 3181 "Gmsh.y"
+#line 3182 "Gmsh.y"
{
(yyval.l) = List_Create(2, 1, sizeof(double));
List_Add((yyval.l), &((yyvsp[0].d)));
@@ -7210,21 +7211,21 @@ yyreduce:
break;
case 320:
-#line 3186 "Gmsh.y"
+#line 3187 "Gmsh.y"
{
(yyval.l) = (yyvsp[0].l);
;}
break;
case 321:
-#line 3190 "Gmsh.y"
+#line 3191 "Gmsh.y"
{
List_Add((yyval.l), &((yyvsp[0].d)));
;}
break;
case 322:
-#line 3194 "Gmsh.y"
+#line 3195 "Gmsh.y"
{
for(int i = 0; i < List_Nbr((yyvsp[0].l)); i++){
double d;
@@ -7236,21 +7237,21 @@ yyreduce:
break;
case 323:
-#line 3207 "Gmsh.y"
+#line 3208 "Gmsh.y"
{
(yyval.u) = CTX.PACK_COLOR((int)(yyvsp[-7].d), (int)(yyvsp[-5].d), (int)(yyvsp[-3].d), (int)(yyvsp[-1].d));
;}
break;
case 324:
-#line 3211 "Gmsh.y"
+#line 3212 "Gmsh.y"
{
(yyval.u) = CTX.PACK_COLOR((int)(yyvsp[-5].d), (int)(yyvsp[-3].d), (int)(yyvsp[-1].d), 255);
;}
break;
case 325:
-#line 3223 "Gmsh.y"
+#line 3224 "Gmsh.y"
{
int flag;
(yyval.u) = Get_ColorForString(ColorString, -1, (yyvsp[0].c), &flag);
@@ -7260,7 +7261,7 @@ yyreduce:
break;
case 326:
-#line 3230 "Gmsh.y"
+#line 3231 "Gmsh.y"
{
unsigned int (*pColOpt)(int num, int action, unsigned int value);
StringXColor *pColCat;
@@ -7282,14 +7283,14 @@ yyreduce:
break;
case 327:
-#line 3252 "Gmsh.y"
+#line 3253 "Gmsh.y"
{
(yyval.l) = (yyvsp[-1].l);
;}
break;
case 328:
-#line 3256 "Gmsh.y"
+#line 3257 "Gmsh.y"
{
(yyval.l) = List_Create(256, 10, sizeof(unsigned int));
GmshColorTable *ct = Get_ColorTable((int)(yyvsp[-3].d));
@@ -7304,7 +7305,7 @@ yyreduce:
break;
case 329:
-#line 3271 "Gmsh.y"
+#line 3272 "Gmsh.y"
{
(yyval.l) = List_Create(256, 10, sizeof(unsigned int));
List_Add((yyval.l), &((yyvsp[0].u)));
@@ -7312,21 +7313,21 @@ yyreduce:
break;
case 330:
-#line 3276 "Gmsh.y"
+#line 3277 "Gmsh.y"
{
List_Add((yyval.l), &((yyvsp[0].u)));
;}
break;
case 331:
-#line 3283 "Gmsh.y"
+#line 3284 "Gmsh.y"
{
(yyval.c) = (yyvsp[0].c);
;}
break;
case 332:
-#line 3287 "Gmsh.y"
+#line 3288 "Gmsh.y"
{
(yyval.c) = (char *)Malloc(32*sizeof(char));
time_t now;
@@ -7337,7 +7338,7 @@ yyreduce:
break;
case 333:
-#line 3295 "Gmsh.y"
+#line 3296 "Gmsh.y"
{
(yyval.c) = (char *)Malloc((strlen((yyvsp[-3].c))+strlen((yyvsp[-1].c))+1)*sizeof(char));
strcpy((yyval.c), (yyvsp[-3].c));
@@ -7348,7 +7349,7 @@ yyreduce:
break;
case 334:
-#line 3303 "Gmsh.y"
+#line 3304 "Gmsh.y"
{
(yyval.c) = (char *)Malloc((strlen((yyvsp[-1].c))+1)*sizeof(char));
int i;
@@ -7365,7 +7366,7 @@ yyreduce:
break;
case 335:
-#line 3317 "Gmsh.y"
+#line 3318 "Gmsh.y"
{
(yyval.c) = (char *)Malloc((strlen((yyvsp[-1].c))+1)*sizeof(char));
int i;
@@ -7382,14 +7383,14 @@ yyreduce:
break;
case 336:
-#line 3331 "Gmsh.y"
+#line 3332 "Gmsh.y"
{
(yyval.c) = (yyvsp[-1].c);
;}
break;
case 337:
-#line 3335 "Gmsh.y"
+#line 3336 "Gmsh.y"
{
char tmpstring[1024];
int i = PrintListOfDouble((yyvsp[-3].c), (yyvsp[-1].l), tmpstring);
@@ -7411,7 +7412,7 @@ yyreduce:
break;
case 338:
-#line 3354 "Gmsh.y"
+#line 3355 "Gmsh.y"
{
char* (*pStrOpt)(int num, int action, char *value);
StringXString *pStrCat;
@@ -7436,7 +7437,7 @@ yyreduce:
break;
case 339:
-#line 3376 "Gmsh.y"
+#line 3377 "Gmsh.y"
{
char* (*pStrOpt)(int num, int action, char *value);
StringXString *pStrCat;
@@ -7464,7 +7465,7 @@ yyreduce:
}
/* Line 1037 of yacc.c. */
-#line 7468 "Gmsh.tab.cpp"
+#line 7469 "Gmsh.tab.cpp"
yyvsp -= yylen;
yyssp -= yylen;
@@ -7692,7 +7693,7 @@ yyreturn:
}
-#line 3399 "Gmsh.y"
+#line 3400 "Gmsh.y"
void DeleteSymbol(void *a, void *b){
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index 573f55a54c5a7dac3fed71ca72c99d7cf5f42b45..29a205d9ab15bcf4e0305437491aa22638e08f68 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -1,5 +1,5 @@
%{
-// $Id: Gmsh.y,v 1.234 2006-08-19 08:26:47 remacle Exp $
+// $Id: Gmsh.y,v 1.235 2006-08-29 10:39:54 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -1887,7 +1887,8 @@ Command :
SleepInSeconds($2);
}
else if(!strcmp($1, "Remesh")){
- ReMesh();
+ Msg(GERROR, "Surface ReMeshing must be reinterfaced");
+ // ReMesh();
}
else if(!strcmp($1, "Mesh")){
yymsg(GERROR, "Mesh directives are not (yet) allowed in scripts");
diff --git a/Parser/Gmsh.yy.cpp b/Parser/Gmsh.yy.cpp
index f8f0a5499a4ea5784e8b841794b9c7d27f7451c6..5ef163627c103c468f4ede23edf1576a88da0267 100644
--- a/Parser/Gmsh.yy.cpp
+++ b/Parser/Gmsh.yy.cpp
@@ -2,7 +2,7 @@
/* A lexical scanner generated by flex */
/* Scanner skeleton version:
- * $Header: /cvsroot/gmsh/Parser/Gmsh.yy.cpp,v 1.272 2006-08-19 08:26:47 remacle Exp $
+ * $Header: /cvsroot/gmsh/Parser/Gmsh.yy.cpp,v 1.273 2006-08-29 10:39:54 remacle Exp $
*/
#define FLEX_SCANNER
@@ -727,7 +727,7 @@ char *yytext;
#line 1 "Gmsh.l"
#define INITIAL 0
#line 2 "Gmsh.l"
-// $Id: Gmsh.yy.cpp,v 1.272 2006-08-19 08:26:47 remacle Exp $
+// $Id: Gmsh.yy.cpp,v 1.273 2006-08-29 10:39:54 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
diff --git a/Plugin/ExtractEdges.cpp b/Plugin/ExtractEdges.cpp
index 522b343ce23d53235b6ab84efe344cd9f3ad1702..a64f745ccbd9adb1b0d170cadbee1e0f73815b3b 100644
--- a/Plugin/ExtractEdges.cpp
+++ b/Plugin/ExtractEdges.cpp
@@ -1,4 +1,4 @@
-// $Id: ExtractEdges.cpp,v 1.2 2006-03-15 19:00:21 geuzaine Exp $
+// $Id: ExtractEdges.cpp,v 1.3 2006-08-29 10:39:54 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -99,7 +99,10 @@ Post_View *GMSH_ExtractEdgesPlugin::execute(Post_View * v)
BDS_Mesh bds;
bds.import_view(v1, CTX.lc * 1.e-12);
- bds.classify(angle * M_PI / 180.);
+ // bds.classify(angle * M_PI / 180.);
+
+ Msg(GERROR, "BDS->classify(angle, edge_prolongation) must be reinterfaced");
+
std::list<BDS_Edge*>::iterator it = bds.edges.begin();
std::list<BDS_Edge*>::iterator ite = bds.edges.end();
while (it != ite){
diff --git a/benchmarks/2d/many_holes.geo b/benchmarks/2d/many_holes.geo
index 2b8c9c9c941105e8d23e87de9e9cd9da6add1773..839d2c13aafaca8a4605c545db0129494f12e86c 100644
--- a/benchmarks/2d/many_holes.geo
+++ b/benchmarks/2d/many_holes.geo
@@ -12,10 +12,10 @@ Line Loop(10) = {1,-3,2};
Function TubeHole
- p1 = newp; Point(p1) = {x,y,0,8};
- p2 = newp; Point(p2) = {x+r,y,0,8};
- p3 = newp; Point(p3) = {x,y+r,0,8};
- p4 = newp; Point(p4) = {x-0.71*r,y-0.71*r,0,8};
+ p1 = newp; Point(p1) = {x,y,0,1};
+ p2 = newp; Point(p2) = {x+r,y,0,1};
+ p3 = newp; Point(p3) = {x,y+r,0,1};
+ p4 = newp; Point(p4) = {x-0.71*r,y-0.71*r,0,1};
c1 = newreg; Circle(c1) = {p2,p1,p3};
c2 = newreg; Circle(c2) = {p3,p1,p4};