Tormod Landet

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)

    # Add the vertices (nodes)
    editor.init_vertices(nvert)
    i_vert = 0
    for x in xpos:
        for y in ypos:
            editor.add_vertex(i_vert, x, y)
            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
            editor.add_cell(i_cell, i_vert0, i_vert1, i_vert2)
            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
            editor.add_cell(i_cell, i_vert0, i_vert1, i_vert2)
            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 ...