Problem in evaluating function at point in parallel with dolfinx

dolfinx
1

Hi everyone,
Following the great Dolfinx Tutorials I tried to evaluate a function at one point in parallel. However, I noticed that points which are in the border of the mesh partition may have a different value depending on the process.

Consider this:

import dolfinx
from mpi4py import MPI
import numpy as np

#
comm = MPI.COMM_WORLD
rank = comm.Get_rank()

# define mesh
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
# define function space
V = dolfinx.fem.FunctionSpace(mesh, ("CG", 1))
# function
foo = dolfinx.fem.Function(V)
foo.interpolate(lambda x: rank * np.ones(x[0].shape))
# xdmf file
xdmf = dolfinx.io.XDMFFile(comm, "scratches/foo.xdmf", "w")
xdmf.write_mesh(mesh)
xdmf.write_function(foo)

# select a point in the mesh in the border between the two processes
rand_global_point = [0.5, 0.5,  0.]

# get bbt for the mesh
mesh_bbt = dolfinx.geometry.BoundingBoxTree(mesh, mesh.topology.dim)

# convert point in array with one element
points_list_array = np.array([rand_global_point, ])
# for each point, compute a colliding cells and append to the lists
points_on_proc = []
cells = []
cell_candidates = dolfinx.geometry.compute_collisions(mesh_bbt, points_list_array)  # get candidates
colliding_cells = dolfinx.geometry.compute_colliding_cells(mesh, cell_candidates, points_list_array)  # get actual
for i, point in enumerate(points_list_array):
    if len(colliding_cells.links(i)) > 0:
        cc = colliding_cells.links(i)[0]
        points_on_proc.append(point)
        cells.append(cc)
# convert to numpy array
points_on_proc = np.array(points_on_proc)
cells = np.array(cells)

if len(points_on_proc) > 0:
    print(f"r{rank}: \t * foo.eval(points = {points_on_proc}, cells = {cells}) = "
          f"{foo.eval(points_on_proc, cells)}")

What I get running with mpirun -n 2 is:

r1:      * foo.eval(points = [[0.5 0.5 0. ]], cells = [104]) = [1.]
r0:      * foo.eval(points = [[0.5 0.5 0. ]], cells = [24]) = [0.]

However, looking at my mesh in Paraview the correct value should be 0., regardless the process (the point Iā€™m testing is [0.5, 0.5]):

Imgur
Imgur1491Ɨ830 82.3 KB

Am I doing something wrong?

2

As interpolation is a local function (no parallel communication), it simply computes the interpolation on each process (where the dofs at 0.5, 0.5 is on both processes, and the solution is not unique).
Adding

foo.x.scatter_forward()

after interpolation will send the data from the process owning the DOF to the other process, and you will get your expected result:

r1:      * foo.eval(points = [[0.5 0.5 0. ]], cells = [104]) = [-5.55111512e-17]
r0:      * foo.eval(points = [[0.5 0.5 0. ]], cells = [24]) = [-5.55111512e-17]
4 Likes