1
$\begingroup$

I have an irregular grid of points describing this surface (a large subduction fault in South America). The color is depth. Anyway I have 3D coordinates (lon,lat,depth) at irregular intervals. I'm trying to generate a triangular mesh using gmsh but I'm struggling with how. I can make every irregular grid point be a "Point" in gmsh with a little python function:

def xyz2gmsh(fout,x,y,z):
    f=open(fout,'w')
    for k in range(len(x)):
        line='Point('+str(k+1)+') = {%.6f, %.6f, %.6f, 0.01};\n' %(x[k],y[k],z[k])
        f.write(line)
    f.close()

And this file (fout) gnerated loads into gmsh but nothing shows! As I understand it I need to somehow tell gmsh that these points collectively represent a surface to be meshed. How? And can I tell gmsh to shoot for making elements of a certain size?

Thanks!

Improve this question
$\endgroup$
5
  • $\begingroup$ This is a question best asked on the gmsh mailing lists. $\endgroup$
    – Wolfgang Bangerth
    Jan 17 '15 at 20:30
  • $\begingroup$ How complex is the surface you're looking at? Does it fold under itself at any point? If you can (temporarily) toss out the z-values, generate a 2d triangular mesh based on the $x$, $y$-values of the sample points, and then "lift" that up to a surface mesh, it would be easier than generating the surface mesh from the start. $\endgroup$
    – Daniel Shapero
    Jan 17 '15 at 21:00
  • $\begingroup$ As long as the gradients of depth aren't too steep, @DanielShapero's method should work pretty well. You could end up with some pretty distorted triangles if they are large since points that are a separated by a large distance get projected only a short distance away from one another in the (x,y) plane. As a side note, the "0.01" that you've hard-coded into your python function is the suggested element size. $\endgroup$
    – Tyler Olsen
    Jan 17 '15 at 21:07
  • $\begingroup$ The surface is very smooth and close to horizontal, so yeah, ignoring the z is not too bad. So how do I do this in gmsh? I need to specify the perimeter then? The tutorial has something about a Line Loop? I just started with gmsh yesterday... $\endgroup$
    – user3684658
    Jan 17 '15 at 21:12
  • 2
    $\begingroup$ Since you already have the points, gmsh may not be the right tool. Gmsh uses a Boundary Representation (BREP) to represent the geometry, onto which it places points that will become the mesh. Since you already have the mesh points, half of that work is already done. What you need to do is to identify the exterior boundary (if not convex) and just use a constrained Delaunay triangulation tool on the points. There are many such tools ready-to-use in the wild. Here's some bonus reading: en.wikipedia.org/wiki/Chew%27s_second_algorithm $\endgroup$
    – Tyler Olsen
    Jan 17 '15 at 21:20
4
$\begingroup$

I'll expand my comment to an answer. Since your surface is fairly smooth, rather than generating a surface mesh, you can generate a 2D mesh of just the $(x, y)$-points that have been sampled, and then create a surface mesh by adding in the $z$-values later.

This might suffice or it might not. Triangulation algorithms, like the one Tyler Olsen linked, are optimized for certain criteria (maximize the minimum angle). When you lift the 2D triangulation to a surface mesh, the outcome may be less than ideal depending on the slopes of the surface. Chapter 8 in Hjelle's book on triangulations covers scattered data interpolation and may be of some use to you.

You've mentioned gmsh, but I actually prefer using the program Triangle for most meshing tasks. Its input/output file formats are much simpler, so you can quite easily write scripts to either create or parse them. I have code for this in C, Fortran and Python if you'd rather be spared the trouble.

In gmsh, you'd have to specify a line loop that parameterizes the convex hull of your input points. With Triangle, you can just give it a point cloud; it will compute the convex hull for you and then triangulate the interior. You only have to describe the boundary if the domain isn't convex. To get triangles of given area, you can either give it a command line switch, or you can write a special .area file if the desired sizes are depend on where you are. I have example code for doing just that in Python for a case where I needed a mesh that was finer in some areas than others.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Does Triangle do constrained triangulations? His domain is not convex, so the Delaunay triangulation may not preserve the original boundary, which may or may not be important to him. $\endgroup$
    – Tyler Olsen
    Jan 17 '15 at 22:04
  • $\begingroup$ @TylerOlsen, yes Triangle does constrained triangulations. In that case, you have to provide an input representing the boundary, but that's no more difficult than in gmsh and all the other output file formats are so much easier to parse. I realize it's kind of obnoxious to answer the question "How to use gmsh?" with "Don't use gmsh," but I just find it to be overkill a lot of the time. $\endgroup$
    – Daniel Shapero
    Jan 17 '15 at 22:13
  • $\begingroup$ This is helpful. I'll give it a go. And the boundary doesn't have to be too strictly obeyed. I don't care if the solution is not gmsh btw. I just want any solution. Cheers. $\endgroup$
    – user3684658
    Jan 17 '15 at 22:18
  • 2
    $\begingroup$ Andreas Kloeckner has written a Python interface to Triangle as part of meshpy; you may want to check that out. $\endgroup$
    – Geoff Oxberry
    Jan 18 '15 at 6:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.