diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h new file mode 100644 index 0000000000000000000000000000000000000000..51dd647f5c41423120019feb5da9b2240c6c3c06 --- /dev/null +++ b/Common/GmshDefines.h @@ -0,0 +1,58 @@ +#ifndef _GMSH_DEFINES_H_ +#define _GMSH_DEFINES_H_ +#define FORMAT_MSH 1 +#define FORMAT_UNV 2 +#define FORMAT_GREF 3 +#define FORMAT_XPM 4 +#define FORMAT_PS 5 +#define FORMAT_BMP 6 +#define FORMAT_GIF 7 +#define FORMAT_GEO 8 +#define FORMAT_JPEG 9 +#define FORMAT_AUTO 10 +#define FORMAT_PPM 11 +#define FORMAT_YUV 12 +#define FORMAT_DMG 13 +#define FORMAT_SMS 14 +#define FORMAT_OPT 15 +#define FORMAT_VTK 16 +#define FORMAT_JPEGTEX 17 +#define FORMAT_TEX 18 +#define FORMAT_VRML 19 +#define FORMAT_EPS 20 +#define FORMAT_EPSTEX 21 +#define FORMAT_PNG 22 +#define FORMAT_PNGTEX 23 +#define FORMAT_PDF 24 +#define FORMAT_PDFTEX 25 +#define FORMAT_LC 26 +#define FORMAT_STL 27 +#define FORMAT_P3D 28 + +#define CONV_VALUE 0.8 + +#define VIS_GEOM (1<<0) +#define VIS_MESH (1<<1) + +#define NOTTOLINK 1 +#define TOLINK 2 + +#define BOF 1 +#define A_TOUT_PRIX 2 + +#define CENTER_CIRCCIRC 1 +#define BARYCENTER 2 + +#define EXTERN 1 +#define INTERN 2 + +#define ONFILE 2 +#define WITHPOINTS 3 + +#define TRANSFINI 1 +#define LIBRE 2 +#define ELLIPTIC 3 + +#define NB_HISTOGRAM 100 + +#endif diff --git a/Fltk/Main.cpp b/Fltk/Main.cpp index e89a5d05baeac653e25a15b9a708b9e8dae875c3..3838c8a6a22f9d0914d9a150e8f889919a894f89 100644 --- a/Fltk/Main.cpp +++ b/Fltk/Main.cpp @@ -1,4 +1,4 @@ -// $Id: Main.cpp,v 1.88 2006-03-17 21:04:34 geuzaine Exp $ +// $Id: Main.cpp,v 1.89 2006-07-11 13:41:22 remacle Exp $ // // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle // @@ -39,10 +39,13 @@ #include "CommandLine.h" #include "Numeric.h" #include "Solvers.h" +#include "GModel.h" Context_T CTX; Mesh M, *THEM = NULL; GUI *WID = NULL; +GModel *GMODEL = 0; + int main(int argc, char *argv[]) { diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp index dbd1550110bbb6db683fa1c99d36e591bbf42a33..1f69045e0fa2e77112540f6a851c0630a1bd16ad 100644 --- a/Geo/GEdge.cpp +++ b/Geo/GEdge.cpp @@ -1,4 +1,5 @@ #include "GEdge.h" +#include "GmshDefines.h" #include <algorithm> void GEdge::addFace ( GFace *e ) { @@ -13,10 +14,11 @@ GEdge::GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1) - : GEntity (model,tag),v0(_v0),v1(_v1) + : GEntity (model,tag),v0(_v0),v1(_v1) { v0->addEdge (this); v1->addEdge (this); + meshGEdgeAttributes.Method = LIBRE; } GEdge::~GEdge() diff --git a/Geo/GEdge.h b/Geo/GEdge.h index aef62e3084b16739d3f62bc4eb326ade6b333829..897be18801419e18a935657d05a30c392d188b00 100644 --- a/Geo/GEdge.h +++ b/Geo/GEdge.h @@ -44,18 +44,29 @@ public: void delFace ( GFace *f ); GVertex * getBeginVertex () const - { - return v0; - } + { + return v0; + } GVertex * getEndVertex () const - { - return v1; - } - + { + return v1; + } + + struct + { + char Method; + double coeffTransfinite; + int nbPointsTransfinite; + int typeTransfinite; + } meshGEdgeAttributes ; + + virtual int minimumEdgeSegments () const {return 1;} + protected: GVertex *v0,*v1; std::list<GFace *> l_faces; + }; diff --git a/Geo/GEntity.h b/Geo/GEntity.h index b7cf855b59ebfe21d80fb41d85c3431bbfe8cadd..5fd4ea7ad7268edd6d149b005bf69493216a5ce4 100644 --- a/Geo/GEntity.h +++ b/Geo/GEntity.h @@ -4,7 +4,9 @@ #include "Range.h" #include "SPoint3.h" #include "SBoundingBox3d.h" +#include "MVertex.h" #include <list> +#include <vector> class GModel; class GVertex; @@ -94,6 +96,8 @@ public: // The bounding box SBoundingBox3d bounds() const{throw;} + std::vector<MVertex*> mesh_vertices; + }; diff --git a/Geo/GModel.h b/Geo/GModel.h index 17c7bee1ed342048ad6b4dc274b9380b3aad09f0..7907daaa9d1e7a431a9b88a9e31bf62b119f4e24 100644 --- a/Geo/GModel.h +++ b/Geo/GModel.h @@ -15,8 +15,8 @@ public: virtual ~GModel() {} typedef std::list<GRegion*>::iterator riter; - typedef std::list<GFace*>::iterator fiter; - typedef std::list<GEdge*>::iterator eiter; + typedef std::list<GFace*> ::iterator fiter; + typedef std::list<GEdge*> ::iterator eiter; typedef std::list<GVertex*>::iterator viter; /** Returns the geometric tolerance for the entire model. */ diff --git a/Geo/GPoint.h b/Geo/GPoint.h index 31deb8b6df5040939b6b1986c76552981869462d..2d35451d5ecb77e611e226e80c418ee8de14f058 100644 --- a/Geo/GPoint.h +++ b/Geo/GPoint.h @@ -2,11 +2,15 @@ #define _GPOINT_H_ class GEntity; -struct GPoint +class GPoint { double X,Y,Z; const GEntity *e; double par[2]; +public: + inline double x() const {return X;} + inline double y() const {return Y;} + inline double z() const {return Z;} GPoint (double _x, double _y, double _z, const GEntity *onwhat) : X(_x), Y(_y), Z(_z), e(onwhat){ } diff --git a/Geo/GVertex.h b/Geo/GVertex.h index a54fc54d30c2d1cb4ab2367f987383543aa71edf..b9ec60ef32df3e04f3d786339f1bd95f5496e12a 100644 --- a/Geo/GVertex.h +++ b/Geo/GVertex.h @@ -15,6 +15,7 @@ public: void delEdge ( GEdge *e ); virtual int dim() const {return 0;} virtual GeomType geomType() const {return Point;} + virtual double prescribedMeshSizeAtVertex() const {return 0;} protected: std::list<GEdge*> l_edges; diff --git a/Geo/MVertex.h b/Geo/MVertex.h new file mode 100644 index 0000000000000000000000000000000000000000..fd35b92cfe09f603773fbd9295a2195887faeb85 --- /dev/null +++ b/Geo/MVertex.h @@ -0,0 +1,24 @@ +#ifndef _MVERTEX_H_ +#define _MVERTEX_H_ +class GEntity; +class MVertex +{ + double x,y,z; + GEntity *ge; +public : + MVertex ( double _x, double _y, double _z , GEntity * _ge ) : + x(_x),y(_y),z(_z),ge(_ge) + { + } +}; + +class MEdgeVertex : public MVertex +{ + double u; +public : + MEdgeVertex ( double _x, double _y, double _z , GEntity * _ge , double _u) : + MVertex (_x,_y,_z,_ge),u(_u) + { + } +}; +#endif diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp index 55803409800d07665bb322d477ac68ef4baa39da..6ceda771e2fdad4ae03117bc35af89e45b911cd0 100644 --- a/Geo/gmshEdge.cpp +++ b/Geo/gmshEdge.cpp @@ -3,6 +3,9 @@ #include "Interpolation.h" #include "CAD.h" #include "Geo.h" +#include "Context.h" + +extern Context_T CTX; gmshEdge::gmshEdge(GModel *model,Curve *edge,GVertex *v1,GVertex *v2) : GEdge ( model, edge->Num, v1, v2 ), c(edge) @@ -134,3 +137,16 @@ int gmshEdge::containsPoint(const SPoint3 &pt) const { throw; } + +int gmshEdge::minimumEdgeSegments () const +{ + if(c->Typ == MSH_SEGM_CIRC || + c->Typ == MSH_SEGM_CIRC_INV || + c->Typ == MSH_SEGM_ELLI || + c->Typ == MSH_SEGM_ELLI_INV) { + return (int)(fabs(c->Circle.t1 - c->Circle.t2) * + (double)CTX.mesh.min_circ_points / Pi) - 1; + } + else + return GEdge::minimumEdgeSegments () ; +} diff --git a/Geo/gmshEdge.h b/Geo/gmshEdge.h index 5b58c7cc57f9e28b064af68c239d44fce81adc59..40df2172a1a0d257532e69288ea65dfbcf49a1fe 100644 --- a/Geo/gmshEdge.h +++ b/Geo/gmshEdge.h @@ -29,6 +29,8 @@ public: } void * getNativePtr() const; virtual double parFromPoint(const SPoint3 &pt) const; + virtual int minimumEdgeSegments () const; + protected: Curve *c; }; diff --git a/Geo/gmshVertex.h b/Geo/gmshVertex.h index e581fb7a1b81bdf0caa8592701e73436ea42fe70..db69b6f32cd187091614411ea9d2752eef7a764b 100644 --- a/Geo/gmshVertex.h +++ b/Geo/gmshVertex.h @@ -9,15 +9,12 @@ class gmshVertex : public GVertex { public: gmshVertex(GModel *m, Vertex *_v) : GVertex(m, _v->Num), v(_v){} virtual ~gmshVertex() {} - - virtual GPoint point() const - { - return GPoint ( v->Pos.X,v->Pos.Y,v->Pos.Z,this); - } + virtual GPoint point() const{return GPoint ( v->Pos.X,v->Pos.Y,v->Pos.Z,this);} virtual double tolerance() const {return 1.e-14;} void * getNativePtr() const {return v;} - Vertex *v; + virtual double prescribedMeshSizeAtVertex () const {return v->lc;} protected: + Vertex *v; }; #endif diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp index 0ae02dca4869dfc87a13c90eb66b45c32a23d105..55e4c4f34add30871f355dad18cbc05a6698092e 100644 --- a/Mesh/Generator.cpp +++ b/Mesh/Generator.cpp @@ -1,4 +1,4 @@ -// $Id: Generator.cpp,v 1.82 2006-03-08 17:04:59 remacle Exp $ +// $Id: Generator.cpp,v 1.83 2006-07-11 13:41:22 remacle Exp $ // // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle // @@ -29,9 +29,12 @@ #include "Views.h" #include "PartitionMesh.h" #include "OS.h" +#include "meshGEdge.h" +#include "GModel.h" extern Mesh *THEM; extern Context_T CTX; +extern GModel *GMODEL; static int nbOrder2 = 0; @@ -265,6 +268,8 @@ void Maillage_Dimension_1(Mesh * M) Tree_Action(M->Curves, Maillage_Curve); + std::for_each (GMODEL->firstEdge(),GMODEL->lastEdge(), meshGEdge() ); + double t2 = Cpu(); M->timing[0] = t2 - t1; } diff --git a/Mesh/Makefile b/Mesh/Makefile index c6ccbb97ebae2c74a726c837da31689bc7b60828..470ebcbd10f5c35efbb3659e869b7a5ad93cac3d 100644 --- a/Mesh/Makefile +++ b/Mesh/Makefile @@ -1,4 +1,4 @@ -# $Id: Makefile,v 1.108 2006-07-10 12:16:35 remacle Exp $ +# $Id: Makefile,v 1.109 2006-07-11 13:41:22 remacle Exp $ # # Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle # @@ -66,6 +66,7 @@ SRC = 1D_Mesh.cpp \ SwapEdge.cpp \ Utils.cpp \ Metric.cpp \ + meshGEdge.cpp \ Nurbs.cpp \ Interpolation.cpp \ SecondOrder.cpp \ diff --git a/Mesh/Mesh.h b/Mesh/Mesh.h index 5b8a07c6476c21df42d31b59d9f628a517f58039..412021f2aac5f8b30fc41b98a49089ac99aa6d82 100644 --- a/Mesh/Mesh.h +++ b/Mesh/Mesh.h @@ -29,61 +29,7 @@ #include "ExtrudeParams.h" #include "VertexArray.h" #include "SmoothNormals.h" - -#define FORMAT_MSH 1 -#define FORMAT_UNV 2 -#define FORMAT_GREF 3 -#define FORMAT_XPM 4 -#define FORMAT_PS 5 -#define FORMAT_BMP 6 -#define FORMAT_GIF 7 -#define FORMAT_GEO 8 -#define FORMAT_JPEG 9 -#define FORMAT_AUTO 10 -#define FORMAT_PPM 11 -#define FORMAT_YUV 12 -#define FORMAT_DMG 13 -#define FORMAT_SMS 14 -#define FORMAT_OPT 15 -#define FORMAT_VTK 16 -#define FORMAT_JPEGTEX 17 -#define FORMAT_TEX 18 -#define FORMAT_VRML 19 -#define FORMAT_EPS 20 -#define FORMAT_EPSTEX 21 -#define FORMAT_PNG 22 -#define FORMAT_PNGTEX 23 -#define FORMAT_PDF 24 -#define FORMAT_PDFTEX 25 -#define FORMAT_LC 26 -#define FORMAT_STL 27 -#define FORMAT_P3D 28 - -#define CONV_VALUE 0.8 - -#define VIS_GEOM (1<<0) -#define VIS_MESH (1<<1) - -#define NOTTOLINK 1 -#define TOLINK 2 - -#define BOF 1 -#define A_TOUT_PRIX 2 - -#define CENTER_CIRCCIRC 1 -#define BARYCENTER 2 - -#define EXTERN 1 -#define INTERN 2 - -#define ONFILE 2 -#define WITHPOINTS 3 - -#define TRANSFINI 1 -#define LIBRE 2 -#define ELLIPTIC 3 - -#define NB_HISTOGRAM 100 +#include "GmshDefines.h" class BDS_Mesh; diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp new file mode 100644 index 0000000000000000000000000000000000000000..28d91e5e5b226c450b2d5cc2e7f5ce5c5709ec4f --- /dev/null +++ b/Mesh/meshGEdge.cpp @@ -0,0 +1,171 @@ +#include "meshGEdge.h" +#include "GEdge.h" +#include "Gmsh.h" +#include "Utils.h" +#include "Mesh.h" +#include "Context.h" +#include "Message.h" + +static GEdge * _myGEdge; +static double _myGEdgeLength, t_begin,t_end,lc_begin,lc_end; +static Range<double> _myGEdgeBounds; +// boooooooh !!!!!! +extern Mesh *THEM; +extern Context_T CTX; + +double F_Lc_bis(double t) +{ + const double fact = (t-t_begin)/(t_end-t_begin); + const double lc_here = lc_begin + fact * (lc_end-lc_begin); + SVector3 der = _myGEdge -> firstDer(t) ; + const double d = norm(der); + + if(THEM->BackgroundMeshType == 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); + else + return d / Lc; + } + else + return d/lc_here; +} + +double F_Transfini_bis(double t) +{ + double val,r; + int i; + + SVector3 der = _myGEdge -> firstDer(t) ; + double d = norm(der); + + double coef = _myGEdge->meshGEdgeAttributes.coeffTransfinite; + int type = _myGEdge->meshGEdgeAttributes.typeTransfinite; + int nbpt = _myGEdge->meshGEdgeAttributes.nbPointsTransfinite; + + if(coef <= 0.0 || coef == 1.0) { + // coef < 0 should never happen + val = d * coef / _myGEdgeLength; + } + else { + switch (abs(type)) { + + case 1: // Geometric progression ar^i; Sum of n terms = THEC->l = a (r^n-1)/(r-1) + if(sign(type) >= 0) + r = coef; + else + r = 1. / coef; + double a = _myGEdgeLength * (r - 1.) / (pow(r, nbpt - 1.) - 1.); + int i = (int)(log(t * _myGEdgeLength / a * (r - 1.) + 1.) / log(r)); + val = d / (a * pow(r, (double)i)); + break; + + case 2: // Bump + if(coef > 1.0) { + a = -4. * sqrt(coef - 1.) * + atan2(1., sqrt(coef - 1.)) / + ((double)nbpt * _myGEdgeLength); + } + else { + a = 2. * sqrt(1. - coef) * + log(fabs((1. + 1. / sqrt(1. - coef)) + / (1. - 1. / sqrt(1. - coef)))) + / ((double)nbpt * _myGEdgeLength); + } + double b = -a * _myGEdgeLength * _myGEdgeLength / (4. * (coef - 1.)); + val = d / (-a * DSQR(t * _myGEdgeLength - (_myGEdgeLength) * 0.5) + b); + break; + + default: + Msg(WARNING, "Unknown case in Transfinite Line mesh"); + val = 1.; + } + } + return val; +} + +double F_One_bis(double t) +{ + SVector3 der = _myGEdge -> firstDer(t) ; + return norm(der); +} + + +void deMeshGEdge :: operator() (GEdge *ge) +{ + for (int i=0;i<ge->mesh_vertices.size();i++) delete ge->mesh_vertices[i]; + ge->mesh_vertices.clear(); +} + +void meshGEdge :: operator() (GEdge *ge) +{ + + deMeshGEdge dem; + dem(ge); + + // Send a messsage to the GMSH environment + Msg(STATUS3, "Meshing curve %d", ge->tag()); + + // Create a list of integration points + List_T *Points = List_Create(10, 10, sizeof(IntPoint)); + + // For avoiding the global variable : + // We have to change the Integration function in order + // to pass an extra argument... + _myGEdge = ge; + + // first compute the length of the curve by integrating one + _myGEdgeBounds = ge->parBounds(0) ; + _myGEdgeLength = Integration(_myGEdgeBounds.low(), _myGEdgeBounds.high(), F_One_bis, Points, 1.e-4); + List_Reset(Points); + + + lc_begin = _myGEdge->getBeginVertex()->prescribedMeshSizeAtVertex(); + lc_end = _myGEdge->getEndVertex()->prescribedMeshSizeAtVertex(); + + t_begin = _myGEdgeBounds.low(); + t_end = _myGEdgeBounds.high(); + + // Integrate ∫ detJ/lc du + double a; + int N; + if(ge->meshGEdgeAttributes.Method == TRANSFINI) + { + a = Integration(_myGEdgeBounds.low(), _myGEdgeBounds.high(), F_Transfini_bis, Points, 1.e-7); + N = ge->meshGEdgeAttributes.nbPointsTransfinite; + } + else + { + a = Integration(_myGEdgeBounds.low(), _myGEdgeBounds.high(), F_Lc_bis, Points, 1.e-4); + N = std::max (ge->minimumEdgeSegments()+1, (int)(a + 1.)); + } + const double b = a / (double)(N - 1); + + int count = 1, NUMP = 1; + IntPoint P1,P2; + + // do not consider the first and the last vertex + // those are not classified on this mesh edge + ge->mesh_vertices.resize(N-2); + + while(NUMP < N - 1) { + List_Read(Points, count - 1, &P1); + List_Read(Points, count, &P2); + const double d = (double)NUMP *b; + + if((fabs(P2.p) >= fabs(d)) && (fabs(P1.p) < fabs(d))) { + double dt = P2.t - P1.t; + double dp = P2.p - P1.p; + double t = P1.t + dt / dp * (d - P1.p); + GPoint V = ge -> point(t) ; + MEdgeVertex *ev = new MEdgeVertex ( V.x(), V.y(), V.z(), ge, t ); + ge->mesh_vertices [NUMP-1] = ev; + NUMP++; + } + else { + count++; + } + } + List_Delete(Points); +}