psmile_geometry.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
! Copyright 2010, DKRZ, Hamburg, Germany.
! All rights reserved. Use is subject to OASIS4 license terms.
!-----------------------------------------------------------------------
!
! !DESCRIPTION:
!
! Module psmile_geometry contains function for testing whether a point
! is on a line, triangle, or quadrangle
!
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 16.09.10    M. Hanke     created
!
!----------------------------------------------------------------------
!
!  $Id: psmile_geometry.F90 2727 2010-11-12 13:37:17Z hanke $
!  $Author: hanke $
!
!----------------------------------------------------------------------


   module psmile_geometry

      use psmile, only : point_dble, point_real

      implicit none

      double precision, parameter :: epsilon_dble = 1.0d-13
      real, parameter             :: epsilon_real = 1.0d-13

      private ! by default everything is private in this module

      interface point_is_on_line
         module procedure point_is_on_line_dble, &
                          point_is_on_line_real
      end interface

      interface point_is_in_triangle
         module procedure point_is_in_triangle_dble, &
                          point_is_in_triangle_real
      end interface

      interface point_is_in_quadrangle
         module procedure point_is_in_quadrangle_dble, &
                          point_is_in_quadrangle_real
      end interface

      interface intersect
         module procedure rectangular_shapes_intersect_int
      end interface

      public :: point_is_on_line
      public :: point_is_in_triangle
      public :: point_is_in_quadrangle
      public :: epsilon_dble
      public :: epsilon_real
      public :: intersect

      contains

         ! returns true if point is on line
         logical function point_is_on_line_dble(point, line)
            type (point_dble), intent(in) :: point
            type (point_dble), intent(in) :: line(2)

            double precision :: u, norm_err

            ! if A and B have the same x coordinate
            if (line(2)%x == line(1)%x) then

               ! if A, B and P have the same x coordinate and
               ! P is between A and B
               if ( point%x == line(1)%x .and. &
                  (point%y >= line(1)%y .and. point%y <= line(2)%y .or. &
                  point%y <= line(1)%y .and. point%y >= line(2)%y)) then
                  point_is_on_line_dble = .true.
               else
                  point_is_on_line_dble = .false.
               endif
            else

               !formula: A + (B - A) * u = P
               ! u = (Px - Ax) / (Bx - Ax)
               u = (point%x - line(1)%x) / (line(2)%x - line(1)%x)

               if (abs(u) > 1 + epsilon_dble) then
                  ! point P might be on a line which goes through AB
                  ! but is outside of the range of AB
                  point_is_on_line_dble = .false.
               else
                  ! computes normalised error of solution u
                  ! error = |(Ay + (By - Ay) * u - Py) / Py| = 0
                  norm_err = line(1)%y + (line(2)%y - line(1)%y) * u - point%y
                  if (point%y /= 0.0d0) norm_err = norm_err / point%y

                  if (abs(norm_err) > epsilon_dble) then
                     ! point is not on line AB
                     point_is_on_line_dble = .false.
                  else
                     ! point is on line AB
                     point_is_on_line_dble = .true.
                  endif
               endif
            endif

         end function point_is_on_line_dble

         logical function point_is_on_line_real(point, line)
            type (point_real), intent(in) :: point
            type (point_real), intent(in) :: line(2)

            real :: u, norm_err

            ! if A and B have the same x coordinate
            if (line(2)%x == line(1)%x) then

               ! if A, B and P have the same x coordinate and
               ! P is between A and B
               if ( point%x == line(1)%x .and. &
                  (point%y >= line(1)%y .and. point%y <= line(2)%y .or. &
                  point%y <= line(1)%y .and. point%y >= line(2)%y)) then
                  point_is_on_line_real = .true.
               else
                  point_is_on_line_real = .false.
               endif
            else

               !formula: A + (B - A) * u = P
               ! u = (Px - Ax) / (Bx - Ax)
               u = (point%x - line(1)%x) / (line(2)%x - line(1)%x)

               if (abs(u) > 1 + epsilon_real) then
                  ! point P might be on a line which goes through AB
                  ! but is outside of the range of AB
                  point_is_on_line_real = .false.
               else
                  ! computes normalised error of solution u
                  ! error = |(Ay + (By - Ay) * u - Py) / Py| = 0
                  norm_err = line(1)%y + (line(2)%y - line(1)%y) * u - point%y
                  if (point%y /= 0.0d0) norm_err = norm_err / point%y

                  if (abs(norm_err) > epsilon_real) then
                     ! point is not on line AB
                     point_is_on_line_real = .false.
                  else
                     ! point is on line AB
                     point_is_on_line_real = .true.
                  endif
               endif
            endif

         end function point_is_on_line_real

         !algorithm is taken from: http://www.drdobbs.com/184404201?pgno=2
         !it is based on barycentric coordinates
         !elemental
         logical function point_is_in_triangle_dble(point, triangle)
            type (point_dble), intent(in) :: point
            type (point_dble), intent(in) :: triangle(3)

            type (point_dble) :: e1, e2, e3
            double precision  :: u, v, d

            ! formula:
            ! (P-A) = (C-A)*u + (B-A)*v
            ! e1    = e3   *u + e2   *v

            e1%x = point%x - triangle(1)%x
            e2%x = triangle(2)%x - triangle(1)%x
            e3%x = triangle(3)%x - triangle(1)%x

            e1%y = point%y - triangle(1)%y
            e2%y = triangle(2)%y - triangle(1)%y
            e3%y = triangle(3)%y - triangle(1)%y

            ! if A and B have the same x coordinate
            if (e2%x == 0) then
               ! if A, B and C have the same x coordinate
               if (e3%x == 0) then
                  point_is_in_triangle_dble = point_is_on_line_dble(point, triangle(1:2)) .or. &
                                             point_is_on_line_dble(point, (/triangle(1), triangle(3)/))
                  return
               endif
               u = e1%x / e3%x
               if (u < - epsilon_dble .or. u > 1 + epsilon_dble) then
                  point_is_in_triangle_dble = .false.
                  return
               endif
               ! if A and B have the same x and y coordinates
               if (e2%y == 0) then
                  point_is_in_triangle_dble = point_is_on_line_dble(point, (/triangle(1), triangle(3)/))
                  return
               endif
               v = (e1%y - e3%y * u) / e2%y
               if (v < - epsilon_dble) then
                  point_is_in_triangle_dble = .false.
                  return
               endif
            else ! A and B have different x coordinates
               d = e3%y * e2%x - e3%x * e2%y
               if (d == 0) then
                  point_is_in_triangle_dble = .false.
                  return
               endif;
               u = (e1%y * e2%x - e1%x * e2%y) / d
               if (u < - epsilon_dble .or. u > 1 + epsilon_dble) then
                  point_is_in_triangle_dble = .false.
                  return
               endif
               v = (e1%x - e3%x * u) / e2%x
               if (v <  - epsilon_dble) then
                  point_is_in_triangle_dble = .false.
                  return
               endif
            endif

            point_is_in_triangle_dble = (u + v <= 1 + epsilon_dble)

         end function point_is_in_triangle_dble

         logical function point_is_in_triangle_real(point, triangle)
            type (point_real), intent(in) :: point
            type (point_real), intent(in) :: triangle(3)

            type (point_real) :: e1, e2, e3
            real              :: u, v, d

            ! formula:
            ! (P-A) = (C-A)*u + (B-A)*v
            ! e1    = e3   *u + e2   *v

            e1%x = point%x - triangle(1)%x
            e2%x = triangle(2)%x - triangle(1)%x
            e3%x = triangle(3)%x - triangle(1)%x

            e1%y = point%y - triangle(1)%y
            e2%y = triangle(2)%y - triangle(1)%y
            e3%y = triangle(3)%y - triangle(1)%y

            ! if A and B have the same x coordinate
            if (e2%x == 0) then
               ! if A, B and C have the same x coordinate
               if (e3%x == 0) then
                  point_is_in_triangle_real = point_is_on_line_real(point, triangle(1:2)) .or. &
                                             point_is_on_line_real(point, (/triangle(1), triangle(3)/))
                  return
               endif
               u = e1%x / e3%x
               if (u < - epsilon_real .or. u > 1 + epsilon_real) then
                  point_is_in_triangle_real = .false.
                  return
               endif
               ! if A and B have the same x and y coordinates
               if (e2%y == 0) then
                  point_is_in_triangle_real = point_is_on_line_real(point, (/triangle(1), triangle(3)/))
                  return
               endif
               v = (e1%y - e3%y * u) / e2%y
               if (v < - epsilon_real) then
                  point_is_in_triangle_real = .false.
                  return
               endif
            else ! A and B have different x coordinates
               d = e3%y * e2%x - e3%x * e2%y
               if (d == 0) then
                  point_is_in_triangle_real = .false.
                  return
               endif;
               u = (e1%y * e2%x - e1%x * e2%y) / d
               if (u < - epsilon_real .or. u > 1 + epsilon_real) then
                  point_is_in_triangle_real = .false.
                  return
               endif
               v = (e1%x - e3%x * u) / e2%x
               if (v <  - epsilon_real) then
                  point_is_in_triangle_real = .false.
                  return
               endif
            endif

            point_is_in_triangle_real = (u + v <= 1 + epsilon_real)

         end function point_is_in_triangle_real

         !elemental
         logical function point_is_in_quadrangle_dble(point, quadrangle)
            type (point_dble), intent(in) :: point
            type (point_dble), intent(in) :: quadrangle(4)

            logical :: inside_ABC, inside_ACD

            inside_ABC = point_is_in_triangle_dble(point, quadrangle(1:3))
            inside_ACD = point_is_in_triangle_dble(point, &
               (/quadrangle(1), quadrangle(3), quadrangle(4)/))

            !neqv == xor
            if (inside_ABC .neqv. inside_ACD) then
               ! standard in-case:
               ! point is either in triangle ABC or in triangle ACD
               point_is_in_quadrangle_dble = .true.
            else if (.not.(inside_ABC .or. inside_ACD)) then
               ! standard out-case:
               ! point is neither in triangle ABC or in triangle ACD
               point_is_in_quadrangle_dble = .false.
            else
               ! point is in both triangles
               ! the point lies directly on line AC and/or the quadrangle is concave
               if (point_is_on_line_dble(point, (/quadrangle(1), quadrangle(3)/))) then
                  ! the point is on the line AC
                  if (point%x == quadrangle(1)%x .and. point%y == quadrangle(1)%y .or. &
                     point%x == quadrangle(2)%x .and. point%y == quadrangle(2)%y .or. &
                     point%x == quadrangle(3)%x .and. point%y == quadrangle(3)%y .or. &
                     point%x == quadrangle(4)%x .and. point%y == quadrangle(4)%y) then
                     ! point is directly on either A, B, C or D
                     point_is_in_quadrangle_dble = .true.
                  else
                     if (point_is_in_triangle_dble(quadrangle(4), quadrangle(1:3)) .or. &
                        point_is_in_triangle_dble(quadrangle(2), &
                           (/quadrangle(1), quadrangle(3), quadrangle(4)/))) then
                        ! the quadrangle is concave (either D lies in ABC or B lies in ACD)
                        ! cases were A is in BCD or C is in ABD are no problem for this algorithm
                        ! and are therefore not especially handled here
                        point_is_in_quadrangle_dble = .false.
                     else
                        ! the quadrangle is convex
                        point_is_in_quadrangle_dble = .true.
                     endif
                  endif
               else
                  ! the quadrangle has to be concave
                  if (point_is_in_triangle_dble(quadrangle(4), quadrangle(1:3))) then
                     ! point D is in ABC
                     if (point_is_on_line_dble(point, (/quadrangle(1), quadrangle(4)/)) .or. &
                        point_is_on_line_dble(point, quadrangle(3:4))) then
                        ! point is either on line AD or CD
                        point_is_in_quadrangle_dble = .true.
                     else
                        point_is_in_quadrangle_dble = .false.
                     endif
                  else
                     ! point B is in ACD
                     if (point_is_on_line_dble(point, quadrangle(2:3)) .or. &
                        point_is_on_line_dble(point, quadrangle(1:2))) then
                        ! point is either on line AB or BC
                        point_is_in_quadrangle_dble = .true.
                     else
                        point_is_in_quadrangle_dble = .false.
                     endif
                  endif
               endif
            endif

         end function point_is_in_quadrangle_dble

         logical function point_is_in_quadrangle_real(point, quadrangle)
            type (point_real), intent(in) :: point
            type (point_real), intent(in) :: quadrangle(4)

            logical :: inside_ABC, inside_ACD

            inside_ABC = point_is_in_triangle_real(point, quadrangle(1:3))
            inside_ACD = point_is_in_triangle_real(point, &
               (/quadrangle(1), quadrangle(3), quadrangle(4)/))

            !neqv == xor
            if (inside_ABC .neqv. inside_ACD) then
               ! standard in-case:
               ! point is either in triangle ABC or in triangle ACD
               point_is_in_quadrangle_real = .true.
            else if (.not.(inside_ABC .or. inside_ACD)) then
               ! standard out-case:
               ! point is neither in triangle ABC or in triangle ACD
               point_is_in_quadrangle_real = .false.
            else
               ! point is in both triangles
               ! the point lies directly on line AC and/or the quadrangle is concave
               if (point_is_on_line_real(point, (/quadrangle(1), quadrangle(3)/))) then
                  ! the point is on the line AC
                  if (point%x == quadrangle(1)%x .and. point%y == quadrangle(1)%y .or. &
                      point%x == quadrangle(2)%x .and. point%y == quadrangle(2)%y .or. &
                      point%x == quadrangle(3)%x .and. point%y == quadrangle(3)%y .or. &
                      point%x == quadrangle(4)%x .and. point%y == quadrangle(4)%y) then
                     ! point is directly on either A, B, C or D
                     point_is_in_quadrangle_real = .true.
                  else
                     if (point_is_in_triangle_real(quadrangle(4), quadrangle(1:3)) .or. &
                         point_is_in_triangle_real(quadrangle(2), &
                           (/quadrangle(1), quadrangle(3), quadrangle(4)/))) then
                        ! the quadrangle is concave (either D lies in ABC or B lies in ACD)
                        ! cases were A is in BCD or C is in ABD are no problem for this algorithm
                        ! and are therefore not especially handled here
                        point_is_in_quadrangle_real = .false.
                     else
                        ! the quadrangle is convex
                        point_is_in_quadrangle_real = .true.
                     endif
                  endif
               else
                  ! the quadrangle has to be concave
                  if (point_is_in_triangle_real(quadrangle(4), quadrangle(1:3))) then
                     ! point D is in ABC
                     if (point_is_on_line_real(point, (/quadrangle(1), quadrangle(4)/)) .or. &
                         point_is_on_line_real(point, quadrangle(3:4))) then
                        ! point is either on line AD or CD
                        point_is_in_quadrangle_real = .true.
                     else
                        point_is_in_quadrangle_real = .false.
                     endif
                  else
                     ! point B is in ACD
                     if (point_is_on_line_real(point, quadrangle(2:3)) .or. &
                         point_is_on_line_real(point, quadrangle(1:2))) then
                        ! point is either on line AB or BC
                        point_is_in_quadrangle_real = .true.
                     else
                        point_is_in_quadrangle_real = .false.
                     endif
                  endif
               endif
            endif

         end function point_is_in_quadrangle_real

         logical function rectangular_shapes_intersect_int (a, b)
            integer, intent(in) :: a(:,:), b(:,:)

            integer :: i

#ifdef PRISM_ASSERTION
            if (size (a, 1) /= 2 .or. size (b, 1) /= 2) then
               call psmile_assert (__FILE__, __LINE__, &
                  "wrong first dimension of argument a or b in rectangles_cuboids_intersect_int")
            endif
            if (size (a,2) /= size (b,2)) then
               call psmile_assert (__FILE__, __LINE__, &
                  "none matching arguments in function rectangles_cuboids_intersect_int")
            endif
#endif
            rectangular_shapes_intersect_int = .not. any (a(1,:) > b(2,:) .or. a(2,:) < b(1,:))
         end function rectangular_shapes_intersect_int

   end module psmile_geometry