diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 99466dc20..ea42073a0 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -282,10 +282,10 @@ subroutine encounter_check_reduce_broadphase(n, nenc, index1, index2, lencounter logical, dimension(:), allocatable, intent(inout) :: lencounter !! Logical flag indicating which of the pairs identified in the broad phase was selected in the narrow phase logical, dimension(:), allocatable, intent(inout) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals - integer(I4B) :: i, i0, j, k + integer(I4B) :: i, i0, j, k, nenci, klo, khi integer(I4B), dimension(:), allocatable :: itmp - logical, dimension(n) :: lfresh integer(I4B), dimension(:), allocatable :: ind + integer(I4B), dimension(n) :: ibeg, iend logical, dimension(:), allocatable :: ltmp nenc = count(lencounter(:)) ! Count the true number of encounters @@ -307,25 +307,44 @@ subroutine encounter_check_reduce_broadphase(n, nenc, index1, index2, lencounter lencounter(:) = .true. call util_sort(index1, ind) - lfresh(:) = .true. - - ! Remove duplicates - do k = 1, nenc - i = index1(ind(k)) - j = index2(ind(k)) - if (k == 1) then - lfresh(j) = .false. - else - i0 = index1(ind(k-1)) - if (i /= i0) lfresh(:) = .true. - if (.not.lfresh(j)) lencounter(ind(k)) = .false. - lfresh(j) = .false. + call util_sort_rearrange(index1, ind, nenc) + call util_sort_rearrange(index2, ind, nenc) + call util_sort_rearrange(lvdotr, ind, nenc) + + ! Get the bounds on the bodies in the first index + ibeg(:) = n + iend(:) = 1 + i0 = index1(1) + ibeg(i0) = 1 + do k = 2, nenc + i = index1(k) + if (i /= i0) then + iend(i0) = k - 1 + ibeg(i) = k + i0 = i end if + if (k == nenc) iend(i) = k + end do + + ! Sort on the second index and remove duplicates + if (allocated(itmp)) deallocate(itmp) + allocate(itmp, source=index2) + do concurrent(i = 1:n, iend(i) - ibeg(i) > 0) + klo = ibeg(i) + khi = iend(i) + nenci = khi - klo + 1 + if (allocated(ind)) deallocate(ind) + call util_sort(index2(klo:khi), ind) + index2(klo:khi) = itmp(klo - 1 + ind(:)) + do j = klo + 1, khi + if (index2(j) == index2(j - 1)) lencounter(j) = .false. + end do end do if (count(lencounter(:)) == nenc) return nenc = count(lencounter(:)) ! Count the true number of encounters + if (allocated(itmp)) deallocate(itmp) allocate(itmp(nenc)) itmp(:) = pack(index1(:), lencounter(:)) call move_alloc(itmp, index1) @@ -875,6 +894,7 @@ pure subroutine encounter_check_make_ragged_list(lencounteri, ind_arr, nenc,inde integer(I4B), dimension(:), intent(in) :: ind_arr integer(I4B), intent(out) :: nenc integer(I4B), dimension(:), allocatable, intent(out) :: index2 + ! Internals nenc = count(lencounteri(:)) if (nenc > 0) then