[Gmsh] Simple crack: gmsh is fast, c++ is slow

Christophe Geuzaine cgeuzaine at ulg.ac.be
Thu Apr 23 21:47:01 CEST 2015


> On 23 Apr 2015, at 15:57, Franco Milicchio <fmilicchio at me.com> wrote:
> 
> Dear all,
> 
> I just started using gmsh, and I am finding difficulties. My situation is easy right now: write a code that meshes a cracked domain.
> 
> As far as I’ve read, I can use a C++ code to create and mesh a domain. However, it seems that gmsh is way faster than my code, even if compiled in release. So, let’s say I have this .geo file (modified from t1.geo):
> 
> lc = 1;
> 
> Point(1) = {  0.0,  0.0,  0.0, lc } ;
> Point(2) = { 10.0,  0.0,  0.0, lc } ;
> Point(3) = { 10.0, 20.0,  0.0, lc } ;
> Point(4) = {  0.0, 20.0,  0.0, lc } ;
> Point(5) = {  0.0, 11.1,  0.0, lc } ;
> Point(6) = {  5.0, 10.0,  0.0, lc } ;
> Point(7) = {  0.0,  8.9,  0.0, lc } ;
> 
> Line(1) = {1, 2} ;
> Line(2) = {2, 3} ;
> Line(3) = {3, 4} ;
> Line(4) = {4, 5} ;
> Line(5) = {5, 6} ;
> Line(6) = {6, 7} ;
> Line(7) = {7, 1} ;
> 
> Line Loop(5) = { 1, 2, 3, 4, 5, 6, 7 } ;
> 
> Plane Surface(6) = { 5 } ;
> 
> 
> The gmsh tool outputs a very nice mesh that I can visualize in paraview. Cool, it also takes very few seconds:
> 
> Info    : Running 'gmsh -2 a.geo -format vtk' [Gmsh 2.9.2, 1 node, max. 1 thread]
> Info    : Started on Thu Apr 23 15:26:57 2015
> Info    : Reading 'a.geo'...
> Info    : Done reading 'a.geo'
> Info    : Meshing 1D...
> Info    : Meshing curve 1 (Line)
> Info    : Meshing curve 2 (Line)
> Info    : Meshing curve 3 (Line)
> Info    : Meshing curve 4 (Line)
> Info    : Meshing curve 5 (Line)
> Info    : Meshing curve 6 (Line)
> Info    : Meshing curve 7 (Line)
> Info    : Done meshing 1D (0.001055 s)
> Info    : Meshing 2D...
> Info    : Meshing surface 6 (Plane, Delaunay)
> Info    : Done meshing 2D (0.022495 s)
> Info    : 300 vertices 605 elements
> Info    : Writing 'a.vtk'...
> Info    : Done writing 'a.vtk'
> Info    : Stopped on Thu Apr 23 15:26:58 2015
> 
> 
> Ok, now I’ve developed my code, with two domains. The non-compiled one as you can see is a square, and it works like a charm. On the other hand, my cracked domain is a complete failure. So, here’s my test code:
> 

Have you tried by simply using the gmshVertex, gmshEdge, gmshFace classes ?


>     auto t = std::time(nullptr);
>     auto tm = *std::localtime(&t);
>     
>     GmshInitialize();
>     
>     std::cout << "start @ " << std::put_time(&tm, "%d-%m-%Y %H-%M-%S") << std::endl;
>     
>     GModel *m = new GModel("plate");
>     
> #if 0
>     
>     // Add vertices
>     vertex* v0 = new vertex(m, 0, 0.0, 0.0, 0.0);
>     vertex* v1 = new vertex(m, 1, 1.0, 0.0, 0.0);
>     vertex* v2 = new vertex(m, 2, 1.0, 1.0, 0.0);
>     vertex* v3 = new vertex(m, 3, 0.0, 1.0, 0.0);
>     m->add(v0);
>     m->add(v1);
>     m->add(v2);
>     m->add(v3);
> 
>     // Add edges counterclockwise
>     edge* e0 = new edge(m, 0, v0, v1);
>     edge* e1 = new edge(m, 1, v1, v2);
>     edge* e2 = new edge(m, 2, v2, v3);
>     edge* e3 = new edge(m, 3, v3, v0);
>     m->add(e0);
>     m->add(e1);
>     m->add(e2);
>     m->add(e3);
>     
>     // Create a list of edges
>     std::list<GEdge*> e { e0, e1, e2, e3 };
>     
>     // Create the only face
>     face* f = new face(m, 0, e);
>     m->add(f);
> #else
>     
>     std::list<vertex::base*> vertices;
>     std::list<edge::base*>   edges;
>     std::list<face::base*>   faces;
>     
>     vertex* bl = new vertex(m, 0,  0.0,  0.0,  0.0, 1.0);
>     vertex* br = new vertex(m, 1, 10.0,  0.0,  0.0, 1.0);
>     vertex* tr = new vertex(m, 2, 10.0, 20.0,  0.0, 1.0);
>     vertex* tl = new vertex(m, 3,  0.0, 20.0,  0.0, 1.0);
>     vertex* cu = new vertex(m, 4,  0.0, 11.1,  0.0, 1.0);
>     vertex* ct = new vertex(m, 5,  5.0, 10.0,  0.0, 1.0);
>     vertex* cd = new vertex(m, 6,  0.0,  8.9,  0.0, 1.0);
>     
>     // All vertices
>     vertices.push_back(bl);
>     vertices.push_back(br);
>     vertices.push_back(tr);
>     vertices.push_back(tl);
>     vertices.push_back(cu);
>     vertices.push_back(ct);
>     vertices.push_back(cd);
>     for (auto &p : vertices) m->add(p);
> 
>     edge* be = new edge(m, 0, bl, br);
>     edge* re = new edge(m, 1, br, tr);
>     edge* te = new edge(m, 2, tr, tl);
>     edge* l1 = new edge(m, 3, tl, cu);
>     edge* eu = new edge(m, 4, cu, ct);
>     edge* ed = new edge(m, 5, ct, cd);
>     edge* l2 = new edge(m, 6, cd, bl);
>     
>     // All edges
>     edges.push_back(be);
>     edges.push_back(re);
>     edges.push_back(te);
>     edges.push_back(l1);
>     edges.push_back(eu);
>     edges.push_back(ed);
>     edges.push_back(l2);
>     for (auto &p : edges) m->add(p);
> 
>     face* f = new face(m, 0, edges);
>     
>     // Add the only face
>     faces.push_back(f);
>     m->add(f);
> #endif
>     
>     t = std::time(nullptr);
>     tm = *std::localtime(&t);
>     std::cout << "meshing @ " << std::put_time(&tm, "%d-%m-%Y %H-%M-%S") << std::endl;
>     
>     // Dimension is 2
>     m->mesh(2);
> 
>     t = std::time(nullptr);
>     tm = *std::localtime(&t);
>     std::cout << "vtk @ " << std::put_time(&tm, "%d-%m-%Y %H-%M-%S") << std::endl;
> 
>     // Write a VTK file to test
>     m->writeVTK("/Users/sensei/Desktop/a.vtk");
>     
>     // Deallocate or else gmsh fails miserably
>     delete m;
>     
>     GmshFinalize();
>     
>     t = std::time(nullptr);
>     tm = *std::localtime(&t);
>     std::cout << "end @ " << std::put_time(&tm, "%d-%m-%Y %H-%M-%S") << std::endl;
> 
> 
> In release mode, it takes about five minutes! 
> 
> start @ 23-04-2015 15-44-57
> meshing @ 23-04-2015 15-44-57
> vtk @ 23-04-2015 15-49-49
> end @ 23-04-2015 15-49-49
> 
> I must me doing something VERY wrong with my code, so I am resorting to you. Being a newbie with gmsh, I expect I am doing something very stupid! 
> 
> Moreover, when I open the output vtk in paraview, it seems it has NO FACES at all. 
> 
> At the end of the mail you will find my custom classes for vertices, edges, and faces.
> 
> Thanks for any hints you can give me!
> 
> Cheers & Thanks!
> 	Franco
> 	/fm
> -- 
> Franco Milicchio <fmilicchio at me.com>
> 
> DIA - Dept. of Engineering
> University Roma Tre
> http://plm.dia.uniroma3.it/milicchio/
> 
> 
> #include "GVertex.h"
> #include "GEdge.h"
> #include "GFace.h"
> 
> // Forward class definition
> class GModel;
> 
> // Vertex class
> class vertex : public GVertex
> {
>     GPoint p_;
>     
>     
> public:
>     
>     // Base class (useful)
>     typedef GVertex base;
>     
>     // Constructor
>     vertex(GModel *m, int idx, double x, double y, double z, double length = 10e-1) : GVertex(m, idx, length), p_(x, y, z)
>     {
>         // NOP
>     }
>     
>     // Destructor
>     virtual ~vertex()
>     {
>         // NOP
>     }
> 
>     
>     // Return the vertex coordinates
>     virtual GPoint point() const
>     {
>         return p_;
>     }
> 
>     // Return the vertex X coordinate
>     virtual double x() const
>     {
>         return p_.x();
>     }
> 
>     // Return the vertex Y coordinate
>     virtual double y() const
>     {
>         return p_.y();
>     }
> 
>     // Return the vertex Z coordinate
>     virtual double z() const
>     {
>         return p_.z();
>     }
> 
>     // Set the position of the vertex
>     virtual void setPosition(GPoint &p)
>     {
>         p_ = p;
>     }
> };
> 
> 
> // Vertex class
> class edge : public GEdge
> {
> public:
>     
>     // Base class (useful)
>     typedef GEdge base;
>     
>     // Constructor
>     edge(GModel *m, int idx, vertex::base *v0, vertex::base *v1) : GEdge(m, idx, v0, v1)
>     {
>         // NOP
>     }
>     
>     // Destructor
>     virtual ~edge()
>     {
>         // NOP
>     }
> 
>     
>     // Parametric bounds (?)
>     virtual Range<double> parBounds(int i) const
>     {
>         return Range<double>(0.0, 1.0);
>     }
>     
>     // Return the parametric point, with 0 <= p <= 1
>     virtual GPoint point(double p) const
>     {
>         double x = (1 - p) * v0->x() + p * v1->x();
>         double y = (1 - p) * v0->y() + p * v1->y();
>         double z = (1 - p) * v0->z() + p * v1->z();
>         
>         return GPoint(x, y, z, this, p);
> 
>     }
>     
>     // Return the first derivative on the edge (parametrized with 0 <= p <= 1)
>     virtual SVector3 firstDer(double p) const
>     {
>         double x = - v0->x() + v1->x();
>         double y = - v0->y() + v1->y();
>         double z = - v0->z() + v1->z();
>         
>         return SVector3(x, y, z);
>     }
> };
> 
> // Face class
> class face : public GFace
> {
> public:
>     
>     // Base class (useful)
>     typedef GFace base;
> 
>     // Constructor
>     face(GModel *m, int idx, const std::list<edge::base*> &edges) : GFace(m, idx)
>     {
>         edgeLoops.push_back(GEdgeLoop(edges));
>                 
>         for (auto it = edges.begin(); it != edges.end(); ++it)
>         {
>             GEdge *e = *it;
>             
>             l_edges.push_back(e);
>             e->addFace(this);
>             l_dirs.push_back(1);
>         }
>         
>         computeMeanPlane();
> 
>     }
> 
>     // Destructor
>     virtual ~face()
>     {
>         // NOP
>     }
>     
>     // Return the parametric bounds (?)
>     Range<double> parBounds(int i) const
>     {
>         return Range<double>(0.0, 1.0);
>     }
>     
>     // Return the parametric point, with 0 <= p_i <= 1
>     virtual GPoint point(double p_0, double p_1) const
>     {
>         double pp[2] = { p_0, p_1 };
>         
>         double x, y, z, VX[3], VY[3];
>         
>         getMeanPlaneData(VX, VY, x, y, z);
>         
>         return GPoint(x + VX[0] * p_0 + VY[0] * p_1,
>                       y + VX[1] * p_0 + VY[1] * p_1,
>                       z + VX[2] * p_0 + VY[2] * p_1, this, pp);
> 
>     }
>     
>     // Return the first derivative on the edge (parametrized with 0 <= p_i <= 1)
>     virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &p) const
>     {
>         double x, y, z, VX[3], VY[3];
>         
>         getMeanPlaneData(VX, VY, x, y, z);
>         
>         return Pair<SVector3, SVector3>(SVector3(VX[0], VX[1], VX[2]),
>                                         SVector3(VY[0], VY[1], VY[2]));
>     }
>     
>     // Return the second derivative on the face (parametrized with 0 <= p_i <= 1)
>     virtual void secondDer(const SPoint2 &param, SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
>     {
>         // No derivative: this is a PLANAR FACE, no curvature
>         *dudu = SVector3(0.0, 0.0, 0.0);
>         *dvdv = SVector3(0.0, 0.0, 0.0);
>         *dudv = SVector3(0.0, 0.0, 0.0);
>     }
>     
>     // Return the geometric type of the face, i.e., a plane face
>     virtual GeomType geomType() const
>     {
>         return Plane;
>     }
> };
> 
> 
> 
> _______________________________________________
> gmsh mailing list
> gmsh at geuz.org
> http://www.geuz.org/mailman/listinfo/gmsh

-- 
Prof. Christophe Geuzaine
University of Liege, Electrical Engineering and Computer Science 
http://www.montefiore.ulg.ac.be/~geuzaine