Evaluating Function in Parallel

Errors
1

Dear all,

I’m trying to assign values to a function based on the coordinates using another function and in order to do so, I simply loop over the coordinates. The small function I wrote get a function T_0 as an argument and a simplified code snippet is as such:

    function_space = T_0.function_space()
    mesh = function_space.mesh()
    y = mesh.coordinates()[dof_to_vertex_map(function_space)]
    z_max = np.max(mesh.coordinates()[:,2])
    array = np.zeros_like(T_0.vector()[:])
    
    for i in range(len(array)):
        if near(y[i][2], z_max, 0.000001):
            array[i] = T_old(y[i][0], y[i][1], z_max-0.1)
        else:
            array[i] = T_old(y[i,0], y[i,1], y[i,2])

    T_0.vector().set_local(array)
    as_backend_type(T_0.vector()).vec().ghostUpdate()

    return T_0

where T_0 is the function that I would like to assign values to and T_old is another function that takes three coordinates. This snippet runs correctly in series but I get an error when running in parallel.

*** Error:   Unable to evaluate function at point.
*** Reason:  The point is not inside the domain. Consider calling "Function::set_allow_extrapolation(true)" on this Function to allow extrapolation.

which is totally expected since the function T_old gets points that are outside its domain. Using

    parameters['allow_extrapolation'] = True

gets rid of the error but it is not exactly the solution since then the function would be extrapolating outside its domain. That’s not what I want because I have the actual function available, only on another partition. What would be the right way to achieve this?

The proposed solution in https://fenicsproject.discourse.group/t/problem-with-evaluation-at-a-point-in-parallel/1188 doesn’t work, the script hangs when run in parallel.

It feels like a common problem but I wasn’t able to find a working solution. Any help would be much appreciated.

Thanks in advance and cheers,

rek

2

Have you considered MPI - Gather and distribute a function to all processes (That I can evaluate at a point) - #4 by dokken
I would also suggest making a reproducible example; use a built in mesh, add all required imports and make a code that reproduce the error you show above.
It would also be of help to others if you add your implementation of: Problem with evaluation at a point in parallel - #6 by MiroK
which hangs, such that people can verify that it happens for them as well.