The answer by dokkenâ€™s is ideal for projects using dolfinx, but it does not work for me because the project Iâ€™m working on is in dolfin, not dolfinx.

Here I describe the solution that I came up with, and post the code for anyone else to use. The method is inspired by, but different than, the outline given by kamensky. It only works for simplicial meshes (triangles/tetrahedra).

The method is as follows:

- Find which simplex (e.g., tetrahedron) is closest to the point with a bounding box tree.
- Project the point onto all affine subspaces generated by â€śsub-simplicesâ€ť of the simplex (e.g., for a tetrahedron the sub-simplices are the simplex itself, each of the triangular faces, each of the line segment edges, and each of the vertices).
- Reject all projections that do not lie within the corresponding sub-simplex.
- Of all the remaining (not rejected) projections, choose the closest one.

The reason why I did it this way, instead of the way suggested by kameksky, is because this way was more easy to vectorize. I.e., the code can find the closest points for many points at once using bulk numpy operations (fast), rather than finding the projections one at a time within a for loop (slow).

Here is the vectorized code, which is fast, but somewhat involved:

**project_point_onto_affine_subspace.py:**

```
import numpy as np
def project_point_onto_affine_subspace(p, V):
"""Projects a point onto an affine subspace
p.shape = (N,) OR (num_pts, N) for vectorization over many points / affine subspaces
V.shape = (k, N) OR (num_pts, k, N) for vectorization
N = dimension of ambient space
k-1 = dimension of affine subspace
For a single point / affine subspace:
- p is the point in R^N to be projected onto the affine subspace
- The affine subspace is the set of all affine combinations
of rows, V[i,:], of the matrix V
Vectorization:
For many points / affine subspaces, each point is
projected onto its corresponding affine subspace
p[i,:] is the ith point
V[i,:,:] is the matrix defining the ith affine subspace
Example usage:
import numpy as np
spatial_dim = 3
p1 = np.random.randn(spatial_dim)
V1 = np.array([[2., 0., 0.],
[2., 1., 0.],
[2., 0., 1.]])
p2 = np.random.randn(spatial_dim)
V2 = np.array([[0., 3., 0.],
[1., 3., 0.],
[0., 3., 1.]])
p = np.stack([p1, p2])
V = np.stack([V1, V2])
projected_p, affine_coordinates = project_point_onto_affine_subspace(p, V)
projected_p1 = projected_p[0,:]
projected_p2 = projected_p[1,:]
print('p1=', p1)
print('projected_p1=', projected_p1)
print('p2=', p2)
print('projected_p2=', projected_p2)
"""
if len(p.shape) == 1:
PP = p.reshape((1, p.shape[0]))
VV = V.reshape((1, V.shape[0], V.shape[1]))
else:
PP = p
VV = V
num_pts, k, N = VV.shape
if k == 1:
affine_coordinates = np.ones((num_pts, 1))
else:
VV0 = VV[:, 0, :].reshape((num_pts, 1, N))
dVV = VV[:, 1: :].reshape((num_pts, k-1, N)) - VV0
dPP = PP.reshape((num_pts, 1, N)) - VV0
if k == 2:
PHI = np.sum(dVV * dVV, axis=-1).reshape((num_pts))
RHS = np.sum(dVV * dPP, axis=-1).reshape((num_pts))
cc_rest = (RHS / PHI).reshape((num_pts, 1))
else:
PHI = np.einsum('xiz,xjz->xij', dVV, dVV) # shape = (num_pts, k-1, k-1)
iPHI = np.linalg.inv(PHI) # shape = (num_pts, k-1, k-1)
RHS = np.sum(dVV * dPP, axis=-1) # shape = (num_pts, k-1)
cc_rest = np.einsum('pij,pj->pi', iPHI, RHS) # shape = (num_pts, k-1)
cc_first = (1. - np.sum(cc_rest, axis=1)).reshape((num_pts, 1))
affine_coordinates = np.concatenate([cc_first, cc_rest], axis=1) # shape = (num_pts, k)
PP_projected = np.einsum('pi,pij->pj', affine_coordinates, VV) # shape = (num_pts, N)
if len(p.shape) == 1:
PP_projected = PP_projected.reshape(-1)
affine_coordinates = affine_coordinates.reshape(-1)
return PP_projected, affine_coordinates
```

**powerset.py:**

```
def powerset(s):
# NOT MY CODE. FROM USER "hughdbrown" ON STACKOVERFLOW HERE:
# https://stackoverflow.com/a/1482320/484944
x = len(s)
masks = [1 << i for i in range(x)]
for i in range(1 << x):
yield [ss for mask, ss in zip(masks, s) if i & mask]
```

**closest_point_on_simplex.py:**

```
import numpy as np
from powerset import powerset
from project_point_onto_affine_subspace import project_point_onto_affine_subspace
def closest_point_on_simplex(p, V):
"""Projects a point onto a simplex (triangle, tetrahedron, etc)
p.shape = (N,) OR (num_pts, N) for vectorization over many points/simplices
V.shape = (k, N) OR (num_pts, k, N) for vectorization
N = dimension of ambient space
k-1 = dimension of simplex
For a single point/simplex:
- p is the point in R^N to be projected onto the simplex
- The simplex is the set of all convex combinations
of rows, V[i,:], of the matrix V
Vectorization:
For many points/simplices, each point is
projected onto its corresponding simplex
p[i,:] is the ith point
V[i,:,:] is the matrix defining the ith simplex
Example usage:
import numpy as np
import matplotlib.pyplot as plt
p1 = np.array([1.1, 0.4])
V1 = np.array([[0., 0.],
[0., 1.],
[1., 0.]])
p2 = np.array([-0.3, 1.1])
V2 = np.array([[-1.0, 0.],
[0., 0.],
[-0.5, 0.5]])
p = np.stack([p1, p2])
V = np.stack([V1, V2])
projected_p = closest_point_on_simplex(p, V)
projected_p1 = projected_p[0,:]
projected_p2 = projected_p[1,:]
plt.figure()
t1 = plt.Polygon(V1)
plt.gca().add_patch(t1)
plt.plot([p1[0], projected_p1[0]],
[p1[1], projected_p1[1]], 'r')
t2 = plt.Polygon(V2)
plt.gca().add_patch(t2)
plt.plot([p2[0], projected_p2[0]],
[p2[1], projected_p2[1]], 'r')
plt.gca().set_aspect('equal')
plt.show()
"""
if len(p.shape) == 1:
PP = p.reshape((1, p.shape[0]))
VV = V.reshape((1, V.shape[0], V.shape[1]))
else:
PP = p
VV = V
num_pts, k, N = VV.shape
subsets = list(powerset(list(range(k)))) # e.g., [[], [0], [1], [2], [0,1], [0,2], [1,2], [0,1,2]]
QQ = list()
CC = list()
for s in subsets:
if s:
simplicial_facet = VV[:,s,:]
Q, C = project_point_onto_affine_subspace(PP, simplicial_facet)
QQ.append(Q)
CC.append(C)
distances = np.stack([np.linalg.norm(Q - PP, axis=-1) for Q in QQ]) # shape=(num_facets, num_pts)
good_inds = np.stack([(np.all(0. <= C, axis=1) & np.all(C <= 1., axis=1)) for C in CC]) # shape=(num_facets, num_pts)
bad_inds = np.logical_not(good_inds)
distances[bad_inds] = np.inf
closest_inds = np.expand_dims(np.argmin(distances[:, :, None], axis=0), axis=0)
QQ_stack = np.stack(QQ) # shape=(num_facets, num_pts, N)
PP_projected = np.take_along_axis(QQ_stack, closest_inds, axis=0)[0,:,:]
if len(p.shape) == 1:
PP_projected = PP_projected.reshape(-1)
return PP_projected
```

**closest_point_in_mesh.py:**

```
import numpy as np
import dolfin as dl
from closest_point_on_simplex import closest_point_on_simplex
def closest_point_in_mesh(p, mesh):
"""Finds the nearest point in a mesh to a given point
p is the point or points (numpy array)
mesh is the fenics/dolfin mesh
p.shape = (N,) OR (num_pts, N) for vectorization over many points
N = dimension of ambient space
Example usage:
import numpy as np
import dolfin as dl
import matplotlib.pyplot as plt
mesh = dl.UnitSquareMesh(13,9)
num_pts = 20
p = np.random.randn(num_pts, 2) + np.array([0.5, 0.5])
closest_p = closest_point_in_mesh(p, mesh)
plt.figure()
dl.plot(mesh)
for ii in range(num_pts):
plt.plot([p[ii,0], closest_p[ii,0]], [p[ii,1], closest_p[ii,1]], 'b')
plt.plot(p[ii,0], p[ii,1], '*k')
plt.plot(closest_p[ii,0], closest_p[ii,1], '.r')
plt.show()
"""
if len(p.shape) == 1:
PP = p[None,:]
else:
PP = p
num_pts, N = PP.shape
tdim = mesh.topology().dim()
k = tdim + 1
VV = np.zeros((num_pts, k, N))
bbt = mesh.bounding_box_tree()
for ii in range(num_pts):
pi = PP[ii,:]
closest_entity, closest_distance = bbt.compute_closest_entity(dl.Point(pi))
closest_cell = mesh.cells()[closest_entity]
vertices_of_closest_cell = mesh.coordinates()[closest_cell, :]
VV[ii, :, :] = vertices_of_closest_cell
closest_PP = closest_point_on_simplex(PP, VV)
if len(p.shape) == 1:
closest_PP = closest_PP.reshape(-1)
return closest_PP
```

Example usage:

```
import numpy as np
import dolfin as dl
import matplotlib.pyplot as plt
from closest_point_in_mesh import closest_point_in_mesh
mesh = dl.UnitSquareMesh(13,9)
num_pts = 20
p = np.random.randn(num_pts, 2) + np.array([0.5, 0.5])
closest_p = closest_point_in_mesh(p, mesh)
plt.figure()
dl.plot(mesh)
for ii in range(num_pts):
plt.plot([p[ii,0], closest_p[ii,0]], [p[ii,1], closest_p[ii,1]], 'b')
plt.plot(p[ii,0], p[ii,1], '*k')
plt.plot(closest_p[ii,0], closest_p[ii,1], '.r')
plt.show()
```

Example output: