! See copyright notice in the COPYRIGHT file.
! ****************************************************************************** !
!> Geometric methods for the TreElM module to act on a octree based mesh.
!! Make sure to read the introduction to the distributed octree in the
!! Documentation.
!!
!! This module contains methods that work with the actual sparse mesh, and
!! therefore depends on the [[Treelmesh_module]].
module tem_geometry_module

  ! include treelm modules
  use mpi
  use env_module,              only: rk, long_k, output_unit, rk_mpi,          &
    &                                globalMaxLevels, labelLen
  use tem_float_module,        only: operator(.feq.)
  use tem_topology_module,     only: tem_LevelOf, tem_PathOf, tem_path_type,   &
    &                                tem_PathComparison,   &
    &                                tem_coordOfID
  use tem_param_module,        only: q__W, q__E, q__S, q__N, q__B, q__T, q_NE, &
    &                                q_SW, q_NW, q_SE, q_TE, q_BW, q_TW, q_BE, &
    &                                q_TN, q_BS, q_TS, q_BN, qTNE, qBSW, qBSE, &
    &                                qTNW, qTSE, qBNW, qBNE, qTSW,             &
    &                                childPosition
  use tem_tools_module,        only: tem_longList, append
  use treelmesh_module,        only: treelmesh_type
  use tem_property_module,     only: prp_hasBnd, prp_hasQVal
  use tem_subTree_type_module, only: tem_subTree_type
  use tem_debug_module,        only: tem_debug_type, tem_debug, main_debug, &
    &                                dbgUnit
  use tem_logging_module,      only: tem_toStr, logUnit

  implicit none

  private

  public :: tem_CoordOfReal
  public :: tem_tIDinfo
  public :: tem_BaryOfCoord, tem_BaryOfId
  public :: tem_originOfId
  public :: tem_ElemSize, tem_ElemSizeLevel
  public :: tem_vrtxCoordOfId
  public :: tem_eligibleChildren
  public :: tem_findElement, tem_findPath
  public :: tem_PosOfId, tem_PosOfPath
  public :: tem_direction_type
  public :: tem_neighbor_type
  ! public :: tem_GetRealBoundingCube
  public :: tem_GetLocalBoundingCube
  public :: tem_setEffBoundingBox
  public :: tem_determine_discreteVector
  public :: tem_intp_bilinear
  public :: tem_intp_trilinear
  public :: tem_build_treeToProp_pointer

  !> derived type for each direction
  type tem_direction_type
    !> first entry of the element list
    type(tem_longList),pointer :: first => null()
    !> number of entries in element list
    integer :: nElems
    !> which list JZ:(whether it is fluid, ghost or halo, see
    !! tem_levelDesc_type)
    integer :: list
    !> position in list
    integer :: pos
    !> the current neighbors are of a different level
    logical :: otherLevel = .false.
  end type tem_direction_type

  !> type neighbor_type includes the direct neighbors of each tree ID
  type tem_neighbor_type
    !> number of directions in which the neighbors are located
    !! e.g. for cubes and surface neighbors this will be 6
    !! and for cubes with surface, edge and vertex neighbors this will
    !! be 22
    !! because stencils will be mapped by this type as well, we are not able
    !! to use it as a constant
    integer :: nNeighborDirections
    !! directions for this neighbor
    !! JZ: in case of reconstructed DG this has to be extended to all
    !! JZ: cells inside the stencil
    type(tem_direction_type),allocatable :: dir(:)
  end type tem_neighbor_type

  interface tem_intp_bilinear
    module procedure tem_intp_bilinear_scalar
    module procedure tem_intp_bilinear_vec
  end interface
  interface tem_intp_trilinear
    module procedure tem_intp_trilinear_scalar
    module procedure tem_intp_trilinear_vec
  end interface

  interface tem_setEffBoundingBox
    module procedure tem_setEffBoundingBox_fromTree
    module procedure tem_setEffBoundingBox_fromSubTree
  end interface
  interface tem_GetLocalBoundingCube
    module procedure tem_GetLocalBoundingCube_fromTree
    module procedure tem_GetLocalBoundingCube_fromSubTree
  end interface

  contains


! ****************************************************************************** !
  !> This function returns a coordinate in the given treelmesh for a physical
  !! point location on the finest possible level.
  !!
  pure function tem_CoordOfReal(mesh, point, level) result(coord)
    ! ---------------------------------------------------------------------------
    type(treelmesh_type), intent(in) :: mesh !< Mesh to locate the point in
    real(kind=rk), intent(in) :: point(3) !< Point to look up
    integer, intent(in), optional :: level !< optional level to return the coordinate on
    integer :: coord(4) !< x,y,z,level
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: locInCube(3)
    real(kind=rk) :: meshDensity
    integer :: dimLen
    integer :: coordlevel
    ! ---------------------------------------------------------------------------

    if (present(level)) then
      coordlevel = min(level, globalMaxLevels)
    else
      coordlevel = globalMaxLevels
    end if

    locInCube = point - mesh%global%Origin
    dimLen = 2**coordlevel
    meshDensity = real(dimLen, kind=rk) / mesh%global%BoundingCubeLength

    ! Look up the real coordinate on the finest possible resolution.
    coord(4) = coordlevel

    ! Coordinates range from 0 to dimLen-1 in each direction
    ! Do not use periodic domain here, instead use the first
    ! and last element to capture all points outside the
    ! domain to deal with numerical inaccuracies.
    ! That is all elements include their lower boundaries, except
    ! the last one in a given direction, which includes the lower
    ! as well as the upper.
    coord(1:3) = max( min( int(locInCube*meshDensity), dimLen-1 ), 0 )
  end function tem_CoordOfReal
! ****************************************************************************** !


! ****************************************************************************** !
  !> Prints information about a treeID to a file unit.
  !!
  subroutine tem_tIDinfo( me, tree, nUnit)
    ! ---------------------------------------------------------------------------
    !> given level to get the size for
    integer(kind=long_k), intent(in) :: me
    !> Mesh to locate the point in
    type(treelmesh_type), intent(in) :: tree
    !> the file unit to use for printing the information.
    integer, intent(in), optional :: nUnit
    ! ---------------------------------------------------------------------------
    integer :: writeUnit, x(4)
    ! ---------------------------------------------------------------------------
    if( present( nUnit )) then
      writeUnit = nUnit
    else
      writeUnit = output_unit
    endif
    write(writeUnit,*) '  treeID info ', me
    x(:) = tem_coordOfId( me )
    write(writeUnit,*) '       level   ', x(4)
    write(writeUnit,*) '       coord   ', x(1:3)
    write(writeUnit,*) '       elemSize', tem_elemSize( tree, me )
    write(writeUnit,*) '       bary    ', tem_baryOfId( tree, me )


  end subroutine tem_tIDinfo
! ****************************************************************************** !


! ****************************************************************************** !
  !> Return the size of elements on a given levle in the mesh by taking into
  !! account the size of the bounding cube given in the global info of the tree
  !!
  pure function tem_ElemSizeLevel( tree, level ) result( dx )
    ! ---------------------------------------------------------------------------
    !> Mesh to locate the point in
    type(treelmesh_type), intent(in) :: tree
    !> given level to get the size for
    integer, intent(in) :: level
    !> size of Element
    real(kind=rk) :: dx
    ! ---------------------------------------------------------------------------

    dx = tree%global%BoundingCubeLength / real(2**level, kind=rk)

  end function tem_ElemSizeLevel
! ****************************************************************************** !


! ****************************************************************************** !
  !> Return the size of a given treeID in the mesh by taking into account the
  !! size of the bounding cube given in the global info of the tree
  !!
  pure function tem_ElemSize( tree, treeID ) result(dx)
    ! ---------------------------------------------------------------------------
    !> Mesh to locate point in
    type(treelmesh_type), intent(in) :: tree
    !> input elements
    integer(kind=long_k), intent(in) :: treeID
    !> size of element
    real(kind=rk) :: dx
    ! ---------------------------------------------------------------------------
    integer       :: level
    ! ---------------------------------------------------------------------------

    level = tem_levelOf(TreeID)
    dx = tem_elemSizeLevel( tree, level )

  end function tem_ElemSize
! ****************************************************************************** !


! ****************************************************************************** !
  !> Run through all the elements, check the vertices and return the fluid
  !! bounding cube
  !!
  function tem_GetLocalBoundingCube_fromTree( tree ) result(BoundingCube)
    ! ---------------------------------------------------------------------------
    !> global mesh information
    type(treelmesh_type), intent(in) :: tree
    !> xyz coordinate for min and max of bounding cube
    real(kind=rk) :: BoundingCube(3,2)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: vrtxCoord(3,8) ! coordinates for all eight vertices
    real(kind=rk) :: minX(3) ! xyz coordinate for min of bounding cube
    real(kind=rk) :: maxX(3) ! xyz coordinate for max of bounding cube
    integer(kind=long_k) :: tTreeID
    integer :: iElem, iVert
    ! ---------------------------------------------------------------------------

    minX =  huge( minX )
    maxX = -huge( maxX )

    do iElem = 1, tree%nElems
      tTreeID = tree%treeID( iElem )
      ! Calculate coordinates of vertices
      vrtxCoord(:,:) = tem_vrtxCoordOfId( tree, tTreeID )
      do iVert = 1, 8
        ! Compare with min max
        minX(:) = min( minX(:), vrtxCoord(:, iVert ))
        maxX(:) = max( maxX(:), vrtxCoord(:, iVert ))
      enddo
    enddo
    BoundingCube(:,1) = minX(:)
    BoundingCube(:,2) = maxX(:)

  end function tem_GetLocalBoundingCube_fromTree
! ****************************************************************************** !


! ****************************************************************************** !
  !> Run through all the elements, check the vertices and return the fluid
  !! bounding cube
  !!
  function tem_GetLocalBoundingCube_fromSubTree( subTree, globalTree )         &
    &                                                      result(BoundingCube)
    ! ---------------------------------------------------------------------------
    !> subTree to locate point in
    type(tem_subTree_type), intent(in) :: subTree
    !> corresponding global tree
    type(treelmesh_type), intent(in) :: globalTree
    !> xyz coordinate for min and max of bounding cube
    real(kind=rk) :: BoundingCube(3,2)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: vrtxCoord(3,8) ! coordinates for all eight vertices
    real(kind=rk) :: minX(3) ! xyz coordinate for min of bounding cube
    real(kind=rk) :: maxX(3) ! xyz coordinate for max of bounding cube
    integer(kind=long_k) :: tTreeID
    integer :: iElem, iVert
    ! ---------------------------------------------------------------------------

    ! if the subTree equals to the global shape
    if( subTree%useGlobalMesh ) then
      ! ... the effective bounding cube equals to the one of the tree
      BoundingCube = globalTree%global%effboundingCube(3,2)
    else ! if this is not the case
      minX =  huge( minX )
      maxX = -huge( maxX )

      ! loop over all elements
      do iElem = 1, subTree%nElems
        ! ... and distinguish between local shapes and other shapes
        if( subTree%useLocalMesh ) then
          ! ... copy the treeIDs directly
          tTreeID = subTree%treeID( iElem )
        else
          ! ... the treeIDs from the globalTree via the positions in map2global
          tTreeID = globalTree%treeID( subTree%map2global( iElem ))
        end if
        ! Calculate coordinates of vertices
        vrtxCoord(:,:) = tem_vrtxCoordOfId( globalTree, tTreeID )
        do iVert = 1, 8
          ! Compare with min max
          minX(:) = min( minX(:), vrtxCoord(:, iVert ))
          maxX(:) = max( maxX(:), vrtxCoord(:, iVert ))
        enddo
      enddo
      BoundingCube(:,1) = minX(:)
      BoundingCube(:,2) = maxX(:)
    end if ! useGlobalShape

  end function tem_GetLocalBoundingCube_fromSubTree
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate the real bounding box around the fluid domain
  !!
  subroutine tem_setEffBoundingBox_fromTree( tree )
    ! ---------------------------------------------------------------------------
    !> Mesh
    type(treelmesh_type), intent(inout) :: tree
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: boundingBox(3,2)
    ! ---------------------------------------------------------------------------

    boundingBox(:,:) = tem_GetRealBoundingCube( tree )

    ! Set the effective origin and length in the global tree part
    tree%global%effOrigin(:) = boundingBox(:,1)
    tree%global%effLength(:) = boundingBox(:,2) -  boundingBox(:,1)

  end subroutine tem_setEffBoundingBox_fromTree
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate the real bounding box around the fluid domain
  !!
  subroutine tem_setEffBoundingBox_fromSubTree( subTree, globalTree )
    ! ---------------------------------------------------------------------------
    !> subTree to get effective bounding cube from
    type(tem_subTree_type)   :: subTree
    !> corresponding global tree
    type(treelmesh_type), intent(in) :: globalTree
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: boundingBox(3,2)
    real(kind=rk) :: tBounding(3)
    integer :: iErr
    ! ---------------------------------------------------------------------------

    ! if the subTree equals to the global tree take over the settings
    if( subTree%useGlobalMesh )then
      subTree%global%effOrigin = globalTree%global%effOrigin
      subTree%global%effLength = globalTree%global%effLength
    else ! subTree is not equal to the global tree
      ! Calculate process-local bounding cube
      boundingBox = tem_GetLocalBoundingCube( subTree, globalTree )
      ! Exchange with neighbors
      call mpi_allreduce( boundingBox(:,1), tBounding, 3, rk_mpi, mpi_min,     &
        &                 subTree%global%comm, iErr  )
      boundingBox(:,1) = tBounding
      call mpi_allreduce( boundingBox(:,2), tBounding, 3, rk_mpi, mpi_max,     &
        &                 subTree%global%comm, iErr  )
      boundingBox(:,2) = tBounding

      ! Set the effective origin and length in the global subTree part
      subTree%global%effOrigin(:) = boundingBox(:,1)
      subTree%global%effLength(:) = boundingBox(:,2) -  boundingBox(:,1)
    end if

  end subroutine tem_setEffBoundingBox_fromSubTree
! ****************************************************************************** !


! ****************************************************************************** !
  !> The following function provides the coordinates of the barycenter for a
  !! given treeID in the complete mesh.
  !!
  pure function tem_BaryOfId(tree,TreeID) result(bary)
    ! ---------------------------------------------------------------------------
    !> mesh information
    type(treelmesh_type), intent(in) :: tree
    !> input Element ID
    integer(kind=long_k), intent(in) :: TreeID
    !> barycenter return value
    real(kind=rk) :: bary(3)
    ! ---------------------------------------------------------------------------
    integer       :: coord(4) ! spatial index triple for a given ID
    real(kind=rk) :: dx ! size of Element
    ! ---------------------------------------------------------------------------
    coord = tem_CoordOfId(TreeID)
    dx = tem_elemSizeLevel( tree, coord(4) )
    bary = tree%global%origin + (real(coord(:3), kind=rk) + 0.5_rk)*dx

  end function tem_BaryOfId
! *******************************************************************************


! *******************************************************************************
  !> The following function provides the coordinates of the origin for a
  !! given ID in the complete mesh.
  !!
  pure function tem_originOfId(tree,TreeID) result(origin)
    ! ---------------------------------------------------------------------------
    !> mesh information
    type(treelmesh_type), intent(in) :: tree
    !> input element ID
    integer(kind=long_k), intent(in) :: TreeID
    !> origin return value
    real(kind=rk) :: origin(3)
    ! ---------------------------------------------------------------------------
    integer       :: coord(4) !< spatial index triple for a given ID
    real(kind=rk) :: dx !< size of Element
    ! ---------------------------------------------------------------------------
    coord = tem_CoordOfId( TreeID )
    dx = tem_elemSizeLevel( tree, coord(4) )
    origin = tree%global%origin + real(coord(:3), kind=rk)*dx

  end function tem_originOfId
! *******************************************************************************


! *******************************************************************************
  !> This function provides the coordinates of the element barycenters
  !! for a set of given element coordinates on the same refinement level.
  pure function tem_BaryOfCoord(coord, nPoints, origin, dx) result(bary)
    ! ------------------------------------------------------------------------
    integer, intent(in) :: nPoints !< Number of points to evaluate
    integer, intent(in) :: coord(nPoints,3)
    !> spatial index triple for a given ID
    real(kind=rk), intent(in) :: origin(3) !< origin of the universe cube
    real(kind=rk), intent(in) :: dx !< size of the elements
    real(kind=rk) :: bary(nPoints,3) !< barycenter to return
    ! ------------------------------------------------------------------------
    integer :: ii
    real(kind=rk) :: c(3)
    ! ------------------------------------------------------------------------

    ! bary = origin + ( coord + 0.5 ) * dx
    !      = origin + coord * dx + 0.5dx
    !      = c      + coord * dx
    c(:) = origin(:) + 0.5_rk * dx
    do ii = 1, nPoints
      bary(ii,1) = c(1) + real(coord(ii,1), kind=rk) * dx
      bary(ii,2) = c(2) + real(coord(ii,2), kind=rk) * dx
      bary(ii,3) = c(3) + real(coord(ii,3), kind=rk) * dx
    end do

  end function tem_BaryOfCoord
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate all eight vertex coordinates of a given tree ID
  !!
  !! the numbering of the vertices is according to
  !! tem_param_module::childPosition
  !! \author Manuel Hasert
  !!
  pure function tem_vrtxCoordOfId( tree, treeID) result(coord)
    ! ---------------------------------------------------------------------------
    !> complete tree for info about dimensions
    type(treelmesh_type), intent(in) :: tree
    !> input element ID
    integer(kind=long_k), intent(in) :: treeID
    !> all vertices coordinates function return value
    real(kind=rk) :: coord(3,8)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: elemBary(3) !< barycenter of current element
    real(kind=rk) :: length  !< element size
    integer :: x(4), iCoord, iVrtx
    ! ---------------------------------------------------------------------------

    x = tem_coordOfId( treeID )
    ! Get element size
    length = tree%global%BoundingCubeLength / real( 2**x(4), kind=rk)

    ! Get Barycenter
    elemBary = tem_BaryOfId( tree, treeID )

    ! ... and calculate from there the vertex coordinates
    do iVrtx = 1, 8
      do iCoord = 1, 3
        coord( iCoord, iVrtx ) = elemBary( iCoord )                            &
          &                    + length*0.5_rk*childPosition( iVrtx, iCoord )
      enddo
    end do

  end function tem_vrtxCoordOfId
! ****************************************************************************** !


! ****************************************************************************** !
  !> Identify all possible children local ids for each of the 27 direct
  !! neighbors results are saved in the ElemList
  !!
  !! How to use this information:
  !! There are two main ways of doing so.
  !!
  !! 1.) Use eligibleChildren array as a subset array of directChildren
  !! Let's say, you have a neighbor on the right border (=East border),
  !! which is on a coarser level.
  !! But you need the treeIDs of the neighbor on the same level.
  !! You then call tem_directChildren( coarseNEighborTreeID )
  !! to get all the children.
  !! Then, you can access the treeIDs on the east border (that is the
  !! treeIDs of the children on the West!! border of the coarser neighbor)
  !! by using the eligibleChildren as a subset to access the tem_directChildren
  !! childrenIDs( eligible_child( iChild ))
  !!
  !! 2.) Use elgibleChildren as an offset from the lower, left, bottom
  !! child tree ID of the parent.
  !!
  !! This is routine is Morton curve specific
  !! HK: and solver specific?
  !!
  pure subroutine tem_eligibleChildren(eligible_child, Direction)
    ! ---------------------------------------------------------------------------
    !> Candidate children, which might be considered as neighbors
    integer, intent(out),allocatable :: eligible_child(:)
    !> In which direction to search for neighbors
    integer, intent(in) :: Direction
    ! ---------------------------------------------------------------------------

    select case ( Direction )

    ! direct neighbor directions: 4 children
    case( q__W)
      allocate(eligible_child(4))
      eligible_child(1) = 1
      eligible_child(2) = 3
      eligible_child(3) = 5
      eligible_child(4) = 7
    case( q__E)
      allocate(eligible_child(4))
      eligible_child(1) = 2
      eligible_child(2) = 4
      eligible_child(3) = 6
      eligible_child(4) = 8
    case( q__S)
      allocate(eligible_child(4))
      eligible_child(1) = 1
      eligible_child(2) = 2
      eligible_child(3) = 5
      eligible_child(4) = 6
    case( q__N)
      allocate(eligible_child(4))
      eligible_child(1) = 3
      eligible_child(2) = 4
      eligible_child(3) = 7
      eligible_child(4) = 8
    case( q__B)
      allocate(eligible_child(4))
      eligible_child(1) = 1
      eligible_child(2) = 2
      eligible_child(3) = 3
      eligible_child(4) = 4
    case( q__T)
      allocate(eligible_child(4))
      eligible_child(1) = 5
      eligible_child(2) = 6
      eligible_child(3) = 7
      eligible_child(4) = 8

    ! edge neighbor directions: 2 children
    case( q_SW)
      allocate(eligible_child(2))
      eligible_child(1) = 1
      eligible_child(2) = 5
    case( q_NW)
      allocate(eligible_child(2))
      eligible_child(1) = 3
      eligible_child(2) = 7
    case( q_SE)
      allocate(eligible_child(2))
      eligible_child(1) = 2
      eligible_child(2) = 6
    case( q_NE)
      allocate(eligible_child(2))
      eligible_child(1) = 4
      eligible_child(2) = 8

    case( q_BW)
      allocate(eligible_child(2))
      eligible_child(1) = 1
      eligible_child(2) = 3
    case( q_TW)
      allocate(eligible_child(2))
      eligible_child(1) = 5
      eligible_child(2) = 7
    case( q_BE)
      allocate(eligible_child(2))
      eligible_child(1) = 2
      eligible_child(2) = 4
    case( q_TE)
      allocate(eligible_child(2))
      eligible_child(1) = 6
      eligible_child(2) = 8

    case( q_BS)
      allocate(eligible_child(2))
      eligible_child(1) = 1
      eligible_child(2) = 2
    case( q_TS)
      allocate(eligible_child(2))
      eligible_child(1) = 5
      eligible_child(2) = 6
    case( q_BN)
      allocate(eligible_child(2))
      eligible_child(1) = 3
      eligible_child(2) = 4
    case( q_TN)
      allocate(eligible_child(2))
      eligible_child(1) = 7
      eligible_child(2) = 8

   ! corner neighbor directions: only one child
    case( qBSW)
      allocate(eligible_child(1))
      eligible_child(1) = 1
    case( qTNE)
      allocate(eligible_child(1))
      eligible_child(1) = 8
    case( qTNW)
      allocate(eligible_child(1))
      eligible_child(1) = 7
    case( qBSE)
      allocate(eligible_child(1))
      eligible_child(1) = 2
    case( qBNW)
      allocate(eligible_child(1))
      eligible_child(1) = 3
    case( qTSE)
      allocate(eligible_child(1))
      eligible_child(1) = 6
    case( qTSW)
      allocate(eligible_child(1))
      eligible_child(1) = 5
    case( qBNE)
      allocate(eligible_child(1))
      eligible_child(1) = 4

    end select ! select ( Direction )

  end subroutine tem_eligibleChildren
! ****************************************************************************** !


! ****************************************************************************** !
  !> Recursive routine to find all actual (eligible) leave nodes in the
  !! local partition for a given treeID.
  !! Alternatively use tem_findPath, which uses precomputed paths in the tree
  !! and should speed up the search (at the expense of storing the paths
  !! beforehand).
  !!
  recursive subroutine tem_findElement( TreeID, eligible_child, ElemList,      &
    &                                   treeIDlist, nElems, Part_First,        &
    &                                   Part_Last, otherLevel)
    ! ---------------------------------------------------------------------------
    !> TreeID to find in the array of Elements
    integer(kind=long_k) :: TreeID
    !> Candidate childs, which might be considered as neighbors
    integer :: eligible_child(:)
    !> linked list of resulting elements building the neighbor
    type(tem_longList), pointer :: ElemList
    !> number of elements in list
    integer, intent(in)  :: nElems
    !> array of treeIDs
    integer(kind=long_k), intent(in)  :: treeIDlist(nElems)
    !> parts first entry
    integer(kind=long_k), intent(in) :: Part_First(:)
    !> parts last entry
    integer(kind=long_k), intent(in) :: Part_Last(:)
    !> entry is on another level
    logical,optional,intent(inout) :: otherLevel
    ! ---------------------------------------------------------------------------
    integer(kind=long_k) :: pos
    integer :: i
    integer(kind=long_k) :: childID, off
    ! ---------------------------------------------------------------------------

    ! binary search of the ID in the array of actual present elements
    ! Return pos < 0 if
    pos = tem_PosOfId(TreeID, treeIDlist)

    ! If the neighbor is on a level higher than myself, it should be
    ! delivered by binary search

    if (pos > 0 ) then
      ! Element actually exists, append it to the list
       call append(ElemList, pos)
    else
      if (pos < 0) then
        ! Element is a virtual neighbor, look for childs
        if( present( otherLevel ) ) otherLevel = .true.
        off = TreeID*8
        do i=1,size(eligible_child)
          childID = off + eligible_child(i)
          call tem_findElement( childID, eligible_child, ElemList,             &
            &                   treeIDlist, nElems, Part_First, Part_Last )
        end do
      end if
    end if

  end subroutine tem_findElement
! ****************************************************************************** !


! ****************************************************************************** !
  !> Recursive routine to find all actual (eligible) leave nodes in the
  !! local partition for a given treeID.
  !!@todo HK: when doing this for the complete domain it would probably better
  !!          just compute the ID on the finest level along with the information
  !!          on the level of each leaf, in order to speed up things a little.
  !!          This way comparison would be just a simple integer comparison.
  !!
  recursive subroutine tem_findPath( Path, eligible_child, ElemList,           &
    &                                Pathlist, nElems, otherLevel )
    ! ---------------------------------------------------------------------------
    !> Path to the leaf to find in the array of Elements
    type(tem_path_type), intent(in) :: Path
    !> Candidate childs, which might be considered as neighbors
    integer, intent(in) :: eligible_child(:)
    !> linked list of resulting elements building the neighbor
    type(tem_longList), pointer :: ElemList
    !> number of elements in list
    integer, intent(in)  :: nElems
    !> array of paths
    type(tem_path_type), intent(in)  :: pathlist(nElems)
    !> entry is on another level
    logical,optional,intent(inout) :: otherLevel
    ! ---------------------------------------------------------------------------
    integer(kind=long_k) :: pos
    integer :: i
    integer(kind=long_k) :: off
    type(tem_path_type) :: childPath
    ! ---------------------------------------------------------------------------

    ! binary search of the Path in the array of actual present elements
    ! Return pos < 0 if
    pos = tem_PosOfPath(Path, Pathlist)

    ! If the neighbor is on a level higher than myself, it should be
    ! delivered by binary search

    if (pos > 0 ) then
      ! Element actually exists, append it to the list
       call append(ElemList, pos)
    else if (pos < 0) then
      ! Element is a GhostFromFiner, look for childs
      if( present( otherLevel ) ) otherLevel = .true.
      off = Path%Node(1)*8
      do i=1,size(eligible_child)
        childPath = tem_pathOf(off + eligible_child(i))
        call tem_findPath( childPath, eligible_child, ElemList, &
          &                Pathlist, nElems, otherLevel)
      end do
    else ! pos == 0
    end if

  end subroutine tem_findPath
! ****************************************************************************** !


! ****************************************************************************** !
  !> This subroutine does a binary search on a given (sparse) list of elements.
  !! The result is the position of the given tree ID in the list,
  !! 0 if no corresponding node is found,
  !! or the negative of the found ID, if it is a virtual node.
  !!
  pure function tem_PosOfId(sTreeID, treeIDlist, lower, upper) result(IdPos)
    ! ---------------------------------------------------------------------------
    !> tree ID to search for
    integer(kind=long_k), intent(in)  :: sTreeID
    !> List to search in
    integer(kind=long_k), intent(in)  :: treeIDlist(:)
    !> lowerbound of search interval
    integer, intent(in), optional     :: lower
    !> upperbound of search interval
    integer, intent(in), optional     :: upper
    !> position of sTreeID in the list of elements
    integer                           :: IdPos
    ! ---------------------------------------------------------------------------
    integer :: lb, ub
    integer :: middleSearch
    type(tem_path_type) :: searched
    type(tem_path_type) :: current
    integer :: pathRelation
    ! ---------------------------------------------------------------------------

    if (present(lower)) then
      lb = lower
    else
      lb = lbound(treeIDList,1)
    end if
    if (present(upper)) then
      ub = upper
    else
      ub = ubound(treeIDList,1)
    end if

    !> Build the path to the searched TreeID from the leaf to the root.
    searched = tem_PathOf(sTreeID)

    ! Start the Binary search for the neighbor elements
    binSearchLoop: do
      middleSearch = (lb + ub) / 2
      ! Build the path to the currently investigated element from leaf to root.
      current = tem_PathOf(treeIDList(middleSearch))

      pathRelation = tem_PathComparison(searched, current)

      if ((pathRelation == 0) .or. (lb >= ub)) then
        ! Leave the loop, if element has been found, or this
        ! was the last element to investigate.
        exit binSearchLoop
      else
        halves: if (pathRelation == 1) then
          ! Continue the search in the higher half, as the looked up element is
          ! to small.
          lb = min(middleSearch + 1, ub)
        else
          ! Continue search in the lower half, as the looked up element is to
          ! large.
          ub = max(middleSearch - 1, lb)
        end if halves
      end if
    end do binSearchLoop

    if (pathRelation == 0) then
      if (current%Level <= searched%Level) then
        ! The found ID is actually a leaf
        IdPos = middleSearch
      else
        ! The found ID is a parent of the searched
        ! virtual treeID
        IdPos = -middleSearch
      end if
    else
      IdPos = 0 ! no matching element found.
    end if

  end function tem_PosOfId
! ****************************************************************************** !


! ****************************************************************************** !
  !> @todo: Add description
  !!
  pure function tem_PosOfPath(sPath, Pathlist, lower, upper) result(IdPos)
    ! ---------------------------------------------------------------------------
    !>
    type(tem_path_type), intent(in) :: sPath
    !>
    type(tem_path_type), intent(in) :: Pathlist(:)
    !>
    integer, intent(in), optional :: lower
    !>
    integer, intent(in), optional :: upper
    !>
    integer  :: IdPos
    ! ---------------------------------------------------------------------------
    integer :: lb, ub
    integer :: middleSearch
    type(tem_path_type) :: current
    integer(kind=long_k) :: pathRelation
    integer :: maxPathLen
    ! ---------------------------------------------------------------------------

    if (present(lower)) then
      lb = lower
    else
      lb = lbound(PathList,1)
    end if
    if (present(upper)) then
      ub = upper
    else
      ub = ubound(PathList,1)
    end if

    ! Start the Binary search for the neighbor elements
    binSearchLoop: do
      middleSearch = (lb + ub) / 2

      ! Build the path to the currently investigated element from leaf to root.
      current = PathList(middleSearch)

      ! pathRelation = tem_PathComparison(sPath, current)
      ! Inlined tem_PathComparison:
      maxPathLen = min(sPath%level-1, current%level-1)

      pathRelation = sPath%Node(sPath%level - maxPathLen)                      &
        &          - current%Node(current%level - maxPathLen)

      if ((pathRelation == 0_long_k) .or. (lb >= ub)) then
        ! Leave the loop, if element has been found, or this
        ! was the last element to investigate.
        exit binSearchLoop
      else
        halves: if (pathRelation > 0_long_k) then
          ! Continue the search in the higher half, as the looked up element is
          ! to small.
          lb = min(middleSearch + 1, ub)
        else
          ! Continue search in the lower half, as the looked up element is to
          ! large.
          ub = max(middleSearch - 1, lb)
        end if halves
      end if

    end do binSearchLoop

    if (pathRelation == 0_long_k) then
      if (current%Level <= sPath%Level) then
        ! The found ID is actually a leaf
        IdPos = middleSearch
      else
        ! The found ID is a child of the searched
        ! virtual treeID
        IdPos = -middleSearch
      end if
    else
      IdPos = 0 ! no matching element found.
    end if

  end function tem_PosOfPath
! ****************************************************************************** !


! ****************************************************************************** !
  !> Compare the incoming discrete vector against a set of prevailing vectors
  !! and return the found prevailing integer vector
  !!
  !! a set of float vectors must be given
  !! A minimazation of the angle between the vector and the compareVector
  !! entries is performed.
  !! \author Manuel Hasert
  !!
  subroutine tem_determine_discreteVector( vector, compareVector, angle )
    ! ---------------------------------------------------------------------------
    !> The given vector, will be set to the vector in compareVector
    !! which is closest in terms of the angle
    integer, intent(inout) :: vector(:)
    !> set of vectors to compare against.
    !! size is (3, nVectors)
    real( kind=rk ), intent(in) :: compareVector(:,:)
    !> angle between vectors
    real( kind=rk ), intent(out), optional :: angle
    ! ---------------------------------------------------------------------------
    integer :: length, iDir, minIndex
    real(kind=rk) :: rVectIn(3), angleLoc, dotProduct, minLen
    ! ---------------------------------------------------------------------------
    ! normalize vector
    length = vector(1)**2 + vector(2)**2 + vector(3)**2
    angleLoc = 0.0_rk
    if( length > 0 ) then
      ! Only do checks if non-zero vector
      ! Normalize vector
      rVectIn(:) = dble( abs(vector(:)) )/dble(sqrt( dble( length )))
      minLen = 1000._rk
      minIndex = 0
      ! Go over all vectors in the compareVector
      do iDir = 1, size(compareVector, 2)
        ! Go through all prevailing directions and check the minimum
        dotProduct = rVectIn(1)*compareVector( 1, iDir )                       &
          &        + rVectIn(2)*compareVector( 2, iDir )                       &
          &        + rVectIn(3)*compareVector( 3, iDir )

        dotProduct = min( dotProduct , 1.0_rk )
        angleLoc = acos( dotProduct )
        if( angleLoc < minLen ) then
          minLen = angleLoc
          minIndex = iDir
          if( angleLoc .feq. 0._rk ) exit
        endif
      enddo
      ! Copy direction
      ! What about the sign?
      vector(:) = sign(int(0.5+compareVector( :, minIndex)), vector(:))
    endif ! not zero-vector

    if (present(angle)) angle = angleLoc

  end subroutine tem_determine_discreteVector
! ****************************************************************************** !


! ****************************************************************************** !
  !> This function returns the bi-linearly interpolated values from the four
  !! source points to the target position located at targetCoord.
  !! The source points are arranged in a square from (0,0)x(1,1)
  !! The order of the source points are according to the morton curve
  !!  1      2      3      4
  !! (0,0); (1,0); (0,1); (1,1)
  !!
  !!\author Manuel Hasert
  !!
  function tem_intp_bilinear_scalar( srcVal, targetCoord ) result( phi )
    ! ---------------------------------------------------------------------------
    !> source values of the square corners
    real(kind=rk), intent(in) :: srcVal(4)
    !> interpolation location within the square
    real(kind=rk), intent(in) :: targetCoord(2)
    !> interpolated value
    real(kind=rk) :: phi
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: phi_north, phi_south
    ! ---------------------------------------------------------------------------
    ! Linear interpolation on the north nodes
    phi_north = (1._rk - targetCoord(1))* srcVal(3) + targetCoord(1)*srcVal(4)
    ! Linear interpolation on the south nodes
    phi_south = (1._rk - targetCoord(1))* srcVal(1) + targetCoord(1)*srcVal(2)
    ! Linear interpolation on the obtained north and south values
    phi       = (1._rk - targetCoord(2))* phi_south + targetCoord(2)*phi_north

  end function tem_intp_bilinear_scalar
! ****************************************************************************** !


! ****************************************************************************** !
  !> This function returns the bi-linearly interpolated values from the four
  !! source points to the target position located at targetCoord.
  !! The source points are arranged in a square from (0,0)x(1,1)
  !! The order of the source points are according to the morton curve
  !!  1      2      3      4
  !! (0,0); (1,0); (0,1); (1,1)
  !!
  !!\author Manuel Hasert
  !!
  function tem_intp_bilinear_vec( srcVal, targetCoord, nVals ) result( phi )
    ! ---------------------------------------------------------------------------
    !> number of values
    integer, intent(in) :: nVals
    !> source values of the square corners
    real(kind=rk), intent(in) :: srcVal(nVals,4)
    !> interpolation location within the square
    real(kind=rk), intent(in) :: targetCoord(2)
    !> interpolated values
    real(kind=rk) :: phi(nVals)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: phi_north, phi_south
    integer :: iVal
    ! ---------------------------------------------------------------------------
    do iVal = 1, nVals
      ! Linear interpolation on the north nodes
      phi_north   = (1._rk - targetCoord(1))* srcVal(iVal,3)                   &
        &                  + targetCoord(1)*srcVal(iVal,4)
      phi_south   = (1._rk - targetCoord(1))* srcVal(iVal,1)                   &
        &                   + targetCoord(1)*srcVal(iVal,2)
      phi( iVal ) = (1._rk - targetCoord(2))* phi_south                        &
        &                   + targetCoord(2)*phi_north
    enddo

  end function tem_intp_bilinear_vec
! ****************************************************************************** !


! ****************************************************************************** !
  !> This function returns the tri-linearly interpolated values from the eight
  !! source points to the target position located at targetCoord.
  !! The source points are arranged in a square from (0,0,0)x(1,1,1)
  !! The order of the source points are according to the morton curve
  !!  1        2        3        4
  !! (0,0,0); (1,0,0); (0,1,0); (1,1,0)
  !!  5        6        7        8
  !! (0,0,1); (1,0,1); (0,1,1); (1,1,1)
  !!
  !!\author Manuel Hasert
  !!
  function tem_intp_trilinear_scalar( srcVal, targetCoord ) result( phi )
    ! ---------------------------------------------------------------------------
    !> source values of the square corners
    real(kind=rk), intent(in) :: srcVal(8)
    !> interpolation location within the square
    real(kind=rk), intent(in) :: targetCoord(3)
    !> interpolated value
    real(kind=rk) :: phi
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: phi_northFront, phi_southFront
    real(kind=rk) :: phi_northBack, phi_southBack
    real(kind=rk) :: phi_front, phi_back
    ! ---------------------------------------------------------------------------
    ! Linear interpolation on the cube front side (z = 0 )
    phi_northFront = (1._rk - targetCoord(1))* srcVal(3)                       &
      &                     + targetCoord(1) * srcVal(4)
    phi_southFront = (1._rk - targetCoord(1))* srcVal(1)                       &
      &                     + targetCoord(1) * srcVal(2)
    ! Linear interpolation on the cube back side (z = 1 )
    phi_northBack  = (1._rk - targetCoord(1))* srcVal(7)                       &
      &                     + targetCoord(1) * srcVal(8)
    phi_southBack  = (1._rk - targetCoord(1))* srcVal(5)                       &
      &                     + targetCoord(1) * srcVal(6)
    ! Linear interpolation on the cube front side (z = 0 )
    phi_front = (1._rk - targetCoord(2))* phi_southFront                       &
      &                + targetCoord(2) * phi_northFront
    ! Linear interpolation on the cube back side (z = 1 )
    phi_back  = (1._rk - targetCoord(2))* phi_southBack                        &
      &                + targetCoord(2) * phi_northBack
    phi       = (1._rk - targetCoord(3))* phi_front  + targetCoord(3)*phi_back

  end function tem_intp_trilinear_scalar
! ****************************************************************************** !


! ****************************************************************************** !
  !> This function returns the tri-linearly interpolated values from the eight
  !! source points to the target position located at targetCoord.
  !! The source points are arranged in a square from (0,0,0)x(1,1,1)
  !! The order of the source points are according to the morton curve
  !!  1        2        3        4
  !! (0,0,0); (1,0,0); (0,1,0); (1,1,0)
  !!  5        6        7        8
  !! (0,0,1); (1,0,1); (0,1,1); (1,1,1)
  !!
  !!\author Manuel Hasert
  !!
  function tem_intp_trilinear_vec( srcVal, targetCoord, nVals ) result( phi )
    ! ---------------------------------------------------------------------------
    !> number of values
    integer, intent(in) :: nVals
    !> source values of the square corners
    real(kind=rk), intent(in) :: srcVal(nVals,8)
    !> interpolation location within the square
    real(kind=rk), intent(in) :: targetCoord(3)
    !> interpolated value
    real(kind=rk) :: phi(nVals)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: phi_northFront, phi_southFront
    real(kind=rk) :: phi_northBack,  phi_southBack
    real(kind=rk) :: phi_front, phi_back
    integer :: iVal
    ! ---------------------------------------------------------------------------
    do iVal = 1, nVals
      ! Linear interpolation on the north nodes
      phi_northFront = (1._rk - targetCoord(1))* srcVal(iVal,3)                &
        &                     + targetCoord(1) * srcVal(iVal,4)
      phi_southFront = (1._rk - targetCoord(1))* srcVal(iVal,1)                &
        &                     + targetCoord(1) * srcVal(iVal,2)
      ! Linear interpolation on the cube back side (z = 1 )
      phi_northBack  = (1._rk - targetCoord(1))* srcVal(iVal,7)                &
        &                     + targetCoord(1) * srcVal(iVal,8)
      phi_southBack  = (1._rk - targetCoord(1))* srcVal(iVal,5)                &
        &                     + targetCoord(1) * srcVal(iVal,6)
      ! Linear interpolation on the cube front side (z = 0 )
      phi_front = (1._rk - targetCoord(2))* phi_southFront                     &
        &                + targetCoord(2) * phi_northFront
      ! Linear interpolation on the cube back side (z = 1 )
      phi_back  = (1._rk - targetCoord(2))* phi_southBack                      &
        &                + targetCoord(2) * phi_northBack
      phi(iVal) = (1._rk - targetCoord(3))* phi_front  + targetCoord(3)*phi_back
    enddo

  end function tem_intp_trilinear_vec
! ****************************************************************************** !


! ****************************************************************************** !
  !> This routine builds mapping from tree%treeID to given property data array
  !! like boundary_ID or qVal
  subroutine tem_build_treeToProp_pointer( treeToProp, nElems, ElemPropertyBits, &
    &                                      prp_bit)
    ! --------------------------------------------------------------------------
    !> mesh elements
    integer :: nElems
    !> Pointer from tree to property data to be filled by this routine
    integer, intent(out) :: treeToProp(nElems)
    !> Elements Property Bits
    integer(kind=long_k), intent(in) :: ElemPropertyBits(nElems)
    !> property bit
    integer, intent(in) :: prp_bit
    ! --------------------------------------------------------------------------
    character(len=labelLen) :: buffer
    integer :: iElem, counter
    ! --------------------------------------------------------------------------
    select case(prp_bit)
    case (prp_hasBnd)
      buffer = 'boundary_ID'
    case (prp_hasQVal)
      buffer = 'qVal'
    end select

    write(logUnit(7),*) 'Building map from tree to '//trim(buffer)

    counter = 0
    do iElem = 1, nElems
      if ( btest(ElemPropertyBits( iElem ), prp_bit ) ) then
        counter = counter + 1
        treeToProp(iElem) = counter
      else
        treeToProp(iElem) = -1
      end if
    end do

  end subroutine tem_build_treeToProp_pointer
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate the real bounding box around the fluid domain
  !! and return only to root (efficiency reasons)
  !!
  function tem_GetRealBoundingCube( tree ) result( boundingCube )
    ! ---------------------------------------------------------------------------
    !>
    type(treelmesh_type), intent(in) :: tree
    !>
    real(kind=rk) :: boundingCube(3,2)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: tBounding(3)
    integer :: iErr
    ! ---------------------------------------------------------------------------
    ! Calculate process-local bounding cube
    boundingCube = tem_GetLocalBoundingCube( tree )
    ! Exchange with neighbors
    call mpi_reduce( boundingCube(:,1), tBounding, 3, rk_mpi, mpi_min,         &
      &              0, tree%global%comm, iErr  )
    boundingCube(:,1) = tBounding
    call mpi_reduce( boundingCube(:,2), tBounding, 3, rk_mpi, mpi_max,         &
      &              0, tree%global%comm, iErr  )
    boundingCube(:,2) = tBounding

  end function tem_GetRealBoundingCube
! ****************************************************************************** !

end module tem_geometry_module
! ****************************************************************************** !
