Locate_entities returning incorrect points with MPI

dolfinx
,
1

Hi all, I’m using dolfinx 0.8.0 with Python 3.11 and trying to locate all the vertices on a boundary. For a single process this code works ok, but as soon as I use 2 or more I find that locate_entities returns vertices that clearly don’t satisfy the marker function:

import numpy as np
from mpi4py import MPI
from dolfinx.mesh import create_unit_cube, locate_entities

comm = MPI.COMM_WORLD
mesh = create_unit_cube(comm, 5, 5, 5)
idx = locate_entities(mesh, 0, lambda x: np.isclose(x[2], 0))
coords = mesh.geometry.x[idx]

assert np.allclose(coords[:, 2], 0, atol=1e-6), f"rank={comm.rank}: Boundary error: {coords[:, 2]}"

I run this using mpirun -n 2 python error.py and it returns something like:

    assert np.allclose(coords[:, 2], 0, atol=1e-6), f"rank={comm.rank}: Boundary error: {coords[:, 2]}"
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AssertionError: rank=0: Boundary error: [0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
 0.  0.  0.  0.  0.8 1.  0.6 0.8 0.  0.  0.  0.2]
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AssertionError: rank=1: Boundary error: [0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.2 0.4 0.2 0.2]

mpirun -n 1 python error.py runs fine. I assume I am doing something wrong in terms of selecting vertices across different processes?

2

You are mixing the indices of topology (vertices, edges,facets and cells) with geometry (the physical nodes making up the grid).

To map from vertices to geometry nodes (those located in mesh.geometry.x) use
entities_to_geometry
https://docs.fenicsproject.org/dolfinx/v0.8.0/python/generated/dolfinx.mesh.html#dolfinx.mesh.entities_to_geometry
that can map the indices given in idx to indices in mesh.geometry.x

3

I see thanks, that makes sense! For those that come across this later this now works:

import numpy as np
from mpi4py import MPI
from dolfinx.mesh import create_unit_cube, locate_entities, entities_to_geometry

comm = MPI.COMM_WORLD
mesh = create_unit_cube(comm, 5, 5, 5)
entity_idx = locate_entities(mesh, 0, lambda x: np.isclose(x[2], 0))
idx = entities_to_geometry(mesh, 0, entity_idx)
coords = mesh.geometry.x[idx[:, 0]]

assert np.allclose(coords[:, 2], 0, atol=1e-6), f"rank={comm.rank}: coordinates: {coords[:, 2]}"

Also, is there a particular reason that entities_to_geometry returns an (N, 1) shaped array of indices rather than (N,)?

mesh.geometry.x[idx] still works but returns a (N,1,3) shaped array of coordinates which seemed odd.

4

Since entities can be edges, facets or cells, it would return (num entities, num nodes per entity). Since your case is a vertex, there is only one node (1-1 map).
I would call:

idx = entities_to_geometry(mesh, 0, entity_idx).flatten()

Rather than doing it after

1 Like