I may end up doing a large project using FEniCS which is a collection of tools for automated, efficient solution of differential equations using the finite element method. While playing with implementing a simple solver for fluid flow in a tank by use of FEniCS to solve the Laplace equation I needed to create my own mesh in code.

FEniCS is quite well documented, but I had to look at the source code for some of the mesh conversion routines to find out how to build a mesh from scratch. So, for posterity, here is a reimplementation of dolfin.RectangleMesh in a much more cumbersome (and flexible) way:

import numpy
import dolfin as df

def create_mesh(length, height, nx, ny, show=False):
"""
Make a square mesh manually

Should give exactly the same results as using the built in dolfin.RectangleMesh() class
"""
# The number of mesh entities
nvert = nx*ny
ncell = 2*(nx-1)*(ny-1)

# Positions of the vertices
xpos = numpy.linspace(0, length, nx)
ypos = numpy.linspace(0, height, ny)

# Create the mesh and open for editing
mesh = df.Mesh()
editor = df.MeshEditor()
editor.open(mesh, 2, 2)

editor.init_vertices(nvert)
i_vert = 0
for x in xpos:
for y in ypos:
i_vert += 1

# Add the cells (triangular elements)
# Loop over the vertices and build two cells for each square
# where the selected vertex is in the lower left corner
editor.init_cells(ncell)
i_cell = 0
for ix in xrange(nx-1):
for iy in xrange(ny-1):
# Upper left triangle in this square
i_vert0 = ix*ny + iy
i_vert1 = ix*ny + iy+1
i_vert2 = (ix+1)*ny + iy + 1
i_cell += 1

# Lower right triangle in this square
i_vert0 = ix*ny + iy
i_vert1 = (ix+1)*ny + iy+1
i_vert2 = (ix+1)*ny + iy
i_cell += 1

# Close the mesh for editing
editor.close()
print 'Created mesh with %d vertices and %d cells' % (nvert, ncell)

if show:
df.plot(mesh)
df.interactive()

return mesh


If anyone finds this by searching the web for how to build a mesh programmatically in FEniCS/dolfin, then I hope the above was understandable. The documentation is quite good when you know that you need to search for the MeshEditor class ...