mesh_index_helpers_m.f90 Source File

module mesh_index_helpers_m
    !! Module containing useful routines for handling mesh indices
    use error_handling_grillix_m, only: handle_error, error_info_t
    use status_codes_grillix_m, only : GRILLIX_ERR_OTHER
    use mesh_cart_m, only: mesh_cart_t
    implicit none

contains 

    subroutine find_neighbor_inds(mesh, ind, diagonal, ind_nb)
        !! Scan points in 3x3 vicinity and return their index only
        type(mesh_cart_t), intent(in) :: mesh
        !! Mesh
        integer, intent(in) :: ind
        !! Index of input boundary point
        logical, intent(in) :: diagonal
        !! Whether to include diagonal neighbors
        integer, dimension(:), intent(out) :: ind_nb
        !! Index of neighboring points, must be already allocated
        !! If diagonal neighbors included, must have 8 entries
        !! If they are excluded, must be of size 4
        
        integer :: i, j, k, l_nb

        ! Assert that output array size is congruent with diagonal setting
        if ( ( diagonal .and. (size(ind_nb) /= 8) ) .or. &
            ( .not. diagonal .and. (size(ind_nb) /= 4 ) ) ) then
            call handle_error('ind_nb field is not consistent with input', &
                                GRILLIX_ERR_OTHER, __LINE__, __FILE__, &
                                error_info_t('size(ind_nb) ', r_info=[size(ind_nb)]))
        endif

        ! Initialize index array to default values
        ind_nb = -1

        k = 0
        ! Search in 3x3 square surrounding input index
        do i = -1, 1
            do j = -1, 1
                if ( (i == 0) .and. (j == 0) ) then ! exclude self
                    cycle
                else if ( .not. diagonal .and. (i*j /= 0) ) then ! exclude diagonals
                    cycle
                else
                    k = k + 1
                    l_nb = mesh%get_index_neighbor(i, j, ind)
                    if ( l_nb > 0 ) then ! Neighbor exists
                        ind_nb(k) = mesh%get_index_neighbor(i, j, ind)
                    endif
                endif
            enddo
        enddo

    end subroutine

end module