Problem with evaluation at a point in parallel

2

The mesh is partitioned between processes when running in parallel. You need to ensure that the point at which you wish to evaluate the FE function is contained in the mesh partition on your current process, otherwise extrapolation is required.

However, if you would like every process to see the same point evaluation of a function, you can do something along these lines which I’ve sketched out below

from dolfin import *
import numpy as np

def evaluate_function(u, x):
    comm = u.function_space().mesh().mpi_comm()
    if comm.size == 1:
        return u(*x)

    # Find whether the point lies on the partition of the mesh local
    # to this process, and evaulate u(x)
    cell, distance = mesh.bounding_box_tree().compute_closest_entity(Point(*x))
    u_eval = u(*x) if distance < DOLFIN_EPS else None

    # Gather the results on process 0
    comm = mesh.mpi_comm()
    computed_u = comm.gather(u_eval, root=0)

    # Verify the results on process 0 to ensure we see the same value
    # on a process boundary
    if comm.rank == 0:
        global_u_evals = np.array([y for y in computed_u if y is not None], dtype=np.double)
        assert np.all(np.abs(global_u_evals[0] - global_u_evals) < 1e-9)
    
        computed_u = global_u_evals[0]
    else:
        computed_u = None

    # Broadcast the verified result to all processes
    computed_u = comm.bcast(computed_u, root=0)

    return computed_u


mesh = UnitSquareMesh(32, 32)
f = interpolate(Expression("sin(pi*x[0])*sin(pi*x[1])", degree=4),
                FunctionSpace(mesh, "CG", 1))

x = (0.5, 0.5)
f_x = evaluate_function(f, x)
info("f(%f, %f) = %.3e" % (x[0], x[1], f_x))
3 Likes