submodule (polygon_in_mesh_m) polygon_in_mesh_build_s !! Routines related with building penalisation implicit none contains module subroutine init_polygon_in_mesh(self, equi, mesh, inds, dir_counterclockwise, dbgout) class(polygon_in_mesh_t), intent(inout) :: self type(mesh_cart_t), intent(in) :: mesh class(equilibrium_t) :: equi integer, dimension(:), intent(in) :: inds logical, intent(in), optional :: dir_counterclockwise integer, intent(in), optional :: dbgout logical :: print_warning integer :: outi, nz, k, kp, nsingles, it integer, allocatable, dimension(:) :: inds_wrk, nadj, kends, rnk, np_segs, inds_singles integer, allocatable, dimension(:,:) :: adjacency ! set output level outi = 0 if (present(dbgout)) then if (is_rank_info_writer) then outi = dbgout endif if (dbgout >= 3) then outi = dbgout ! every rank writes endif endif if (present(dir_counterclockwise)) then self%is_cc_wise = dir_counterclockwise else self%is_cc_wise = .true. endif if (outi > 0) then write(get_stdout(),*)'' write(get_stdout(),*)'Building a polygon' endif ! Get unique list and fill inds nz = size(inds) allocate(inds_wrk, source=inds) call unique_tuples(ndim=1,nz=nz, arr=inds_wrk) self%np = nz allocate(self%inds(self%np)) self%inds(1:self%np) = inds_wrk(1:self%np) deallocate(inds_wrk) do it = 1, 2 ! For it = 1 the adjacency is built for the given dataset ! and number of single points (nsingles) (no adjacency) is determined ! Only if nsingles > 0, the single points are removed from the dataset ! and the adjacency is built again (it=2) for the clean dataset ! Adjacency allocate(nadj(self%np)) allocate(adjacency(self%np,2)) call self%build_adjacency(mesh, nadj, adjacency) ! Get number of single points having no adjacency nsingles = 0 do k = 1, self%np if (nadj(k) == 0) then nsingles = nsingles + 1 endif enddo if ((it == 2) .and. (nsingles > 0)) then call handle_error('Removal of single points not succesful', & GRILLIX_ERR_POLYGON, __LINE__, __FILE__) endif if (nsingles == 0) then ! Adjacency is ok, exit it loop exit else ! Remove single points and rerun adjacency allocate(inds_wrk(self%np-nsingles)) kp = 0 do k = 1, self%np if (nadj(k) > 0) then kp = kp + 1 inds_wrk(kp) = self%inds(k) endif enddo deallocate(self%inds) self%np = self%np - nsingles allocate(self%inds, source=inds_wrk) deallocate(adjacency) deallocate(nadj) endif enddo ! Catch resulting empty list if (self%np <= 0) then self%is_closed = .false. self%nsegs = 0 allocate(self%ki_seg(self%nsegs)) if (outi > 0) then call self%display() endif return endif ! Segmentation call self%build_segmentation(nadj, kends) ! Ranking allocate(rnk(self%np)) allocate(np_segs(self%nsegs)) call self%build_ranking(nadj, adjacency, kends, rnk, np_segs) deallocate(adjacency) deallocate(nadj) deallocate(kends) ! Execute sorting call self%execute_sort(rnk, np_segs) deallocate(np_segs) deallocate(rnk) ! Order segments counterclockwise call self%direction_order(equi, mesh) ! Print warning message if it could not be determined doubtlessly ! if polygon or any segment is ordered counter-clockwise print_warning = .false. if (self%warning_self_intersection) then if (equi%is_axisymmetric()) then if (is_rank_info_writer) then print_warning = .true. endif else print_warning = .true. endif endif if (print_warning) then call handle_error('Could not doubtlessly determine if polygon & or one of its segment is ordered & clock or counterclockwise. You may continue, & but should check direction of polygon manually', & GRILLIX_WRN_POLYGON, __LINE__, __FILE__) endif if (outi > 0) then call self%display() endif end subroutine module subroutine build_adjacency(self, mesh, nadj, adjacency) class(polygon_in_mesh_t), intent(in) :: self type(mesh_cart_t), intent(in) :: mesh integer, dimension(self%np), intent(out) :: nadj integer, dimension(self%np, 2), intent(out) :: adjacency integer :: k, l, ic, jc, kcand, lcand, icand, jcand, dsti, dstj adjacency = 0 do k = 1, self%np l = self%inds(k) ic = mesh%get_cart_i(l) jc = mesh%get_cart_j(l) nadj(k) = 0 ! Firstly look for horizontal or vertical neighbours as candidates do kcand = 1, self%np lcand = self%inds(kcand) icand = mesh%get_cart_i(lcand) jcand = mesh%get_cart_j(lcand) dsti = abs(ic - icand) dstj = abs(jc - jcand) if ( ((dsti==1) .and. (dstj==0)) .or. & ((dsti==0) .and. (dstj==1)) ) then nadj(k) = nadj(k) + 1 if (nadj(k) > 2) then call handle_error('Polygon cannot be build uniquely', & GRILLIX_ERR_POLYGON, __LINE__, __FILE__, & error_info_t('for (mesh,polygon) point:', [l,k]) ) endif adjacency(k, nadj(k)) = kcand endif enddo ! If less than 2 neighbours where found, search for diagonal neighbours as candidates if (nadj(k) < 2) then do kcand = 1, self%np lcand = self%inds(kcand) icand = mesh%get_cart_i(lcand) jcand = mesh%get_cart_j(lcand) dsti = abs(ic - icand) dstj = abs(jc - jcand) if ((dsti==1) .and. (dstj==1)) then nadj(k) = nadj(k) + 1 if (nadj(k) > 2) then call handle_error('Polygon cannot be build uniquely', & GRILLIX_ERR_POLYGON, __LINE__, __FILE__, & error_info_t('for (mesh,polygon) point:', [l,k]) ) endif adjacency(k, nadj(k)) = kcand endif enddo endif enddo end subroutine module subroutine build_segmentation(self, nadj, kends) class(polygon_in_mesh_t), intent(inout) :: self integer, dimension(self%np), intent(in) :: nadj integer, allocatable, dimension(:), intent(out) :: kends integer :: nends, k, iends ! Check sanity of polygon and count number of ends nends = 0 do k = 1, self%np if (nadj(k) > 2) then ! Note: This is redundant with build_adjacency, but kept here for further safety. call handle_error('Polygon cannot be build uniquely', & GRILLIX_ERR_POLYGON, __LINE__, __FILE__, & error_info_t('for polygon point:', [k]) ) endif if (nadj(k) < 1) then call handle_error('No adjacent point for polygon available', & GRILLIX_ERR_POLYGON, __LINE__, __FILE__, & error_info_t('for polygon point:', [k]) ) endif if (nadj(k) == 1) then nends = nends + 1 endif enddo if (mod(nends,2) /= 0) then call handle_error('Polygon has odd number of ends', & GRILLIX_ERR_POLYGON, __LINE__, __FILE__) endif ! get number of segments if (nends == 0) then self%is_closed = .true. self%nsegs = 1 else self%is_closed = .false. self%nsegs = nends/2 endif ! Get end points of polygon allocate(kends(nends)) iends = 0 if (.not. self%is_closed) then do k = 1, self%np if (nadj(k) == 1) then iends = iends + 1 kends(iends) = k endif enddo endif end subroutine module subroutine build_ranking(self, nadj, adjacency, kends, rnk, np_segs) class(polygon_in_mesh_t), intent(inout) :: self integer, dimension(self%np), intent(in) :: nadj integer, dimension(self%np, 2), intent(in) :: adjacency integer, dimension(2*self%nsegs), intent(in) :: kends integer, dimension(self%np), intent(out) :: rnk integer, dimension(self%nsegs), intent(out) :: np_segs integer :: nends, k, iseg, iends, ip, kn, kadj integer, dimension(2*self%nsegs) :: kends_copy, wrk_kends, klast_seg nends = 2*self%nsegs kends_copy = kends k = 0 ! This is actually the running index over he whole polygon do iseg = 1, self%nsegs ! Determine starting point of segment if (self%is_closed) then k = k + 1 rnk(k) = 1 else do iends = 1, nends if (kends_copy(iends) /= 0) then k = k + 1 rnk(k) = kends_copy(iends) wrk_kends(2*(iseg-1)+1) = kends_copy(iends) kends_copy(iends) = 0 exit endif enddo endif ! Choose any adjacent point as second point of segment ! Direction will be determined later k = k + 1 rnk(k) = adjacency(rnk(k-1), 1) ! Put subsequently next adjacent point into ranking do ip = 2, self%np-1 do kn = 1, nadj(rnk(k)) kadj = adjacency(rnk(k), kn) if (kadj /= rnk(k-1)) then k = k + 1 rnk(k) = kadj exit endif enddo enddo ! Running index of last point of segment klast_seg(iseg) = k ! Update kends wrk_kends(2*(iseg-1)+2) = rnk(k) do iends = 1, nends if ( wrk_kends(2*(iseg-1)+2) == kends_copy(iends) ) then kends_copy(iends) = 0 exit endif enddo enddo ! Number of points per segment np_segs(1) = klast_seg(1) do iseg = 2, self%nsegs np_segs(iseg) = klast_seg(iseg) - klast_seg(iseg-1) enddo end subroutine module subroutine execute_sort(self, rnk, np_segs) class(polygon_in_mesh_t), intent(inout) :: self integer, intent(in), dimension(self%np) :: rnk integer, dimension(self%nsegs), intent(in) :: np_segs integer :: k, iseg integer, dimension(self%np) :: wrk do k = 1, self%np wrk(k) = self%inds(rnk(k)) self%inds(rnk(k)) = 0 enddo ! Check if all points have been accessed if (maxval(self%inds) > 0) then call handle_error('Not all polygon points were accessed', & GRILLIX_ERR_POLYGON, __LINE__, __FILE__) endif self%inds = wrk ! Build ki_seg allocate(self%ki_seg(self%nsegs+1)) self%ki_seg(1) = 1 do iseg = 1, self%nsegs self%ki_seg(iseg+1) = self%ki_seg(iseg) + np_segs(iseg) enddo end subroutine module subroutine direction_order(self, equi, mesh) class(polygon_in_mesh_t), intent(inout) :: self class(equilibrium_t) :: equi type(mesh_cart_t), intent(in) :: mesh integer :: iseg, np_seg, k, ip, ks, kf real(GP), allocatable, dimension(:) :: xp, yp integer, dimension(self%np) :: wrk logical :: is_counter_clockwise, self_intersection_occured do iseg = 1, self%nsegs np_seg = self%ki_seg(iseg+1) - self%ki_seg(iseg) if (self%is_closed) then allocate(xp(np_seg)) allocate(yp(np_seg)) else allocate(xp(np_seg+1)) allocate(yp(np_seg+1)) endif ip = 0 do k = self%ki_seg(iseg), self%ki_seg(iseg+1) - 1 ip = ip + 1 xp(ip) = mesh%get_x(self%inds(k)) yp(ip) = mesh%get_y(self%inds(k)) enddo ! Add magnetic axis to open segment if (.not. self%is_closed) then ip = ip + 1 select type(equi) type is(slab_t) ! For slab use center of mesh xp(ip) = 0.5_GP * (equi%xmin + equi%xmax) yp(ip) = 0.5_GP * (equi%ymin + equi%ymax) class default xp(ip) = equi%x0 yp(ip) = equi%y0 end select endif call determine_direction(xp, yp, is_counter_clockwise, self_intersection_occured) if (self_intersection_occured) then self%warning_self_intersection = .true. endif if (is_counter_clockwise .neqv. self%is_cc_wise) then ! Flip indices of segment ks = self%ki_seg(iseg) kf = self%ki_seg(iseg+1) - 1 wrk(ks:kf) = self%inds(kf:ks:-1) self%inds(ks:kf) = wrk(ks:kf) endif deallocate(yp) deallocate(xp) enddo end subroutine module subroutine destructor(self) type(polygon_in_mesh_t), intent(inout) :: self if (allocated(self%inds)) deallocate(self%inds) if (allocated(self%ki_seg)) deallocate(self%ki_seg) end subroutine subroutine determine_direction(xp, yp, is_counter_clockwise, self_intersection_occured) !! Determines if closed (!) polygon prescribed via points xp, yp is ordered clockwise !! Method based on computing signed area of polygon !! (see https://mathworld.wolfram.com/PolygonArea.html) !! The algorithm assumes that the polygon does not have any self-intersections, !! This is checked and returned in the additional variable self_intersection_occured. real(GP), dimension(:), intent(in) :: xp !! x-coordinates of polygon real(GP), dimension(:), intent(in) :: yp !! y-coordinates of polygon logical, intent(out) :: is_counter_clockwise !! If true (false) and self_intersection_occured=.false.: !! The polygon is counter-clockwise (clockwise) !! !! If self_intersection_occured=.true.: !! There is still a good chance that the result is usable !! for the given purpose, but the user should check manually. logical, intent(out) :: self_intersection_occured !! If true, the polygon is self-intersecting, !! direction of polygon cannot doubtlessly be determined integer :: np, ip real(GP) :: area np = size(xp) if (size(xp) /= size(yp)) then call handle_error('Sizes of polygon coordinates do not match', & GRILLIX_ERR_POLYGON, __LINE__, __FILE__) endif ! Check for self intersection of polygon self_intersection_occured = is_self_intersecting(xp, yp) area = 0.0_GP do ip = 1, np - 1 area = area + xp(ip) * yp(ip+1) - xp(ip+1) * yp(ip) enddo area = area + xp(np) * yp(1) - xp(1) * yp(np) area = 0.5_GP * area if (area < 0.0_GP) then is_counter_clockwise = .false. else is_counter_clockwise = .true. endif end subroutine logical function is_self_intersecting(xp, yp) !! Returns true if closed (!) polygon prescribed via points xp, yp has self-intersection real(GP), dimension(:), intent(in) :: xp !! x-coordinates of polygon real(GP), dimension(:), intent(in) :: yp !! y-coordinates of polygon integer :: np, ip, jp, info real(GP) :: xi, yi, ta, tb real(GP), dimension(2) :: pa, qa, pb, qb np = size(xp) if (size(xp) /= size(yp)) then call handle_error('Sizes of polygon coordinates do not match', & GRILLIX_ERR_POLYGON, __LINE__, __FILE__) endif is_self_intersecting = .false. do ip = 1, np pa(1) = xp(ip) pa(2) = yp(ip) if (ip < np) then qa(1) = xp(ip+1) qa(2) = yp(ip+1) else qa(1) = xp(1) qa(2) = yp(1) endif do jp = 1, np pb(1) = xp(jp) pb(2) = yp(jp) if (jp < np) then qb(1) = xp(jp+1) qb(2) = yp(jp+1) else qb(1) = xp(1) qb(2) = yp(1) endif call intersection_lines(pa, qa, pb, qb, xi, yi, info, ta, tb) if (info == 0) then if ( (ta > GP_EPS) .and. (ta < (1.0_GP-GP_EPS)) .and. & (tb > GP_EPS) .and. (tb < (1.0_GP-GP_EPS)) ) then is_self_intersecting = .true. return endif endif enddo enddo end function end submodule