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

Franco Milicchio fmilicchio at me.com
Thu Apr 23 15:57:33 CEST 2015


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:

    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;
    }
};



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.geuz.org/pipermail/gmsh/attachments/20150423/2dd513ae/attachment.html>