Error with compute_closest_entity

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…


RuntimeError: GJK error: max iteration limit reached

my code

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

Thank you!

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

Sorry… I’ve updated the code.

NOte that your example can be reduced to:

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

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

would avoid large python loops

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

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]