! Copyright (c) 2012-2014 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2012-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2012, 2014 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2013-2014 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2013-2016, 2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2013 Monika Harlacher <monika.harlacher@uni-siegen.de>
! Copyright (c) 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2016 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ****************************************************************************** !
!> This module contains the definition of geometry type and routines to 
!! geometry information like mesh, boundary, immersed_boundary and restart
module mus_geom_module

  ! include treelm modules
  use env_module,              only: rk, long_k, globalMaxLevels
  use tem_global_module,       only: tem_global_mesh_read
  use tem_property_module,     only: prp_hasBnd, prp_hasQVal
  use treelmesh_module,        only: treelmesh_type, load_tem
  use tem_bc_prop_module,      only: tem_bc_prop_type, init_tem_bc_prop, &
    &                                load_tem_bc_qVal
  use tem_geometry_module,     only: tem_build_treeToProp_pointer
  use tem_solveHead_module,    only: tem_solveHead_type
  use tem_timeControl_module,  only: tem_timeControl_start_at_sim
  use tem_restart_module,      only: tem_restart_type, tem_load_restart
  use tem_comm_env_module,     only: tem_comm_env_type
  use tem_simControl_module,   only: tem_simControl_type
  use tem_tools_module,        only: tem_horizontalSpacer
  use tem_logging_module,      only: logUnit

  ! include musubi modules
  use mus_geomIncrHead_module, only: mus_geomIncrHead_type, &
    &                                mus_geomIncrHead_load
  use mus_IBM_module,          only: mus_IBM_globType, mus_load_IBM

  implicit none

  private

  public :: mus_geom_type
  public :: mus_build_posInProp
  public :: mus_load_geom
  public :: mus_load_bc_data

  !> Geometric information and definitions
  type mus_geom_type

    !> tree data type
    type( treelmesh_type ) :: tree

    !> boundary information as stored on disk
    type( tem_bc_prop_type ) :: boundary

    !> The header type containing all the geometry increase information
    type( mus_geomIncrHead_type ), allocatable :: geomIncr(:)

    !> Logical to define whether geometry increase is active or not
    logical :: dynamicGeom = .false.

    !> Tree element position in the boundary_ID( nDir, nElems) in bc_prop_type
    !! it has a size of tree%nElems
    !! How to use:
    !! do iElem = 1, tree%nElems
    !!   posInBndID = posInBndID( iElem )
    !!   ! current element has boundary only if posInBndID>0 
    !!   ! else posInBndID = -1
    !!   if (posInBnd > 0 )
    !!     bnd_ID(1:nDir) = bc_prop%boundary_ID( 1:nDir, posInBndID )
    !!   end if
    !! end do
    integer, allocatable :: posInBndID(:)

    !> Tree element position in the qVal( nDir, nElems) in bc_prop_type
    !! it has a size of tree%nElems
    !! How to use:
    !! do iElem = 1, tree%nElems
    !!   posInQVal = posInQVal( iElem )
    !!   ! current element has qVal if posInQVal>0 else posInQVal = -1
    !!   if (posInQVal > 0 )
    !!     qVal(1:nDir) = bc_prop%qVal( 1:nDir, posInQVal )
    !!   end if
    !! end do
    integer, allocatable :: posInQVal(:)

    !> tree element position in level descriptor total list
    !! it has a size of tree%nElems
    !! How to use:
    !! do iElem = 1, tree%nElems
    !!   treeID = tree%treeID( iElem )
    !!   level = tem_levelOf( treeID )
    !!   posInTotal = levelPointer( iElem )
    !!   treeID = LevelDesc( iLevel )%total( posInTotal )
    !! end do
    integer, allocatable :: levelPointer(:)

    !> Boundary element poisition in the levelwise globBC%elemLvl(:)%elem%val
    !! It has a size of geometry%boundary%property%nElems.
    !! It is used in tracking to extract value stored in boundary types.
    !! Hot to use this access normal direction of boundary element:
    !! do iElem = 1, tree%nElems
    !!   level = tem_levelOf( treeID )
    !!   posInBndID = posInBndID(iElem)
    !!   if (posInBndID > 0) then
    !!     BCIDs = bc_prop%boundary_ID(:, posInBndID)
    !!     minBCID = minval(BCIDs, BCIDs > 0)
    !!     posInBcElem = bcLevelPointer(posInBndID)
    !!     normal = globBC%elemLvl(iLevel)%normal(posInBcElem)
    !!   end if
    !! end do
    integer, allocatable :: bcLevelPointer(:)

    !> Minimum bcID for each boundary element.
    !! if a element has more than one boundary then use minBcID which depends
    !! on boundary order in seeder configuration.
    integer, allocatable :: minBcID(:)

    !> immersed boundary data
    type(mus_IBM_globType) :: globIBM

    !> Contains Forces on boundary elements computed using momentum exchange
    !! method. This will be used to derive_bndForce routine to compute force
    !! of certain boundaries.
    !! Forces are stored in level-independent fashion as geometry%boundaryID
    !! loaded from mesh files.
    !! Dim1: geometry%boundaryi%property%nElems
    !! Dim2: 3
    real(kind=rk), allocatable :: bndForce(:,:)
  end type mus_geom_type

contains

  ! ************************************************************************** !
  !> This routine load all geometry related datas like mesh, boundary
  !! and immersed_boundary. Restart is also loaded here because mesh is loaded
  !! in tem_load_restart if restart read is defined.
  subroutine mus_load_geom(me, restart, solverHead, simControl, proc, &
    &                      scaleFactor, initial_balance)
    ! --------------------------------------------------------------------------
    !< contains geometry information which are loaded in this routine
    type(mus_geom_type), intent(out) :: me 
    !> contains restart information
    type(tem_restart_type), intent(out) :: restart
    !> contains general description of the solver including flu_state 
    type(tem_solveHead_type), intent(inout) :: solverHead
    !> contains simulation time control information
    type(tem_simControl_type), intent(inout) :: simControl
    !> contains MPI communication environment 
    type(tem_comm_env_type), intent(in) :: proc
    !> Temporal scaling factor for multilevel mesh 
    integer, intent(in) :: scaleFactor
    !> If true, do initial balancing using level_weights 
    logical, intent(in) :: initial_balance
    ! --------------------------------------------------------------------------
    integer :: iLevel, minLevel, maxLevel
    real(kind=rk) :: level_weights(globalMaxLevels)
    ! --------------------------------------------------------------------------

    ! -------------------------------------------------------------------------
    !  Load mesh                                                             !
    ! -------------------------------------------------------------------------
    ! First check, if we are starting from a restart
    ! KJ: Now the restart is read from initial conditions, which are loaded
    ! much later than the config in the process flow
    call tem_load_restart( me       = restart,            &
      &                    conf     = solverHead%conf(1), &
      &                    tree     = me%tree,            &
      &                    timing   = simControl%now,     &
      &                    globProc = proc                )

    if( restart%controller%readRestart ) then
      ! If there is a restart, the timings in the params type have to be
      ! updated to those read from the restart
      call tem_timeControl_start_at_sim( me  = simControl%timeControl, &
        &                                now = simControl%now          )
    else

      write(logUnit(1),*) 'No read restart given. Loading mesh file'

      ! First load global info
      call tem_global_mesh_read( me     = me%tree%global,     &
        &                        conf   = solverHead%conf(1), &
        &                        myPart = proc%rank,          &
        &                        nParts = proc%comm_size,     &
        &                        comm   = proc%comm           )

      minLevel = me%tree%global%minLevel
      maxLevel = me%tree%global%maxLevel

      if( (minLevel /= maxLevel) .and. initial_balance ) then

        call tem_horizontalSpacer(fUnit = logUnit(1))
        write(logUnit(1),"(A)") 'Do initial level-wise load balancing.'

        ! Calculate level wise weights to scaleFactor^( level - minlevel )
        !                 acoustic      diffusive
        !  ------------------------------------------------
        !  fac =         |   2              4
        !  ------------------------------------------------
        !  minLevel      |   1              1
        !  minLevel + 1  |   2              4
        !  minLevel + 2  |   4              16
        !  ...
        !  ------------------------------------------------
        do iLevel = minLevel, maxLevel
          level_weights( iLevel ) = dble( scaleFactor ** (iLevel-minLevel) )
        end do

        call load_tem( me          = me%tree,              &
          &            conf        = solverHead%conf(1),   &
          &            myPart      = proc%rank,            &
          &            nParts      = proc%comm_size,       &
          &            comm        = proc%comm,            &
          &            levelWeight = level_weights,        &
          &            meshDir     = solverHead%meshFolder )
      else
        ! load the tree from the mesh = '' definition in case no restartRead is
        ! given in this case the mesh from the restart header is read 
        ! (in tem_load_restart)
        call load_tem( me      = me%tree,              &
          &            conf    = solverHead%conf(1),   &
          &            myPart  = proc%rank,            &
          &            nParts  = proc%comm_size,       &
          &            comm    = proc%comm,            &
          &            meshDir = solverHead%meshFolder )
      end if ! minLevel /= maxLevel .and. initial_balance
    end if ! readRestart
    ! ---------------------------------------------------------------------------
    ! Done loading mesh.                                                       !
    ! ---------------------------------------------------------------------------

    ! Load boundary and qval
    call mus_load_bc_data( geometry = me,        &
      &                    rank     = proc%rank, &
      &                    comm     = proc%comm  )

    ! Load solidification/fluidification settings
    call mus_geomIncrHead_load( me          = me%geomIncr,        &
      &                         conf        = solverHead%conf(1), &
      &                         dynamicGeom = me%dynamicGeom      )

    ! load IBM data
    call mus_load_IBM( me     = me%globIBM,         &
      &                conf   = solverHead%conf(1), &
      &                rank   = proc%rank           )

  end subroutine mus_load_geom
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> This routine invokes the treelm routines to load the boundary conditions
  !!
  subroutine mus_load_bc_data( geometry, rank, comm )
    ! --------------------------------------------------------------------------
    type( mus_geom_type ),  intent(inout)  :: geometry !< Treelmesh data
    integer, intent(in) :: rank, comm
    ! --------------------------------------------------------------------------
    integer :: iProp
    ! --------------------------------------------------------------------------

    ! ----------------------Load the boundary conditions------------------------
    ! Boundary conditions are loaded even in the case of dynamic load balancing
    ! when the restart is done at requisite load balancing intervals
    call init_tem_bc_prop( geometry%tree, rank, &
      &                    comm, geometry%boundary )

    do iProp = 1, geometry%tree%global%nProperties
      if( geometry%tree%global%property( iProp )%nElems > 0 ) then
        select case( geometry%tree%global%property( iProp )%bitpos )
        case( prp_hasBnd )
          ! Already loaded
        case( prp_hasQVal )
          ! Load qVal from disk
          ! prp_hasQVal is the 2nd property in mesh
          write(logUnit(1),"(A)") 'Loading qVal data from directory: '&
            &                     //trim(geometry%tree%global%dirname)
          call load_tem_BC_qVal(                                               &
            &           me       = geometry%boundary,                          &
            &           offset   = geometry%tree%Property(iProp)%Offset,       &
            &           nElems   = geometry%tree%Property(iProp)%nElems,       &
            &           basename = trim(geometry%tree%global%dirname)//'qval', &
            &           mypart   = rank,                   &
            &           comm     = comm                    )
          write(logUnit(1),*)'Done, reading the qVal!'

        end select ! property( iProp )%bitpos
      endif ! property( iProp )%nElems > 0
    enddo ! iProp

    call mus_build_posInProp( geometry )

    ! when qVal exist, it is allocated inside load_tem_BC_qVal
    ! otherwise allocate it with size 0
    if ( .not. allocated( geometry%boundary%qVal ) ) then
      allocate( geometry%boundary%qVAl(0,0) )
    end if

  end subroutine mus_load_bc_data
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This routine builds mapping between elements in tree to to propery list
  subroutine mus_build_posInProp( me )
    ! --------------------------------------------------------------------------
    type( mus_geom_type ), intent(inout) :: me !< Treelmesh data
    ! --------------------------------------------------------------------------
    integer :: iProp
    ! --------------------------------------------------------------------------

    do iProp = 1, me%tree%global%nProperties
      if( me%tree%global%property( iProp )%nElems > 0 ) then
        select case( me%tree%global%property( iProp )%bitpos )
        case( prp_hasBnd )
          if (  allocated( me%posInBndID )) deallocate( me%posInBndID )
          allocate( me%posInBndID(me%tree%nElems) )
          ! build me%posInBnd which maps tree to boundary%boundary_ID
          call tem_build_treeToProp_pointer( &
            &        treeToProp = me%posInBndID, &
            &        nElems     = me%tree%nElems,&
            &        ElemPropertyBits = me%tree%ElemPropertyBits,&
            &        prp_bit    = prp_hasBnd           )

          call mus_build_minBcID( minBcID    = me%minBcID,   &
            &                     bc_prop    = me%boundary,  &
            &                     posInBndID = me%posInBndID )

        case( prp_hasQVal )
          if ( allocated( me%posInQval )) deallocate( me%posInQval )
          allocate( me%posInQval(me%tree%nElems) )
          ! build me%posInQVal which maps tree to boundary%qVal
          call tem_build_treeToProp_pointer( &
            &        treeToProp = me%posInQVal, &
            &        nElems     = me%tree%nElems,&
            &        ElemPropertyBits = me%tree%ElemPropertyBits,&
            &        prp_bit    = prp_hasQVal         )

        end select ! property( iProp )%bitpos
      endif ! property( iProp )%nElems > 0
    enddo ! iProp

  end subroutine mus_build_posInProp
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This routine build minBcID for boundary elements, it is required if a
  !! element has more than one boundary in its directions.
  !! if a element has more than one boundary then use minBcID which depends
  !! on boundary order in seeder configuration
  subroutine mus_build_minBcID(minBcID, bc_prop, posInBndID)
    ! --------------------------------------------------------------------------
    integer, allocatable, intent(out) :: minBcID(:)
    !> boundary information from mesh
    type(tem_bc_prop_type), intent(in) :: bc_prop
    !> tree element position in boundaryID
    integer, intent(in) :: posInBndID(:)
    ! --------------------------------------------------------------------------
    integer :: iElem
    integer(kind=long_k) :: bcIDs(bc_prop%nSides)
    ! --------------------------------------------------------------------------
    allocate(minBcID(bc_prop%property%nElems))

    do iElem = 1, bc_prop%property%nElems
      bcIDs = bc_prop%boundary_ID(:, iElem)
      minBcID(iElem) = int(minval(bcIDs, bcIDs > 0_long_k))
    end do
  end subroutine mus_build_minBcID
  ! ************************************************************************** !

end module mus_geom_module
! ****************************************************************************** !
