module polygon_in_mesh_m !! Definition of polygon running within mesh use precision_grillix_m, only : GP, GP_EPS use error_handling_grillix_m, only: handle_error, error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_POLYGON, GRILLIX_WRN_POLYGON use equilibrium_m, only : equilibrium_t use slab_equilibrium_m, only : slab_t use mesh_cart_m, only : mesh_cart_t use list_operations_m, only : unique_tuples use euclidean_geo_m, only : intersection_lines use screen_io_m, only : get_stdout use parallelisation_setup_m, only : is_rank_info_writer implicit none type, public :: polygon_in_mesh_t !! Datatype for polygon running on mesh points !! This can either be a single closed polygon, !! or any positive number of open polygon segments (but not both simultaneously). !! The polgon must be sufficiently well defined, !! such that only up to two adjacent points can be uniquely determined for each point. integer, private :: np = 0 !! Total number of polygon points integer, private, allocatable, dimension(:) :: inds !! Indices of polygon on mesh logical, private :: is_closed !! True if polygon is closed (periodic) integer, private :: nsegs = 0 !! Number of segments integer, allocatable, dimension(:) :: ki_seg !! Index of initial point of polygon segment logical, private :: is_cc_wise !! If true polygon is ordered counter-clocwise, otherwise clock-wise logical, private :: warning_self_intersection = .false. !! If true indicates that during build phase some polygon(-segment) !! had self_intersections. Indicates that polygon(-segment) might not be ordered clockwise contains procedure, public :: init => init_polygon_in_mesh procedure, public :: is_periodic procedure, public :: get_np_total procedure, public :: get_nsegs procedure, public :: get_nps procedure, public :: get_ind_on_mesh procedure, public :: display => display_polygon_in_mesh procedure, public :: write_netcdf => write_netcdf_polygon_in_mesh procedure, public :: read_netcdf => read_netcdf_polygon_in_mesh procedure, private :: build_adjacency procedure, private :: build_segmentation procedure, private :: build_ranking procedure, private :: execute_sort procedure, private :: direction_order final :: destructor end type interface module subroutine init_polygon_in_mesh(self, equi, mesh, inds, dir_counterclockwise, dbgout) !! Initialises polygon in mesh class(polygon_in_mesh_t), intent(inout) :: self !! Instance of type type(mesh_cart_t), intent(in) :: mesh !! Mesh class(equilibrium_t) :: equi !! Equilibrium integer, dimension(:), intent(in) :: inds !! Indices (on mesh) of polygon logical, intent(in), optional :: dir_counterclockwise !! If true (=default) the direction of the polygons is set counter-clockwise !! otherwise it is clockwise integer, intent(in), optional :: dbgout !! Debug output level end subroutine module subroutine build_adjacency(self, mesh, nadj, adjacency) !! Returns adjacency for given point list !! Note: The polygon must be sufficiently uniquely defined, !! i.e. there must be not more than 2 candidates as adjacent points for any point class(polygon_in_mesh_t), intent(in) :: self !! Instance of type type(mesh_cart_t), intent(in) :: mesh !! Mesh integer, dimension(self%np), intent(out) :: nadj !! Number of adjacent points for given point integer, dimension(self%np, 2), intent(out) :: adjacency !! Adjacency information, i.e. indices (of polygon) of up to two adjacent points !! If value is zero, no adjacent point could be determined end subroutine module subroutine build_segmentation(self, nadj, kends) !! Determines segmentation of polygon, i.e. if !! polygon is closed, number of segments, end points of segments class(polygon_in_mesh_t), intent(inout) :: self !! Instance of type integer, dimension(self%np), intent(in) :: nadj !! Number of adjacent points for given point (from build_adjacency) integer, allocatable, dimension(:), intent(out) :: kends !! Polygon indices of end points of polygon segments end subroutine module subroutine build_ranking(self, nadj, adjacency, kends, rnk, np_segs) !! Builds ranking of polygon indices for sorting class(polygon_in_mesh_t), intent(inout) :: self !! Instance of type integer, dimension(self%np), intent(in) :: nadj !! Number of adjacent points for given point (from build_adjacency) integer, dimension(self%np, 2), intent(in) :: adjacency !! Adjacency information (from build_adjacency) integer, dimension(2*self%nsegs), intent(in) :: kends !! Polygon indices of end points of polygon segments. !! On succesfull return, will be sorted according to segments. integer, dimension(self%np), intent(out) :: rnk !! Ranking array, i.e. indices of polygon to sort integer, dimension(self%nsegs), intent(out) :: np_segs !! Number of points of segment end subroutine module subroutine execute_sort(self, rnk, np_segs) !! Executes sorting of indices based on given ranking class(polygon_in_mesh_t), intent(inout) :: self !! Instance of type integer, intent(in), dimension(self%np) :: rnk !! Ranking array (from build_ranking) integer, dimension(self%nsegs), intent(in) :: np_segs !! Number of points of segment (from build_ranking) end subroutine module subroutine direction_order(self, equi, mesh) !! Sorts polgon in co- or counter-clockwise order, depending on is_cc_wise !! For open segments, direction of each segment is defined w.r.t. magnetic axis class(polygon_in_mesh_t), intent(inout) :: self !! Instance of type class(equilibrium_t) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh !! Mesh end subroutine module subroutine write_netcdf_polygon_in_mesh(self, fgid) !! Writes polygon_in_mesh to netcdf file class(polygon_in_mesh_t), intent(in) :: self !! Instance of the type integer, intent(in) :: fgid !! Netcdf file or group id end subroutine module subroutine read_netcdf_polygon_in_mesh(self, fgid) !! Reads polygon_in_mesh from netcdf file class(polygon_in_mesh_t), intent(inout) :: self !! Instance of the type integer, intent(in) :: fgid !! Netcdf file or group id end subroutine module subroutine destructor(self) !! Frees memory associated with polygon_in_mesh_t type(polygon_in_mesh_t), intent(inout) :: self end subroutine end interface contains pure function is_periodic(self) result(res) !! Returns tru if polygon is closed (periodic) class(polygon_in_mesh_t), intent(in) :: self !! Instance of the type logical :: res res = self%is_closed end function pure function get_np_total(self) result(res) !! Returns total number of polygon points class(polygon_in_mesh_t), intent(in) :: self !! Instance of the type integer :: res res = self%np end function pure function get_nsegs(self) result(res) !! Returns number of segements class(polygon_in_mesh_t), intent(in) :: self !! Instance of the type integer :: res res = self%nsegs end function pure function get_nps(self, iseg) result(res) !! Returns number of points of segment class(polygon_in_mesh_t), intent(in) :: self !! Instance of the type integer, intent(in), optional :: iseg !! Number of segment (default 1) integer :: res integer :: iseg_actual if (present(iseg)) then iseg_actual = iseg else iseg_actual = 1 endif res = self%ki_seg(iseg_actual+1) - self%ki_seg(iseg_actual) end function pure function get_ind_on_mesh(self, ks, iseg) result(res) !! Returns mesh index of point of segement class(polygon_in_mesh_t), intent(in) :: self !! Instance of the type integer, intent(in) :: ks !! Index on segment !! For closed polygon this is treated periodically !! For open segments zero is returned integer, intent(in), optional :: iseg !! Number of segment (default 1) integer :: res integer :: iseg_actual, kg, kseg_last if (self%nsegs == 0) then ! Catch case of empty polygon res = 0 return endif if (present(iseg)) then iseg_actual = iseg else iseg_actual = 1 endif kg = self%ki_seg(iseg_actual) - 1 + ks if ( (kg >= self%ki_seg(iseg_actual)) .and. (kg < self%ki_seg(iseg_actual+1)) ) then res = self%inds(kg) return else if (self%is_closed) then ! Apply periodicity kseg_last= self%ki_seg(iseg_actual+1) - 1 kg = modulo(ks, kseg_last) if (kg == 0) then kg = kseg_last endif res = self%inds(kg) return else res = 0 return endif endif end function subroutine display_polygon_in_mesh(self) !! Displays information on polygon_in_mesh class(polygon_in_mesh_t), intent(in) :: self !! Instance of the type if (.not. is_rank_info_writer) then return endif write(get_stdout(),*)'Information on polygon_in_mesh -------------------------' write(get_stdout(),101)self%is_closed, self%nsegs, self%np, self%is_cc_wise, & self%warning_self_intersection 101 FORMAT( 3X,'is_closed = ',L1 /, & 3X,'nsegs = ',I9 /, & 3X,'np = ',I9 /, & 3X,'is_cc_wise = ',L1 /, & 3X,'wrn_slfisect = ',L1 ) write(get_stdout(),*)'-------------------------------------------------------' end subroutine end module