Compute the point of the triangle closest to the given point.
This function computes the point of the triangle with the smallest distance to the provided point. Input parameters:
me - the triangle to checkpoint - the point for which the closest neighbor point within the
triangle is to be foundOutput parameters:
closest - point of the triangle which is the closest to the given
pointclosekind - there are three different kinds of closest point:
- it may be inside the triangle
- it may be on one of the edges
- it may be one of the three vertices
A fourth case arises with degenerate triangles for which
no proper distance can be computed, if this is set, the
other values have no proper values assigned to them.
This parameter returns which kind of closest point we have.distance - the smallest distance between point and the triangle (me)normal - depending on the closekind this returns:
- the normal vector on the triangle pointing from closest
to point (has the length distance) if closest is withing
the triangle
- the unit normal vector to the side of point if closest
is on an edge or vertex of the triangle| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(tem_triangle_type), | intent(in) | :: | me | |||
| real(kind=rk), | intent(in) | :: | point(3) | |||
| real(kind=rk), | intent(out) | :: | closest(3) | |||
| integer, | intent(out) | :: | closekind | |||
| real(kind=rk), | intent(out) | :: | distance | |||
| real(kind=rk), | intent(out) | :: | normal(3) |
subroutine tem_triangle_normal_proximity(me, point, closest, closekind, &
& distance, normal )
! ---------------------------------------------------------------------- !
type(tem_triangle_type), intent(in) :: me
real(kind=rk), intent(in) :: point(3)
real(kind=rk), intent(out) :: closest(3)
integer, intent(out) :: closekind
real(kind=rk), intent(out) :: distance
real(kind=rk), intent(out) :: normal(3)
! ---------------------------------------------------------------------- !
integer :: min_edge, min_vertex
real(kind=rk) :: edge(3,3) ! The three edges of the triangle
real(kind=rk) :: edgelen_squared(3) ! Squared length the edges
real(kind=rk) :: edgescale(3)
real(kind=rk) :: edgeprod1_2 ! Dot Product of first two edges
real(kind=rk) :: tri2p(3) ! Vector pointing from the triangles first vertex
! to the given point.
real(kind=rk) :: plane_point(3) ! Point projected onto plane of triangle
real(kind=rk) :: edge_point(3,3) ! Projection of point onto the three edges
real(kind=rk) :: plane_vector(3) ! Vector from first node to projected point
real(kind=rk) :: tri_prod(3) ! Point projections onto edges
real(kind=rk) :: tri_coord(2) ! Point in coordinates of first and second
! edge
real(kind=rk) :: vertex_distance(3) ! Distances to each of the three
! vertices
real(kind=rk) :: edge_distance(3) ! Distances to each of the three edges
real(kind=rk) :: normlen_squared
real(kind=rk) :: normlen
real(kind=rk) :: det
! ---------------------------------------------------------------------- !
closekind = tem_triangle_close_invalid
closest = 0.0_rk
distance = huge(distance)
! First: project the point onto the plane of the triangle
! if the projected point is inside the triangle the projected point
! is the `closest`.
edge(:,1) = me%nodes(:,2) - me%nodes(:,1)
edge(:,2) = me%nodes(:,3) - me%nodes(:,1)
edge(:,3) = me%nodes(:,3) - me%nodes(:,2)
normal = cross_product3D(edge(:,1), edge(:,2))
normlen_squared = normal(1)**2 + normal(2)**2 + normal(3)**2
validtri: if (normlen_squared > epsilon(normlen_squared)) then
normlen = sqrt(normlen_squared)
normal = normal / normlen
tri2p = point - me%nodes(:,1)
edgelen_squared(1) = dot_product(edge(:,1), edge(:,1))
edgelen_squared(2) = dot_product(edge(:,2), edge(:,2))
edgeprod1_2 = dot_product(edge(:,1), edge(:,2))
plane_point = point - dot_product(tri2p, normal)
plane_vector = plane_point - me%nodes(:,1)
tri_prod(1) = dot_product(plane_vector, edge(:,1))
tri_prod(2) = dot_product(plane_vector, edge(:,2))
det = edgeprod1_2**2 - (edgelen_squared(1)*edgelen_squared(2))
tri_coord(1) = ( (edgeprod1_2 * tri_prod(2)) &
& - (edgelen_squared(2) * tri_prod(1)) )
tri_coord(2) = ( (edgeprod1_2 * tri_prod(1)) &
& - (edgelen_squared(1) * tri_prod(2)) )
if (tri_coord(1) > tiny(0.0_rk) .and. tri_coord(2) > tiny(0.0_rk) &
& .and. (tri_coord(1) + tri_coord(2)) < 1.0_rk ) then
! The projected point lies within the triangle
closekind = tem_triangle_close_face
closest = plane_point
normal = point - plane_point
distance = sqrt(dot_product(normal, normal))
else
edgelen_squared(3) = dot_product(edge(:,3), edge(:,3))
tri_prod(3) = dot_product(plane_point - me%nodes(:,2), edge(:,3))
! The projected point lies outside the triangle (or on its boundary).
! Thus we check for the closest point on the boundary
! (vertices and edges).
vertex_distance(1) = sqrt(sum((point - me%nodes(:,1))**2))
vertex_distance(2) = sqrt(sum((point - me%nodes(:,2))**2))
vertex_distance(3) = sqrt(sum((point - me%nodes(:,3))**2))
min_vertex = minloc(vertex_distance, 1)
edgescale = tri_prod / edgelen_squared
! Limit the points to the length of the edge
edgescale = min(edgescale, 1.0_rk)
edgescale = max(edgescale, 0.0_rk)
edge_point(:,1) = me%nodes(:,1) + edgescale(1) * edge(:,1)
edge_point(:,2) = me%nodes(:,1) + edgescale(2) * edge(:,2)
edge_point(:,3) = me%nodes(:,2) + edgescale(3) * edge(:,3)
edge_distance(1) = sqrt(sum((point - edge_point(:,1))**2))
edge_distance(2) = sqrt(sum((point - edge_point(:,2))**2))
edge_distance(3) = sqrt(sum((point - edge_point(:,3))**2))
min_edge = minloc(vertex_distance, 1)
closekind = tem_triangle_close_edge
closest = edge_point(:, min_edge)
distance = edge_distance(min_edge)
if (vertex_distance(min_vertex) < distance + epsilon(distance)) then
closekind = tem_triangle_close_vertex
closest = me%nodes(:, min_vertex)
distance = vertex_distance(min_vertex)
end if
if (dot_product(normal, point-closest) < 0.0_rk) then
! flip the normal to point to the side of the given point
normal = -normal
end if
end if
end if validtri
end subroutine tem_triangle_normal_proximity