! See copyright notice in the COPYRIGHT file.
! ****************************************************************************** !
!> author: Harald Klimach
!! This module provides informations and operations on
!! the topology of the complete linearized octree.
!! It is completely independent of the actual sparse
!! octree and therefore this module is independent of the Treelmesh_module.
!!
!! Node identification in the octree
!! ---
!! Due to the breadth first numbering of the octree
!! it is straight forward to move in vertical direction through the tree.
!!
!! The parent \(i_p\) of a node with a certain
!! [[treelmesh_type:treeID]] \(i_n\) is given by:
!! \[  i_p = \left\lfloor \frac{i_n - 1}{8} \right\rfloor \]
!! Where the bracket \(\left\lfloor\cdot\right\rfloor\) denotes the integer floor
!! rounding to the next smaller integer number.
!! While the children \(i_{c1..8}\) of a given node \(i_n\) can be obtained by:
!! \[   i_{cj} = 8 \cdot i_n + j \]
!! The numbering scheme begins at the root node with a
!! [[treelmesh_type:treeID]] of 0 for
!! this element, that covers the complete enclosed cubic domain.
!! All further [[treelmesh_type:treeID]]s are then
!! given by the rules stated above.
!! What remains to be defined is the actual spatial placement of the nodes in
!! the tree as elements in the mesh.
!! This is achieved by a space filling curve, that describes the numbering
!! scheme of the 8 children for a given node.
!!
!! Our choice here is the Morton ordering [6] (bibliography of treelm) , that
!! allows an efficient computation of the index from the coordinates with the
!! same orientation on each refinement level.
!! For this ordering the coordinates of the children within their parent element
!! are first noted as a three dimensional tuple, where each index might be 0 or
!! 1.
!! Now this tuple is interpreted as a single number of the basis 2 with three
!! digits.
!! The recursive in depth traversal results then in a concatenation of these
!! numbers visited during the traversal.
!! This concatenated number provides the sorting key for all elements on a given
!! level according to the space filling curve.
!!
module tem_topology_module

  ! include treelm modules
  use env_module,       only: long_k, globalMaxLevels

  implicit none

  private

  public :: tem_ParentOf
  public :: tem_childNumber
  public :: tem_SiblingsOfId
  public :: tem_PathOf
  public :: tem_LevelOf
  public :: tem_FirstIdAtLevel, tem_LastIdAtLevel
  public :: tem_directChildren
  public :: tem_TreeIDComparison
  public :: tem_PathComparison
  public :: tem_path_type
  public :: tem_coordOfID, tem_IDofCoord

  integer(kind=long_k), parameter, public :: firstIdAtLevel(20) =  [&
    &                         1_long_k, & ! level =  1
    &                         9_long_k, & ! level =  2
    &                        73_long_k, & ! level =  3
    &                       585_long_k, & ! level =  4
    &                      4681_long_k, & ! level =  5
    &                     37449_long_k, & ! level =  6
    &                    299593_long_k, & ! level =  7
    &                   2396745_long_k, & ! level =  8
    &                  19173961_long_k, & ! level =  9
    &                 153391689_long_k, & ! level = 10
    &                1227133513_long_k, & ! level = 11
    &                9817068105_long_k, & ! level = 12
    &               78536544841_long_k, & ! level = 13
    &              628292358729_long_k, & ! level = 14
    &             5026338869833_long_k, & ! level = 15
    &            40210710958665_long_k, & ! level = 16
    &           321685687669321_long_k, & ! level = 17
    &          2573485501354569_long_k, & ! level = 18
    &         20587884010836553_long_k, & ! level = 19
    &        164703072086692425_long_k  ] ! level = 20

  !> Paths of elements from the root node to themselves going through the
  !! hierarchy of the tree. Is used to compare to elements
  !!
  !! For all
  !! [[tem_construction_module::identify_local_element]] "local elements",
  !! their actual position in the sparse mesh has to be identified for a given
  !! \ref treelmesh_module::treelmesh_type::treeid "treeID".
  !! Due to the sparsity of the mesh, the position of a certain element in the
  !! total list of elements can not be directly deduced from its
  !! \ref treelmesh_module::treelmesh_type::treeid "treeID".
  !! As already pointed out, we can rely on binary searches in the sorted list
  !! of elements for this task, as long as a comparison between any two elements
  !! is possible.
  !! A comparison of two nodes to decide their relative position in the ordered
  !! list of elements has to take into account the hierarchy of the mesh.
  !! The comparison operator is therefore defined as follows.
  !!
  !! - Build the path from each leaf to the root of the tree, by repeatedly
  !!   computing the parents with the equation used in the
  !!   [tem_ParentOf function] (@ref tem_parentatlevel), until the root with
  !!   \ref treelmesh_module::treelmesh_type::treeid "treeID" \f$0\f$ is reached.
  !! - After the path through the tree is known for both leaves,
  !!   comparison can be done on the greatest common level by a simple integer
  !!   subtraction of the treeIDs in the path on this level.
  !! - The result of this subtraction indicates the relation between the two
  !!   paths, if it is zero they are considered to be equal.
  !!   Otherwise the sign indicates the ordering of the two compared leaves.
  !!
  !! In a tree with leaves on different levels, this results in children to be
  !! considered to be equal to their parents.
  type tem_path_type
    !> Levels counted from 1
    !! This level is 1 higher than the actual level.
    integer :: Level
    !> treeIDs from current leaf to root
    integer(kind=long_k) :: Node(globalMaxLevels+1)
  end type

  ! interface tem_LevelOf
  !   module procedure tem_levelOf_singleTreeID
  !   module procedure tem_levelOf_arrayOfTreeIDs
  ! end interface

  interface tem_ParentOf
    module procedure tem_directParent
    module procedure tem_ParentAtLevel
  end interface

  contains

! ****************************************************************************** !
  !> This function delivers the refinement level of a TreeID
  !!
  !! The refinement level of a
  !! \ref treelmesh_module::treelmesh_type::treeid "treeID" is found by
  !! repeatedly computing the parents with the equation used in the
  !! [tem_ParentOf function] (@ref tem_parentatlevel) and
  !! counting the iterations, until the universal root is reached.
  !!
  ! elemental function tem_LevelOf_singleTreeID( TreeID ) result(res)
  elemental function tem_LevelOf( TreeID ) result(res)
    ! ---------------------------------------------------------------------------
    !> ID in the complete tree to look up
    integer(kind=long_k), intent(in) :: TreeID
    !> Level of the given TreeID
    integer :: res
    ! ---------------------------------------------------------------------------
    integer(kind=long_k) :: tElem
    integer :: level
    ! ---------------------------------------------------------------------------

    level = 0
    tElem = TreeID
    do while (tElem /= 0)
      ! bit operation is more efficient than division
      ! tElem = (tElem-1) / 8
      tElem = ishft( (tElem-1), -3 )
      level = level + 1
    end do
    res = level

  ! end function tem_LevelOf_singleTreeID
  end function tem_LevelOf
! ****************************************************************************** !


! ****************************************************************************** !
  !> This function delivers the refinement level of a TreeID
  !!
  !! The refinement level of a
  !! \ref treelmesh_module::treelmesh_type::treeid "treeID" is found by
  !! repeatedly computing the parents with the equation used in the
  !! [tem_ParentOf function] (@ref tem_parentatlevel) and
  !! counting the iterations, until the universal root is reached.
  !!
  ! function tem_LevelOf_arrayOfTreeIDs(TreeID, nElems) result(res)
  !   ! ---------------------------------------------------------------------------
  !   !> number of elements in the array
  !   integer, intent(in) :: nElems
  !   !> Levels of the given TreeIDs
  !   integer :: res(nElems)
  !   !> ID in the complete tree to look up
  !   integer(kind=long_k), intent(in) :: TreeID(nElems)
  !   ! ---------------------------------------------------------------------------
  !   integer(kind=long_k) :: tElem
  !   integer :: level, iElem
  !   ! ---------------------------------------------------------------------------

  !   do iElem = 1, size(treeID)
  !     level = 0
  !     tElem = TreeID(iElem)
  !     do while (tElem /= 0)
  !       ! bit operation is more efficient than division
  !       ! division by 8 is equvalent to right shift 3 bits
  !       ! tElem = (tElem-1) / 8
  !       tElem = ishft( (tElem-1), -3 )
  !       level = level + 1
  !     end do
  !     res(iElem ) = level
  !   end do

  ! end function tem_LevelOf_arrayOfTreeIDs
! ****************************************************************************** !


! ****************************************************************************** !
  !> First TreeID in the complete tree on a given level
  !!
  !! e.g. level 1:  1 ..  8   -> offset = 1
  !!      level 2:  9 .. 72   -> offset = 9
  !!      level 3: 73 ..      -> offset = 73
  !!
  !! A level offset is calculated by summing all the nodes of coarser levels.
  !! With a level count starting at \f$0\f$ on the coarsest level (the universal
  !! root node), this offset \f$s\f$ for a given level \f$L\f$ can be computed by
  !! the follwing equation:
  !! \f[
  !!   s = \sum\limits_{i=0}^{L-1} 8^i = \frac{8^{L} - 1}{7}
  !! \f]
  !!
  !! @todo: calling to this function can be replaced by accessing array
  !!        firstIdAtLevel which explicitly stores the results
  elemental function tem_FirstIdAtLevel( level ) result(res)
    ! ---------------------------------------------------------------------------
    !> level to check
    integer, intent(in) :: level
    !> resulting first treeID
    integer(kind=long_k) :: res
    ! ---------------------------------------------------------------------------

    res = ( 8_long_k**level - 1_long_k ) / 7_long_k

  end function tem_FirstIdAtLevel
! ****************************************************************************** !


! ****************************************************************************** !
  !> Last ID in the complete tree on a given level
  !!
  elemental function tem_LastIdAtLevel( level ) result(res)
    ! ---------------------------------------------------------------------------
    !> level to check
    integer, intent(in) :: level
    !> resulting last treeID
    integer(kind=long_k) :: res
    ! ---------------------------------------------------------------------------
    integer :: i
    integer(kind=long_k) :: monomonial
    ! ---------------------------------------------------------------------------

    res = -1
    monomonial = 1
    do i = 0, level
      res = res + monomonial
      monomonial = monomonial * 8
    end do

  end function tem_LastIdAtLevel
! ****************************************************************************** !


! ****************************************************************************** !
  !> This function delivers all sibling treeIDs of TreeID in an array
  !!
  function tem_SiblingsOfId( TreeID ) result(siblings)
    ! ---------------------------------------------------------------------------
    !> current treeID
    integer(kind=long_k), intent(in) :: TreeID
    !> treeIDs siblings
    integer(kind=long_k) :: siblings(7)
    ! ---------------------------------------------------------------------------
    integer :: myNumber, iElem, i
    integer(kind=long_k) :: offset
    ! ---------------------------------------------------------------------------

    offset   = ishft((TreeID-1_long_k), -3) * 8_long_k
    myNumber = int( TreeID - offset )
    iElem = 0
    do i = 1,8
      if( i /= myNumber ) then
        iElem = iElem + 1
        siblings( iElem ) = offset + int(i,long_k)
      endif
    enddo

  end function tem_SiblingsOfId
! ****************************************************************************** !


! ****************************************************************************** !
  !> This function delivers of TreeID, which child number it is from its parent
  !!
  elemental function tem_childNumber( TreeID ) result(res)
    ! ---------------------------------------------------------------------------
    !> current treeID
    integer(kind=long_k), intent(in) :: TreeID
    !> result of function containing child number
    integer :: res
    ! ---------------------------------------------------------------------------
    integer(kind=long_k) :: offset
    ! ---------------------------------------------------------------------------

    ! offset = ((treeID-1)/8) * 8
    offset = ishft((TreeID-1_long_k), -3) * 8_long_k
    res    = int( TreeID - offset )

  end function tem_childNumber
! ****************************************************************************** !


! ****************************************************************************** !
  !> This function delivers the parent ID of a given TreeID
  !!
  elemental function tem_directParent(TreeID) result(res)
    ! ---------------------------------------------------------------------------
    !> current treeID
    integer(kind=long_k), intent(in) :: TreeID
    !> result of function containing parent ID
    integer(kind=long_k) :: res
    ! ---------------------------------------------------------------------------

    ! res = ( treeID - 1 ) / 8
    res = ishft( (treeID-1_long_k), -3 )

  end function tem_directParent
! ****************************************************************************** !


! ****************************************************************************** !
  !> This function returns the complete path through the tree from
  !! a given treeID to the root (all parents).
  !!
  elemental function tem_PathOf( TreeID ) result(res)
    ! ---------------------------------------------------------------------------
    !> current treeID
    integer(kind=long_k), intent(in) :: TreeID
    !> resulting path
    type(tem_path_type) :: res
    ! ---------------------------------------------------------------------------

    res%level = 1
    res%Node(res%Level) = treeID
    do while (res%Node(res%Level) > 0)
      res%Level = res%Level + 1
      ! res%Node(res%Level) = (res%Node(res%Level-1)-1) / 8
      res%Node(res%Level) = ishft( (res%Node(res%Level-1)-1), -3 )
    end do
  end function tem_PathOf
! ****************************************************************************** !


! ****************************************************************************** !
  !> This function provides the parent ID of a given tree ID on a given level.
  !!
  !! The parent ID on a given level is calculated using the equation
  !! \f[
  !!    id_{l-1}=\frac{id_{l}-1}{8}
  !! \f]
  !! in a loop.
  !!
  function tem_ParentAtLevel(TreeID, level) result(parentID)
    ! ---------------------------------------------------------------------------
    !> treeID of which the pID is requested
    integer(kind=long_k),intent(in) :: TreeID
    !> the level on which the pID is requested
    integer,intent(in) :: level
    !> resulting parent ID
    integer(kind=long_k) :: parentID
    ! ---------------------------------------------------------------------------
    integer :: cLevel ! the current level of the tree ID
    integer :: iLevel ! level iterator
    ! ---------------------------------------------------------------------------

    parentID = TreeID
    cLevel = tem_LevelOf(TreeID)

    do iLevel = 1, cLevel-level
      ! parentID = (parentID - 1) / 8
      parentID = ishft( (parentID-1_long_k), -3 )
    end do

  end function tem_ParentAtLevel
! ****************************************************************************** !


! ****************************************************************************** !
  !> This function delivers the direct children in the full tree for
  !! a given tree ID
  !!
  function tem_directChildren( TreeID ) result(childrenIDs)
    ! ---------------------------------------------------------------------------
    !> given treeID
    integer(kind=long_k), intent(in) :: TreeID
    !> Array for the treeIDs of the 8 direct children
    integer(kind=long_k) :: childrenIDs(8)
    ! ---------------------------------------------------------------------------
    integer :: childCount
    integer(kind=long_k) :: off
    ! ---------------------------------------------------------------------------

    childrenIDs(:) = 0

    off = treeID * 8_long_k
    do childCount = 1,8
      childrenIDs(childCount) = off + int(childCount, kind=long_k)
    end do

  end function tem_directChildren
! ****************************************************************************** !


! ****************************************************************************** !
  !> Compare Position of two treeIDs in the linearized tree
  !! This is relatively expansive, therefore the result should
  !! be stored, if more than one comparison of the same treeIDs
  !! has to be done.
  !! Result:
  !! -1: left is lower  than right
  !!  0: left is same   than right
  !!  1: left is higher than right
  !!
  elemental function tem_TreeIDComparison( left, right ) result(relation)
    ! ---------------------------------------------------------------------------
    !> candidate treeID
    integer(kind=long_k), intent(in) :: left
    !> candidate treeID
    integer(kind=long_k), intent(in) :: right
    !> relation between the treeIDs
    integer :: relation
    ! ---------------------------------------------------------------------------
    type(tem_path_type) :: left_path
    type(tem_path_type) :: right_path
    ! ---------------------------------------------------------------------------

    ! Build the path to the left from the leaf to the root (0).
    left_path = tem_PathOf(left)

    ! Build the path to the right from the leaf to the root (0).
    right_path = tem_PathOf(right)

    relation = tem_PathComparison(left_path, right_path)

  end function tem_TreeIDComparison
! ****************************************************************************** !


! ****************************************************************************** !
  !> Compare two Paths in the linearized tree
  !! Result:
  !! -1: left is lower  than right
  !!  0: left is same   than right
  !!  1: left is higher than right
  !!
  elemental function tem_PathComparison( left, right ) result(relation)
    ! ---------------------------------------------------------------------------
    !> candidate path
    type(tem_path_type), intent(in) :: left
    !> candidate path
    type(tem_path_type), intent(in) :: right
    !> relation between the paths
    integer :: relation
    ! ---------------------------------------------------------------------------
    integer :: maxPathLen
    integer(kind=long_k) :: diff
    ! ---------------------------------------------------------------------------

    ! Find the greatest common level.
    maxPathLen = min(left%level-1, right%level-1)

    ! Difference of treeIDs on that level indicates the relation of the two
    ! paths.
    diff =  left%Node( left%Level - maxPathLen)                                &
      &  - right%Node(right%Level - maxPathLen)

    ! Normalize the relation.
    relation = int(diff / max(1_long_k, abs(diff)))

  end function tem_PathComparison
! ****************************************************************************** !


! ****************************************************************************** !
  !> The following function provides the spatial index triple for a given ID in
  !! the complete mesh, on its refinement level.
  !! The returned coordinates start at 0.
  !! The fourth entry is the level on which the coordinates are located in the
  !! tree.
  !! The steps are following:
  !! 1. calculate the refinement level, store to coord(4).
  !! 2. Calcualte the level offset.
  !! 3. the Morton index is obtained by subtracting the offset from treeID.
  !! 4. Apply bit interleaving rule to generate coordinates.
  !!
  pure function tem_CoordOfId(TreeID, offset) result(coord)
    ! ---------------------------------------------------------------------------
    !> input element ID
    integer(kind=long_k), intent(in) :: TreeID
    !> First Elem of current level, if known, to avoid recomputations
    integer(kind=long_k), intent(in), optional :: offset
    !> coordinate of element return value
    integer :: coord(4)
    ! ---------------------------------------------------------------------------
    integer(kind=long_k) :: tElem ! temporary element
    integer(kind=long_k) :: fak(3)
    integer :: bitlevel
    integer :: i
    ! ---------------------------------------------------------------------------

    coord(:) = 0
    coord(4) = tem_LevelOf(treeID)
    fak(1) = 1
    fak(2) = 2
    fak(3) = 4
    bitlevel = 1

    ! offset has to be subtracted to get the "real" TreeID
    if (present(offset)) then
      tElem = TreeID - offset
    else
      tElem = TreeID - FirstIdAtLevel(coord(4))
    end if

    ! get x coordinate from
    !
    ! x = sum(iLevel=0, inf) { 2**iLevel * mod(tElem / (8**iLevel) ),2)   }
    ! y = sum(iLevel=0, inf) { 2**iLevel * mod(tElem /(2 * (8**iLevel) ),2)   }
    ! z = sum(iLevel=0, inf) { 2**iLevel * mod(tElem /(4 * (8**iLevel) ),2)   }
    !
    do while ((tElem / fak(1)) > 0)
      do i=1,3
        coord(i) = coord(i) + bitlevel * int(mod(tElem / fak(i), 2_long_k))
      end do
      bitlevel = bitlevel * 2
      fak = fak * 8
    end do
  end function tem_CoordOfId
! ****************************************************************************** !


! ****************************************************************************** !
  !> This function calculates the tree ID for a given x,y and z on the given
  !! refinement level. This ID in the complete tree does not have to be in the
  !! list of leafs (treeIDlist)
  !! An example of this procedure is following:
  !! 1. Convert coordinates into binary representation:
  !!   (x,y,z) = (5,9,1) = (0101,1001,0001)
  !! 2. Interleaving the bits by 3 digits for each direction:
  !!        x(0101):   **0 **1 **0 **1
  !!        y(1001):   *1* *0* *0* *1*
  !!        z(0001):   0** 0** 0** 1**
  !!    Combining these bits results in a binary number: 010 001 000 111,
  !!    which is 1095 when seen as a 10-base number.
  !! 3. This number is the cell position in a single level sorted element list.
  !!    Its treeID can be obtained by adding the correspoding level offset.
  !!
  pure function tem_IdOfCoord(coord, offset) result(nTreeID)
    ! ---------------------------------------------------------------------------
    !> 3D Coordinates to translate
    integer, intent(in) :: coord(4)
    !>
    integer(kind=long_k), intent(in), optional :: offset
    !> resulting treeID
    integer(kind=long_k) :: nTreeID
    ! ---------------------------------------------------------------------------
    integer :: iLevel    ! level iterator
    ! this variable defines the upper threshold for x, y and z coordinates
    ! for a given refinement level e.g. for level = 1 only x, y and z
    ! from 0 to 1 are allowed, otherwise this neighbor is outside of the
    ! geometry
    integer :: upperBoundary
    integer :: myCoord(3)
    integer :: fak2
    integer(kind=long_k) :: fak8
    ! ---------------------------------------------------------------------------

    ! the offset has to be added to the level summation
    ! @todo: level offset is hard coded in array firstIdAtLevel,
    !        thus no need to pass offset as an argument.
    if (present(offset)) then
      nTreeID = offset
    else
      nTreeID = firstIdAtLevel(coord(4))
    end if

    !! Periodic displacement for requested elements outside the cube
    upperBoundary = 2**coord(4)
    myCoord = mod(coord(1:3) + upperBoundary, upperBoundary)

    ! To compute the tree ID of a neighbor on the same refinement level the
    ! following formula is used:
    ! newTreeID = sum(i=0, level) {8**i * (mod( (x / (2**level) ) , 2)
    !                                      + 2 * mod( ( y / (2**level) ) , 2)
    !                                      + 4 * mod( (z / (2**level) ) , 2) ) }
    fak8 = 1
    fak2 = 1
    do iLevel = 0, coord(4) - 1
      nTreeID = nTreeID + fak8 * (      mod( (myCoord(1) / fak2), 2) &
        &                         + 2 * mod( (myCoord(2) / fak2), 2) &
        &                         + 4 * mod( (myCoord(3) / fak2), 2) )
      fak2 = fak2 * 2
      fak8 = fak8 * 8
    end do

  end function tem_IdOfCoord
! ****************************************************************************** !


end module tem_topology_module
! ****************************************************************************** !
