FEniCS: The mesh workflow

FEniCS: The mesh workflow

This is my final post for the GSoC2019 program. The primary goal of the project was to ensure that the meshing package of choice gmsh, DOLFIN, and the preferred visualization package, Paraview work seamlessly together. The intention was to make improvements to the process of preserving the information about tagged regions of the mesh when importing it to DOLFIN. Two approaches were finalized for the project.

  • Preserve the string tags when converting from .msh to .xdmf.
  • Add a constructor to the class MeshValueCollection to support its creation from primitive arrays.

The targets were achieved with the following pull requests.


  • PR #425-Add methods to read and write field_data to XDMF – merged on Aug 9
    • This PR adds methods to preserve the mapping between the string tags and int tags when converting from .msh to .xdmf. The idea is to store this mapping inside the <Information> element of XDMF.


  • PR #439 – Add function to read Information tag data in XDMF – open
    • This PR extends the functionality of XDMFFile interface with methods read_information() and write_information(). The methods are designed to read and write <Information> tag data of XDMF file format. The syntax is inline with XDMF standard and works well with PARAVIEW.
  • PR #467 – Method to construct MeshValueCollection from arrays – open
    • This PR extends the functionality of MeshValueCollection class with constructor that supports its creation from primitive arrays.


We can solve many different forms of PDE on simple as well as complex domains in FEniCS. We have many different prebuilt geometries in FEniCS that helps new users to get up and running with simple FEniCS classes and methods. Even though these built-in meshes provide various methods for their construction and refinement, they are limited to simple shapes. To work with complex geometrical structures it is recommended that the user follows the following mesh workflow.

Create a mesh using software Gmsh

Convert the mesh to XDMF using meshio

Read the XDMF file into FEniCS

Visualize the results in PARAVIEW

In this post, I am going to walk you through the process of solving a PDE in FEniCS. I will demonstrate the whole mesh workflow in detail. The domain under consideration in a unit square as presented below. The domain is made up of two different materials and we need to apply different boundary conditions on all of its edges.

Creating mesh in Gmsh

There are excellent tutorials and documentation available on the official website of Gmsh. In this small screencast, I have described the process of creation and tagging of the above mesh using gmsh.

Convert mesh to XDMF using meshio

The .msh file created by gmsh could be converted to .xdmf by using the package meshio. The package could be easily installed by the following command:

pip install meshio

Once you have the package installed you can use the following command to convert the mesh to xdmf. Right now FEniCS does not support mixed topologies so you have to individually export mesh entities of different dimension to different XDMF files. Thus for the current mesh, we need to export one mesh of 2D `triangle` elements and the other of 1D ‘line’ elements.

mesh_of_triangles = meshio.Mesh(points=points[:, :2],
                                cells={'triangle': cells['triangle']},

meshio.write("poisson_subdomain_triangle.xdmf", mesh_of_triangles )

mesh_of_lines = meshio.Mesh(points=points[:, :2],
                                cells={'line': cells['line']},
                                cell_data={'line': {'subdomain':

meshio.write("poisson_subdomain_line.xdmf", mesh_of_lines )

Import XDMF to FEniCS

Once we have the XDMF files we can import them to FEniCS as follows:

with XDMFFile(MPI.comm_world,
              "poisson_subdomain_triangle.xdmf") as xdmf_infile:
    mesh = xdmf_infile.read_mesh(cpp.mesh.GhostMode.none)
    tag_info = xdmf_infile.read_information_int()

FEniCS provide us with two classes ‘MeshFunction’ and ‘MeshValueCollection’ that help us to mark different mesh entities (vertex, line, facet, cell). The way in which the ‘MeshFunction’ class works is that it assigns a marker to every element of a particular dimension (vertex = 0, line = 1, etc.) and then use that marker to differentiate between different regions. ‘MeshValueCollection’ class differs from the ‘MeshFunction’ class in two ways. First, data does not need to be associated with all entities (only a subset). Second, data is associated with entities through the corresponding cell index and local entity number (relative to the cell), not by global entity index, which means that data may be stored robustly to file. You can get a much better understanding of them by exploring them individually as I have done here.

Visualize the solution

We could visualize the solution using opensource package paraview. Once again the filetype of choice is XDMF and we can easily write the solution to the file using the write method of class XDMFFile.

with XDMFFile(dolfin.MPI.comm_world, "output.xdmf") as xdmf_outfile:

You can find the complete code here.


The mesh workflow of FEniCS is not as straight forward as that of some commercial packages like Abaqus or Ansys. But, FEniCS provides you with greater flexibility in creating a computational model and would also help you understand the mathematics behind FEM. I believe that if understood and used properly, FEniCS could help researchers code and test out their mathematical models with a speed and ease that is just simply not possible with other packages.

Categories: FEniCS

9 Comments on “FEniCS: The mesh workflow

  1. Hi,
    I am trying out your method to create 2D meshes for FEniCS – but I encountered an issue.
    When trying to create the “mesh of triangles” I get the following error:

    TypeError: list indices must be integers or slices, not str

    Do you know a solution to this?
    Thanks for your help.

      1. Hi Abhinav,

        Thank you for this great post on workflow. I am experiencing the same slicing error “TypeError: list indices must be integers or slices, not str” but I was not able to figure out how to resolve this given the page you linked. It seems like we are unable to ever use strings to slice lists (i.e. “cells={‘triangle’: cells[‘triangle’]}” so I am at a loss of how you were able to run these lines.

        Any insight you could provide, or a glimpse into the format of your .msh file, would be very appreciated.

        Thanks for your help.

        1. Hello Tyler,

          When you call the mesh.field_data method, it should return you a dictionary and not a list. This is why you are getting the error message.

          1. meshio is constantly under development. Thus the syntax changes every now and then. I personally use meshiov3.2.7. The above workflow works with it. Use the following command to install
            pip3 install --no-cache-dir meshio==3.2.7

  2. This error is related to the version of meshio that you are used. you should to choose the meshio 3.3.0

  3. I am using meshio 4.0.16 (gmsh 4.4.1, python 3.8.5) and there the cells are structured in cell blocks, so my files read
    meshio.Mesh(points=points, cells=selected_cells, cell_data={selection_data_string: selected_cell_data}
    where selected cells is a list of cellblocks
    selected_cells.append(meshio.CellBlock(cell_type, selected_data))
    cell_type=’line’ or ‘triangle’
    selected_data=array(node_numbers belonging to element).

  4. Hi Abhinav. Your blog is very helpful. I am a novice user of FEniCS and I can use some help from you in troubleshooting a problem I am facing when I am trying to increase the size of system of PDEs from 4 to 5 or more. I used PYGMSH to generate the mesh. Here is a part of the code.

    from dolfin import XDMFFile, Mesh
    import pygmsh
    import meshio
    import numpy as np
    import os
    import datetime as dtm
    import math as mt

    def Generate_PBRpygmsh(xlim, ylim, f_ids, meszf=0.003, folder=’pygmeshio_data’):

    # meszf = mesh element size factor
    # rectangle geometry
    xmin, xmax = xlim
    ymin, ymax = ylim
    WAid, INid = f_ids # Wall and inlet boundaries ids.

    # Pygmsh files
    pygmsh_geo = os.path.join(folder, "PBR.geo")
    pygmsh_msh = os.path.join(folder, "PBR.msh")
    pygmsh_mesh = os.path.join(folder, "PBR_mesh.xdmf")
    pygmsh_facets = os.path.join(folder, "PBR_facets.xdmf")

    geom = pygmsh.built_in.Geometry()
    rect = geom.add_rectangle(xmin, xmax, ymin, ymax, 0.0, lcar=meszf)
    geom.add_physical([rect.line_loop.lines[0]], 10)
    geom.add_physical([rect.line_loop.lines[1]], 11)
    geom.add_physical([rect.line_loop.lines[2]], WAid) # Wall boundary
    geom.add_physical([rect.line_loop.lines[3]], INid) # Inlet boundary
    geom.add_physical([rect.surface], 20)

    mesh = pygmsh.generate_mesh(geom, dim=2, prune_z_0=True,

    cells = np.vstack(np.array([cells.data for cells in mesh.cells
    if cells.type == "triangle"]))
    triangle_mesh = meshio.Mesh(points=mesh.points, cells=[("triangle", cells)])

    facet_cells = np.vstack(np.array([cells.data for cells in mesh.cells
    if cells.type == "line"]))
    facet_data = mesh.cell_data_dict["gmsh:physical"]["line"]
    facet_mesh = meshio.Mesh(points=mesh.points,
    cells=[("line", facet_cells)],
    cell_data={"name_to_read": [facet_data]})

    # Write mesh
    meshio.xdmf.write(pygmsh_mesh, triangle_mesh)
    meshio.xdmf.write(pygmsh_facets, facet_mesh)

    return pygmsh_mesh, pygmsh_facets

    from fenics import *
    from dolfin import XDMFFile, Mesh, FiniteElement, MixedElement, \
    FunctionSpace, Measure, Constant, TestFunctions, \
    Function, split, Expression, project, dot, grad, dx, \
    solve, as_vector, MeshValueCollection, cpp, MPI, \
    TimeSeries, as_backend_type, assemble \

    from ufl.operators import exp

    t = 0.0
    tf = 1.0 # Final time, sec
    num_steps = 100 # Number of time steps
    delta_t = tf / num_steps # Time step size, sec
    dt = Constant(delta_t)

    R_path = ‘Results’ # Path for data storage
    Data_postprocessing = True # Data storage for later visualization & FEniCS computations.
    PBR_L = 5.0 # Reactor length, m
    PBR_R = 0.01 # Reactor radius, m
    meszf = 0.001 # mesh element size factor

    Mesh generation and formats conversion

    Wallboundary_id = 12
    Inletboundary_id = 13
    Boundary_ids = Wallboundary_id, Inletboundary_id
    pygmsh_mesh, pygmsh_facets = Generate_PBRpygmsh((0., PBR_L), (0., PBR_R),
    Boundary_ids, meszf)

    Reading mesh data stored in .xdmf files.

    mesh = Mesh()
    with XDMFFile(pygmsh_mesh) as infile:
    mvc = MeshValueCollection(“size_t”, mesh, 1)
    with XDMFFile(pygmsh_facets) as infile:
    infile.read(mvc, “name_to_read”)

    print(‘num_cells :’, mesh.num_cells())
    print(‘num_vertices:’, mesh.num_vertices())
    print(‘cell_type :’, mesh.ufl_cell())
    #print(‘num_facets:’, mvc.num_facets())
    #print(‘num_edges:’, mvc.num_edges())

    mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

    P1 = FiniteElement(‘P’, mesh.ufl_cell(), 1) # Lagrange 1-order polynomials family
    #element = MixedElement([P1, P1, P1, P1]) #this works fine
    element = MixedElement([P1, P1, P1, P1, P1])
    V = FunctionSpace(mesh, element) # starts throwing error

    Thanks for your help.

    P.S. By the way, I am an alumni from IITR 2003-08 Chemical.

Leave a Reply