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.

   Reply

   1. This error is related to slicing in python. You can understand more about it from here.

      Reply

      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.

         Reply

         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.

            Reply

            1. How to deal in that case

            2. 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

   Reply

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))
   and
   cell_type=’line’ or ‘triangle’
   selected_data=array(node_numbers belonging to element).

   Reply

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,
   geo_filename=pygmsh_geo,
   msh_filename=pygmsh_msh)

   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:
   infile.read(mesh)
   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.

   Reply