derive_density Subroutine

  subroutine derive_density(fun, varsys, elempos, time, tree, n, 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

    !> TreeID of the element to get the variable for.
    integer, intent(in) :: elempos(:)

    !> 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) :: n

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

    !> Resulting values for the requested variable.
    !!
    !! The first dimension are the n requested entries and the second
    !! dimension are the components of this variable.
    !! Third dimension are the degrees of freedom.
    real(kind=rk), intent(out) :: res(:)
    ! --------------------------------------------------------------------------!
    integer :: iElem, iComp, iDof, pos, input_varPos, inVar_nComp
    real(kind=rk) :: dens
    type(solver_type), pointer :: fPtr
    ! --------------------------------------------------------------------------!
    write(*,*) 'myPos ', fun%myPos
    write(*,*) 'derive variable label: ', trim(varSys%varname%val(fun%myPos))
    call C_F_POINTER( fun%method_Data, fPtr )

    input_varPos = fun%input_varPos(1)
    inVar_nComp = varSys%method%val(input_varPos)%nComponents
    write(*,*) 'input_varPos ', input_varPos
    res = 0.0_rk

    do iElem = 1, n
      ! if state array is defined level wise then use levelPointer(pos)
      ! to access state array
      pos = elemPos(iElem)
      ! use state_varPos(iComp) if state has more than one variable
      do iDof = 1, nDofs
        dens = 0.0_rk
        do iComp = 1, inVar_nComp
          dens = dens + fPtr%state( (pos-1)*varSys%nScalars*nDofs            &
            &             + (iDof-1)*varSys%nScalars                         &
            &             + varSys%method%val(input_varPos)%state_varPos(iComp) )
        end do
        res( (iElem-1)*nDofs + iDof )  = dens
      end do
    end do

  end subroutine derive_density