# Nearest point on boundary to exterior point

How can I determine the point on the boundary of a mesh that is closest to a given point that is outside of the mesh?

I can determine the closest mesh vertex, but this is not what I want. I want to find the coordinates of the closest point on the boundary, which may be on an edge, or face, rather than a vertex.

Below I made an image showing this:

Are there any built-in features in Fenics that can do this? Thanks for any help!

I donâ€™t know of any built-in functionality for this. I did think of a somewhat-involved workaround:

• Define a CG1 scalar space, `V = FunctionSpace(mesh,"CG",1)`.
• Find the index of the closest cell to `x` in `mesh` using the `compute_first_entity_collision` method of `mesh.bounding_box_tree()`.
• Use some DoF maps for `V` and `V.element().evaluate_basis` to evaluate all shape functions on that cell. The shape function that is negative is the one whose node is the vertex â€śoppositeâ€ť `x`.
• Get that shape functionâ€™s gradient with `V.element().evaluate_basis_derivatives`.
• Given `x`, the shape function evaluation at `x`, and the shape function gradient, you can do some algebra to calculate the closest point to `x` on the cell (where the shape function goes to zero along a line through `x` in the direction of the gradient).

For an example of how to use the basis function evaluation functions, you might take a look at this discussion.

1 Like

There is functionality for these features in dolfinx

MWE:

``````import dolfinx
import dolfinx.geometry
from mpi4py import MPI
import numpy as np

mesh = dolfinx.UnitCubeMesh(MPI.COMM_WORLD, 10, 10, 10)

def closest_point_in_mesh(mesh, point):
points = np.reshape(point, (1,3))
tdim = mesh.topology.dim
bounding_box = dolfinx.geometry.BoundingBoxTree(mesh, tdim)
entity, distance = dolfinx.geometry.compute_closest_entity(bounding_box, points[0], mesh)
mesh_geom = mesh.geometry.x
geom_dofs = dolfinx.cpp.mesh.entities_to_geometry(mesh, tdim, [entity], False)
mesh_nodes = mesh_geom[geom_dofs][0]
displacement =  dolfinx.geometry.compute_distance_gjk(points, mesh_nodes)
return points[0] - displacement

point = [1.25, 0.3, 0]
print(closest_point_in_mesh(mesh, point))

point = [1.25, 1.3, 0.1]
print(closest_point_in_mesh(mesh, point))

point = [-1.2, 0.273, -0.3]
print(closest_point_in_mesh(mesh, point))

point = [0.31, 0.57, -0.2]
print(closest_point_in_mesh(mesh, point))
``````

returning

``````[1.  0.3 0. ]
[1.  1.  0.1]
[0.    0.273 0.   ]
[3.10000000e-01 5.70000000e-01 2.77555756e-17]
``````
3 Likes

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:

1. Find which simplex (e.g., tetrahedron) is closest to the point with a bounding box tree.
2. 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).
3. Reject all projections that do not lie within the corresponding sub-simplex.
4. 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):
``````

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.plot([p1[0], projected_p1[0]],
[p1[1], projected_p1[1]], 'r')
t2 = plt.Polygon(V2)
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)
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:

1 Like

Thanks for this answer. Surely this is the way to go for dolfinx. Unfortunately, I canâ€™t use it because I am working on a project that uses dolfin, not dolfinx. I was able to come up with a solution that works for dolfin, and posted it below.

Thanks for this answer. I had some issues vectorizing this method. But, inspired by this, I was able to come up with another somewhat related method that was easier to vectorize, which I posted below.