[BUG]: Wrong collision detected in `dolfinx.geometry.determine_point_ownership`
- Dominant language
- C++
- Stars
- 1.2k
- Forks
- 261
- Avg merge
- 1d 15h
- Merged PRs (30d)
- 65
Description
### Summarize the issue
Hi,
I found a bug using `dolfinx.geometry.determine_point_ownership`, since a collision is detected for a point that lies outside the mesh. I think this is possible, if a point lies outside the mesh, but still inside the bounding box of the cell.
### How to reproduce the bug
I ran the following code below using `dolfinx 0.11.0` on a single MPI-Process. It creates a mesh consisting of a single tetrahedron and check collision of 3 different points using `dolfinx.geometry.determine_point_ownership`. The first point is contained inside the tetrahedron,
the second point lies outside the tetrahedron but inside its bounding box and the third point lies outside this bounding box.
As expected a collision is found for the first point and no collision is found for the third point, but for the second point there's also a collision detected, which is wrong. See output messages below.
I believe this is a bug in the implementation of `dolfinx.geometry.determine_point_ownership` (probably on the C++-side), where the collisions with bounding boxes were checked correctly, but not the actual collision with the cells.
### Minimal Example (Python)
```Python
from mpi4py import MPI
import numpy as np
import ufl
import basix
import dolfinx
import pyvista as pv
comm = MPI.COMM_WORLD
mpi_size = comm.size
# Define mesh consiting of a single tetrahedron
rank = comm.rank
points = np.array([
[0.0, 0.0, 0.0],
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 1.0]
], dtype=np.float64)
cells = np.array([[0, 1, 2, 3]], dtype=np.int64)
coord_element = basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(3,))
ufl_mesh = ufl.Mesh(coord_element)
mesh = dolfinx.mesh.create_mesh(comm, cells=cells, e=ufl_mesh, x=points)
tdim = mesh.topology.dim
# Define points, for which collisions with the mesh will be detected
if comm.rank == 0:
points_ext = np.array([
[0.25, 0.25, 0.25], # Point inside the mesh
[0.5, 0.5, 0.5], # Point outside the mesh, but inside the bounding box
[1.5, 0.5, 0.5], # Point outside the bounding box
], dtype=np.float64)
else:
points_ext = np.empty((0, 3), dtype=np.float64)
# Check, which rank contains the point, if any does
point_owner_ship_data = dolfinx.geometry.determine_point_ownership(mesh, points_ext, padding=0.0)
dest_points = point_owner_ship_data.dest_points
dest_cells = point_owner_ship_data.dest_cells
dest_owner = point_owner_ship_data.dest_owner
src_owner = point_owner_ship_data.src_owner
num_cells = mesh.geometry.dofmaps[0].shape[0]
if num_cells > 0:
# Only consider ranks that owns a cell
print(f"rank = {comm.rank}")
print(f"dest_points = \n{dest_points}")
print(f"dest_cells = {dest_cells}")
print(f"dest_owner = {dest_owner}")
print(f"src_owner = {src_owner}")
# Check, if the point is indeed contained in the mesh
if len(dest_cells) > 0:
cell_candidates = dolfinx.graph.adjacencylist(dest_cells.reshape(-1, 1))
colliding_cells = dolfinx.geometry.compute_colliding_cells(mesh, cell_candidates, dest_points)
# Only the first point should be found
print(f"colliding_cells.array = {colliding_cells.array}")
print(f"colliding_cells.offsets = {colliding_cells.offsets}")
print(flush=True)
# Visualization of the mesh and the colliding points
pv_mesh = pv.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh, tdim))
pl = pv.Plotter()
pl.add_mesh(pv_mesh, show_edges=True, opacity=0.5)
if dest_points.shape[0] > 0:
pl.add_points(dest_points, render_points_as_spheres=True, point_size=100.0, color="blue")
pl.show_axes()
pl.show()
```
### Output (Python)
```bash
rank = 0
dest_points =
[[0.25 0.25 0.25]
[0.5 0.5 0.5 ]]
dest_cells = [0 0]
dest_owner = [0 0]
src_owner = [ 0 0 -1]
colliding_cells.array = [0]
colliding_cells.offsets = [0 1 1]
```
### Version
0.11.0
### DOLFINx git commit
_No response_
### Installation
I installed dolfinx via conda on Linux (openSUSE Leap 15.6)
### Additional information
_No response_
Contributor guide
Research direction
Start by running the minimal Python example with dolfinx.geometry.determine_point_ownership and inspect the subsequent compute_colliding_cells call. Trace these entry points into the C++ geometry implementation; done means the point inside the tetrahedron remains a collision, the point inside only its bounding box does not, and the reported output matches the issue's expectation.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- cpp, python
- Domain
- hpc
- Issue type
- Bug
- Difficulty
- 3/5
- Estimated time
- 1-2 days
- Activity status
- Active
- Clarity
- Mostly clear
- Newbie friendliness
- 68/100