[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 ¶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;
}
};
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.geuz.org/pipermail/gmsh/attachments/20150423/2dd513ae/attachment.html>