FEniCS / FEniCS/dolfinx

[BUG]: Wrong collision detected in `dolfinx.geometry.determine_point_ownership`

Open
#4,450 3 comments 0 reactions 0 assignees View on GitHub
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

Open the contributing 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

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.