Mesh generation

Dear all,

I would like to know if there exists a generic mesh generator (2D/3D) that allows for:

  1. in 2D: the use either triangles or quadrangles
  2. controlling order of elements
  3. including regional attributes to elements
  4. returning node to element and/or element to element connectivities

So far, I have experience with Triangulate which address point 2 and 3. Both Gmsh and Meshes look very complete. Anyone know whether they can address the 4 points I mentioned above?

thanks for the hints!

Hey,

I have a bit of experience with gmsh and you can do all of the listed things with it. I don’t know about Meshes.jl but it’s probably able to do all of it as well.

Hi, thanks for the reply. I have tried running gmsh. I’ve managed to run the first tutorial by commenting gmsh.fltk.run() which could not run on my machine.
Would you happen to know and example which shows how to extract node to element connectivity list? I only see element to node being stored in the output file. I could make a script that creates a node to element list but I guess this should be directly accessible from the meshing procedure (?).

MeshCore can do this for you: GitHub - PetrKryslUCSD/MeshCore.jl: Curation and manipulation of general unstructured meshes for the Finite Element Methods (FEM).

The functionality with respect to fltk and occ is bundled in the new gmsh_jll.jl you can use it by adding gmsh_jll and

import gmsh_jll
include(gmsh_jll.gmsh_api)

you can retrive the element to node list by

    elementtypes, elementtags, nodetags = gmsh.model.mesh.getElements(dim, -1)

Node to element is probably not directly accessible (at least I don’t know the function call for it) but you can search for stuff by typing gmsh. and use the autocomplete function. You probably need to search for gmsh.model.mesh. and there some function like getNode...

edit: Until no package bundles gmsh_jll you should switch to use this binarybundle package directly, since its well maintained and feature rich. That’s why I archived the gmsh.jl repo

thanks @koehlerson for the detailed explanations, I was until now not using the jll, I will give that one a try and add a post if I find a way to directly get vertex to element and/or element to element lists out of gmsh (I’ve quickly searched the documentation and could not find something like this - I may have searched for wrong keywords).

If I understand you well @PetrKryslUCSD, you’d suggest to (1) generate meshes (e.g. gmsh) and (2) convert the mesh data structure to that used by MeshCore.jl and then (3) use the functions exposed by MeshCore.jl (which look indeed useful). Or you’d rather suggest a more “all-in-one” procedure?

I just mean that in order to derive the inverse relation (nodes → elements), one could use MeshCore. And what to do with it has many options: use MeshCore for everything, or just derive that one relation and store it in your own datastructure…

I am assuming that you are interested in a specific type of mesh generation: tesselation from a set of points. Notice that Gmsh.jl and Triangulate.jl are wrappers to well-established libraries written in other languages that do this type of tesselation. On the other hand, Meshes.jl is a more ambitious project, written in Julia, which aims to include tesselation algorithms, but also other types of mesh generation such as triangulation of areas based on their boundary.

Meshes.jl currently provides the following discretization methods:

https://juliageometry.github.io/Meshes.jl/stable/algorithms/discretization.html

You can also take the result of any of the methods above and refine it into triangles or quadrangles:

https://juliageometry.github.io/Meshes.jl/stable/algorithms/refinement.html

We do not provide tesselation algorithms yet. We have plans to try to revive VoronoiDelaunay.jl and absorb the algorithm in an API that is compatible with the rest of the project, but that will take time.

Feel free to reach out if you have questions. We have a Zulip channel where the conversation can continue if you prefer.

1 Like

you can use the gmsh_jll as I wrote it and translate it yourself or the effortless solution would be using FerriteGmsh

using Ferrite
using FerriteGmsh

grid = saved_file_to_grid("meshfile.msh")
top = Ferrite.ExclusiveTopology(grid.cells)
top.vertex_to_cell

However, it encodes only the vertex to cell and not the node to cell, which could be problematic for the higher order stuff your are planning to do, but that depends on the application

FYI a small note. “Triangle” the underlying C package of Triangulate has an option (-n) that generates neighbor element information. It looks like you can pass that option through Triangulate as well.

@Abhilash, thanks, yes, I do use this option. Shewchuk’s triangle is indeed super useful and I really like to use it. I want to try some stuff on quads and also 3D, hence I’m checking alternatives to triangle.

@juliohm, yes, I’d like to use a tesselation from a set of points. I’d certainly be interested to try meshes.jl once it also provides this. I’d be happy to serve as guinea pig. I’ll meet on Zulip once I’ve also figured out how this one works :smiley:

@koehlerson, I’ve almost managed to do what I want with gmsh. It’s great to have the gmsh popup window with a data picker and possibility to rotate the mesh (without latency!). In fact, much better than any other alternatives I know so far (like writing vtk files and reading them in paraview…), a very good surprise. However, one more question, would you know why the data returned by

element_types, element_ids, element_connectivity = gmsh.model.mesh.getElements(2,1)

is returned in such an unusual (to me) format:

julia> julia> typeof(element_connectivity)
Vector{Vector{UInt64}} (alias for Array{Array{UInt64, 1}, 1})
julia> size(element_connectivity[1])
(192,)

while gmesh indicates ‘Info : 135 nodes 222 elements’

julia> element_connectivity[1][1]
0x000000000000004f

which is does not really read like the connectivity list of element 1 for a human like me.

Every once in a while I perform an FEA analysis. Lately I have been using Salome-Meca as an open source FEA program. It uses Salome as a pre-processor (geometry and meshing), Code-Aster for performing the FEA and Paraview as a post-processor. This may also be worth looking at if you want another option. If you do go in that direction I find the tutorials by Cyprien Rusu to be very helpful.

I don’t know if it meets your criteria, but it will output the nodes and mesh in universal file format.

it only looks like that because of the printing mehtods of UInt64.

julia> a = UInt64(1)
0x0000000000000001

julia> b = [[a]]
1-element Vector{Vector{UInt64}}:
 [0x0000000000000001]

julia> convert(Vector{Vector{Int}},b)
1-element Vector{Vector{Int64}}:
 [1]

Essentially it is just encoding in the number type that your nodeid (they call it node tag in gmsh I think) never can be smaller than 0 if you start to label your nodes at 0. But it’s no big deal, if it distracts you in a REPL based workflow, just convert it to a Vector{Vector{Int}}

thanks @koehlerson, I would have not realised that the data needed a cast to become readable at print (weird, I thougth that the difference between Uint64 and Int64 were about the sign).

What still gives problems is that resulting arrays are not of size of number of elements nor of size of number of nodes (in my previous experiences with triangle, element to node arrays are size number of elements * number of node per element).
Here, if I want to recover the element to node list based, I should do something like this (picked from another discourse post):

# Extract elements
element_types, element_ids, element_connectivity = gmsh.model.mesh.getElements(2,1)
nel = length(element_ids[1])  # number of elements reads from length of number element ids
e2n = zeros(Int16,nel,3)      # element to node numbering

for el in 1:nel
    node1_id = element_connectivity[1][3*(el-1)+1]
    node2_id = element_connectivity[1][3*(el-1)+2]
    node3_id = element_connectivity[1][3*(el-1)+3]
    e2n[el,1]= Int16(node1_id)     
    e2n[el,2]= Int16(node2_id)
    e2n[el,3]= Int16(node3_id)
end

but what if, a mesh contains both quads and triangles? How could one recover the element to node list?

Below I’ve pasted a MWE that attempts to describe the problem. I try to construct a mesh where the west part contains triangles and the right contains quadrangles. Then, I would like to compute element to node list (node to element list would be next on the to do list :smiley: ). Here it’s clear that the length of element_ids[1] does not correspond to the total number of elements (maybe one side?) and I have no idea how to recover the element to node list without known the total number of elment and their type (tri/quad). If anyone would have a hint it would be very welcome! cheers

import gmsh_jll
include(gmsh_jll.gmsh_api)
# Start
gmsh.initialize()
gmsh.model.add("MWE_tri_quads")
# outer box
gmsh.model.geo.addPoint(0.0, 0.0, 0.0) # SW
gmsh.model.geo.addPoint(1.5, 0.0, 0.0) # SE
gmsh.model.geo.addPoint(1.5, 1.0, 0.0) # NE
gmsh.model.geo.addPoint(0.0, 1.0, 0.0) # NW
# Inner limit
gmsh.model.geo.addPoint(0.75, 0.0, 0.0) # S
gmsh.model.geo.addPoint(0.75, 1.0, 0.0) # N
# West 
gmsh.model.geo.addLine(4, 1)
gmsh.model.geo.addLine(1, 5)
gmsh.model.geo.addLine(6, 4)
gmsh.model.geo.addLine(5, 6)
gmsh.model.geo.addCurveLoop([3, 1, 2, 4], 1)
pl1 = gmsh.model.geo.addPlaneSurface([1])
# East
gmsh.model.geo.addLine(5, 2)
gmsh.model.geo.addLine(2, 3)
gmsh.model.geo.addLine(3, 6)
gmsh.model.geo.addCurveLoop([-4, 5, 6, 7], 2) 
pl2 = gmsh.model.geo.addPlaneSurface([2])
# Synchro
gmsh.model.geo.synchronize()
# Make quads in the western part of the domain
gmsh.model.mesh.setRecombine(2, pl2)
gmsh.model.mesh.generate(2)
gmsh.write("MWE_tri_quads.msh")
if !("-nopopup" in ARGS)
    gmsh.fltk.run()
end
########## Try to get things out 
element_types, element_ids, element_connectivity = gmsh.model.mesh.getElements(2,1)
nel = length(element_ids[1])
println(nel)
# for el in 1:nel
#     node1_id = element_connectivity[1][3*(el-1)+1]
#     node2_id = element_connectivity[1][3*(el-1)+2]
#     node3_id = element_connectivity[1][3*(el-1)+3]
#     e2v[el,1]= Int16(node1_id)
#     e2v[el,2]= Int16(node2_id)
#     e2v[el,3]= Int16(node3_id)
# end