Building 2D function g(x,y) from a 3D function f(x,y; z=constant)

Question

Hi,

If I have a function on a 3D space u=Function(V), can I use it to build a function on a 2D space that would represent a slice of the 3D function? For example, I can easily build a 2D grid by doing something like

for x in range(nx):
    for y in range(ny):
        print h*x, h*y, u(h*x,h*y,5.)

but it would be much handier if I had a dolfin Function to manipulate. The fact that some interpolation occurs is not a problem.

Thanks

[edit] - In addition to the answer given below, the following solution works:

from dolfin import *
parameters["reorder_dofs_serial"] = False


#3D FE space
mesh = UnitCubeMesh(20,20,20)
V = FunctionSpace(mesh,"Lagrange",1)


#solve 3D poisson equation on cube
def u0_boundary(x,on_boundary):
    return on_boundary

bc = DirichletBC(V,Constant(0.0) , u0_boundary)

u = TestFunction(V)
v = TrialFunction(V)
a = inner(grad(u),grad(v))*dx
L = Constant("1.0")*v*dx

u = Function(V)
solve(a==L,u,bc)

#Assign slice at z=.5 to 2D function

#2D FE space w/ function
mesh2 = UnitSquareMesh(20,20)
W = FunctionSpace(mesh2,"Lagrange",1)
w = Function(W)


coor = mesh2.coordinates()
w_array = w.vector().array()
if mesh2.num_vertices() == len(w_array):
    for i in range(mesh2.num_vertices()):
        w_array[i] = u(coor[i][0],coor[i][1],.5)

w.vector().set_local(w_array)

plot(w,interactive=True)

Answer 1

The fenicstools package was originally developed with this exact purpose. You can use it to slice your 3D mesh as follows

from dolfin import *
from fenicstools import StructuredGrid, interpolate_nonmatching_mesh

mesh = UnitCubeMesh(30, 30, 30)
V = FunctionSpace(mesh, "CG", 1)

# Create some Function on 3D mesh to be sliced
u = interpolate(Expression("x[0]*x[1]"), V)

# Create a structured grid in x-y plane at z=0.5
dims = (20, 20)
directions = [[1, 0, 0], [0, 1, 0]]
lengths = (1.0, 1.0)
origin = (0, 0, 0.5)
sl = StructuredGrid(V, dims, origin, directions, lengths)

# interpolate u on slice
sl(u)

# Dump result to vtk
sl.tovtk(0, "slice.vtk")

I also have an alternative approach that works without fenicstools. It works by creating a dolfin mesh of topology 2 in 3 dimensions. It goes like this

# Create boundary mesh consisting of the 6 sides of the cube
bmesh = BoundaryMesh(mesh, "exterior")

# Create SubMesh for side at z=0
# This will be a UnitSquareMesh with topology dimension 2 in 3 space dimensions
cc = CellFunction('size_t', bmesh, 0)
xyplane = AutoSubDomain(lambda x: x[2] < 1e-8)
xyplane.mark(cc, 1)
submesh = SubMesh(bmesh, cc, 1)

# Move slice/submesh to z=0.5
x = submesh.coordinates()
x[:, 2] += 0.5

# Create a FunctionSpace on the submesh
Vs = FunctionSpace(submesh, "CG", 1)

# interpolate_nonmatching_mesh required in parallel, 
# interpolate works in series
#us = interpolate_nonmatching_mesh(u, Vs)
us = interpolate(u, Vs)

plot(us)
interactive()