Can't identify cell using `find_cell`

[ldouteau]

Hi,

We are currently trying to use Vedo to figure out where a Point is positioned in a mesh. We identified the find_cell method, but it gives weird results. On the piece of code below we build a 1-cell mesh, position a point inside it, then find out in which cell the point is positioned. This code returns -1, as if no corresponding cell were identified.

We played a bit with ordering of the points to create the cell, but this had no effect. FWIW we tried the same piece of code with PyVista, using find_containing_cell, which returned expected results. We didn't look into Vedo's internals yet.

Let me know if you need more details

import vedo

# 1-cell mesh
verts = [(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)]
cells = [(0, 1, 2, 3)]
type = [10 for _ in cells]
mesh = vedo.UnstructuredGrid([verts, cells, type])

# Point to track
point_to_find = (0.2, 0.2, 0.2)
pt = vedo.Point(point_to_find)

# Cell id is returned as -1, which is unexpected
cell_index = mesh.find_cell(point_to_find)
print(cell_index)

# Visualization
mesh.opacity(0.1)
vedo.show(mesh, pt, viewup="z", axes=3)

[ldouteau]

Hi @marcomusy. Sorry for pinging, but are you aware of this one ?

[ldouteau]

I investigated and have more details on this. Weird behavior, the function find_cell works properly when it's not targeting the first cell of the mesh. I extended the minimal example provided above to compare the results of Vedo's find_cell (uses method FindCell from vtk.vtkUnstructuredGrid if i'm correct) with a VTK implementation using FindCell from vtk.vtkCellLocator. The latter correspond to what is used in Pyvista.

Test case

Cube composed of 6 tetra cells. Identification of the cell containing each point from the list. Each point will be located in the cell that corresponds to it's id: P0 -> cell 0, ...

Works properly for cells with id > 0, fails when id is 0.

Results

point [0.1, 0.5, 0.9]
find_cell -1
vtkCellLocator 0
point [0.1, 0.1, 0.5]
find_cell 1
vtkCellLocator 1
point [0.5, 0.1, 0.1]
find_cell 2
vtkCellLocator 2
point [0.5, 0.9, 0.9]
find_cell 3
vtkCellLocator 3
point [0.9, 0.9, 0.5]
find_cell 4
vtkCellLocator 4
point [0.9, 0.5, 0.1]
find_cell 5
vtkCellLocator 5

Code

import vedo
import vtk

verts = [
    (0, 0, 0),
    (1, 0, 0),
    (0, 1, 0),
    (0, 0, 1),
    (0, 1, 1),
    (1, 0, 1),
    (1, 1, 0),
    (1, 1, 1),
]
cells = [
    (2, 3, 4, 5),
    (0, 2, 3, 5),
    (0, 1, 2, 5),
    (2, 4, 5, 7),
    (2, 5, 6, 7),
    (1, 2, 5, 6),
]
type = [10 for _ in cells]
mesh = vedo.UnstructuredGrid([verts, cells, type])

points_to_find = [
    [0.1, 0.5, 0.9],
    [0.1, 0.1, 0.5],
    [0.5, 0.1, 0.1],
    [0.5, 0.9, 0.9],
    [0.9, 0.9, 0.5],
    [0.9, 0.5, 0.1],
]

for point_to_find in points_to_find:
    print(f"point {point_to_find}")
    pt = vedo.Point(point_to_find)
    # Vedo's find cell. Missing the first cell of the mesh.
    cell_index = mesh.find_cell(point_to_find)
    print(cell_index)
    # Find cell through VTK locator. Works properly.
    locator = vtk.vtkCellLocator()
    locator.SetDataSet(mesh.dataset)
    locator.BuildLocator()
    val = locator.FindCell(point_to_find)
    print(val)

# Diplay the example
mesh.opacity(0.5)
vedo.show(
    mesh.cmap(
        "jet", range(mesh.ncells), on="cells", n_colors=mesh.ncells
    ).add_scalarbar(),
    *[vedo.Point(p) for p in points_to_find],
    viewup="z",
    axes=3,
)