nami
October 19, 2021, 5:51pm
1
Hi,

This might be a very simple question but can’t find an easy solution to it. Suppose I have a point which lies strictly inside of a mesh cell.
1- What is the fastest way to extract the ID and vertices of that cell?
2- What if that point is exactly a mesh point? Which cell ID will I get as the answer?

Here is a minimal code with a mesh and the x,y,z coordinates of a point which I am sure lies strictly inside of a cell. How can I answer the above questions for this code?

```
from dolfin import *
meshlevel =10
degree = 1
dim = 2
mesh = UnitDiscMesh.create(MPI.comm_world,meshlevel, degree, dim)
Px = -0.519
Py = 0.33
Pz = 0.0
```

dokken
October 19, 2021, 8:26pm
2
The following code illustrates all your questions:

```
from dolfin import *
meshlevel =10
degree = 1
dim = 2
mesh = UnitDiscMesh.create(MPI.comm_world,meshlevel, degree, dim)
Px = -0.519
Py = 0.33
Pz = 0.0
p = Point(Px, Py, Pz)
tree = BoundingBoxTree()
tree.build(mesh)
first = tree.compute_first_entity_collision(p)
print(first)
c_to_v = mesh.topology()(mesh.topology().dim(), 0)
print(c_to_v(first), mesh.coordinates()[c_to_v(first)])
p0 = Point(mesh.coordinates()[0])
all_cells = tree.compute_entity_collisions(p0)
print(mesh.coordinates()[0])
for cell in all_cells:
print(cell, c_to_v(cell), mesh.coordinates()[c_to_v(cell)])
```

You can choose to either get the first entity collision (as traversed through the boundingboxtree) or get all of them and then choose one of them.

1 Like

nami
October 19, 2021, 8:33pm
3

dokken:

```
from dolfin import *
meshlevel =10
degree = 1
dim = 2
mesh = UnitDiscMesh.create(MPI.comm_world,meshlevel, degree, dim)
Px = -0.519
Py = 0.33
Pz = 0.0
p = Point(Px, Py, Pz)
tree = BoundingBoxTree()
tree.build(mesh)
first = tree.compute_first_entity_collision(p)
print(first)
c_to_v = mesh.topology()(mesh.topology().dim(), 0)
print(c_to_v(first), mesh.coordinates()[c_to_v(first)])
p0 = Point(mesh.coordinates()[0])
all_cells = tree.compute_entity_collisions(p0)
print(mesh.coordinates()[0])
for cell in all_cells:
print(cell, c_to_v(cell), mesh.coordinates()[c_to_v(cell)])
```

Wow. This is exactly what I was looking for. So the key is BoundingBoxTree. Thank you so much.