! Copyright (c) 2012, 2014-2017 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2012-2014 Simon Zimny <s.zimny@grs-sim.de>
! 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-2016 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2012-2013 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016-2017 Raphael Haupt <raphael.haupt@uni-siegen.de>
! Copyright (c) 2018 Daniel Fleischer <daniel.fleischer@student.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.
?? include 'header/lbm_macros.inc'
! ****************************************************************************** !
!> author: Kartik Jain
!! author: Jiaxing Qi
!! This module provides functionality to perform geometry change  
!! depending on the flow properties. Fluidify/Solidify, proximity
!! condition variable are supplied in the lua
!! file. Checks are performed against those parameters and geometry is 
!! accordingly changed (an element solidified or fluidified)
!!
!> Geometry Increment module
!! 
module mus_geomIncr_module

  ! include treelm modules
  use env_module,              only: rk, long_k, io_buffer_size, LabelLen
  use tem_topology_module,     only: tem_LevelOf
  use tem_property_module,     only: prp_fluidify, prp_solid, prp_fluid, &
    &                                prp_noSolidification
  use tem_comm_module,         only: tem_commPattern_type
  use tem_timeControl_module,  only: tem_timeControl_check
  use tem_logging_module,      only: logUnit
  use tem_general_module,      only: tem_general_type
  use tem_dyn_array_module,    only: PositionOfVal
  use tem_varSys_module,       only: tem_varSys_type
  use tem_construction_module, only: tem_levelDesc_type
  use treelmesh_module,        only: treelmesh_type
  use tem_aux_module,          only: tem_abort

  ! include musubi modules
  use mus_geomIncrHead_module, only: mus_geomIncrHead_type
  use mus_geom_module,         only: mus_geom_type
  use mus_scheme_type_module,  only: mus_scheme_type
  use mus_time_module,         only: mus_timeControl_homogenize
      


  implicit none

  private

  public :: mus_geomIncr, proximity
  public :: mus_init_geomIncr
  public :: prepare_target_elem_list
  ! public :: update_connectivity

contains

! ****************************************************************************** !
  !> This subroutine initializes the geometry increment.
  !!
  subroutine mus_init_geomIncr(me, varSys, dt, reqInt)
    ! ---------------------------------------------------------------------------
    !> geometry increment type
    type( mus_geomIncrHead_type ), intent(inout) :: me(:)
    !> Global variable system
    type( tem_varSys_type ), intent(in)          :: varSys
    !> dt of maxlevel or smallest dt
    real(kind=rk), intent(in) :: dt
    !> Required interval, in which the update MUST occur.
    !! This is required for the musubi multilevel, where the time step should
    !! only be determined active, when the end of the largest cycle is reached.
    integer, intent(in)                          :: reqInt
    ! ---------------------------------------------------------------------------
    integer :: iGInc
    ! ---------------------------------------------------------------------------
    do iGInc = 1, size(me)
      ! Initialize the timer
      call mus_timeControl_homogenize( me     = me( iGInc )%timeControl, &
        &                              dt     = dt,                      &
        &                              reqInt = reqInt                   )

      me(iGInc)%cond_varPos = PositionOfVal(me  = varSys%varname,             &
        &                                   val = trim(me(iGInc)%cond_varName))
      if (me(iGInc)%cond_varPos < 1) then 
        write(logUnit(1),*) 'Error: variable '//trim(me(iGInc)%cond_varName) &
          &              // ' is not added to variable system'
        call tem_abort()  
      end if

    end do

  end subroutine mus_init_geomIncr
! ******************************************************************************!


! ******************************************************************************!
  !> This subroutine checks for various conditions defined in the geomIncr table
  !! within the lua file, calculates the requested macroscpoic variables and 
  !! then compares them against the specified threshold. 
  !! Accordingly then solidification or fluidification of elements is performed.
  !!
  subroutine mus_geomIncr( geometry, scheme, commPattern, general)
    ! ---------------------------------------------------------------------------
    !>
    type( mus_geom_type )                     :: geometry
    !>
    type( mus_scheme_type ), intent(inout)    :: scheme
    !>
    type( tem_commPattern_type ), intent(in)  :: commPattern
    !> global parameters
    type( tem_general_type ),intent(in)       :: general
    ! ---------------------------------------------------------------------------
    integer :: pot_nElems ! potential elements for gi temporarily
    integer :: nElems_gi  ! number of elements which participate in geomincr
    integer :: nElems ! Total number of elements in the tree
    integer :: varPos
    ! levelPointer list of element to be solidify
    integer, allocatable :: LPlist(:), pntTreeID(:)
    ! List of target iElems
    integer, allocatable :: target_ielem_list(:)
    real(kind=rk), allocatable :: chunk(:)
    integer :: iGInc, iChunk
    integer :: nComp, chunkSize, nChunks
    integer :: nChunkElems, iLevel
    integer :: elemoff ! starting index of iElem for different chunks
    real(kind=rk), allocatable :: chunk_res(:)
    integer :: buf_start, buf_end
    integer :: res_size
    logical :: triggered
    ! ---------------------------------------------------------------------------

    ! The number of elements in the tree
    nElems = geometry%tree%nElems

    chunkSize = io_buffer_size

    allocate( LPlist( size(geometry%tree%treeID)) )
    allocate( pntTreeID( size(geometry%tree%treeID)) )
    allocate( target_ielem_list( size(geometry%tree%treeID)) )

    allocate(chunk(io_buffer_size))

    ! Loop over the number of geomIncrs defined in the Lua file
    do iGInc = 1, size( geometry%geomIncr )
      call tem_timeControl_check(                                         &
        &               me        = geometry%geomIncr(iGInc)%timeControl, &
        &               now       = general%simControl%now,               &
        &               comm      = general%proc%comm,                    &
        &               triggered = triggered                             )

      ! time_modulo is to do only when one complete cyle is done with multilevel
      ! when read from restart trigger is active at 1st time step
      if( triggered ) then

        ! Initialize the total number of target elements to zero
        nElems_gi = 0
        ! Call the proximity routine to generate
        ! a list of treeIDs which pass the proximity check
        call proximity( scheme       = scheme,      &
          &             geometry     = geometry,    &
          &             iGInc        = iGInc,       &
          &             pot_nElems   = pot_nElems,  &
          &             pot_lp       = LPlist(:),   &
          &             pot_lp_tree  = pntTreeID(:) )

        nChunks = ceiling( real(pot_nElems, kind=rk) &
          &              / real(chunkSize, kind=rk) )

        do iChunk = 1, nChunks
          elemOff = ( (iChunk-1)*chunkSize )

          ! number of elements written to THIS chunk
          nChunkElems = min(chunkSize, pot_nElems-elemOff)

          ! Compute the element lower and upper bound for the current chunk
          buf_start = (iChunk-1)*chunkSize+1
          buf_end   = (iChunk-1)*chunkSize+nChunkElems

          ! position of condition variable in varSys
          varPos = geometry%geomIncr(iGInc)%cond_varPos
          ! get the no. of components for each variable system
          nComp = scheme%varSys%method%val(varPos)%nComponents

          ! allocate local array where the results are placed
          res_size = nComp*nChunkElems
          allocate( chunk_res(res_size) )

          ! Use the function pointers to calculate variables identified 
          ! from Lua results are placed in array chunk_res.
          ! Variable returns 1 if condition satisfies else 0
          call scheme%varSys%method%val(varpos)%get_element(         &
            &               varSys  = scheme%varSys,                 &
            &               elemPos = pntTreeID(buf_start:buf_end),  &
            &               time    = general%simControl%now, &
            &               tree    = geometry%tree,                 &
            &               nElems  = nChunkElems,                   &
            &               nDofs   = 1,                             &
            &               res     = chunk_res(1:res_size)          )

          ! Prepare the list of target elements and set the property bits 
          ! of these elements
          call prepare_target_elem_list(                    &
            & nComp             = nComp,                    &
            & elemoff           = elemoff,                  &
            & nChunkElems       = nChunkElems,              &
            & res               = chunk_res,                &
            & LPList            = LPList,                   &
            & nElems_gi         = nElems_gi,                &
            & pntTreeID         = pntTreeID,                &
            & target_ielem_list = target_ielem_list,        &
            & scheme            = scheme,                   &
            & tree              = geometry%tree,            &
            & geomIncr          = geometry%geomIncr(iGInc)  )

          ! Deallocate the array containing computed results
          deallocate( chunk_res )
        end do ! iChunk

        if(nElems_gi > 0) then
          if (geometry%geomIncr(iGInc)%solidify) then
            write(logUnit(1),*)'Elements found for solidification:', nElems_gi
          elseif (geometry%geomIncr(iGInc)%fluidify) then
            write(logUnit(1),*)'Elements found for fluidification:', nElems_gi
          end if
          nElems_gi = 0
          ! and update the flag meshChange because something changed
          geometry%tree%global%meshChange = .true.
        end if

        do iLevel = geometry%tree%global%minLevel, geometry%tree%global%maxLevel
          ! Communicate the property bits of halo elements
          call commPattern%exchange_long(                            &
            &  send         = scheme%levelDesc( iLevel )%sendbuffer, &
            &  recv         = scheme%levelDesc( iLevel )%recvbuffer, &
            &  state        = scheme%levelDesc( iLevel )%property,   &
            &  message_flag = iLevel,                                &
            &  comm         = general%proc%comm                    )

          ! Invoke the update connectivity routine which updates the 
          ! neighbor array of target elements and the neighbor array of 
          ! their neighbors. 
          call update_connectivity( scheme, iLevel)
        end do
      end if ! Timing check
    end do ! iGInc
    ! Deallocate the chunk array
    deallocate( chunk )
    deallocate( LPlist )
    deallocate( target_ielem_list )
  end subroutine mus_geomIncr
! ****************************************************************************** !


! ****************************************************************************** !
  !> A subroutine which checks the proximity condition and generates a new 
  !! treeID list of potential target elements which can then be checked against
  !! various thresholds as defined in the lua configuration file.
  !!
  !! Proximity condition means that an element can be solidified only if it has 
  !! at least one wall as a neighboring element.
  !! In biological clotting models it ensures that no isolated clot is produced
  !! in the fluid domain.
  !! Furthermore the noSolidification property is checked for every element.
  !! It is used to define spaces where solidification is not allowed.
  !!
  subroutine proximity( scheme, geometry, iGInc, pot_nELems, pot_lp,      &
    &                   pot_lp_tree)
    ! ---------------------------------------------------------------------------
    !> scheme information
    type(mus_scheme_type), intent(inout)  :: scheme
    !> scheme information
    type(mus_geom_type), intent(in)       :: geometry
    !> Geomtry icrement
    integer, intent(in)                   :: iGInc
    !> number of potential elements found
    integer, intent(inout)                :: pot_nElems
    !> level pointer of potential elements
    integer, intent(inout)                :: pot_lp(:)
    !> position of potential elements in the treeID list
    integer, intent(inout)                :: pot_lp_tree(:)
    ! ---------------------------------------------------------------------------
    ! local variables
    integer :: iElem, iLevel, iField, iDir
    integer :: neighPos, elemPos 
    integer :: QQN, QQ, neighbor, nScalars
    integer(kind=long_k) :: elemProp
    ! ---------------------------------------------------------------------------

    pot_nElems = 0
    pot_lp = 0
    pot_lp_tree = 0
    QQN = scheme%layout%fStencil%QQN
    QQ  = scheme%layout%fStencil%QQ
    nScalars = scheme%varSys%nScalars

    if ( geometry%geomIncr(iGInc)%proximity ) then
      do iElem = 1, size( geometry%tree%treeID )
        elemProp = geometry%tree%ElemPropertyBits( iElem )
        ! Only consider element if solidification is allowed
        if ( .NOT. btest ( elemProp, prp_noSolidification ) ) then
          do iField = 1, scheme%nFields
            iLevel = tem_levelOf( geometry%tree%treeID( iElem ))
            elemPos = geometry%levelPointer( iElem )
            do iDir = 1, QQN
              ! Get the neighbor 
              neighbor = scheme%pdf( iLevel )%neigh(               &
                & ?NGPOS?(iDir, elemPos, scheme%pdf(iLevel)%nSize) )
              ! position in total list
              neighPos = ?ElemIDX?( neighbor, nScalars, scheme%pdf( iLevel )%nSize )

              !> Acc the proximity condition, as soon as its found that any of the
              !! neighbors of current element is solid, that element is stored in
              !! the new list of potential elements. Loop is terminated for this
              !! element and check is performed for other element.
              if (neighPos == elemPos) then
                pot_nElems              = pot_nElems + 1
                pot_lp(pot_nElems)      = elemPos
                pot_lp_tree(pot_nElems) = iElem
                exit
              end if
            end do ! iDir
          end do ! iField
        end if
      end do ! iElem
    else
      do iElem = 1, size( geometry%tree%treeID )
        elemProp = geometry%tree%ElemPropertyBits( iElem )
        if ( .NOT. btest ( elemProp, prp_noSolidification ) ) then
          elemPos = geometry%levelPointer( iElem )
          pot_nElems              = pot_nElems + 1
          pot_lp(pot_nElems)      = elemPos
          pot_lp_tree(pot_nElems) = iElem
        end if
      end do
    end if

  end subroutine proximity
! ****************************************************************************** !


! ****************************************************************************** !
  !> This routine compares the macroscopic quantity obtained for an element
  !! against specified threshold and populates the list of target elements 
  !! which participate in geometry change. In addition it sets the property
  !! bit(s) of target elements to fluidify in case of solidification and clears
  !! this bit in case of fluidification
  !!
  subroutine prepare_target_elem_list( nComp, elemoff, nChunkElems, res, &
    &                                  LPList, nElems_gi, pntTreeID,     &
    &                                  target_ielem_list, scheme, tree,  &
    &                                  geomIncr                          )
    ! ---------------------------------------------------------------------------
    !> total number of components
    integer, intent(in)    :: nComp
    !> offset of elements in chunk
    integer, intent(in)    :: elemoff
    !> number of elements in chunk
    integer, intent(in)    :: nChunkElems
    !> result from condition variable
    real(kind=rk), intent(in) :: res(:)
    !> list of level pointers
    integer, intent(in)    :: LPList(:)
    !> number of final elements
    integer, intent(inout) :: nElems_gi
    !> element position in global treeID list tree%treeID
    integer, intent(in)    :: pntTreeID(:) 
    !> target element list for second iGInc to reduce computation effort
    integer, intent(inout) :: target_ielem_list(:) 
    !> scheme information
    type(mus_scheme_type), intent(inout)  :: scheme
    !> fluid tree from mesh, element property bits are set to solidify
    !! if condition variable is 1
    type( treelmesh_type ), intent(inout)  :: tree
    !> current geometry increment 
    type( mus_geomIncrHead_type ), intent(inout) :: geomIncr
    ! ---------------------------------------------------------------------------
    ! local variables
    integer                :: local_elem
    integer                :: iElem, iLevel, elemPos, neighPos, iField, iDir
    integer(kind=long_k)   :: neighProp
    logical                :: isSolid, cond
    ! ---------------------------------------------------------------------------

    ! Loop over all elements in chunks for the condition variable in var system
    ! For solidify (solidify=true), fluid elements will be solidified if the 
    ! condition is satisfied and for fluidify (fluidify=true) solid elements
    ! will be fluidified.
      local_elem = 0
      do iElem = elemoff+1, elemoff+nChunkElems
        local_elem = local_elem + 1
        elemPos = LPlist( iElem ) ! elem position in total list
        iLevel = tem_LevelOf( tree%treeID( elemPos ))
        if( geomIncr%solidify ) then
          ! if condition variable is a vector then all index must be satisfied 
          cond = all(res( (local_elem-1)*nComp + 1 : local_elem*nComp) > 0.0_rk)
          if ( cond ) then
            ! Update the property bit 
            call mus_setProp( property = scheme%levelDesc(iLevel) &
              &                          %property(elemPos),      &
              &               geomIncr = geomIncr                 )
            ! and the ones in the global tree
            call mus_setProp( property = tree%ElemPropertyBits &
              &                          (elemPos),            &
              &               geomIncr = geomIncr              )
            nElems_gi = nElems_gi + 1
          end if ! condition satisfies
        elseif ( geomIncr%fluidify ) then
          ! if condition variable is a vector then all index must be satisfied 
          cond = all(res( (local_elem-1)*nComp + 1 : local_elem*nComp) > 0.0_rk)
          if ( cond ) then
            ! Get the direct neighbor elements for fluidification
            do iField = 1, scheme%nFields
              do iDir = 1, 6
                neighPos = scheme%levelDesc( iLevel )%neigh(1)%               &
                  &        nghElems( ?NgDir?( iDir, scheme%layout%fStencil ), &
                  &                  elemPos )
                ! neighPos < 0 means boundary in that direction
                if ( neighPos > 0 ) then
                  neighProp = scheme%levelDesc(iLevel)%property( neighPos )
                  
                  isSolid = btest( neighProp, prp_Solid )

                  ! Only need to fluidify if neighbor is solid
                  if( isSolid ) then
                    ! Update the property bit 
                    call mus_setProp( property = scheme%levelDesc(iLevel) &
                      &                          %property(neighPos),     &
                      &               geomIncr = geomIncr                 )
                    ! and the ones in the global tree
                    call mus_setProp( property = tree%ElemPropertyBits  &
                      &                          (neighPos),            &
                      &               geomIncr = geomIncr               )
                    nElems_gi = nElems_gi + 1
                    ! As soon as one solid neighbor is fluidifed 
                    ! dir-loop is terminated. This will keep a balance
                    ! between soldification and fluidification.
          !!          EXIT
                  end if
                end if
              end do ! iDir
            end do ! iField
          end if ! condition satisfies
        end if ! is my target
      end do ! iElem

! DF: Following code is only used to improve performance.
!     Still needs a couple of changes to run with 
!     solidification AND fluidification.

!!    else ! iVar != 1 or iDepend != 1
!!      ! loop over the tree ID list which was found for the first geomIncr loop. 
!!      ! Doing this we ensure the AND condition among 
!!      ! different variables within a depend table and as the number of target 
!!      ! elements is hardly 5% of the total geometry, it improves computational
!!      ! efficiency of the code.
!!      local_elem = 0
!!      do iElem = 1, nElems_gi
!!        elemPos = LPlist( target_ielem_list( iElem ))
!!        iLevel = tem_LevelOf( tree%treeID( &
!!          &                   pntTreeID( target_ielem_list(iElem) ) ) )
!!        elemProp = levelDesc( iLevel )%property( elemPos )
!!        isSolid = btest( elemProp, prp_fluidify )
!!        if( isSolid .neqv. geomIncr%solidify ) then
!!          ! if condition variable is a vector then all index must be satisfied 
!!          cond = all(res( (target_ielem_list(iElem)-1)*nComp + 1 &
!!            &            : target_ielem_list(iElem)*nComp) > 0.0_rk)
!!          
!!          ! Update the local tree ID list
!!          if ( cond ) then
!!            ! Update the property bit
!!            call mus_setSolidProp(                                    &
!!              &       property = levelDesc(iLevel)%property(elemPos), &
!!              &       solidify = geomIncr%solidify                    )
!!            ! and the ones in the global tree
!!            call mus_setSolidProp(                             &
!!              &   property = tree%ElemPropertyBits(pntTreeID(iElem)), &
!!              &   solidify = geomIncr%solidify                        )
!!            ! and update the flag meshChange because something changed
!!            ! geometry%tree%global%meshChange = .true.
!!            local_elem = local_elem + 1
!!            target_ielem_list(local_elem) = target_ielem_list(iElem)
!!          end if ! condition satisfies 
!!        end if ! is my target
!!      end do ! iElem
!!      ! number of target elements is updated
!!      nElems_gi = local_elem
!!    end if ! nElems_gi == 0

  end subroutine prepare_target_elem_list
! ****************************************************************************** !


! ****************************************************************************** !
  !> Construct the propagation stencil for each element
  !!
  !! the bounce back rules have to be applied here
  !!
  subroutine update_connectivity( scheme, iLevel)
    ! -------------------------------------------------------------------------- 
    !> scheme information
    type( mus_scheme_type ), intent(inout) :: scheme
    !> current level
    integer, intent(in)     :: iLevel
    ! -------------------------------------------------------------------------- 
    ! defining local variables
    integer :: iElem  ! element counter for current element in treeID list
    integer :: iField ! field counter
    integer :: idx_idir ! current direction in global varsys
    integer :: QQ
    integer :: iDir ! counter for current direction in iElem
    integer :: invDir ! inverse direction for iDir
    ! result: from which direction to take pdf from (inverse for bounceback) 
    integer :: GetFromDir
    ! result: from which element to take the pdf from (local for bounceback) 
    integer :: GetFromPos
    integer :: neighPos      ! position of current neighElem in treeID list
    integer(kind=long_k)  :: neighProp, elemProp
    integer :: nElems, nSize
    ! ---------------------------------------------------------------------------
    QQ = scheme%layout%fStencil%QQ

    nElems = scheme%pdf( iLevel )%nElems_local
    nSize  = scheme%pdf( iLevel )%nSize
    do iElem = 1, nElems
      elemProp = scheme%levelDesc( iLevel )%property( iElem )

      ! Loop over every link
      do iDir  = 1, QQ
        if ( iDir == QQ  ) then
          ! Rest density takes from itself, iElem and iDir
          GetFromPos = iElem
          GetFromDir = iDir
        else ! iDir not rest density
          ! For push we need to find different neighbors than for pull
          ! -> ngDir
          neighPos = scheme%levelDesc( iLevel )%neigh(1)%               &
            &        nghElems( ?NgDir?( iDir, scheme%layout%fStencil ),    &
            &                  iElem )

          neighProp = 0_long_k
          ! get neighbor's property
          if( neighPos > 0 )                                                 &
            & neighProp = scheme%levelDesc( iLevel )%property( neighPos )

          if( (neighPos <= 0) .or. (btest(neighProp, prp_fluidify))        &
              &  .or. (btest(elemProp, prp_fluidify))) then
            ! 1. Fluid element has boundary in current link direction
            ! 2. neighbor element has fluidify property (still solid now)
            ! Get inverse direction of LOCAL element
            invDir = scheme%layout%fStencil%cxDirInv( iDir )
            GetFromDir = invDir
            GetFromPos = iElem
            ! I set the source to the target link for halo cells 
            ! with no neighbor in that direction 
          else
            ! neighElem = scheme%levelDesc( iLevel )%total( neighPos )
            GetFromDir = iDir 
            GetFromPos = neighPos
          endif ! has no boundary -> regular treatment
        endif ! not the 0 velocity
 
        scheme%pdf( iLevel )%neigh( ?NGPOS?( iDir, iElem, nElems )) =    &
          &  ?IDX?( GetFromDir, GetFromPos, scheme%varSys%nScalars, nSize )

        do iField = 1, scheme%nFields
          idx_idir = scheme%varSys%method%val(                      &
            &               scheme%stateVarMap%varPos%val(iField) ) &
            &                     %state_varPos(iDir)
!!          ! Set all the pdf links of target element to an infinitesimally
!!          ! small value. This will prevent density from becoming NaN if pdfs
!!          ! are set to zero and besides changed elements will have a very
!!          ! small density which will help in visualization and differentiate
!!          ! them from other fluid elements
!!          if (btest(elemProp, prp_fluidify) ) then
!!             scheme%state(iLevel)%val(                                         &
!!               &   ?IDX?(idx_idir, iElem, scheme%varSys%nScalars, nSize),:) =  &
!!               &   0.000000001_rk
!!          end if
        end do

      enddo ! neighLoop
    enddo ! elemLoop

  end subroutine update_connectivity
! ****************************************************************************** !

! ****************************************************************************** !
  !> This routine updates the propertybits of an element
  !! If solidify == true, set prp_fluidify and set prp_solid
  !! If fluidify == treu, clean prp_fluidify and set prp_fluid
  subroutine mus_setProp( property, geomIncr )
    ! ---------------------------------------------------------------------------
    integer( kind=long_k ) :: property
    type( mus_geomIncrHead_type ), intent(in) :: geomIncr
    ! ---------------------------------------------------------------------------

    if( geomIncr%solidify ) then
      ! set solid property and clean fluid property
      property = ibset( property, prp_fluidify )
      property = ibset( property, prp_solid )
      property = ibclr( property, prp_fluid )
    elseif( geomIncr%fluidify ) then
      ! set fluid property and clean solid property
      property = ibclr( property, prp_fluidify )
      property = ibset( property, prp_fluid )
      property = ibclr( property, prp_solid )
    endif

  end subroutine mus_setProp
! ****************************************************************************** !


end module mus_geomIncr_module
! ****************************************************************************** !
