Pygmsh Tutorial?

mesh
1

Hi all,

Is there a tutorial on how to use pygmsh with FEniCS? For example, I have something like:

import pygmsh
import meshio

-snip-

def generate_box(geom, h1, h2, w):
    lcar = 0.001
    poly = geom.add_polygon(
        [[0.0, 0.0, h1],
         [w, 0.0, h1],
         [w, w, h1],
         [0.0, w, h1]],
        lcar)

    top, volume, lat = geom.extrude(poly.surface, [0, 0, h2-h1])
    bottom = poly.surface

    return [volume,top,bottom]

geom = pygmsh.built_in.Geometry()

-snip-

vol_0,_,_ = generate_box(geom, 0.0, h0, w0)
geom.add_physical(vol_0, label=0)

vol_1,_,_ = generate_box(geom, h0, h1, w0)
geom.add_physical(vol_1, label=1)

-snip-

msh = pygmsh.generate_mesh(geom)

This creates 3D box:es (which I can view in Paraview if I use meshio to write a vtk-file to disk, for example).

How do I make FEniCS understand the mesh with domain/subdomains and how do I set the subdomain properties?

2

See this post Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio for many examples. If you downloading the latest version of meshio, see: Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio (second to last post in this thread)

3

Thanks for the info! I have looked these links before and tried, for example,

-snip-

msh = pygmsh.generate_mesh(geom)


meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]}))
meshio.write("mf.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": msh.cells["triangle"]},
                                    cell_data={"triangle": {"name_to_read": msh.cell_data["triangle"]["gmsh:physical"]}}))

But I only get errors, like

-snip-

Info    : Done optimizing 3D mesh (0.019183 s)
Info    : 657 vertices 4009 elements
Info    : Writing '/tmp/tmp_l_a2ud8.msh'...
Info    : Done writing '/tmp/tmp_l_a2ud8.msh'
Info    : Stopped on Mon Feb 24 10:58:33 2020
Traceback (most recent call last):
  File "run_single_element.py", line 190, in <module>
    meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]}))
TypeError: list indices must be integers or slices, not str

Also, is it possible to convert the “msh” data structure to a FEniCS compatible format directly or do one need to write it to disk and use external tools to convert it?

4

FEniCS does not support direct import from msh, so you have to rewrite the mesh data. As I pointed out in the last post, if you are using the latest version of meshio, you need to use the following syntax:

import meshio
msh = meshio.read("mesh.msh")
for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data
    elif  cell.type == "tetra":
        tetra_cells = cell.data

for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]
    elif key == "tetra":
        tetra_data = msh.cell_data_dict["gmsh:physical"][key]
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells})
triangle_mesh =meshio.Mesh(points=msh.points,
                           cells=[("triangle", triangle_cells)],
                           cell_data={"name_to_read":[triangle_data]})
meshio.write("mesh.xdmf", tetra_mesh)

meshio.write("mf.xdmf", triangle_mesh)

as meshio changed its user interface recently.

2 Likes
5

Thank you! Now I can write the two xdmf-files. I still struggle with reading the data back, though. By looking at the link above one should do something like (?):

mesh_file = XDMFFile(self.comm, "mesh.xdmf")
mesh = Mesh()
mesh_file.read(mesh);

mvc = MeshValueCollection("size_t", mesh, 2)
mf_file = XDMFFile(self.comm, "mf.xdmf")
mf_file.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

to read the data back which then could be used as if I had used mshr, that is,

  -snip-

  dx = Measure('dx', domain=mesh, subdomain_data=mf)

 boundary_parts = FacetFunction("size_t", mesh)
 boundary_parts.set_all(0)
 -snip-

However, I get an error:

    mf_file.read(mvc, "name_to_read")
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to find entity in map.
*** Reason:  Error reading MeshValueCollection.
*** Where:   This error was encountered inside HDF5File.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  86617f912d4035520fa988d6af07a19d6069d5ee
*** -------------------------------------------------------------------------
6

To clarify, I don’t generate or read any msh-files directly, I run the pygmsh command:

msh = pygmsh.generate_mesh(geom)

which returns an “msh” object. I don’t know if that is the same thing as reading an msh-file generated by gmsh (pygmsh generates temporary .geo->.msh files, in /tmp, but the user never see them).

7

I tried to analyze the contents of the msh, tetra_mesh, and, triangle_mesh objects using the python debugger where I got:

(Pdb) p msh       
<meshio mesh object>
  Number of points: 657
  Number of cells:
    triangle: 68
    tetra: 916
    tetra: 204
    tetra: 214
    tetra: 207
    tetra: 213
    tetra: 205
    tetra: 480
  Cell data: gmsh:physical

(Pdb) p triangle_mesh
<meshio mesh object>
  Number of points: 657
  Number of cells:
    triangle: 68
  Cell data: name_to_read

(Pdb) p tetra_mesh
<meshio mesh object>
  Number of points: 657
  Number of cells:
    tetra: 480

and it looks like only the last Box (corresponding to tetra: 480) can be found in tetra_mesh (I also verified this by opening the mesh.xdmf file in Paraview). I’m also not really sure that this is the way to implement subdomains when using pygmsh.

8

That is because you are not making your mesh properly in pygmsh. You need to use boolean fragments in 3D (using opencascade and not built_in). Additionaly, some handling of these unions are required when reading in the mesh.
For generating a simple geometry consisting of a box inside another, with unique markers for each box:

import pygmsh
import meshio

geom = pygmsh.opencascade.Geometry()
box0 =  geom.add_box([0.0, 0, 0], [1, 1, 1], char_length=0.05)
box1 =  geom.add_box([0.25, 0.25, 0.25], [0.4, 0.3, 0.2], char_length=0.05)
frags = geom.boolean_fragments([box0],[box1])

geom.add_raw_code("Physical Volume(1) = {13};")
geom.add_raw_code("Physical Volume(2) = {14};")
msh = pygmsh.generate_mesh(geom, geo_filename = "mesh.msh")
import numpy as np
tetra_cells = None
for cell in msh.cells:
    if cell.type == "tetra":
        # Collect the individual meshes
        if tetra_cells is None:
            tetra_cells = cell.data
        else:
            tetra_cells = np.vstack((tetra_cells, cell.data))

tetra_data = None
for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "tetra":
        tetra_data = msh.cell_data_dict["gmsh:physical"][key]
# Create tetra mesh with cell data attached
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells},
                           cell_data={"name_to_read":[tetra_data]})

meshio.write("mesh.xdmf", tetra_mesh)

To load the mesh and corresponding cell_data into dolfin, you can do the following:


from dolfin import MPI, MeshValueCollection, cpp, Mesh, XDMFFile
mesh_file = XDMFFile(MPI.comm_world, "mesh.xdmf")
mesh = Mesh()
mesh_file.read(mesh);

mvc = MeshValueCollection("size_t", mesh, 3)
mesh_file.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
9

Thanks again for the info! I’m starting to have something that runs now and I will try to sum this up when I understand this better here. Just a few more questions:

  1. Why do I need the geom.boolean_fragments([box0],[box1]) call? Is there documentation somewhere where I can get more info on how this really works?

  2. Why do I need the low level

geom.add_raw_code("Physical Volume(1) = {13};")
geom.add_raw_code("Physical Volume(2) = {14};")

gmsh calls? Doesn’t

geom.add_physical(vol_1, label=1)
geom.add_physical(vol_2, label=2)

do the same thing?

  1. For 2D meshing do I simply replace “tetra” with “triangle” when I loop over the cells and cell_data_dict:s?
10
  1. So, you need the boolean fragments to make a mesh that respects the boundaries of each of the domains that you are trying to make a union of. See the gmsh-tutorials or Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio for more info on this subject.
  2. Since you have to use opencascade to be able to make unions of 3D geometries, a lot of the interfacing in pygmsh does not really work. For 3D geometries, I usually use a combination of the Gmsh GUI and pygmsh to script physical labels.

For 2D meshing you can simply replace “tetra” with “triangle”.