I have 1D finite-volume code written in python for a cell-centred mesh,

First generate a sequence which is the location of the faces, $
\{ x_{j-1/2} \}$. For example, for uniform spacing over the domain [0,1] this is a simple as,
a = 0
b = 1
J = 50
faces = numpy.linspace(a, b, J}
After you have the faces the rest is easy because the faces uniquely determine the cell centres,
$$
x_j = \frac{1}{2}\left(x_{j-1/2} + x_{j+1/2} \right)
$$
and the cell volumes (only a true volume in 3D but we stick with the wording),
$$
h_j = x_{j+1/2} - x_{j-1/2}.
$$
For connivence you can define the centroid spacing, this is useful when interpolating values to/from centres/faces. Following Hundsorfer's approach,
$$
h_{-} = x_j - x_{j-1} = \frac{1}{2}\left( h_{j-1} + h_j \right) \\
h_{+} = x_{j+1} - x_j = \frac{1}{2}\left( h_{j} + h_{j+1} \right)
$$
The following project on github might be useful, https://github.com/danieljfarrell/FVM to see how this can be implemented in practice. I use a Mesh
and CellVariable
objects (inspired by Fipy). You can instantiate the mesh and the query it for centroid position, mesh volume etc. The CellVariable
class using the mesh to do transparent linear interpolation of solution variable.
Here is the mesh class and an example of how to use it. It might be a useful basis for your implementation,
class Mesh(object):
"""A 1D cell centered mesh defined by faces for the finite volume method."""
def __init__(self, faces):
super(Mesh, self).__init__()
# Check for duplicated points
if len(faces) != len(set(faces)):
raise ValueError("The faces array contains duplicated positions. No cell can have zero volume so please update with unique face positions.")
self.faces = np.array(faces)
self.cells = 0.5 * (self.faces[0:-1] + self.faces[1:])
self.J = len(self.cells)
self.cell_widths = (self.faces[1:] - self.faces[0:-1])
def h(self, i):
"""Returns the width of the cell at the specified index."""
return self.cell_widths[i]
def hm(self, i):
"""Distance between centroids in the backwards direction."""
if not check_index_within_bounds(i,1,self.J-1):
raise ValueError("hm index runs out of bounds")
return (self.cells[i] - self.cells[i-1])
def hp(self, i):
"""Distance between centroids in the forward direction."""
if not check_index_within_bounds(i,0,self.J-2):
raise ValueError("hp index runs out of bounds")
return (self.cells[i+1] - self.cells[i])
import numpy
faces = numpy.linspace(0,1,6)
mesh = Mesh(faces)
mesh.cells # i.e. cell centres
# array([ 0.1, 0.3, 0.5, 0.7, 0.9])
mesh.cell_widths # i.e. cell volumes
# array([ 0.2, 0.2, 0.2, 0.2, 0.2])