How can I obtain plane averages?

Question

Hi,
I am a newbie, and I am trying to write a simple poisson solver in 3D with PBC.

I would like to access plane averages of the quantities I calculate, in the sense that I would like to calculate $\frac{1}{A}\int_{z=z_0} u\ ds$ where ds represents an xy plane at $z=z_0$ for a set of $z_0$s

However, the scitools.boxfield method suggested in the tutorial does not work (in the sense that dolfin_function2BoxField is unable to convert the solution to a boxfield).

In short, could anyone please let me know how to access the information contained in U in a structured sense so that I can calculate plane averages?

For the record, the mesh is a BoxMesh, and I am following the pure Neumann example to handle singular matrix (null space example).

Answer 1

Hi

Here is a solution that creates slices as meshes of topology 2 in 3 dimensions and that interpolates from the 3D solution to the slice. You can also use the StructuredGrid routine in the fenicstools package

from dolfin import *

# Create original mesh with solution to be sliced
mesh = UnitCubeMesh(10, 10, 10)
V = FunctionSpace(mesh, "CG", 1)
u0 = interpolate(Expression("x[2]"), V)

# Create a UnitSquareMesh of topology 2 in 3 dimensions
bmesh = BoundaryMesh(mesh, "exterior")
# Choose a side with normal in z-direction
z0 = AutoSubDomain(lambda x, on_bnd: x[2] < DOLFIN_EPS)
ff = CellFunction("size_t", bmesh, 0)
z0.mark(ff, 1)
smesh = SubMesh(bmesh, ff, 1)

# Move the slice to the z-location you want
xyz = smesh.coordinates()
xyz[:, 2] = 0.5

# interpolate from u0 to the slice
Q = FunctionSpace(smesh, 'CG', 1)
lp = LagrangeInterpolator()
u = Function(Q)
lp.interpolate(u, u0)
plot(u)
print "Integral over slice = ", assemble(u*dx)