Gmsh, Meshio, Dolfin

mesh
1

I had a couple of questions regarding the whole meshing workflow when using gmsh and meshio to get a mesh suitable for FEniCS.

1. What kind of mesh format is suitable for FEniCS? Is it XML or XDMF? What is the difference between the two? What does the H5 file have to do with the XDMF format?

2. Are there any severe differences when using the gmsh API vs using the pygmsh module to create the geometry and the corresponding mesh?

3. I have used the following code to convert a 2D .msh file to .xdmf format and have a couple of remarks on what I think is happening and please correct me if I am wrong :

a) The cell_data_dict() function returns the cell data of the values under the key ‘triangle’ which is again value under the key ‘gmsh:physical’. So is it basically a nested dictionary? I am a beginner in python so please correct me if my understanding is wrong.

b) Inside meshio.write, in cell_data={“surf_marker” : [tri_data]}, is it basically just writing a dictionary with the key as “surf_marker”? Where is this key useful to us later on? Is it during the post processing of the mesh? Or as I read it during the input of mesh using Dolfin? as seen in question 4)

c) I have seen in a code that @Dokken provided, uses a prune z command to remove the z cordinates before using the points for a 2d case? What was the need and cause of that?

import meshio

mesh_name = "<filename>"
msh = meshio.read(mesh_name+".msh")

tri_data = msh.cell_data_dict["gmsh:physical"]["triangle"]
meshio.write(mesh_name+"_surf.xdmf",
    meshio.Mesh(points=msh.points,
        cells={"triangle": msh.cells_dict["triangle"]},
        cell_data={"surf_marker": [tri_data]}
    )
)

line_data = msh.cell_data_dict["gmsh:physical"]["line"]
meshio.write(mesh_name+"_line.xdmf",
    meshio.Mesh(points=msh.points,
        cells={"line": msh.cells_dict["line"]},
        cell_data={"line_marker": [line_data]}
    )
)

4. My final question, if I have obtained the .xdmf files for both the surface and the lines, how should be inputting them into a fenics program? Please find my MWE based on an answer by Dokken in another question I came across…

from dolfin import *

mesh = Mesh()
with XDMFFile("<surfacemesh>.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 1) 
with XDMFFile("<linemesh>.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

V = FunctionSpace(mesh, CG, 1)

u_D = constant(0.0)
bc_left = DiricheltBC(V, u_D, mf, 1) 
bc_right = DirichletBC(V, u0, mf, 3)

bcs = [bc_left, bc_right]

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx

# Constant forcing function and no Nuemann conditions
f = Constant(-6.0)
L = f*v*dx

# Compute solution
uh = Function(V)

solve(a == L, uh, bcs, solver_parameters={'linear_solver': 'gmres','preconditioner': 'ilu'})

# Save solution in VTK format for Paraview
file = File("poisson.pvd")
file << uh

where the “name_to_read” would be the “surf_marker” or “line_marker” depending on a 2d or 3d mesh.?
Can someone explain what MeshValueCollection and MeshFunctionSize do here?
Also after this if I wanted to set boundary conditions for the left and right boundary would I do that just using as shown in the last line? (Its physical_group tag is 1 and 3 as per the script given below)

Please consider also my script that generates a unit square mesh:

import gmsh
import sys 

gmsh.initialize()

lc = 0.25
gmsh.model.geo.add_point(-0.5,-0.5,0,lc,1)
gmsh.model.geo.add_point(-0.5,0.5,0,lc,2)
gmsh.model.geo.add_point(0.5,0.5,0,lc,3)
gmsh.model.geo.add_point(0.5,-0.5,0,lc,4)

gmsh.model.geo.add_line(1,2,1)
gmsh.model.geo.add_line(2,3,2)
gmsh.model.geo.add_line(3,4,3)
gmsh.model.geo.add_line(4,1,4)

gmsh.model.geo.add_curve_loop([1,2,3,4],1)

gmsh.model.geo.add_plane_surface([1],1)

gmsh.model.geo.synchronize()


gmsh.model.add_physical_group(1,[1],1,name="left_boundary")
gmsh.model.add_physical_group(1,[2],2,name="bottom_boundary")
gmsh.model.add_physical_group(1,[3],3,name="right_boundary")
gmsh.model.add_physical_group(1,[4],4,name="top_boundary")
gmsh.model.add_physical_group(2,[1],1,name="domain")

gmsh.model.mesh.generate(2)


gmsh.write("unitsquare.msh")

if '-nopopup' not in sys.argv:
    gmsh.fltk.run()


gmsh.finalize()

Thanks in advance!!

Wanted to add that this workflow is currently working and I was able to produce the following result for the 2d Poisson problem: :smile:

Screenshot from 2023-06-03 03-05-51
Screenshot from 2023-06-03 03-05-51823×582 157 KB

1 Like
2

XML is plaintext, Which makes it hard to read in parallel.
To quote the xdmf site: XDMF Model and Format - XdmfWeb

XDMF categorizes data by two main attributes; size and function. Data can be Light (typically less than about a thousand values) or Heavy (megabytes, terabytes, etc.). In addition to raw values, data can refer to Format (rank and dimensions of an array) or Model (how that data is to be used. i.e. XYZ coordinates vs. Vector components).
XDMF uses XML to store Light data and to describe the data Model. Either HDF5[3] or binary files can be used to store Heavy data. The data Format is stored redundantly in both XML and HDF5. This allows tools to parse XML to determine the resources that will be required to access the Heavy data. For the binary Heavy data option, the xml must list a filename where the binary data is stored.

This means that xdmf is suitable for large files and parallel computing.

pygmsh is a convenience wrapper around the Python API of GMSH, and offers a subset of the commands in the Gmsh Python api. Historically pygmsh has had a “neater” interface than the GMSH Python api (thus its popularity). Personally, I use the GMSH Python API, as the user interface is stable, and allows me to do anything that is part of GMSH.

That said, you can combine the two (as pygmsh calls gmsh), But it might be confusing.

Yes

This writes information abou physical mesh groups that has received a tag/marker in Gmsh to a format that can be read by Dolfin/Paraview.

Meshes generated with Gmsh always has 3D coordinates. For Finite element applications, it is important to know if your mesh is 2D, or a 2D manifold in a 3D space (as gradients, normals and general vectors differ).

Yes

See Difference between meshvaluecollection and meshfunction - #2 by dokken

As Im not at a computer this script is a bit to complex to read for me (as you are not using the Gmsh API with return values, but hardcoded tags, so I cant comment on this at the moment

3

Thanks for the detailed reply!

I had a couple more questions…I used the same mesh converter and read the mesh for solving the NS equations for flow past a cylinder and I got the following error:
Illegal value dimension (2), expecting (3). FEniCS
I think it was due to the absence of pruning function as mentioned earlier, since it worked when I used an example done by you in the forum which pruned the z-direction coordinates. So is the difference coming from solving for the scalar function space vs vector function space?

4

Well, Yes and no.

If you dont prune your third dimension of the mesh, everything VectorFunctionSpaces, grad (gradients), divergence (div), and any spatial operator is assumed to be 3D. Thus if you create a vector space on a manifold, it had three components.

5

I haven’t completely understood your answer, but that is probably because I need to brush up on some math jargon better…also I don’t really understand what a manifold is…
Could you link any articles that explain this better?
Sorry for the trouble…
Thanks again!

6

See for instance this paper by @meg et al.

1 Like