! Copyright (c) 2012-2014 Jens Zudrop <j.zudrop@grs-sim.de>
! Copyright (c) 2012, 2016-2017 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2012 Laura Didinger <l.didinger@grs-sim.de>
! Copyright (c) 2012-2013 Melven Zoellner <yameta@freenet.de>
! Copyright (c) 2012-2016, 2018 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2012 Vyacheslav Korchagin <v.korchagin@grs-sim.de>
! Copyright (c) 2012 Jan Hueckelheim <j.hueckelheim@grs-sim.de>
! Copyright (c) 2012 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2013-2014, 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2013-2016 Nikhil Anand <nikhil.anand@uni-siegen.de>
! Copyright (c) 2013-2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2016-2017 Tobias Girresser <tobias.girresser@student.uni-siegen.de>
! Copyright (c) 2018 Robin Weihe <robin.weihe@student.uni-siegen.de>
! Copyright (c) 2018 Neda Ebrahimi Pour <neda.epour@uni-siegen.de>
!
! Permission to use, copy, modify, and distribute this software for any
! purpose with or without fee is hereby granted, provided that the above
! copyright notice and this permission notice appear in all copies.
!
! THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHORS DISCLAIM ALL WARRANTIES
! WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
! MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
! ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
! WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
! ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
! OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
! **************************************************************************** !

!> author: Harald Klimach
!! author: Peter Vitt
!! author: Melven Zoellner
!! author: Jens Zudrop
!!
!! Provides routines to initialize and update source terms (accumulate the RHS
!! of the PDE)
!!
module atl_source_module
  use, intrinsic :: iso_c_binding,   only: c_loc,       &
    &                                      c_f_pointer, &
    &                                      c_ptr
  use env_module,                    only: rk, labellen
  use aotus_module,                  only: flu_state

  ! include treelm modules
  use tem_varSys_module,             only: tem_varSys_type,              &
    &                                      tem_varSys_append_derVar,     &
    &                                      tem_varSys_proc_point,        &
    &                                      tem_varSys_proc_element,      &
    &                                      tem_varSys_proc_setParams,    &
    &                                      tem_varSys_proc_getParams,    &
    &                                      tem_varSys_proc_setupIndices, &
    &                                      tem_varSys_proc_getValOfIndex
  use tem_varMap_module,             only: tem_possible_variable_type, &
    &                                      tem_variable_loadMapping
  use tem_stringKeyValuePair_module, only: grw_stringKeyValuePairArray_type, &
    &                                      append
  use tem_dyn_array_module,          only: PositionOfVal
  use treelmesh_module,              only: treelmesh_type
  use tem_element_module,            only: eT_fluid
  use tem_time_module,               only: tem_time_type
  use tem_aux_module,                only: tem_abort
  use tem_logging_module,            only: logUnit
  use tem_spacetime_fun_module,      only: tem_st_fun_listElem_type
  use tem_grow_array_module,         only: append
  use tem_varSys_module,             only: tem_varSys_solverData_evalElem_type
  use tem_stringKeyValuePair_module, only: tem_stringKeyValuePair_type, &
    &                                      init

  use ply_poly_project_module,       only: ply_poly_project_type, &
    &                                      ply_prj_body_type,     &
    &                                      assignment(=)

  ! include ateles modules
  use atl_equation_module,           only: atl_equations_type
  use atl_scheme_module,             only: atl_modg_scheme_prp,  &
    &                                      atl_scheme_type,      &
    &                                      atl_modg_2d_scheme_prp
  use atl_cube_elem_module,          only: atl_cube_elem_type
  use atl_materialPrp_module,        only: atl_material_type, &
    &                                      atl_ConstMatIdx
  use atl_varSys_module,             only: atl_varSys_solverData_type, &
    &                                      atl_set_stFun_getElement,   &
    &                                      atl_get_new_varSys_data_ptr
  use atl_operator_module,           only: atl_set_opVar_getElement
  use atl_source_types_module,       only: atl_source_type,        &
    &                                      atl_eqn_sourceMap_type, &
    &                                      atl_init_source_type
  use atl_reference_element_module,  only: atl_refToPhysCoord

  implicit none
  private


  public :: atl_update_sourcedata
  public :: atl_allocate_sourceData
  public :: atl_deallocate_sourceData
  public :: atl_initialize_sources
  public :: atl_append_newSourceVars
  public :: atl_fill_sourceIndex
  public :: atl_source_prim2cons


contains


! ******************************************************************************
  subroutine atl_initialize_sources( source, initSource, conf, equation, &
    &                                poly_proj_list, mesh_list, tree,    &
    &                                varSys_data                         )
    ! --------------------------------------------------------------------------
    !> Instance of atl_source_type to be initialized.
    !! This instance will contain the source definitions from lua as well as
    !! the corresponding variables in the global variable system.
    type(atl_source_type), intent(inout) :: source

    !> Initialize source type contains possible source terms and function
    !! pointers to update those source terms
    type(atl_init_source_type), intent(in) :: initSource

    !> lua state
    type(flu_state), intent(inout) :: conf

    !> Description on the equation system to solve.
    type(atl_Equations_type), intent(inout) :: equation

    !> unique list for projection methods
    type(ply_poly_project_type), intent(in) :: poly_proj_list(:)

    !> Mesh data in treelmesh format.
    type(treelmesh_type), intent(in) :: tree

    !> Mesh list to access the level descriptors
    type(atl_cube_elem_type), intent(in) :: mesh_list(tree%global%minLevel:)

    type(atl_varSys_solverData_type), intent(in), target :: varSys_data
    ! --------------------------------------------------------------------------
    type(tem_stringKeyValuePair_type) :: predef_src
    ! --------------------------------------------------------------------------

    ! Initilized keyValuePair dict for sources
    call init( me = source%varDict )

    call tem_variable_loadMapping(                     &
      & possVars            = initSource%poss_srcVars, &
      & conf                = conf,                    &
      & key                 = "source",                &
      & varSys              = equation%varSys,         &
      & varDict             = source%varDict           )

    ! In case the equation has a pre-defined/permament source.
    ! The key and the value is appended to the sourceDict here.
    ! Then, this source term would get added if it was added as a
    ! poss_srcVars
    ! @todo NA: This needs to fixed and done elegantly
    select case(equation%eq_kind)
      case('filtered_navier_stokes')
        predef_src%key = 'ransSource'
        predef_src%value = ''
        call append(me = source%varDict, val = predef_src)
      case('filtered_navier_stokes_2d')
        predef_src%key = 'rans2dsource'
        predef_src%value = ''
        call append(me = source%varDict, val = predef_src)
    end select

    ! As we need some information about mesh and projection during source term
    ! evaluation, we have to reference them here.
    call atl_append_newSourceVars( me           = source,                  &
      &                            varSys       = equation%varSys,         &
      &                            varSys_data  = varSys_data,             &
      &                            poss_srcVars = initSource%poss_srcVars, &
      &                            eval_source  = initSource%eval_source   )
    ! convert the sourceVars from prim to cons
    ! At the moment this is only for SpongeLayer which is checked inside this fn
    ! KM: Do conversion from prim 2 cons before setupIndices because for
    ! spatial function values are stored in st_fun setupIndices
    if (equation%hasPrimitiveVariables) then
      call atl_source_prim2cons( equation = equation,        &
        &                        source   = source,          &
        &                        varSys   = equation%varSys  )
    end if
    ! Fill source idx returned from setupIndices of variable to retreive
    ! a value later on in update_sourceTerm using get_valOfIndex.
    ! Source index are not filled for predefined or permanent source terms.
    call atl_fill_sourceIndex( source         = source,                 &
      &                        varSys         = equation%varSys,        &
      &                        nDim           = equation%nDimensions,   &
      &                        tree           = tree,                   &
      &                        poly_proj_list = poly_proj_list,         &
      &                        mesh_list      = mesh_list               )

  end subroutine atl_initialize_sources
! ******************************************************************************


! *******************************************************************************
  subroutine atl_append_newSourceVars( me, varSys, varSys_data, poss_srcVars,  &
    &                                  eval_source )
    ! ---------------------------------------------------------------------------
    !> Instance of atl_source_type to be initialized.
    !! This instance will contain the source definitions from lua as well as
    !! the corresponding variables in the global variable system.
    type(atl_source_type), intent(inout) :: me
    !> The variable system to which the souce variables have to be added
    type(tem_varSys_type), intent(inout) :: varSys
    !> Data for the variable System
    type(atl_varSys_solverData_type), intent(in), target :: varSys_data
    !> The list of possible source variables. This is used to determine the
    !! index of the eval-source_routine in eval_source.
    type(tem_possible_variable_type) :: poss_srcVars
    !> List of function pointers for each possible source to apply the source
    !! term to the state.
    type(atl_eqn_sourceMap_type), allocatable :: eval_source(:)
    ! ---------------------------------------------------------------------------
    integer :: addedSourcePosition, iSource, nComponents, iState, poss_srcVarPos
    integer :: nComp_data
    logical :: wasAdded
    character(labellen), allocatable :: inputVarNames(:)
    character(labellen) :: varname
    procedure(tem_varSys_proc_point), pointer :: get_point => null()
    procedure(tem_varSys_proc_element), pointer :: get_element => null()
    procedure(tem_varSys_proc_setParams), pointer :: set_params => Null()
    procedure(tem_varSys_proc_getParams), pointer :: get_params => Null()
    procedure(tem_varSys_proc_setupIndices), pointer ::setup_indices => NULL()
    procedure(tem_varSys_proc_getValOfIndex), pointer :: get_ValOfIndex => NULL()
    ! ---------------------------------------------------------------------------

    ! List all state variables to use them as input variables for the routines
    ! that apply the source variables to the state.
    ! Additionally, the space-time-function for this source variable is also
    ! needed as an input variable. Thus we allocatate one additional slot for
    ! the name of the space-time-function.
    allocate(inputVarNames(varSys%nStateVars+1))
    do iState = 1, varSys%nStateVars
      write(logUnit(5),*) 'Adding state variable ' &
        & // trim(varSys%varname%val(iState))            &
        & // ' to the list of input variables'
      inputVarNames(iState) = varSys%varname%val(iState)
    end do

    ! Now we have to add the source variables. These are the 'real' source
    ! variables that can make use of the space-time-functions added earlier.
    ! Therefore we first have to allocate the slots in atl_source_type to store
    ! information like affection elements or functions to evaluate the source
    ! term.
    write(logUnit(9),'(A,I1,a)') 'Allocating space for ', me%varDict%nVals, &
      & ' sources'
    allocate( me%method( me%varDict%nVals ) )
    do iSource = 1, me%varDict%nVals

      ! get the name of the current variable
      varname = trim(me%varDict%val(iSource)%key)

      ! The function to apply the source term to the state has already been
      ! added to the eval_source list by the equation systems during their
      ! initialization of the source terms. So we need to find the position of
      ! the variable in the list of possible variables, as this is also the
      ! position of the function pointer in the eval_source list.
      poss_srcVarPos = PositionOfVal( me = poss_srcVars%varName, &
        &                             val = trim(varname)        )
      me%method(iSource)%updateSrc => eval_source(poss_srcVarPos)%do
      nComponents = poss_srcVars%nComponents%val(poss_srcVarPos)

      if (len_trim(me%varDict%val(iSource)%value) == 0 ) then
        ! In case of the predefined sources the input_varname should just contain
        ! the state variables
        ! We use the global atl_varSys_data here. This is possible because the
        ! first access to atl_varSys_data will happen after it is initialized.
        call tem_varSys_append_derVar(                                 &
          & me             = varSys,                                   &
          & varName        = varname,                                  &
          & operType       = 'predefined_source',                      &
          & nComponents    = nComponents,                              &
          & input_varname  = inputVarNames(1:varSys%nStateVars),       &
          & method_data    = atl_get_new_varSys_data_ptr(varSys_data), &
          & get_point      = get_point,                                &
          & get_element    = get_element,                              &
          & set_params     = set_params,                               &
          & get_params     = get_params,                               &
          & setup_indices  = setup_indices,                            &
          & get_valOfIndex = get_valOfIndex,                           &
          & pos            = addedSourcePosition,                      &
          & wasAdded       = wasAdded                                  )

      else
        ! Here we add the source variable specific space-time-function
        ! as the last input variable.
        ! The state variables in the first slots stay unchanged.
        inputVarNames(varSys%nStateVars+1) = me%varDict%val(iSource)%value

        ! We use the global atl_varSys_data here. This is possible because the
        ! first access to atl_varSys_data will happen after it is initialized.
        call tem_varSys_append_derVar(                                 &
          & me             = varSys,                                   &
          & varName        = varname,                                  &
          & operType       = 'source',                                 &
          & nComponents    = nComponents,                              &
          & input_varname  = inputVarNames,                            &
          & method_data    = atl_get_new_varSys_data_ptr(varSys_data), &
          & get_point      = get_point,                                &
          & get_element    = get_element,                              &
          & set_params     = set_params,                               &
          & get_params     = get_params,                               &
          & setup_indices  = setup_indices,                            &
          & get_valOfIndex = get_valOfIndex,                           &
          & pos            = addedSourcePosition,                      &
          & wasAdded       = wasAdded                                  )

      endif

      if (wasAdded) then

        write(logUnit(5),*) 'Appended variable ' // trim(varname) &
          & // ' to the variable system'
        me%method(iSource)%srcTerm_varPos = addedSourcePosition
        ! Take the name of the variable that provides access to the space time
        ! function for this source term, find it's position in the
        ! varSys%varname array as this is also the position in the
        ! varSys%method array and assign it to the source. Thus we can access
        ! the space time variable immediately when evaluating the source.
        ! it makes sense to do this only for the sources which are not pre-defined
        if (len_trim(me%varDict%val(iSource)%value) > 0 ) then

          me%method(iSource)%data_varPos = positionOfVal( &
            & me  = varSys%varname,                       &
            & val = me%varDict%val(iSource)%value )

          nComp_data = varSys%method%val(me%method(iSource)%data_varPos) &
            &                %nComponents
          if ( nComp_data /= nComponents ) then
            write(logUnit(1),'(a)') 'Error: Appending source variable'
            write(logUnit(1),'(a,i0)') 'nComponent of defined variable: "' &
              & // trim(me%varDict%val(iSource)%value)//'"= ', nComp_data
            write(logUnit(1),'(a,i0)') '/= nComponent of expected variable: "' &
              & // trim(varname)//'"= ', nComponents
            call tem_abort()
          end if
        else
          ! permanent source terms
          me%method(iSource)%data_varPos = -1
          me%method(iSource)%isPermanent = .true.
        end if

      else

        write(logUnit(1),*) 'Error: variable ' // trim(varname)          &
          & // ' is not added to the variable system. Adding it to the ' &
          & // ' list for the next turn'

      end if

    end do

  end subroutine atl_append_newSourceVars
! *******************************************************************************


! *******************************************************************************
  !> Create source elements list for given source variable
  !!
  !! This routine was taken from musubi's mus_source_module. This was
  !! necessary, because we can't move the routine to treelm and work with
  !! type extensions, as there are differences between atl_source_elems_type and
  !! mus_source_elems_type, which prevents us from using a shared and extendable
  !! tem_source_op_type.
  !!
  !! If this routine is changed, please also keep the coresponding routine in
  !! mus_source_module up to date.
  subroutine atl_fill_sourceIndex( source, varSys, nDim, tree, poly_proj_list, &
    &                              mesh_list )
    ! ---------------------------------------------------------------------------
    !> Instance of atl_source_type to be initialized.
    !! This instance will contain the source definitions from lua as well as
    !! the corresponding variables in the global variable system.
    type(atl_source_type), intent(inout) :: source

    !> global variable system
    type(tem_varSys_type), intent(in) :: varSys

    !> equation nDimensions
    integer, intent(in) :: nDim

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

    !> unique list for projection methods
    type(ply_poly_project_type), target, intent(in) :: poly_proj_list(:)

    !> Mesh list to access the level descriptors
    type(atl_cube_elem_type), intent(in) :: mesh_list(tree%global%minLevel:)
    ! ---------------------------------------------------------------------------
    integer :: iElem, iLevel, iSource, nElems, nQuadPnts
    integer :: data_varPos, minLevel, maxLevel
    type(ply_prj_body_type), pointer :: polyProjBody
    real(kind=rk), allocatable, dimension(:,:) :: refPoints, srcPhysPnts
    integer, allocatable :: idx(:)
    ! ---------------------------------------------------------------------------
    write(logUnit(1),*) 'Fill source index'
    ! return if there are not active source terms
    if ( source%varDict%nVals == 0 ) then
      write(logUnit(1),*) 'No active source terms'
      return
    end if

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

    do iSource = 1, source%varDict%nVals
      data_varPos = source%method(iSource)%data_varPos
      if ( data_varPos > 0 ) then
        ! Set params for source variables
        call varSys%method%val(data_varPos)%set_params( &
          & varSys   = varSys,                          &
          & instring = 'isSurface = false'              )
      end if

      ! allocate source elems level array
      allocate(source%method(iSource)%elems(minLevel:maxLevel))
    end do


    write(logUnit(3),*) ' Setup indices for source terms'
    ! run over all levels
    do iLevel = minLevel, maxLevel

      ! First get the amount of points to allocate arrays accordingly
      select case(nDim)
      case (1)
        ! set the correct projection type
        polyProjBody => poly_proj_list(source%poly_proj_pos(iLevel))%body_1d
      case (2)
        polyProjBody => poly_proj_list(source%poly_proj_pos(iLevel))%body_2d
      case (3)
        polyProjBody => poly_proj_list(source%poly_proj_pos(iLevel))%body_3d
      end select

      nQuadPnts = polyProjBody%nQuadPoints
      allocate(refPoints(nQuadPnts,3))
      refPoints = polyProjBody%nodes

      allocate(srcPhysPnts(nQuadPnts,3))
      allocate(idx(nQuadPnts))

      ! total number of fluid elements on this level
      nElems = mesh_list(iLevel)%descriptor%elem%nElems(eT_fluid)

      do iElem = 1, nElems

        !  Get the cheb. phys. coordinate for the element
        call atl_refToPhysCoord(                                   &
          & refPoints  = refPoints,                                &
          & nPoints    = nQuadPnts,                                &
          & baryCoord  = mesh_list(iLevel)%bary_coord( iElem, :),  &
          & elemLength = mesh_list(iLevel)%length,                 &
          & physPoints = srcPhysPnts                               )

        do iSource = 1, source%varDict%nVals
          data_varPos = source%method(iSource)%data_varPos
          ! we do it only for the list of sources defined in the lua files
          if ( data_varPos > 0 ) then
            idx = 0
            call varSys%method%val(data_varPos)%setup_indices( &
              &  varSys     = varSys,                          &
              &  point      = srcPhysPnts,                     &
              &  iLevel     = iLevel,                          &
              &  tree       = tree,                            &
              &  nPnts      = nQuadPnts,                       &
              &  idx        = idx                              )

            ! store index only if idx>0 i.e point is active on data variable
            ! shape
            if ( all(idx > 0) ) then
              call append( me  = source%method(iSource)%elems(iLevel) &
                &                                      %posInTotal,   &
                &          val = iElem                                )
              call append( me  = source%method(iSource)%elems(iLevel)%idx, &
                &          val = idx                                       )
            end if
          end if
        end do ! iSource
      end do ! iElem

      ! set number of elements active per source term.
      ! For predefined or permanent sources, set nFluids on this level
      do iSource = 1, source%varDict%nVals
        data_varPos = source%method(iSource)%data_varPos
        if ( data_varPos > 0 ) then
          source%method(iSource)%elems(iLevel)%nElems               &
            & = source%method(iSource)%elems(iLevel)%posInTotal%nVals
        else
          source%method(iSource)%elems(iLevel)%nElems = nElems
        end if
      end do

      ! Deallocate the allocated arrays before the end of ilevel loop
      deallocate(refPoints)
      deallocate(srcPhysPnts)
      deallocate(idx)

    end do ! iLevel

    do iSource = 1, source%varDict%nVals
      write(logUnit(10),*) ' source: (iLevel , nElem) for srcTerm: ' &
        &                // trim(source%varDict%val(iSource)%key)
      do iLevel = minLevel, maxLevel
        write(logUnit(10),*) '(',iLevel,' , ', source%method(iSource)   &
          &                                    %elems(iLevel)%nElems, ')'
      end do
    end do
  end subroutine atl_fill_sourceIndex
! *******************************************************************************

! ******************************************************************************
  !> This routine converts primitive variables in source terms to convervative
  subroutine atl_source_prim2cons(equation, source, varSys)
    ! --------------------------------------------------------------------------
    type(atl_equations_type), intent(in) :: equation
    type(atl_source_type), intent(in) :: source
    type(tem_varSys_type), intent(in) :: varSys
    ! --------------------------------------------------------------------------
    character(len=labelLen) :: vartype
    integer :: iSrc, iFun, data_varPos
    real(kind=rk) :: cons(1,1,varSys%nScalars)
    type(tem_st_fun_listElem_type), pointer :: stFun => NULL()
    logical :: isSpatialSponge
    ! --------------------------------------------------------------------------
    do iSrc = 1, source%varDict%nVals

      data_varPos = source%method(iSrc)%data_varPos

      select case ( trim(source%varDict%val(iSrc)%key) )
      case ('spongelayer')
        call varSys%method%val(data_varPos)%get_params( &
          & varSys    = varSys,                         &
          & instring  = 'vartype',                      &
          & outstring = vartype                         )
        if (trim(vartype) == 'st_fun') then
          call c_f_pointer(                               &
            & varSys%method%val(data_varPos)%method_data, &
            & stFun                                       )

          isSpatialSponge = .true.
          do iFun = 1,size(stFun%val)
            isSpatialSponge = isSpatialSponge                                 &
              & .and. ( stFun%val(iFun)%spatial%kind(1:11) == 'spongelayer' )
          end do
          if (.not. isSpatialSponge) then
            call tem_abort('Error: spongelayer source variable st_fun is ' &
              &          //'not spatial spongelayer')
          end if

          do iFun = 1,size(stFun%val)
            ! write the primitive values in target state to cons
            cons(1,1,:) = stFun%val(iFun)%spatial%sponge%targetstate

            call equation%prim2cons( instate  = cons,     &
              &                      nElems   = 1         )

            ! write it back to targetState
            stFun%val(iFun)%spatial%sponge%targetstate = cons(1,1,:)
          enddo
        else
          call tem_abort('Error: spongelayer source variable is not st_fun')
        end if
      end select
    enddo
  end subroutine atl_source_prim2cons
! ******************************************************************************


! ******************************************************************************
  !> summary: subroutine to calculate the RHS of the PDE from the sum of all
  !! source terms
  subroutine atl_update_sourcedata( equation, time, mesh, poly_proj,       &
    &                               currentLevel, state, material, source, &
    &                               scheme                                 )
    ! --------------------------------------------------------------------------
    !> The equation with source term data
    type(atl_equations_type), intent(in) :: equation
    !> current time
    type(tem_time_type), intent(in)             :: time
    !> Current level mesh information
    type(atl_cube_elem_type), intent(in) :: mesh
    !> The projection used for update the source terms for.
    type(ply_poly_project_type), intent(inout)  :: poly_proj
    !> The current Level
    integer, intent(in)                         :: currentLevel
    !> state vector (sources may depend on the state)
    real(kind=rk), intent(in)                   :: state(:,:,:)
    !> The material description.
    type(atl_material_type), intent(inout)      :: material
    !> sources for this level
    type(atl_source_type), intent(inout)        :: source
    !> The scheme you update the source terms for.
    type(atl_scheme_type), intent(in)               :: scheme
    ! --------------------------------------------------------------------------
    integer :: iSource
    ! --------------------------------------------------------------------------

    ! loop over all sources
    do iSource = 1, source%varDict%nVals

      ! init sourcedata with zeros
      source%method(iSource)%val = 0.0_rk

      select case(scheme%scheme)
      case(atl_modg_scheme_prp, atl_modg_2d_scheme_prp)

        ! build the resulting source term ( call the do function pointer)
        call source%method(iSource)%updateSrc(                        &
          & varSys       = equation%varSys,                           &
          & time         = time,                                      &
          & mesh         = mesh,                                      &
          & poly_proj    = poly_proj,                                 &
          & currentLevel = currentLevel,                              &
          & state        = state,                                     &
          & material     = material%material_dat                      &
          &                        %elemMaterialData(atl_ConstMatIdx) &
          &                        %materialDat(1,1,:),               &
          & sourcedata   = source%method(iSource)%val                 )

        !> @todo PV 20150925 The call above gives the material value for the
        !!                   first point of the first element to evaluate the
        !!                   source term. As this is not necessarily the correct
        !!                   value, we thought about handing over the whole
        !!                   material_dat. Out of luck, this is not necessary to
        !!                   get material up and running, so I didn't do it yet.

      case default
        write(logUnit(1),*) 'ERROR in update_sourcedata: not able to ' &
          & // 'evaluate source terms for this scheme type, stopping...'
        call tem_abort()
      end select

    end do

  end subroutine atl_update_sourcedata
! *******************************************************************************

! *******************************************************************************
  subroutine atl_allocate_sourceData( source, nDofs, nComponents )
    ! ------------------------------------------------------------------- !
    type(atl_source_type), intent(inout) :: source
    integer, intent(in) :: nDofs
    integer, intent(in) :: nComponents
    ! ------------------------------------------------------------------- !
    integer :: iSource
    ! ------------------------------------------------------------------- !

    ! Allocate the array for storing the StateData
    do iSource = 1, size(source%method,1)
      ! In case of multi-leveles it should only be allocated once
      if (.not. allocated (source%method(iSource)%val) ) then
        allocate( source%method(iSource)%val(                       &
          &         maxVal(source%method(iSource)%elems(:)%nElems), &
          &         nDofs,                                          &
          &         nComponents )                                   )
      end if
    end do

  end subroutine atl_allocate_sourceData
! *******************************************************************************

! *******************************************************************************
  !> Deallocates the array for storing the sourceData for the currentLevel
  !  after it is no longer needed
  subroutine atl_deallocate_sourceData( source )
    !> Levelwise list of sources
    type(atl_source_type), intent(inout) :: source
    ! --------------------------------------------------------------------------
    integer :: iSource
    ! --------------------------------------------------------------------------

    do iSource = 1, size(source%method)
      deallocate(source%method(iSource)%val)
    end do

  end subroutine atl_deallocate_sourceData
! *******************************************************************************


end module atl_source_module
