# Error with compute\_closest\_entity

**URL:** <https://fenicsproject.discourse.group/t/error-with-compute-closest-entity/10969>\
**Category:** dolfinx\
**Created:** [April 18, 2023, 6:51pm UTC](https://fenicsproject.discourse.group/t/error-with-compute-closest-entity/10969 "2023-04-18T18:51:32Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Karlosjusus](https://avatars.discourse-cdn.com/v4/letter/k/cc9497/32.png) [@Karlosjusus](https://fenicsproject.discourse.group/u/Karlosjusus)\
**Post date:** [April 18, 2023, 6:51pm UTC](https://fenicsproject.discourse.group/t/error-with-compute-closest-entity/10969/1 "2023-04-18T18:51:32Z")

</div>

Hi everyone!  
I used dolfinx 0.4.0 version and my code runs. But i upgraded ububtu to version 22.04 to install dolfinx 0.5.2.  
And now I have a problem with dolfinxgeometry.compute\_closest\_entity.  
Now I get this error in a vertex…

```auto

RuntimeError: GJK error: max iteration limit reached

```

my code

```auto
import dolfinx
from mpi4py import MPI
import numpy as np

vx_np = np.loadtxt('vx.txt')
vx_aux = vx_np[:,1:4]

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "venturi.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid")
    cf = xdmf.read_meshtags(domain, name = 'Grid')

tdim = domain.topology.dim
tree = dolfinx.geometry.BoundingBoxTree(domain, tdim)
num_entities_local = domain.topology.index_map(tdim).size_local + domain.topology.index_map(tdim).num_ghosts
entities = np.arange(num_entities_local, dtype = np.int32)
mid_tree = dolfinx.geometry.create_midpoint_tree(domain, tdim,entities)
mesh_geom = domain.geometry.x
tol = 1e-5

def closest_point(mesh, point,tol):
    points = point
    entity = dolfinx.geometry.compute_closest_entity(tree, mid_tree, mesh, points)
    geom_dof = dolfinx.cpp.mesh.entities_to_geometry(mesh,tdim,[entity], False)
    #print('geom',geom_dof)
    mesh_nodes = mesh_geom[geom_dof][0]
    #print('mesh_nodes', mesh_nodes)
    dis = dolfinx.geometry.compute_distance_gjk(points, mesh_nodes)
    #print('dis',dis)
    index = -100
    for i in range(len(mesh_nodes)):
        #print(mesh_nodes[i])
        if np.abs(mesh_nodes[i][0]-points[0]+ dis[0]) < tol and np.abs(mesh_nodes[i][1]-points[1]+ dis[1]) < tol and np.abs(mesh_nodes[i][2]-points[2]+ dis[2]) < tol :
            index = i
            
    if (index == -100):
        return point, index
    else:
        return points - dis , geom_dof[0][index]

b = []
indices = []
ind = int(0)
for i in range(len(vx_aux)):
    if(i%20000 == 0):
        print(i)
    a,j = closest_point(domain,vx_aux[i],tol)
    if (j == -100):
        b.append(a)
        indices.append(0)
    else:
        indices.append(j)

```

The flies are here [136.02 MB folder on MEGA](https://mega.nz/folder/lQUnWa4D#WJhvRPoCgiO_0736IRG6IQ)

Thank you!

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 19, 2023, 8:28am UTC](https://fenicsproject.discourse.group/t/error-with-compute-closest-entity/10969/2 "2023-04-19T08:28:22Z")

</div>

This code is not complete, as you are missing definitions of your `domain`. Please make sure that your code is executable by others.

---

<div class="post-metadata">

**Author:** ![Karlosjusus](https://avatars.discourse-cdn.com/v4/letter/k/cc9497/32.png) [@Karlosjusus](https://fenicsproject.discourse.group/u/Karlosjusus)\
**Post date:** [April 19, 2023, 1:38pm UTC](https://fenicsproject.discourse.group/t/error-with-compute-closest-entity/10969/3 "2023-04-19T13:38:24Z")

</div>

Sorry… I’ve updated the code.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 19, 2023, 3:16pm UTC](https://fenicsproject.discourse.group/t/error-with-compute-closest-entity/10969/4 "2023-04-19T15:16:03Z")

</div>

NOte that your example can be reduced to:

```python
import dolfinx
from mpi4py import MPI
import numpy as np

vx_np = np.loadtxt('vx.txt')
vx_aux = vx_np[:,1:4]

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "venturi_coarse.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid")
tdim = domain.topology.dim

num_entities_local = domain.topology.index_map(tdim).size_local + domain.topology.index_map(tdim).num_ghosts
entities = np.arange(num_entities_local, dtype = np.int32)
tree = dolfinx.geometry.BoundingBoxTree(domain, tdim, entities)
mid_tree = dolfinx.geometry.create_midpoint_tree(domain, tdim,entities)

N = 48369 #fails
#N = 48368 #works
assert (N<vx_aux.shape[0])
entity = dolfinx.geometry.compute_closest_entity(tree, mid_tree, domain, vx_aux[:N,:])

```

This is quite strange (i’ll try to look into what’s up).  
However, I would strongly avoid the loops you do in your python code, calling the code as

```python
entities = dolfinx.geometry.compute_closest_entity(tree, mid_tree, domain, vx_aux)

```

would avoid large python loops

---

<div class="post-metadata">

**Author:** ![Karlosjusus](https://avatars.discourse-cdn.com/v4/letter/k/cc9497/32.png) [@Karlosjusus](https://fenicsproject.discourse.group/u/Karlosjusus)\
**Post date:** [April 20, 2023, 3:29pm UTC](https://fenicsproject.discourse.group/t/error-with-compute-closest-entity/10969/5 "2023-04-20T15:29:57Z")

</div>

Thanks for the suggestion @dokken.  
Only 2 points fails in the previous code.  
I’ve modified the code as follow and it works.

```auto
def closest_point(mesh, point,tol):
    points = point
    while True:
        try:
            entity = dolfinx.geometry.compute_closest_entity(tree, mid_tree, mesh, points)
            break
        except RuntimeError:
            print(points)
            for j in range(len(mesh.geometry.x)):
                if (np.abs(mesh_geom[j][0]-points[0])< tol and np.abs(mesh_geom[j][1]-points[1])< tol and np.abs(mesh_geom[j][2]-points[2])< tol):
                     return points, j
    
    geom_dof = dolfinx.cpp.mesh.entities_to_geometry(mesh,tdim,[entity], False)
    #print('geom',geom_dof)
    mesh_nodes = mesh_geom[geom_dof][0]
    #print('mesh_nodes', mesh_nodes)
    dis = dolfinx.geometry.compute_distance_gjk(points, mesh_nodes)
    for i in range(len(mesh_nodes)):
        if np.abs(mesh_nodes[i][0]-points[0]+ dis[0]) < tol and np.abs(mesh_nodes[i][1]-points[1]+ dis[1]) < tol and np.abs(mesh_nodes[i][2]-points[2]+ dis[2]) < tol :
            index = i
    return points - dis , geom_dof[0][index]

```