module polygon_m use, intrinsic :: IEEE_ARITHMETIC use, intrinsic :: iso_fortran_env use precision_m, only : FP, FP_EPS use MPI use netcdf use descriptors_m, only : ERR_UNHANDLED use screen_io_m, only : get_stderr, get_stdout, nf90_handle_err use comm_handling_m, only: is_master, get_communicator implicit none private ! Class-wide private constant integer, parameter, private :: string_buffer_length = 100 real, parameter, private :: threshold_factor = 2.5_FP ! If the polygon2d is marked as disconnected (i.e. disconnected sub-polyhedra) then ignore ! edges if the distance between two points is greater than the average by threshold_factor integer, parameter, private :: twodim = 2 ! Named constant for number of dimensions of polygon2d_t (easier to identify where factors of 2 come from) integer, parameter, private :: x = 1 ! Named constant for x dimension integer, parameter, private :: y = 2 ! Named constant for y dimension type, public, abstract :: polygon2d_t private real(FP), dimension(:, :), allocatable, public :: pts ! Array of polygon points (size N_pts x dim) integer, public :: N_pts ! Number of polygon points integer, public :: dim = twodim ! Dimensionality real(FP), dimension(:, :), allocatable, public :: p2p ! Point-to-point spacing in each dimension (first dimension length N_pts-1, second dimension corresponds to x, y) real(FP), dimension(:), allocatable, public :: p2p_mag ! Magnitude of point-to-point spacing (length N_pts-1) real(FP) :: threshold_p2p = 0.0_FP ! For Disconnectedpolygon2d, exclude edges of length > threshold_p2p in dist_to_polyedge ! For other polygon types, this is simply a placeholder contains procedure, pass(self) :: create_polygon2d ! Shared constructor (derived types extend this method) procedure, pass(self) :: dist_to_polypoint ! Finds the shortest distance between a given (x, y) point and a polygon2d vertex procedure, pass(self) :: dist_to_polyedge ! Finds the shortest distance between a given (x, y) point and any point on the polygon2d procedure, pass(self) :: closest_polyedge ! Finds the closest polygon2d to a given (x, y) point procedure, pass(self) :: polygon_print ! Prints basic information about the polygon2d procedure(polygon_string), deferred, pass(self) :: polygon_string ! For returning a string in another program procedure, pass(self) :: destroy_polygon ! Deallocate all arrays procedure, pass(self) :: signed_area ! Returns the signed_area of the polygon. >0 if counterclockwise (required) ! procedure, pass(self) :: write_polygon_points ! Writes polygon points to trunk directory (helps for post-processing) end type abstract interface character(len=string_buffer_length) function polygon_string(self) import :: polygon2d_t, string_buffer_length class(polygon2d_t), intent(in) :: self end function polygon_string end interface type, public, extends(polygon2d_t) :: closed_polygon2d_t ! Closed (area) polygon. ! First point must equal last - this is achieved in constructor contains procedure, pass(self) :: polygon_string => string_closed_polygon2d_t final :: finalize_closed_polygon2d_t procedure :: pt_inside ! Uses the 'Winding Algorithm' to determine whether a point is within a closed polygon end type closed_polygon2d_t type, public, extends(closed_polygon2d_t) :: limiting_polygon2d_t ! Closed polygon with 'local_max' and 'local_min' values (useful for defining regions with ! custom limits) real(FP), public :: local_min, local_max logical, public :: use_local_min, use_local_max contains procedure, pass(self) :: polygon_string => string_limiting_polygon2d_t final :: finalize_limiting_polygon2d_t end type limiting_polygon2d_t public :: create_closed_polygon2d_t public :: create_limiting_polygon2d_t contains subroutine polygon_print(self) class(polygon2d_t), intent(in) :: self print *, self%polygon_string() end subroutine polygon_print subroutine magnitude(array, N_pts, dim, res) integer, intent(in) :: N_pts, dim real(FP), dimension(N_pts, dim), intent(in) :: array real(FP), dimension(N_pts), intent(out) :: res res = sqrt(sum(array * array, dim=2)) end subroutine subroutine create_polygon2d(self, N_pts, X_pts, Y_pts) ! Shared constructor for polygon2d_t class(polygon2d_t), intent(inout) :: self integer, intent(in) :: N_pts ! Number of polygon2d points real(FP), dimension(N_pts), intent(in) :: X_pts, Y_pts ! arrays of (x, y) polygon2d points self%N_pts = N_pts allocate(self%pts(self%N_pts, twodim)) self%pts(:, x) = X_pts self%pts(:, y) = Y_pts allocate(self%p2p(self%N_pts-1, twodim)) allocate(self%p2p_mag(self%N_pts-1)) ! Vector of distances between adjacent vertices self%p2p(:, :) = self%pts(2:, :) - self%pts(1:self%N_pts-1, :) ! Find the magnitude of distance vector call magnitude(self%p2p, self%N_pts-1, twodim, self%p2p_mag) if (.not.(minval(self%p2p_mag) > 0.0_FP)) then write(get_stderr(), *) & "Error (polygon_m):: & &repeated points give side-length of zero. & &Will result in floating invalid" ERROR STOP ERR_UNHANDLED endif end subroutine function create_closed_polygon2d_t(N_pts, X_pts, Y_pts) result(poly) type(closed_polygon2d_t) :: poly integer, intent(in) :: N_pts real(FP), dimension(N_pts), intent(in) :: X_pts, Y_pts integer :: N_extended real(FP), dimension(:), allocatable :: X_extended, Y_extended if ((abs(X_pts(1) - X_pts(N_pts)) < FP_EPS) .and. (abs(Y_pts(1) - Y_pts(N_pts)) < FP_EPS)) then ! If the first and last point are the same call poly%create_polygon2d(N_pts, X_pts, Y_pts) else ! If the first and last points are different N_extended = N_pts+1 allocate(X_extended(N_extended)) allocate(Y_extended(N_extended)) X_extended(1:N_pts) = X_pts Y_extended(1:N_pts) = Y_pts ! Add an extra point, to have the first and last points be equal X_extended(N_pts+1) = X_pts(1) Y_extended(N_pts+1) = Y_pts(1) call poly%create_polygon2d(N_extended, X_extended, Y_extended) deallocate(X_extended, Y_extended) endif end function create_closed_polygon2d_t function create_limiting_polygon2d_t(N_pts, X_pts, Y_pts, local_min, local_max, use_local_min, use_local_max) result(poly) type(limiting_polygon2d_t) :: poly integer, intent(in) :: N_pts real(FP), dimension(N_pts), intent(in) :: X_pts, Y_pts real(FP), intent(in) :: local_min, local_max logical, intent(in) :: use_local_min, use_local_max call poly%create_polygon2d(N_pts, X_pts, Y_pts) poly%local_min = local_min poly%local_max = local_max poly%use_local_min = use_local_min poly%use_local_max = use_local_max end function create_limiting_polygon2d_t !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Finalisation routines MUST be type-bound - so have to implement for each dynamic type subroutine destroy_polygon(self) class(polygon2d_t), intent(inout) :: self if(allocated(self%pts)) deallocate(self%pts) if(allocated(self%p2p)) deallocate(self%p2p) if(allocated(self%p2p_mag)) deallocate(self%p2p_mag) end subroutine subroutine finalize_closed_polygon2d_t(self) type(closed_polygon2d_t), intent(inout) :: self if(allocated(self%pts)) deallocate(self%pts) if(allocated(self%p2p)) deallocate(self%p2p) if(allocated(self%p2p_mag)) deallocate(self%p2p_mag) end subroutine subroutine finalize_limiting_polygon2d_t(self) type(limiting_polygon2d_t), intent(inout) :: self if(allocated(self%pts)) deallocate(self%pts) if(allocated(self%p2p)) deallocate(self%p2p) if(allocated(self%p2p_mag)) deallocate(self%p2p_mag) end subroutine subroutine write_polygon_points(self, output_folder, polygon_name, invert_polygon) ! writes polygon points to trunk class(polygon2d_t), intent(in) :: self character(len=20), parameter :: srtn_name = 'write_polygon_points' character(len=*), intent(in) :: polygon_name, output_folder logical, intent(in) :: invert_polygon integer :: iostat, ierr integer :: fid, iddim_Npts, id_x, id_y if (is_master()) then ! write meta information of penalisation write(get_stdout(), *) & "Writing "//output_folder//"/"//trim(polygon_name)//'_points.nc' iostat = nf90_create(output_folder//"/"//trim(polygon_name)//'_points.nc', NF90_CLOBBER, fid) if (iostat /= 0) then write(get_stderr(), *) & 'error(write_polygon_points):: error creating file', iostat ERROR STOP ERR_UNHANDLED endif ! Write the name of the polygon for reference call nf90_handle_err(nf90_put_att(fid, NF90_GLOBAL, 'identifier', polygon_name), srtn_name) ! Write the logical 'invert' as an integer equivalent if (invert_polygon) then call nf90_handle_err(nf90_put_att(fid, NF90_GLOBAL, 'invert', 1), srtn_name) else call nf90_handle_err(nf90_put_att(fid, NF90_GLOBAL, 'invert', 0), srtn_name) endif ! define dimensions call nf90_handle_err(nf90_def_dim(fid, 'dim_Npts', self%N_pts, iddim_Npts), srtn_name) ! define variables call nf90_handle_err(nf90_def_var(fid, 'x', NF90_DOUBLE, iddim_Npts, id_x), srtn_name) call nf90_handle_err(nf90_def_var(fid, 'y', NF90_DOUBLE, iddim_Npts, id_y), srtn_name) ! end of definition call nf90_handle_err(nf90_enddef(fid) , srtn_name) ! write variables call nf90_handle_err(nf90_put_var(fid, id_x, self%pts(:, x)), srtn_name) call nf90_handle_err(nf90_put_var(fid, id_y, self%pts(:, y)), srtn_name) ! close file call nf90_handle_err(nf90_close(fid), srtn_name) endif call MPI_Barrier(get_communicator(), ierr) end subroutine write_polygon_points pure real(FP) function test_side(P0x, P0y, P1x, P1y, P2x, P2y) ! Tests whether a point is to the left or right of an infinite line ! ! From http://geomalgorithms.com/ ! ! Returns: >0 for P2 left of the line through P0 and P1 ! =0 for P2 on the line ! <0 for P2 right of the line ! real(FP), intent(in) :: P0x, P0y, P1x, P1y, P2x, P2y test_side = (((P1x - P0x) * (P2y - P0y)) - ((P2x - P0x) * (P1y - P0y))) end function test_side pure logical function pt_inside(self, Px, Py) ! Determine whether a point (Px, Py) is within a closed polygon2d self ! Uses the 'Winding Algorithm' ! ! From http://geomalgorithms.com/a03-_inclusion.html#wn_PnPoly() ! class(closed_polygon2d_t), intent(in) :: self real(FP), intent(in) :: Px, Py integer :: winding, i ! winding is the 'winding number' winding = 0 ! loop through all edges of the polygon2d do i = 1, self%N_pts-1 ! edge from V[i] to V[i+1] if (self%pts(i, 2) <= Py) then ! start y <= Py if (self%pts(i+1, 2) > Py) then ! an upward crossing if (test_side(P0x = self%pts(i, x), & P0y = self%pts(i, y), & P1x = self%pts(i+1, x), & P1y = self%pts(i+1, y), & P2x = Px, & P2y = Py) > 0) then ! P left of edge winding = winding +1 ! have a valid up intersect end if end if else ! start y > Py (no test needed) if (self%pts(i+1, 2) <= Py) then ! a downward crossing if (test_side(P0x = self%pts(i, x), & P0y = self%pts(i, y), & P1x = self%pts(i+1, x), & P1y = self%pts(i+1, y), & P2x = Px, & P2y = Py) < 0) then ! P right of edge winding = winding - 1 ! have a valid down intersect end if end if end if end do if (winding == 0) then pt_inside = .false. else pt_inside = .true. endif end function pt_inside real(FP) function dist_to_polypoint(self, Px, Py) ! Determines shortest distance from a point (Px, Py) to a vertex of the polygon2d class(polygon2d_t), intent(in) :: self real(FP), intent(in) :: Px, Py real(FP), dimension(self%N_pts, twodim) :: vector_distance ! Distance (Px - Pi[x]), (Py - Pi[y]) where Pi is the 'i'th polygon2d point real(FP), dimension(self%N_pts) :: vector_distance_mag vector_distance(:, x) = Px - self%pts(:, x) vector_distance(:, y) = Py - self%pts(:, y) call magnitude(vector_distance, self%N_pts, twodim, vector_distance_mag) dist_to_polypoint = minval(vector_distance_mag) end function dist_to_polypoint subroutine dist_to_all_polyedges(self, Px, Py, dist) ! Uses vector projection to determine distance from a point (Px, Py) to all the polygon2d edges ! See https://en.wikipedia.org/wiki/Vector_projection class(polygon2d_t), intent(in) :: self real(FP), intent(in) :: Px, Py real(FP), dimension(self%N_pts), intent(out) :: dist ! Note: for a closed_polygon2d_t number of vertices = number of edges ! for an Openpolygon2d number of vertices = number of edges + 1 real(FP), dimension(self%N_pts, twodim) :: vector_distance ! Distance (Px - Pi[x]), (Py - Pi[y]) where Pi is the 'i'th polygon2d point real(FP), dimension(self%N_pts - 1) :: vector_projection, a2_mag real(FP), dimension(self%N_pts - 1, twodim) :: a1, a2 real(FP), dimension(self%N_pts) :: vector_distance_mag integer :: i, j ! Vector of distances from the point to the vertices vector_distance(:, x) = Px - self%pts(:, x) vector_distance(:, y) = Py - self%pts(:, y) call magnitude(vector_distance, self%N_pts, twodim, vector_distance_mag) ! (inner product) vector_projection = sum(vector_distance(1:self%N_pts-1, :) * self%p2p(:, :), dim=2)/self%p2p_mag ! Find projection of a onto b a1(:, x) = self%p2p(:, x) * vector_projection / self%p2p_mag a1(:, y) = self%p2p(:, y) * vector_projection / self%p2p_mag ! Find the rejection of a from b a2(:, :) = vector_distance(1:self%N_pts-1, :) - a1(:, :) call magnitude(a2, self%N_pts -1, twodim, a2_mag) do i = 1, self%N_pts - 1 if (vector_projection(i) < 0) then ! If dot product is less than zero, use the P1 radial distance j = mod(i, self%N_pts) dist(i) = vector_distance_mag(j) else if (vector_projection(i) > self%p2p_mag(i)) then ! If dot product is greater than the length of P1_P2, use the P2 radial distance dist(i) = vector_distance_mag(i) else ! Else use the distance to the line dist(i) = a2_mag(i) endif end do ! The last point must be treated exceptionally, since it doesn't have a vector to project onto -- so just use radial distance dist(self%N_pts) = vector_distance_mag(self%N_pts) end subroutine dist_to_all_polyedges real(FP) function dist_to_polyedge(self, Px, Py) ! Uses vector projection to determine distance from a point (Px, Py) to the closest polygon2d edge class(polygon2d_t), intent(in) :: self real(FP), intent(in) :: Px, Py real(FP), dimension(self%N_pts) :: dist call dist_to_all_polyedges(self, Px, Py, dist) dist_to_polyedge = minval(dist) end function dist_to_polyedge integer function closest_polyedge(self, Px, Py) ! Uses vector projection to determine the closest polygon2d edge class(polygon2d_t), intent(in) :: self real(FP), intent(in) :: Px, Py real(FP), dimension(self%N_pts) :: dist call dist_to_all_polyedges(self, Px, Py, dist) closest_polyedge = minloc(dist, dim=1) end function closest_polyedge real(FP) function signed_area(self) class(polygon2d_t), intent(in) :: self integer :: n_points, index real(FP) :: running_sum ! http://mathworld.wolfram.com/PolygonArea.html ! Returns the signed area of a non-intersecting polygon ! If the sign is positive, points are counterclockwise ! Otherwise, points are clockwise n_points = self%N_pts running_sum = 0 do index = 1, self%N_pts running_sum = running_sum & + self%pts(index, x) * self%pts(mod(index, n_points)+1, y) & - self%pts(index, y) * self%pts(mod(index, n_points)+1, x) enddo signed_area = running_sum/2.0_FP end function signed_area character(len=string_buffer_length) function string_closed_polygon2d_t(self) class(closed_polygon2d_t), intent(in) :: self character(len=string_buffer_length) :: temporary_string write(temporary_string, 102) self%N_pts 102 FORMAT('Closed 2D Polygon with ', I3, ' points') string_closed_polygon2d_t = temporary_string end function string_closed_polygon2d_t character(len=string_buffer_length) function string_limiting_polygon2d_t(self) class(limiting_polygon2d_t), intent(in) :: self character(len=string_buffer_length) :: temporary_string if (.not.(self%use_local_min) .and. .not.(self%use_local_max)) then write(temporary_string, 104) self%N_pts 104 FORMAT('Limiting, closed 2D polygon with ', I3, ' points') else if (.not.(self%use_local_max)) then write(temporary_string, 105) self%N_pts, self%local_min 105 FORMAT('Limiting, closed 2D polygon with ', I3, ' points, local_min: ', ES10.3E2) else if (.not.(self%use_local_min)) then write(temporary_string, 106) self%N_pts, self%local_max 106 FORMAT('Limiting, closed 2D polygon with ', I3, ' points, local_max: ', ES10.3E2) else write(temporary_string, 107) self%N_pts, self%local_min, self%local_max 107 FORMAT('Limiting, closed 2D polygon with ', I3, ' points, local_min: ', ES10.3E2, ', local_max: ', ES10.3E2) endif string_limiting_polygon2d_t = temporary_string end function string_limiting_polygon2d_t end module polygon_m