visualization of 2d cross-sections from 3d mesh solution

Question

Hey, I have the following code:

import from dolfin import *

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

# Define boundary conditions
u0 = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]')

def u0_boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u0, u0_boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = inner(nabla_grad(u), nabla_grad(v))*dx
L = f*v*dx

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

# Plot solution and mesh
plot(u)
plot(mesh)
interactive()

Now, I have to visualize the solution to this problem on a cross section of the cube, for example on plane " x - y = 0" or "x - z =0".
How can I do that..?

Answer 1

Hi, if you don't insist on visualizing in FEniCS then perhaps the easiest way is to use Paraview. Dump the solution in your code by

File("u.pvd") << u

Then open the file and use the slice filter. The result looks like this.

Comments

Thanks a lot MiroK

Answer 2

The fenicstools Python package was created to do just this. You create the slice in FEniCS and dump the slice to a vtk file that can be visualised in, e.g., ParaView. This is cheaper and faster than saving the whole solution, which can be very memory demanding. I wrote fenicstools since my laptop didn't have enough memory to view a single snapshot of the CSF flow you see on the wiki.

You could put this at the end of your program. You can visualise the solution directly in Python if you have scitools installed. Otherwise dump to vtk.

from fenicstools import StructuredGrid
origin = [0., 0., 0.5]           # origin of slice
vectors = [[1, 0, 0], [0, 1, 0]] # directional tangent directions (scaled in StructuredGrid)
dL = [1., 1.]                    # extent of slice in both directions
N  = [50, 50]                    # number of points in each direction

sl = StructuredGrid(V, N, origin, vectors, dL)
sl(u)
sl.surf(0)  # Visualize slice using gnuplot (if scitools is installed)
sl.tovtk(0, filename='slice.vtk')