! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016-2017, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2016, 2019-2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016-2017 Raphael Haupt <raphael.haupt@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.
! **************************************************************************** !
!> author: Kannan Masilamani, Verena Krupp
!! This module contains solver data type with pointers to musubi
!! [[mus_scheme_type]], [[mus_param_type]], etc which are required for
!! variable operation method data.
!! Also contains all general routines for the variable system.
!!
! Copyright (c) 2013 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2013-2014 Nikhil Anand <nikhil.anand@uni-siegen.de>
! Copyright (c) 2014, 2016 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2015, 2018, 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@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 COPYRIGHT HOLDERS AND CONTRIBUTORS "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 THE COPYRIGHT HOLDER 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.
!--------------------------------------------
!    A O S - Array of structures layout new
!-------------------------------------------
! Access to get_point value output
! Access to get_element value output
module mus_varSys_module

  use, intrinsic :: iso_c_binding,  only: c_ptr, c_f_pointer, c_loc
  use env_module,                   only: rk

  use tem_logging_module,           only: logUnit
  use tem_debug_module,             only: dbgUnit
  use tem_aux_module,               only: tem_abort
  use tem_float_module,             only: operator(.fne.)
  use tem_varSys_module,            only: tem_varSys_type,                    &
    &                                     tem_varSys_op_type,                 &
    &                                     tem_varSys_solverData_evalElem_type
  use tem_pointData_module,         only: tem_pointData_list_type
  use treelmesh_module,             only: treelmesh_type
  use tem_geometry_module,          only: tem_CoordOfReal, &
    &                                     tem_PosofId, tem_BaryOfId
  use tem_time_module,              only: tem_time_type
  use tem_topology_module,          only: tem_IdOfCoord
  use tem_topology_module,          only: tem_levelOf
  use tem_operation_var_module,     only: tem_varSys_op_data_type
  use tem_spacetime_fun_module,     only: tem_st_fun_listelem_type
  use tem_element_module,           only: eT_fluid
  use tem_stencil_module,           only: tem_stencilHeader_type

  ! include musubi modules
  use mus_scheme_type_module,  only: mus_scheme_type
  use mus_physics_module,      only: mus_physics_type
  use mus_geom_module,         only: mus_geom_type


  implicit none

  private

  public :: mus_varSys_solverData_type
  public :: mus_varSys_data_type
  public :: mus_get_new_solver_ptr
  public :: mus_init_varSys_solverData
  public :: mus_deriveVar_forPoint
  public :: mus_generic_varFromPDF_fromIndex
  public :: mus_set_stFun_getElement
  public :: mus_derive_fromPDF
  public :: mus_generic_fromPDF_forElement

  !> Method data container for every variable
  type mus_varSys_data_type
    ! type with all solver relevant data for the variable
    type(mus_varSys_solverData_type), pointer :: solverData

    ! the point_datas need to be stored levelwise
    type(tem_pointData_list_type) :: pointData

    !> data array for operation or derived varibales
    !! consists the index arrys for points stored in the
    !! poingtData of input variable
    !! size is number of input variables
    type(tem_varSys_op_data_type) :: opData
  end type mus_varSys_data_type

  !> Contains pointer to musubi data types required for variable operation
  !! method_data
  type mus_varSys_solverData_type
    !> scheme data type
    type(mus_scheme_type), pointer :: scheme => NULL()

    !> contains basic SI units to convert from lattice to physical and
    !! vice versa
    type(mus_physics_type), pointer :: physics => NULL()

    !> Contains geometry information and definitions
    type(mus_geom_type), pointer :: geometry => NULL()

  end type mus_varSys_solverData_type

  abstract interface
    !> This interface describes the arguments to be used for routines that do
    !! the derivation of variables from the variable system.
    subroutine mus_derive_fromPDF( fun, varSys, stencil, iLevel, posInState, &
      &                            pdf, res, nVals )
      import :: tem_varSys_op_type,     &
        &       tem_varSys_type,        &
        &       tem_stencilHeader_type, &
        &       rk
      !> Description of the method to obtain the variables,
      class(tem_varSys_op_type), intent(in) :: fun
      !> The variable system to obtain the variable from.
      type(tem_varSys_type), intent(in) :: varSys
      !> fluid stencil defintion
      type(tem_stencilHeader_type), intent(in)  :: stencil
      !> current level
      integer, intent(in) :: iLevel
      !> Position of element in levelwise state array
      integer, intent(in) :: posInState(:)
      !> pdf
      real(kind=rk), intent(in) :: pdf(:)
      !> results
      real(kind=rk), intent(out) :: res(:)
      !> nVals to get
      integer, intent(in) :: nVals
    end subroutine mus_derive_fromPDF
  end interface

contains


  ! ************************************************************************* !
  !> Routine to get a pointer to a new instance of mus_varSys_solverData_type
  !! to be used as method data for a variable in the variable system.
  !!
  !! A new instance is allocated and a c_ptr to this type is returned.
  !! Be aware that local pointer are not automatically deallocated when leaving
  !! the routine
  function mus_get_new_solver_ptr( solver ) result(resPtr)
    ! --------------------------------------------------------------------- !
    !> The prototype is used to initialize the new instance.
    type(mus_varSys_solverData_type), intent(in), optional, target :: solver
    !> Pointer to the newly created instance.
    type(c_ptr) :: resPtr
    ! --------------------------------------------------------------------- !
    !> Local variable to allocate a new instance.
    type(mus_varSys_data_type), pointer :: res
    ! --------------------------------------------------------------------- !

    allocate(res)
    if (present(solver)) res%solverData => solver

    !>TODO init the point arrays

    resPtr = c_loc(res)

  end function mus_get_new_solver_ptr
  ! ************************************************************************* !


  ! ************************************************************************* !
  ! This routine sets varSys_solverData pointers
  subroutine mus_init_varSys_solverData( me, scheme, physics, geometry )
    ! --------------------------------------------------------------------- !
    type(mus_varSys_solverData_type),  intent(out) :: me
    type(mus_scheme_type),  target, intent(in) :: scheme
    type(mus_physics_type), target, intent(in) :: physics
    type(mus_geom_type),    target, intent(in) :: geometry
    ! --------------------------------------------------------------------- !
    me%scheme => scheme
    me%physics => physics
    me%geometry => geometry
  end subroutine mus_init_varSys_solverData
  ! ************************************************************************* !


  ! ************************************************************************ !
  !> Routine to store musubi varSys Data in stFun variable
  subroutine mus_set_stFun_getElement(solData_evalElem, fun)
    ! -------------------------------------------------------------------- !
    !> Description on how to set the element retrieval function for stfuns.
    class(tem_varSys_solverData_evalElem_type), intent(in) :: solData_evalElem

    !> Description of the method to obtain the variables, here some preset
    !! values might be stored, like the space time function to use or the
    !! required variables.
    type(tem_varSys_op_type), intent(inout) :: fun
    ! -------------------------------------------------------------------- !
    type(tem_st_fun_listElem_type), pointer :: fptr
    type(mus_varSys_solverData_type), pointer :: fSDptr
    ! -------------------------------------------------------------------- !

    write(logunit(10),*) "Setting different solver_bundle and" &
      & // " get_element routine for variable at position ",   &
      & fun%myPos
    call C_F_Pointer(fun%method_data, fptr)
    call c_f_pointer(solData_evalElem%solver_bundle, fSDptr)
    fptr%solver_bundle = mus_get_new_solver_ptr( fSDptr )

  end subroutine mus_set_stFun_getElement
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Derive variable for a given set of points using linear interpolation.
  !! This is a generic routine for any variable.
  !! Limitation: If neighbor is halo element then its not considered for
  !! interpolation, only the fluid (non-ghost) elements in the local process
  !! are used for interpolation.
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_point]].
  !!
  recursive subroutine mus_deriveVar_forPoint(fun, varsys, point, time, tree, &
    &                                         nPnts, res                      )
    !--------------------------------------------------------------------- !
    !> Description of the method to obtain the variables, here some preset
    !! values might be stored, like the space time function to use or the
    !! required variables.
    class(tem_varSys_op_type), intent(in) :: fun

    !> The variable system to obtain the variable from.
    type(tem_varSys_type), intent(in) :: varSys

    !> Three-dimensional coordinates at which the variable should be
    !! evaluated. Only useful for variables provided as space-time functions.
    real(kind=rk), intent(in) :: point(:,:)

    !> Point in time at which to evaluate the variable.
    type(tem_time_type), intent(in)  :: time

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

    !> Number of values to obtain for this variable (vectorized access).
    integer, intent(in) :: nPnts

    !> Resulting values for the requested variable.
    !!
    !! Dimension: n requested entries x nComponents of this variable
    !! Access: (iElem-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! -------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    type(mus_scheme_type), pointer :: scheme
    integer :: statePos, iPnt, iLevel, elemPos
    integer :: coord(4)
    real(kind=rk) :: bary(3)
    integer :: iNeigh, neighPos, neighLvlPos
    integer, allocatable :: srcElemPos(:)
    real(kind=rk), allocatable :: weights(:), srcRes(:), pntVal(:)
    real(kind=rk) :: dist(3), dist_len, sumWeight
    integer :: iSrc, nSrcElems
    ! -------------------------------------------------------------------- !
!write(dbgUnit(1),*) 'Derive for point :'//trim(varSys%varname%val(fun%myPos))
    call C_F_POINTER( fun%method_Data, fPtr )
    scheme => fPtr%solverData%scheme
    allocate(srcElemPos(scheme%layout%fStencil%QQ))
    allocate(weights(scheme%layout%fStencil%QQ))
    allocate(srcRes(scheme%layout%fStencil%QQ*fun%nComponents))
    allocate(pntVal(fun%nComponents))

    res = 0.0_rk
    do iPnt = 1, nPnts
      weights = 0.0_rk
      srcRes = 0.0_rk
      srcElemPos = 0
      coord =  tem_CoordOfReal(tree, point(iPnt,:), tree%global%maxLevel )
      ! returns position of existing element in tree which contains point
      ! 0 if no corresponding node is found,
      ! or the negative of the found ID, if it is a virtual node.
      elemPos = abs(tem_PosofId( tem_IdOfCoord(coord), tree%treeID ))

      ! skip this point if no corresponding element is found
      if (elemPos == 0) cycle

      ! level of existing element
      iLevel = tem_levelOf( tree%treeID( elemPos ) )

      nSrcElems = 0
      ! use this element as 1st source element for interpolation
      ! if restPosition is part of Stencil
      if (scheme%layout%fStencil%QQ /= scheme%layout%fStencil%QQN) then
        nSrcElems = 1
        srcElemPos(nSrcElems)  = elemPos
      end if

      bary = tem_BaryOfID(tree, tree%treeID(elemPos ))
      dist = abs(bary - point(iPnt,:))
      dist_len = sqrt(dot_product(dist, dist))
      ! if point is exact bary center of current element then
      ! no need to do interpolation
      if ( dist_len .fne. 0.0_rk ) then
        ! position of element in levelDesc total list
        statePos = fPtr%solverData%geometry%levelPointer( elemPos )

        ! get existing neighbor elements in actual tree
        do iNeigh = 1, scheme%layout%fStencil%QQN
          ! neighbor position in levelwise list
          neighLvlPos = scheme%levelDesc( iLevel )%neigh(1) &
            &                 %nghElems( iNeigh, statePos )
          if (neighLvlPos > 0 .and. neighLvlPos <= scheme                      &
            &                       %levelDesc(iLevel)%offset(2, eT_fluid)) then
            ! position of neighbor in actual tree
            neighPos = scheme%levelDesc(iLevel)%pntTID(neighLvlPos)
            ! neighbor exist in same level
            if (neighPos > 0) then
              nSrcElems = nSrcElems + 1
              srcElemPos(nSrcElems)  = neighPos
            else
              ! if neighbor not in same level then use tem_PosOfId
              neighPos = abs(tem_PosOfId( scheme%levelDesc(iLevel) &
                &                               %total(neighLvlPos),    &
                &            tree%treeID ))
              if (neighPos > 0) then
                nSrcElems = nSrcElems + 1
                srcElemPos(nSrcElems)  = neighPos
              end if
            end if !neighbor in same level
          end if
        end do !iNeigh

        ! calculate weight
        do iSrc = 1, nSrcElems
          bary = tem_BaryOfID(tree, tree%treeID( srcElemPos(iSrc) ))
          dist = abs(bary - point(iPnt,:))
          weights(iSrc) = 1.0_rk/sqrt(dot_product( dist, dist ))
        end do
      else
        ! if point is exact bary center of current element then
        ! no need to do interpolation
        weights(nSrcElems) = 1.0_rk
      end if

      sumWeight = sum(weights)
      weights = weights/sumWeight

      ! get source element values
      call varSys%method%val(fun%myPos)%get_element( &
         & varSys  = varSys,                         &
         & elemPos = srcElemPos(1:nSrcElems),        &
         & time    = time,                           &
         & tree    = tree,                           &
         & nElems  = nSrcElems,                      &
         & nDofs   = 1,                              &
         & res     = srcRes                          )

      ! Linear interpolation res = sum(weight_i*phi_i)
      pntVal = 0.0_rk
      do iSrc = 1, nSrcElems
        pntVal(:) = pntVal(:) + weights(iSrc) &
          & * srcRes( (iSrc-1)*fun%nComponents+1 : iSrc*fun%nComponents )
      end do
!      write(*,*) 'Linear intp pntVal ', pntVal
!      pntVal = mus_interpolate_quadratic2d_leastSq( &
!        &          srcMom       = srcRes,           &
!        &          targetCoord  = point(iPnt,1:2),  &
!        &          nSourceElems = nSrcElems         )
!      write(*,*) 'quadratic intp pntVal ', pntVal
      res( (iPnt-1)*fun%nComponents+1 : iPnt*fun%nComponents ) = pntVal

!write(dbgUnit(1),*) 'treeID ', tree%treeID(elemPos), 'val ', pntVal

    end do !iPnt
  end subroutine mus_deriveVar_forPoint
  ! ************************************************************************* !


  ! ************************************************************************* !
  !> Routine to get the actual value for a given array of indices for
  !! musubi derive variables
  !! The indices belong to the grwarray of points storing levelwise in
  !! Pointdata%pntLvl(iLevel).
  !! Hence this routines takes the indeices as input, can refer to the pointData
  !! and evaluate the variable and returns the values
  subroutine mus_generic_varFromPDF_fromIndex( fun, varSys, time, iLevel, idx, &
    &                                       idxLen, nVals, fnCalcPtr, res )
    ! -------------------------------------------------------------------------- !
    !> Description of the method to obtain the variables,
    class(tem_varSys_op_type), intent(in)   :: fun

    !> The variable system to obtain the variable from.
    type(tem_varSys_type), intent(in)       :: varSys

    !> Point in time at which to evaluate the variable.
    type(tem_time_type), intent(in)         :: time

    !> Level on which values are requested
    integer, intent(in)                     :: iLevel

    !> Index of points in the growing array and variable val array to
    !! return.
    !! Size: n
    integer, intent(in)                     :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: dependes on number of first index for contiguous array,
    !! but the sum of all idxLen is equal to nVals
    integer, optional, intent(in)           :: idxLen(:)

    !> Number of values to obtain for this variable (vectorized access).
    integer, intent(in)                     :: nVals

    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer  :: fnCalcPtr

    !> Resulting values for the requested variable.
    !!
    !! Dimension: n requested entries x nComponents of this variable
    !! Access: (iElem-1)*fun%nComponents + iComp
    real(kind=rk), intent(out)              :: res(:)
    ! -------------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer     :: fPtr, state_fPtr
    type(mus_scheme_type), pointer          :: scheme
    real(kind=rk), allocatable              :: pdf(:)
    integer, allocatable                    :: posInState(:)
    integer                                 :: varPos, pdfElemSize
    integer                                 :: elemPos, iVal
    ! -------------------------------------------------------------------------- !
    write(dbgUnit(4),*) 'get the derived values of indices for variable ',  &
      &                  trim(varSys%varname%val(fun%myPos))
    ! distinguish if we have an array of index or we have contingous memory
    ! access where index are always first entries!
    if (present(idxLen)) then
      call tem_abort('Error: idxLen is not supported in get_valOfIndex for ' &
        &          //'state variable')
    end if

    !convert pointer from C to Fotran
    call C_F_POINTER( fun%method_Data, fPtr )
    scheme => fPtr%solverData%scheme
    ! pdf entries per Element
    pdfElemSize = scheme%layout%fStencil%QQ*scheme%nFields

    allocate( pdf( pdfElemSize*nVals ) )
    allocate( posInState(nVals) )

    varPos = fun%input_varPos(1)

    ! get pdf values for IDX
    call varSys%method%val( varPos )%get_ValOfIndex(  &
      &     varSys = varSys,                          &
      &     time   = time,                            &
      &     iLevel = iLevel,                          &
      &     idx    = fPtr%opData%input_pntIndex(1)    &
      &              %indexLvl(iLevel)%val( idx(:) ), &
      &     nVals  = nVals,                           &
      &     res    = pdf                              )

    ! the pointData is stored only in state variable function pointer
    call C_F_POINTER( varSys%method%val(varPos)%method_Data, state_fPtr )
    do iVal = 1, nVals
      ! Position in original tree
      elemPos = state_fPtr%pointData%pntLvl(iLevel)%elemPos%val(idx(iVal))
      ! position in levelwise state list
      posInState(iVal) = fPtr%solverData%geometry%levelPointer(elemPos)
    end do

    ! Call the procedure that does the calculation
    call fnCalcPtr(fun        = fun,                    &
      &            varSys     = varSys,                 &
      &            iLevel     = iLevel,                 &
      &            posInState = posInState,             &
      &            stencil    = scheme%layout%fstencil, &
      &            pdf        = pdf,                    &
      &            nVals      = nVals,                  &
      &            res        = res                     )

    deallocate( pdf )
    deallocate( posInState )

  end subroutine mus_generic_varFromPDF_fromIndex
  ! ************************************************************************* !

  ! ************************************************************************* !
  !> This routine prepares the data for variable derivation or operators. It
  !! gathers all input variables from the variable system, calls the function
  !! with the actual calculation.
  recursive subroutine mus_generic_fromPDF_forElement(fun, varSys, elempos,    &
    &  tree, time, nVals, fnCalcPtr, nDofs, res )
    ! --------------------------------------------------------------------------
    ! -------------------------------------------------------------------------- !
    !> description of the method to obtain the variables, here some preset
    !! values might be stored, like the space time function to use or the
    !! required variables.
    class(tem_varsys_op_type), intent(in)   :: fun

    !> the variable system to obtain the variable from.
    type(tem_varsys_type), intent(in)       :: varsys

    !> position of the treeid of the element to get the variable for in the
    !! global treeid list.
    integer, intent(in)                     :: elempos(:)

    !> point in time at which to evaluate the variable.
    type(tem_time_type), intent(in)         :: time

    !> Number of values to obtain for this variable
    integer, intent(in)                     :: nVals

    !> global treelm mesh info
    type(treelmesh_type), intent(in)        :: tree

    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer  :: fnCalcPtr

    !> number of degrees of freedom within an element.
    integer, intent(in)                     :: ndofs

    !> resulting values for the requested variable.
    !!
    !! linearized array dimension:
    !! (n requested entries) x (ncomponents of this variable)
    !! x (ndegrees of freedom)
    !! access: (ielem-1)*fun%ncomponents*ndofs +
    !!         (idof-1)*fun%ncomponents + icomp
    real(kind=rk), intent(out)              :: res(:)
    ! --------------------------------------------------------------------------
    !> contain all pdf's for element
    real(kind=rk), allocatable              :: pdf(:)
    real(kind=rk), allocatable              :: loc_res(:)
    type(mus_varSys_data_type), pointer     :: fPtr
    type(mus_scheme_type), pointer          :: scheme
    integer                                 :: iElem, pdfpos
    integer                                 :: iLevel, posInState
    ! --------------------------------------------------------------------------
    call C_F_POINTER( fun%method_Data, fPtr )
    scheme => fPtr%solverData%scheme
    ! position of PDF in glob System
    pdfpos = fun%input_varPos(1)
    allocate( pdf(scheme%layout%fStencil%QQ*scheme%nFields) )
    allocate( loc_res( fun%nComponents ) )
    ! loop over all elements ( could also be done in chunks or in total )
    do iElem = 1, nVals
     ! get source element value
     call varSys%method%val(pdfpos)%get_element( &
       & varSys  = varSys,                       &
       & elemPos = (/ elempos(iElem) /),         &
       & time    = time,                         &
       & tree    = tree,                         &
       & nElems  = 1,                            &
       & nDofs   = nDofs,                        &
       & res     = pdf                           )

     ! get iLevel for element
     iLevel = tem_levelOf( tree%treeID( elemPos(iElem ) ) )
     ! get posInState
     posInState = fPtr%solverData%geometry%levelPointer( elemPos(iElem) )

     ! Call the procedure that does the calculation
     call fnCalcPtr(fun        = fun,                    &
       &            varSys     = varSys,                 &
       &            iLevel     = iLevel,                 &
       &            posInState = (/ posInState /),       &
       &            stencil    = scheme%layout%fstencil, &
       &            pdf        = pdf,                    &
       &            nVals      = 1,                      &
       &            res        = loc_res                 )
      res((iElem-1)*fun%nComponents + 1: iElem*fun%nComponents ) = loc_res
    end do !iElem
    deallocate(pdf)
    deallocate(loc_res)
  end subroutine mus_generic_fromPDF_forElement
! **************************************************************************** !

end module mus_varSys_module
! **************************************************************************** !

