How to extract 1-d slice from 2-d mesh

Question

Suppose I have some function defined on a 2d mesh, for example:

mesh = BoundaryMesh(SphereMesh(Point(0,0,0), 1, 0.1), 'exterior')
V = FunctionSpace(mesh, 'CG', 1)
f = Function(V)
(xs, ys, zs) = V.dofmap().tabulate_all_coordinates(mesh).reshape((-1,3)).T
thetas = np.arctan2(ys, xs)
f.vector()[:] = np.cos(thetas) * (1-zs**2)

How do I extract the function values along a specific path, say with theta = constant? The straightforward approach using e.g.

phi = np.linspace(-math.pi/2, math.pi/2, 100)
values = [ f(x, 0, z) for (x, z) in zip(np.cos(phi), np.sin(phi)) ]

doesn't work (even with parameters['allow_extrapolation'] = True), I think because the spherical mesh does not describe a perfect circle. Is there a way to project my path on the mesh? Or some other way to do what I want?

Answer 1

Here is a workaround to get the desired data manually using ParaView:

1. Save F in a pvd file: File('F.pvd') << F
2. Open F.pvd in paraview
3. Slice the view in the X-Z plane
4. Cut the view along the Y-Z plane
5. Run the Plot on sorted lines filter
6. Export the data as csv
7. Import back into Python using numpy.loadtxt

Obviously that's not feasible if this has to be carried out repeatedly, but it's good enough for temporary debugging.

Answer 2

You may want to check out Mikael Mortensen's fenicstools package: https://github.com/mikaem/fenicstools/wiki

In particular, take a look at the Probes class: https://github.com/mikaem/fenicstools/wiki/Multiple-Probes

For your application, this might do the job:

phi = np.linspace(-math.pi/2, math.pi/2, 100)
probing_pts = np.array([(x, 0, z) for (x, z) in zip(np.cos(phi), np.sin(phi))])
probes = Probes(probing_pts.flatten(), V)
probes(f)  # evaluate f at all probing points
print probes.array()

Note that if you plot the result, it is apparent that it can't compute all the function values and throws a few spurious zeros in there (which is probably what caused the error in your original code).

Alternatively, you could build your own mesh with the vertices you are interested in (e.g. using dolfin's MeshEditor class; see chapter 10.3.2 of the FEniCS book for an example) and then use fenicstools.interpolate_nonmatching_mesh to compute the function values on that mesh. This is more involved, though, and I don't know if it would give better (or worse) results.