Fast function evaluation on a grid

Question

I am trying to calculate a 3D Function f on a large set of points, say a uniform cubic grid. Besides, the cubic grid points do not coincide with the mesh vertices where the Function f is defined. For example:

from dolfin import *
mesh1 = Mesh( meshfile) # A uniform mesh
V1 = FunctionSpace(mesh1, 'CG', 1)
f1 = Function(V1)
f1 = ... # Construct f1 (e.g. by solving a PDE based on mesh1)

mesh2 = UnitCube(100,100,100) # Geometrically, the unit cube domain is a subset of mesh1's domain

Now I want to calculate the values of f1 at all vertices of mesh2.

I know that I can either use interpolate function:

V2 = FunctionSpace(mesh2, ‘CG', 1)
f2 = interpolate(f1, V2)

Then f2.vector().array() gives me those values (although I understand that the order is not the same as the mesh vertices for the current version).

Or I can loop all vertices of mesh2 explicitly, and at each step invoke f1(x, y, z) to calculate at that vertex.

The problem is that both these two methods are very very costly, and the interpolate method is even more time consuming than the explicit loop. Does anybody know if there is a more efficient way to do this job?

Thanks a lot,

Answer 1

There will be non-negligible overhead when moving between meshes (some current developments will make this faster than it is now). You could first interpolate, then to evaluate vertex values efficiently use the function Function::compute_vertex_values. For spaces that are not continuous it can yield some unexpected results.