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