[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 ¶m, 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