2D slice from 3D solution

You can try defining an expression that you can interpolate to get rid of your third axis.

Something like this:

#Dummy problem

from fenics import *
import matplotlib.pyplot as plt

# Create mesh and define function space
mesh = UnitCubeMesh(8, 8, 8)
V = FunctionSpace(mesh, 'CG', 1)

# Define boundary condition
u_D = Constant(0)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(100)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Plot solution and mesh
plt.colorbar(plot(u))
plt.title('u')

full_sol

Define an expression to evaluate u at some value of z:

class u_slice(UserExpression):
    def eval(self, value, x):

        value[0] = u(x[0],x[1],0.5)

    def value_shape(self):
        return () 

mesh_2D = UnitSquareMesh(8, 8)
V_2D = FunctionSpace(mesh_2D, 'CG', 1)
u_sl = interpolate(u_slice(), V_2D)
plt.colorbar(plot(u_sl))
plt.title('u in the middle plane')

plane_sol