From 173921b9f5ae1bfffeaf40906c5375d41d307821 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Sat, 26 Aug 2006 15:13:22 +0000
Subject: [PATCH] *** empty log message ***
---
Geo/GEdge.cpp | 12 +++-
Geo/GEdge.h | 7 ++
Geo/GFace.cpp | 42 ++++++++++-
Geo/GFace.h | 3 +
Geo/SVector3.cpp | 6 +-
Geo/SVector3.h | 2 +-
Geo/gmshEdge.cpp | 6 +-
Geo/gmshModel.cpp | 4 +-
Mesh/BDS.cpp | 6 +-
Mesh/BDS.h | 6 +-
Mesh/Interpolation.cpp | 5 +-
Mesh/Utils.cpp | 19 ++---
Mesh/meshGEdge.cpp | 31 +++++---
Mesh/meshGFace.cpp | 138 +++++++++++++++++++++++++++---------
Parser/OpenFile.cpp | 6 +-
benchmarks/2d/Square-01.geo | 2 +
16 files changed, 227 insertions(+), 68 deletions(-)
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 6f4af768da..de954a6eb5 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: GEdge.cpp,v 1.14 2006-08-26 13:34:46 geuzaine Exp $
+// $Id: GEdge.cpp,v 1.15 2006-08-26 15:13:22 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -108,3 +108,13 @@ std::string GEdge::getAdditionalInfoString()
sprintf(tmp, "{%d,%d}", v0->tag(), v1->tag());
return std::string(tmp);
}
+
+
+/// use central differences
+SVector3 GEdge::secondDer(double par) const
+{
+ const double eps = 1.e-3;
+ SVector3 x1 = firstDer(par-eps);
+ SVector3 x2 = firstDer(par+eps);
+ return 500*(x2-x1);
+}
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 47850c6255..3b242c5098 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -54,6 +54,9 @@ class GEdge : public GEntity {
// The bounding box
SBoundingBox3d bounds() const;
+ // Faces that bound this entity or that this entity bounds.
+ virtual std::list<GFace*> faces() const{return l_faces;}
+
// Get the parameter location for a point in space on the edge.
virtual double parFromPoint(const SPoint3 &) const = 0;
@@ -69,6 +72,10 @@ class GEdge : public GEntity {
// Get first derivative of edge at the given parameter.
virtual SVector3 firstDer(double par) const = 0;
+ // Get second derivative of edge at the given parameter.
+ // Default implentation using central differences
+ virtual SVector3 secondDer(double par) const ;
+
// Reparmaterize the point onto the given face.
virtual SPoint2 reparamOnFace(GFace *face, double epar,int dir) const = 0;
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index cff8a2dbbf..252c206793 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.14 2006-08-26 13:34:46 geuzaine Exp $
+// $Id: GFace.cpp,v 1.15 2006-08-26 15:13:22 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -341,3 +341,43 @@ void GFace::getMeanPlaneData(double VX[3], double VY[3],
y = meanPlane.y;
z = meanPlane.z;
}
+
+// X=X(u,v) Y=Y(u,v) Z=Z(u,v)
+// curv = div n = dnx/dx + dny/dy + dnz/dz
+
+// dnx/dx = dnx/du du/dx + dnx/dv dv/dx
+
+
+double GFace::curvature (const SPoint2 ¶m) const
+{
+
+ if (geomType() == Plane)return 0;
+
+ const double eps = 1.e-3;
+
+ Pair<SVector3,SVector3> der = firstDer(param) ;
+
+ SVector3 du = der.first();
+ SVector3 dv = der.second();
+ SVector3 nml = crossprod(du,dv);
+
+ double detJ = norm ( nml );
+
+ du.normalize();
+ dv.normalize();
+
+ SVector3 n1 = normal(SPoint2(param.x() - eps , param.y() )) ;
+ SVector3 n2 = normal(SPoint2(param.x() + eps , param.y() )) ;
+ SVector3 n3 = normal(SPoint2(param.x() , param.y() - eps )) ;
+ SVector3 n4 = normal(SPoint2(param.x() , param.y() + eps )) ;
+
+ SVector3 dndu = 500 * ( n2-n1 );
+ SVector3 dndv = 500 * ( n4-n3 );
+
+ double c = fabs(dot(dndu,du) + dot(dndv,dv)) / detJ;
+
+ // Msg (INFO,"c = %g detJ %g",c,detJ);
+
+ return c;
+
+}
diff --git a/Geo/GFace.h b/Geo/GFace.h
index b6ae1a5197..def562358c 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -95,6 +95,9 @@ class GFace : public GEntity
// Return the first derivate of the face at the parameter location.
virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const = 0;
+
+ // Return the curvature i.e. the divergence of the normal
+ virtual double curvature (const SPoint2 ¶m) const;
// True if the surface underlying the face is periodic and we need
// to worry about that.
diff --git a/Geo/SVector3.cpp b/Geo/SVector3.cpp
index c04b392f7e..d25b432232 100644
--- a/Geo/SVector3.cpp
+++ b/Geo/SVector3.cpp
@@ -1,4 +1,4 @@
-// $Id: SVector3.cpp,v 1.3 2006-08-15 06:26:52 geuzaine Exp $
+// $Id: SVector3.cpp,v 1.4 2006-08-26 15:13:22 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -23,7 +23,7 @@
#include "math.h"
-SVector3 cross(const SVector3 &a, const SVector3 &b)
+SVector3 crossprod(const SVector3 &a, const SVector3 &b)
{
return SVector3(a.y()*b.z()-b.y()*a.z(),
-(a.x()*b.z()-b.x()*a.z()),
@@ -44,7 +44,7 @@ double angle(const SVector3 &a, const SVector3 &b, const SVector3 &n)
ddot = -1.0;
double dang = acos(ddot);
// check if we just found the angle or 2*PI-angle
- SVector3 check = cross(a,b);
+ SVector3 check = crossprod(a,b);
// ***** check this below ******
if( norm(check) > 0.0000001*(norm(a)+norm(b))){ // check if a,b are colinear
double dir = dot(check,n);
diff --git a/Geo/SVector3.h b/Geo/SVector3.h
index 9f4777826a..9789905990 100644
--- a/Geo/SVector3.h
+++ b/Geo/SVector3.h
@@ -53,7 +53,7 @@ class SVector3 {
double norm(const SVector3 &v);
double dot(const SVector3 &a, const SVector3 &b);
-SVector3 cross(const SVector3 &a, const SVector3 &b);
+SVector3 crossprod(const SVector3 &a, const SVector3 &b);
double angle(const SVector3 &a, const SVector3 &b, const SVector3 &n);
SVector3 operator*(double m,const SVector3 &v);
SVector3 operator*(const SVector3 &v, double m);
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index 7a948ade6b..7b705f18ca 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshEdge.cpp,v 1.15 2006-08-15 06:26:53 geuzaine Exp $
+// $Id: gmshEdge.cpp,v 1.16 2006-08-26 15:13:22 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -34,6 +34,10 @@ extern Mesh *THEM;
gmshEdge::gmshEdge(GModel *model, Curve *edge, GVertex *v1, GVertex *v2)
: GEdge(model, edge->Num, v1, v2), c(edge)
{
+ meshAttributes.Method = c->Method;
+ meshAttributes.nbPointsTransfinite = c->ipar[0];
+ meshAttributes.coeffTransfinite = c->dpar[0];
+ meshAttributes.typeTransfinite = c->ipar[1];
}
gmshEdge::gmshEdge(GModel *model, int num)
diff --git a/Geo/gmshModel.cpp b/Geo/gmshModel.cpp
index 4ba63d8bc1..f0f48bf848 100644
--- a/Geo/gmshModel.cpp
+++ b/Geo/gmshModel.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshModel.cpp,v 1.13 2006-08-18 17:49:35 geuzaine Exp $
+// $Id: gmshModel.cpp,v 1.14 2006-08-26 15:13:22 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -41,6 +41,8 @@ void gmshModel::import()
GVertex *v = vertexByTag(p->Num);
if(!v){
v = new gmshVertex(this, p);
+ MVertex *mv = *(v->mesh_vertices.begin());
+ // Msg(INFO,"vertex %d %g %g %g -- %g %g %g",v->tag(),v->x(),v->y(),v->z(),mv->x(),mv->y(),mv->z());
add(v);
}
if(!p->Visible) v->setVisibility(0);
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 133d06e060..46a2dbeda8 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.58 2006-08-21 13:32:42 remacle Exp $
+// $Id: BDS.cpp,v 1.59 2006-08-26 15:13:22 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -429,6 +429,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));
+ // pp->radius() = (curvature ==0.0) ? 1.e22:1./curvature;
+ pp->radius() = curvature;
MAXPOINTNUMBER = (MAXPOINTNUMBER < num) ? num : MAXPOINTNUMBER;
return pp;
}
@@ -2105,6 +2108,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));
eit = p->edges.begin();
while(eit != p->edges.end()) {
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 4d1a462381..83da04f8ff 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -223,7 +223,7 @@ public:
class BDS_Point
{
- double radius_of_curvature;
+ double radius_of_curvature,_lc;
public:
double X,Y,Z; // Real COORDINATES
double u,v; // Parametric COORDINATES
@@ -234,7 +234,7 @@ public:
// just a transition
double & radius () {return radius_of_curvature;}
- double & lc () {return radius_of_curvature;}
+ double & lc () {return _lc;}
BDS_Vector N() const;
@@ -258,7 +258,7 @@ public:
void getTriangles(std::list<BDS_Triangle *> &t) const;
void compute_curvature();
BDS_Point(int id, double x=0, double y=0, double z=0)
- : radius_of_curvature(1.e22),X(x),Y(y),Z(z),u(0),v(0),config_modified(true),iD(id),g(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)
{
}
};
diff --git a/Mesh/Interpolation.cpp b/Mesh/Interpolation.cpp
index 16b51a20dd..23fb12bb8b 100644
--- a/Mesh/Interpolation.cpp
+++ b/Mesh/Interpolation.cpp
@@ -1,4 +1,4 @@
-// $Id: Interpolation.cpp,v 1.27 2006-01-06 00:34:26 geuzaine Exp $
+// $Id: Interpolation.cpp,v 1.28 2006-08-26 15:13:22 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -30,8 +30,9 @@
extern Mesh *THEM;
-// Curves
+// X = X(u), Y = Y(u), Z = Z(u)
+// Curves
Vertex InterpolateCurve(Curve * Curve, double u, int derivee)
{
int N, i, j;
diff --git a/Mesh/Utils.cpp b/Mesh/Utils.cpp
index b0cef3eeb6..4a1ededfc1 100644
--- a/Mesh/Utils.cpp
+++ b/Mesh/Utils.cpp
@@ -1,4 +1,4 @@
-// $Id: Utils.cpp,v 1.33 2006-08-21 13:32:42 remacle Exp $
+// $Id: Utils.cpp,v 1.34 2006-08-26 15:13:22 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -246,7 +246,7 @@ end:
}
-#define Precision 1.e-10
+#define Precision 1.e-8
#define MaxIter 25
#define NumInitGuess 11
@@ -369,7 +369,7 @@ void invert_singular_matrix3x3(double MM[3][3], double II[3][3])
void XYZtoUV(Surface * s, double X, double Y, double Z, double *U, double *V,
double relax)
{
- double Unew = 0., Vnew = 0., err;
+ double Unew = 0., Vnew = 0., err,err2;
int iter;
Vertex D_u, D_v, P;
double mat[3][3], jac[3][3];
@@ -383,8 +383,8 @@ void XYZtoUV(Surface * s, double X, double Y, double Z, double *U, double *V,
vmax = s->kv[s->OrderV + s->Nv];
}
else {
- umin = vmin = 0.0-1.e-6;
- umax = vmax = 1.0+1.e-6;
+ umin = vmin = 0.0-1.e-8;
+ umax = vmax = 1.0+1.e-8;
}
for(int i = 0; i < NumInitGuess; i++){
@@ -418,9 +418,9 @@ void XYZtoUV(Surface * s, double X, double Y, double Z, double *U, double *V,
(jac[0][1] * (X - P.Pos.X) + jac[1][1] * (Y - P.Pos.Y) +
jac[2][1] * (Z - P.Pos.Z));
- // err = DSQR(Unew - *U) + DSQR(Vnew - *V);
+ err = DSQR(Unew - *U) + DSQR(Vnew - *V);
// A BETTER TEST !! (JFR/AUG 2006)
- err = DSQR(X - P.Pos.X) + DSQR(Y - P.Pos.Y) + + DSQR(Z - P.Pos.Z);
+ err2 = DSQR(X - P.Pos.X) + DSQR(Y - P.Pos.Y) + DSQR(Z - P.Pos.Z);
iter++;
*U = Unew;
@@ -429,10 +429,11 @@ void XYZtoUV(Surface * s, double X, double Y, double Z, double *U, double *V,
- if(iter < MaxIter && err <= Precision &&
+ if(iter < MaxIter && err <= Precision && err2 <= 1.e-5 &&
Unew <= umax && Vnew <= vmax &&
Unew >= umin && Vnew >= vmin){
- //printf("converged for i=%d j=%d (err=%g iter=%d)\n", i, j, err, iter);
+ if (err2 > Precision)
+ Msg(WARNING,"converged for i=%d j=%d (err=%g iter=%d) BUT err2 = %g", i, j, err, iter,err2);
return;
}
}
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 173618737a..92354a1e3e 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdge.cpp,v 1.13 2006-08-21 13:32:42 remacle Exp $
+// $Id: meshGEdge.cpp,v 1.14 2006-08-26 15:13:22 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -21,6 +21,7 @@
#include "meshGEdge.h"
#include "GEdge.h"
+#include "GFace.h"
#include "Gmsh.h"
#include "Utils.h"
#include "Mesh.h"
@@ -44,20 +45,34 @@ double F_LC_ANALY (double xx, double yy, double zz)
return 0.05 + .1*fabs(xx*yy) ;
}
+double max_surf_curvature ( GPoint & gp )
+{
+ std::list<GFace *> faces = _myGEdge->faces();
+ std::list<GFace *>::iterator it = faces.begin();
+ double curv = 0;
+ while (it != faces.end())
+ {
+ SPoint2 par = (*it)->parFromPoint(SPoint3 (gp.x(),gp.y(),gp.z()));
+ curv = std::max(curv, (*it)->curvature ( par ) );
+ ++it;
+ }
+ return curv;
+}
+
double F_Lc_bis(double t)
{
+ // const double nb_points_per_radius_of_curv = 2;
+ GPoint point = _myGEdge -> point(t) ;
const double fact = (t-t_begin)/(t_end-t_begin);
double lc_here = lc_begin + fact * (lc_end-lc_begin);
- SVector3 der = _myGEdge -> firstDer(t) ;
- const double d = norm(der)/CTX.mesh.lc_factor;
+ SVector3 der = _myGEdge -> firstDer(t) ;
+ const double d = norm(der);
- // TESTTTT
- GPoint points = _myGEdge -> point(t) ;
- // double l_bgm = F_LC_ANALY(points.x(),points.y(),points.z()) ;
- // lc_here = l_bgm;
+ // double curv = max_surf_curvature ( point );
+ // if (curv != 0)
+ // lc_here = std::min( 1./(curv * nb_points_per_radius_of_curv),lc_here);
if(CTX.mesh.bgmesh_type == ONFILE) {
- GPoint point = _myGEdge -> point(t) ;
const double Lc = BGMXYZ(point.x(), point.y(), point.z());
if(CTX.mesh.constrained_bgmesh)
return std::max(d / Lc, d / lc_here);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index f8e14e3ec0..9438f876fa 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.7 2006-08-21 13:32:42 remacle Exp $
+// $Id: meshGFace.cpp,v 1.8 2006-08-26 15:13:22 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -219,7 +219,11 @@ extern double F_LC_ANALY (double xx, double yy, double zz);
double NewGetLc ( BDS_Point *p )
{
- return p->lc();
+ // const double curv = p->radius();
+ // if (curv == 0.0)
+ return p->lc();
+ // else
+ // return std::min(p->lc(),1./(2*curv));
return F_LC_ANALY(p->X,p->Y,p->Z);
if(BGMExists())
{
@@ -311,10 +315,12 @@ void OptimizeMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
}
{
// swap edges that provide a better configuration
- std::list<BDS_Edge *> temp (m.edges);
- std::list<BDS_Edge*>::iterator it = temp.begin();
- while (it != temp.end())
+ int NN1 = m.edges.size();
+ int NN2 = 0;
+ std::list<BDS_Edge*>::iterator it = m.edges.begin();
+ while (1)
{
+ if (NN2++ >= NN1)break;
if (!(*it)->deleted && (*it)->numfaces() == 2)
{
BDS_Point *op[2];
@@ -374,8 +380,6 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
// ++ittt;
// }
- std::list<BDS_Edge *> temp (m.edges);
- std::list<BDS_Edge*>::iterator it = temp.begin();
int nb_split=0;
int nb_smooth=0;
int nb_collaps=0;
@@ -384,8 +388,12 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
// split long edges
double minL=1.E22,maxL=0;
- while (it != temp.end())
+ int NN1 = m.edges.size();
+ int NN2 = 0;
+ std::list<BDS_Edge*>::iterator it = m.edges.begin();
+ while (1)
{
+ if (NN2++ >= NN1)break;
if (!(*it)->deleted)
{
(*it)->p1->config_modified = false;
@@ -398,12 +406,15 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
}
Msg(STATUS2," %d triangles : conv %g -> %g (%g sec)",m.triangles.size(),maxL,1.5,(double)(clock()-t1)/CLOCKS_PER_SEC);
- if ((minL > 0.2 && maxL < 1.5) || IT > NIT)break;
+ if ((minL > 0.4 && maxL < 1.5) || IT > NIT)break;
- it = temp.begin();
- while (it != temp.end())
+ NN1 = m.edges.size();
+ NN2 = 0;
+ it = m.edges.begin();
+ while (1)
{
+ if (NN2++ >= NN1)break;
if (!(*it)->deleted)
{
double lone = NewGetLc ( *it);
@@ -423,10 +434,12 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
}
// swap edges that provide a better configuration
- temp = m.edges;
- it = temp.begin();
- while (it != temp.end())
+ NN1 = m.edges.size();
+ NN2 = 0;
+ it = m.edges.begin();
+ while (1)
{
+ if (NN2++ >= NN1)break;
if (!(*it)->deleted && edgeSwapTest(*it))
if (m.swap_edge ( *it , BDS_SwapEdgeTestParametric()))
nb_swap++;
@@ -434,10 +447,12 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
}
// collapse short edges
- temp = m.edges;
- it = temp.begin();
- while (it != temp.end())
+ NN1 = m.edges.size();
+ NN2 = 0;
+ it = m.edges.begin();
+ while (1)
{
+ if (NN2++ >= NN1)break;
double lone = NewGetLc ( *it);
if (!(*it)->deleted && (*it)->numfaces() == 2 && lone < 0.6 )
{
@@ -453,10 +468,12 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
}
// swap edges that provide a better configuration
- temp = m.edges;
- it = temp.begin();
- while (it != temp.end())
+ NN1 = m.edges.size();
+ NN2 = 0;
+ it = m.edges.begin();
+ while (1)
{
+ if (NN2++ >= NN1)break;
if (!(*it)->deleted && edgeSwapTest(*it))
if (m.swap_edge ( *it , BDS_SwapEdgeTestParametric()))
nb_swap++;
@@ -694,7 +711,7 @@ bool recover_medge ( BDS_Mesh *m, GEdge *ge)
void gmsh2DMeshGenerator ( GFace *gf )
{
- // if(gf->tag() == 138|| gf->tag() == 196|| gf->tag() == 200|| gf->tag() == 262) return;
+ // if(gf->tag() != 136) return;
std::set<MVertex*>all_vertices;
std::map<int, MVertex*>numbered_vertices;
@@ -707,6 +724,12 @@ void gmsh2DMeshGenerator ( GFace *gf )
{
all_vertices.insert ( (*it)->mesh_vertices.begin() , (*it)->mesh_vertices.end() );
all_vertices.insert ( (*it)->getBeginVertex()->mesh_vertices.begin() , (*it)->getBeginVertex()->mesh_vertices.end() );
+// GVertex *v1 = (*it)->getBeginVertex();
+// GVertex *v2 = (*it)->getEndVertex();
+// MVertex *mv1 = *(v1->mesh_vertices.begin());
+// MVertex *mv2 = *(v2->mesh_vertices.begin());
+// Msg(INFO,"%d -> %d = %g %g %g -- %g %g %g",v1->tag(),v2->tag(),v1->x(),v1->y(),v1->z(),v2->x(),v2->y(),v2->z());
+// Msg(INFO,"%g %g %g -- %g %g %g",mv1->x(),mv1->y(),mv1->z(),mv2->x(),mv2->y(),mv2->z());
all_vertices.insert ( (*it)->getEndVertex()->mesh_vertices.begin() , (*it)->getEndVertex()->mesh_vertices.end() );
++it;
}
@@ -719,6 +742,30 @@ void gmsh2DMeshGenerator ( GFace *gf )
++it;
}
+
+ double * X_ = new double [all_vertices.size()];
+ double * Y_ = new double [all_vertices.size()];
+ double * Z_ = new double [all_vertices.size()];
+
+ std::set<MVertex*>::iterator itv = all_vertices.begin();
+
+ // FILE *fdeb = fopen("debug.dat","w");
+ // fprintf(fdeb,"surface %d\n" ,gf->tag());
+ int count = 0;
+ while(itv != all_vertices.end())
+ {
+ MVertex *here = *itv;
+ // fprintf(fdeb,"%d %g %g %g\n" ,here->getNum(),here->x(),here->y(),here->z());
+ X_[count] = here->x();
+ Y_[count] = here->y();
+ Z_[count] = here->z();
+ count ++;
+ ++itv;
+ }
+
+ // fclose (fdeb);
+
+
// mesh the domain in the parametric space ->
// project all points in their parametric space
fromCartesianToParametric c2p ( gf );
@@ -728,7 +775,7 @@ void gmsh2DMeshGenerator ( GFace *gf )
// I do not have SBoundingBox, so I use a 3D one...
// At the same time, number the vertices locally
SBoundingBox3d bbox;
- std::set<MVertex*>::iterator itv = all_vertices.begin();
+ itv = all_vertices.begin();
int NUM = 0;
while(itv != all_vertices.end())
{
@@ -746,6 +793,8 @@ void gmsh2DMeshGenerator ( GFace *gf )
/// Fill the DocRecord Data Structure with the points
DocRecord doc;
+
+
doc.points = (PointRecord*)malloc((all_vertices.size()+4) * sizeof(PointRecord));
itv = all_vertices.begin();
int j = 0;
@@ -764,6 +813,7 @@ void gmsh2DMeshGenerator ( GFace *gf )
++itv;
}
+
/// Increase the size of the bounding box by 20 %
bbox *= 1.5;
/// add 4 points than encloses the domain
@@ -820,13 +870,6 @@ void gmsh2DMeshGenerator ( GFace *gf )
recover_medge ( m, *it);
++it;
}
- it = emb_edges.begin();
- while(it != emb_edges.end())
- {
- recover_medge ( m, *it);
- ++it;
- }
-
// Msg(INFO,"Boundary Edges recovered for surface %d",gf->tag());
// Look for an edge that is on the boundary for which one of the
// two neighbors has a negative number node. The other triangle
@@ -856,6 +899,16 @@ void gmsh2DMeshGenerator ( GFace *gf )
++ite;
}
}
+
+
+ it = emb_edges.begin();
+ while(it != emb_edges.end())
+ {
+ recover_medge ( m, *it);
+ ++it;
+ }
+
+
// delete useless stuff
{
std::list<BDS_Triangle*>::iterator itt = m->triangles.begin();
@@ -894,7 +947,7 @@ void gmsh2DMeshGenerator ( GFace *gf )
// goto hhh;
// start mesh generation
- RefineMesh (gf,*m,27);
+ RefineMesh (gf,*m,20);
OptimizeMesh (gf,*m,2);
// fill the small gmsh structures
@@ -934,15 +987,32 @@ void gmsh2DMeshGenerator ( GFace *gf )
// delete the mesh
- // outputScalarField(m->triangles, "lc.pos");
+ // char name[245];
+ // sprintf(name,"s%d.pos",gf->tag());
+ // outputScalarField(m->triangles, name);
delete m;
-
+
+
+ itv = all_vertices.begin();
+ count = 0;
+ while(itv != all_vertices.end())
+ {
+ MVertex *here = *itv;
+ here->x() = X_[count] ;
+ here->y() = Y_[count] ;
+ here->z() = Z_[count] ;
+ count ++;
+ ++itv;
+ }
+ delete [] X_;
+ delete [] Y_;
+ delete [] Z_;
fromParametricToCartesian p2c ( gf );
- std::for_each(all_vertices.begin(),all_vertices.end(),p2c);
+ // std::for_each(all_vertices.begin(),all_vertices.end(),p2c);
std::for_each(gf->mesh_vertices.begin(),gf->mesh_vertices.end(),p2c);
}
-
+
diff --git a/Parser/OpenFile.cpp b/Parser/OpenFile.cpp
index 098b3155f2..a649313a08 100644
--- a/Parser/OpenFile.cpp
+++ b/Parser/OpenFile.cpp
@@ -1,4 +1,4 @@
-// $Id: OpenFile.cpp,v 1.117 2006-08-26 13:34:46 geuzaine Exp $
+// $Id: OpenFile.cpp,v 1.118 2006-08-26 15:13:22 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -205,8 +205,6 @@ int ParseFile(char *f, int close, int warn_if_missing)
if(close)
fclose(yyin);
- GMODEL->import();
-
strncpy(yyname, yyname_old, 255);
yyin = yyin_old;
yyerrorstate = yyerrorstate_old;
@@ -332,6 +330,8 @@ int MergeProblem(char *name, int warn_if_missing)
}
}
+ GMODEL->import();
+
SetBoundingBox();
CTX.mesh.changed = ENT_ALL;
Msg(STATUS2, "Read '%s'", name);
diff --git a/benchmarks/2d/Square-01.geo b/benchmarks/2d/Square-01.geo
index 9d91879787..0feef348ca 100644
--- a/benchmarks/2d/Square-01.geo
+++ b/benchmarks/2d/Square-01.geo
@@ -14,3 +14,5 @@ Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5};
//Attractor Point{2} = {0.05,0.05,2};
//Mesh.Algorithm = 2;
+Transfinite Line {4,-2} = 130 Using Bump 5;
+Transfinite Line {1,3} = 5 Using Progression 1;
--
GitLab