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:
V = FunctionSpace(mesh,"CG",1).x in mesh using the compute_first_entity_collision method of mesh.bounding_box_tree().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.V.element().evaluate_basis_derivatives.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.
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]
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:
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:
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.