tem_triangle_normal_proximity Subroutine

public subroutine tem_triangle_normal_proximity(me, point, closest, closekind, distance, normal)

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 check
  • point - the point for which the closest neighbor point within the triangle is to be found

Output parameters:

  • closest - point of the triangle which is the closest to the given point
  • closekind - 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

Arguments

Type IntentOptional 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)

Calls

proc~~tem_triangle_normal_proximity~~CallsGraph proc~tem_triangle_normal_proximity tem_triangle_normal_proximity proc~cross_product3d cross_product3D proc~tem_triangle_normal_proximity->proc~cross_product3d

Source Code

  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