! See copyright notice in the COPYRIGHT file.
! ****************************************************************************** !
!> This module provides the description of subresolution for elements.
!!
!! Subresolutions provide polynomial respresentations of boundaries within
!! elements. This property depends on the color property. Only colored elements
!! can provide subresolution data. For each color there can be an independent
!! subresolved representation of the according boundaries.
module tem_subres_prop_module
  use mpi
  use env_module, only: pathlen, long_k, newunit

  use treelmesh_module, only: treelmesh_type
  use tem_prophead_module, only: tem_prophead_type
  use tem_property_module, only: tem_property_type, prp_hasPolynomial
  use tem_color_prop_module, only: tem_color_prop_type, colors_per_char

  implicit none

  type elemid_list_type
    integer, allocatable :: id(:)
  end type elemid_list_type

  type tem_subres_prop_type
    !> Pointer to treelmesh_type%global%property
    type(tem_prophead_type),  pointer :: header => null()

    !> Indication which of the colors contain subresolution information.
    !!
    !! The first index has length nChars from the color property.
    !! The second runs through all elements with with subresolution property.
    character, allocatable :: subresolved_colors(:,:)

    !> Pointer to treelmesh_type%property
    type(tem_property_type),  pointer :: property => null()

    !> Number of subresolved elements in each color.
    integer, allocatable :: nElems(:)

    !> List of indices of elements with subresolution for each color.
    type(elemid_list_type), allocatable :: elem(:)

    !> Offset for the subresolved elements on this partition for each color.
    integer(kind=long_k), allocatable :: offset(:)

  end type tem_subres_prop_type


contains

  ! **************************************************************************** !
  !> Load the subresolution property from disk.
  !!
  !! Before this can be done, the coloring information has to have been loaded.
  subroutine tem_subres_prop_load( me, tree, coloring )
    ! --------------------------------------------------------------------------!
    !> Color definitions to load.
    type(tem_subres_prop_type), intent(out) :: me

    !> Tree to build the polynomial subresolution information for
    type(treelmesh_type), intent(in) :: tree

    !> Information on the colors in the mesh.
    type(tem_color_prop_type), intent(in) :: coloring
    ! --------------------------------------------------------------------------!
    integer :: rl
    integer :: fUnit
    character(len=pathLen) :: datafile
    integer :: iColor
    integer :: iProp
    integer :: colChar, colBit
    integer(kind=long_k), allocatable :: long_counts(:)
    integer :: iError
    integer :: iElem
    integer :: ice
    ! --------------------------------------------------------------------------!

    prp_loop: do iprop=1, tree%global%nProperties
      if (tree%global%Property(iprop)%bitpos == prp_hasPolynomial) then
        me%header => tree%global%Property(iprop)
        me%property => tree%property(iprop)

        datafile   = trim(tree%global%dirname)//'subres.ascii'

        allocate(me%subresolved_colors(coloring%nChars, me%property%nElems))
        allocate(me%nElems(coloring%nColors))
        allocate(me%Offset(coloring%nColors))
        me%Offset = 0_long_k

        ! If there are actually subresolved elements on the local process,
        ! read them now.
        if (me%property%nElems > 0) then

          allocate(me%elem(coloring%nColors))
          inquire(iolength=rl) me%subresolved_colors(:,1)
          fUnit = newUnit()
          open( unit = fUnit, file = datafile, action = 'read', &
            &   access = 'stream', form = 'unformatted',        &
            &   status = 'old'                                  )
          read(fUnit, pos=me%property%offset+1) me%subresolved_colors
          close(fUnit)

          do iColor=1,coloring%nColors
            colChar = (iColor-1)/colors_per_char + 1
            colBit = mod(iColor-1, colors_per_char)
            me%nElems(iColor)                                            &
              &  = count( btest(ichar(me%subresolved_colors(ColChar,:)), &
              &                 ColBit)                                  )
            allocate( me%elem(iColor)%id(me%nElems(iColor)) )

            ! Store the element link for each color.
            if (me%nElems(iColor) > 0) then
              ice = 0
              do iElem=1,me%property%nElems
                if ( btest(ichar(me%subresolved_colors(ColChar,iElem)), &
                  &        colBit) ) then
                  ice = ice + 1
                  me%elem(iColor)%id(ice) = me%property%elemid(iElem)
                end if
              end do
            end if

          end do

        else
          ! No subresolved elements on the local partition:
          me%nElems = 0
        end if

        allocate( long_counts(coloring%nColors) )
        long_counts = me%nElems
        call MPI_Exscan( long_counts, me%Offset, coloring%nColors,       &
          &              MPI_INTEGER8, MPI_SUM, tree%global%comm, iError )

        EXIT prp_loop

      end if
    end do prp_loop

  end subroutine tem_subres_prop_load

end module tem_subres_prop_module
