How to evaluate dolfinx Function when some points aren't on the current processor?

dolfinx
1

I have a 2D mesh where I am trying to extract the values at many sets of points on the mesh. When I run the set of points through the code below, I don’t get all points. If I have 1000 points to find, I only get 950 of them when I check the lenght of the “cells” list.
How do I get all of the points I want and not just those listed on the current proccesor? I followed the example from this
https://jsdokken.com/dolfinx-tutorial/chapter1/membrane_code.html#making-curve-plots-throughout-the-domain

cells_collisions = dolfinx.geometry.compute_collisions(tree,points)
colliding_cells = dolfinx.geometry.compute_colliding_cells(msh,cells_collisions,points)
cells = []
for i, point in enumerate(points):
    if len(colliding_cells.links(i))>0: 
        cells.append(colliding_cells.links(i)[0])
2

You should evalutate the points on the process that They are found on, and then use MPI.gather to gather results to a single (or multiple) processes if you want to plot the result.

3

Thanks for replying! I have not intentionally split them into multiple processes so I am unsure how to find which points are on which process. I am also unfamiliar with MPI.gather. Could I run something like:

cells_collisions = dolfinx.geometry.compute_collisions(tree,points)
colliding_cells = dolfinx.geometry.compute_colliding_cells(msh,cells_collisions,points)
gathered_cells = MPI.COMM_WORLD.gather(colliding_cells)
cells = []
for i, point in enumerate(points):
    if len(gathered_cells.links(i))>0: 
        cells.append(gathered_cells.links(i)[0])
4

No.

You need to consider what information each processor has.

If you use MPI.COMM_WORLD as input communicator to the mesh creation, you get a mesh that is partitioned. This means that the cells, and corresponding dofmaps, functions etc also are distributed.

See for instance
https://jsdokken.com/dolfinx_docs/meshes.html
or
Parallel Computations with Dolfinx using MPI — MPI tutorial (You should Also consider the intro pages to mpi in the latter link).

Ive also covered something similar in Gather solutions in parallel in FEniCSX - #2 by dokken

5

Maybe I misunderstood my own problem. Now I am currently running the code with MPI.COMM_SELF to just have it run on one process (if I understood correctly). I check the rank and size and get 0 and 1 respectively so I believe there is not multiple processes running.

When I run

for i, point in enumerate(points):
    if len(colliding_cells.links(i))>0: 
        cells.append(colliding_cells.links(i)[0])

The length of “cells” is less than the length of “points” I input.

If I run

for i, point in enumerate(points): 
        cells.append(colliding_cells.links(i)[0])

I return the error “index 0 is out of bounds for axis 0 with size 0”
This makes sense from the ‘if’ statement above.

If I run

for i, point in enumerate(points): 
        cells.append(colliding_cells.links(i))
uh.eval(points,cells)

I return the error “setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (1001000,) + inhomogeneous part.”

Is there a reason some “colliding_cells.links(i)” have no length? I appreciate your help!

6

You would have to make a minimal reproducible example. If you are running in serial and you don’t get a colliding cell for a point, are you sure that the point in question is inside the mesh?

7

I realized I was making a mistake with what points are in the mesh vs not. The 2D grid of points I was using is within the bounds of the mesh but I had some small spots cut out of the mesh for boundary conditions. I hadn’t thought of this because I was previously using BoundingBoxTree and create_midpoint_tree with compute_closest_entity and hadn’t run into problems with the small cut out areas.