atl_modg_2d_multiLevel_module.f90 Source File


This file depends on

sourcefile~~atl_modg_2d_multilevel_module.f90~~EfferentGraph sourcefile~atl_modg_2d_multilevel_module.f90 atl_modg_2d_multiLevel_module.f90 sourcefile~atl_scheme_module.f90 atl_scheme_module.f90 sourcefile~atl_modg_2d_multilevel_module.f90->sourcefile~atl_scheme_module.f90 sourcefile~atl_kerneldata_module.f90 atl_kerneldata_module.f90 sourcefile~atl_modg_2d_multilevel_module.f90->sourcefile~atl_kerneldata_module.f90 sourcefile~ply_modg_basis_module.f90 ply_modg_basis_module.f90 sourcefile~atl_modg_2d_multilevel_module.f90->sourcefile~ply_modg_basis_module.f90 sourcefile~atl_timer_module.f90 atl_timer_module.f90 sourcefile~atl_modg_2d_multilevel_module.f90->sourcefile~atl_timer_module.f90 sourcefile~atl_modg_2d_scheme_module.f90 atl_modg_2d_scheme_module.f90 sourcefile~atl_modg_2d_multilevel_module.f90->sourcefile~atl_modg_2d_scheme_module.f90 sourcefile~atl_cube_elem_module.f90 atl_cube_elem_module.f90 sourcefile~atl_modg_2d_multilevel_module.f90->sourcefile~atl_cube_elem_module.f90 sourcefile~atl_facedata_module.f90 atl_facedata_module.f90 sourcefile~atl_modg_2d_multilevel_module.f90->sourcefile~atl_facedata_module.f90 sourcefile~atl_scheme_module.f90->sourcefile~ply_modg_basis_module.f90 sourcefile~atl_scheme_module.f90->sourcefile~atl_modg_2d_scheme_module.f90 sourcefile~atl_modg_1d_scheme_module.f90 atl_modg_1d_scheme_module.f90 sourcefile~atl_scheme_module.f90->sourcefile~atl_modg_1d_scheme_module.f90 sourcefile~atl_stabilization_module.f90 atl_stabilization_module.f90 sourcefile~atl_scheme_module.f90->sourcefile~atl_stabilization_module.f90 sourcefile~atl_modg_scheme_module.f90 atl_modg_scheme_module.f90 sourcefile~atl_scheme_module.f90->sourcefile~atl_modg_scheme_module.f90 sourcefile~ply_dof_module.f90 ply_dof_module.f90 sourcefile~atl_kerneldata_module.f90->sourcefile~ply_dof_module.f90 sourcefile~ply_modg_basis_module.f90->sourcefile~ply_dof_module.f90 sourcefile~ply_space_integration_module.f90 ply_space_integration_module.f90 sourcefile~ply_modg_basis_module.f90->sourcefile~ply_space_integration_module.f90 sourcefile~atl_modg_2d_scheme_module.f90->sourcefile~ply_dof_module.f90 sourcefile~atl_boundary_module.f90 atl_boundary_module.f90 sourcefile~atl_facedata_module.f90->sourcefile~atl_boundary_module.f90 sourcefile~atl_modg_1d_scheme_module.f90->sourcefile~ply_dof_module.f90 sourcefile~atl_boundary_module.f90->sourcefile~atl_scheme_module.f90 sourcefile~atl_boundary_module.f90->sourcefile~atl_cube_elem_module.f90 sourcefile~atl_reference_element_module.f90 atl_reference_element_module.f90 sourcefile~atl_boundary_module.f90->sourcefile~atl_reference_element_module.f90 sourcefile~atl_bc_header_module.f90 atl_bc_header_module.f90 sourcefile~atl_boundary_module.f90->sourcefile~atl_bc_header_module.f90 sourcefile~ply_poly_project_module.f90 ply_poly_project_module.f90 sourcefile~atl_boundary_module.f90->sourcefile~ply_poly_project_module.f90 sourcefile~atl_equation_module.f90 atl_equation_module.f90 sourcefile~atl_boundary_module.f90->sourcefile~atl_equation_module.f90 sourcefile~atl_spectral_viscosity_module.f90 atl_spectral_viscosity_module.f90 sourcefile~atl_stabilization_module.f90->sourcefile~atl_spectral_viscosity_module.f90 sourcefile~atl_cons_positivity_preserv_module.f90 atl_cons_positivity_preserv_module.f90 sourcefile~atl_stabilization_module.f90->sourcefile~atl_cons_positivity_preserv_module.f90 sourcefile~atl_covolume_module.f90 atl_covolume_module.f90 sourcefile~atl_stabilization_module.f90->sourcefile~atl_covolume_module.f90 sourcefile~atl_positivity_preserv_module.f90 atl_positivity_preserv_module.f90 sourcefile~atl_stabilization_module.f90->sourcefile~atl_positivity_preserv_module.f90 sourcefile~atl_modg_scheme_module.f90->sourcefile~ply_dof_module.f90

Files dependent on this one

sourcefile~~atl_modg_2d_multilevel_module.f90~~AfferentGraph sourcefile~atl_modg_2d_multilevel_module.f90 atl_modg_2d_multiLevel_module.f90 sourcefile~atl_stabilize_module.f90 atl_stabilize_module.f90 sourcefile~atl_stabilize_module.f90->sourcefile~atl_modg_2d_multilevel_module.f90 sourcefile~atl_compute_module.f90 atl_compute_module.f90 sourcefile~atl_compute_module.f90->sourcefile~atl_modg_2d_multilevel_module.f90 sourcefile~atl_imexrk_module.f90 atl_imexrk_module.f90 sourcefile~atl_imexrk_module.f90->sourcefile~atl_stabilize_module.f90 sourcefile~atl_imexrk_module.f90->sourcefile~atl_compute_module.f90 sourcefile~atl_ssprk2_module.f90 atl_ssprk2_module.f90 sourcefile~atl_ssprk2_module.f90->sourcefile~atl_stabilize_module.f90 sourcefile~atl_ssprk2_module.f90->sourcefile~atl_compute_module.f90 sourcefile~atl_harvesting.f90 atl_harvesting.f90 sourcefile~atl_harvesting.f90->sourcefile~atl_stabilize_module.f90 sourcefile~atl_fwdeuler_module.f90 atl_fwdEuler_module.f90 sourcefile~atl_fwdeuler_module.f90->sourcefile~atl_stabilize_module.f90 sourcefile~atl_fwdeuler_module.f90->sourcefile~atl_compute_module.f90 sourcefile~atl_predcor_cerk4_module.f90 atl_predcor_cerk4_module.f90 sourcefile~atl_predcor_cerk4_module.f90->sourcefile~atl_compute_module.f90 sourcefile~atl_rk4_module.f90 atl_rk4_module.f90 sourcefile~atl_rk4_module.f90->sourcefile~atl_stabilize_module.f90 sourcefile~atl_rk4_module.f90->sourcefile~atl_compute_module.f90 sourcefile~atl_rktaylor_module.f90 atl_rktaylor_module.f90 sourcefile~atl_rktaylor_module.f90->sourcefile~atl_stabilize_module.f90 sourcefile~atl_rktaylor_module.f90->sourcefile~atl_compute_module.f90 sourcefile~atl_global_time_integration_module.f90 atl_global_time_integration_module.f90 sourcefile~atl_global_time_integration_module.f90->sourcefile~atl_imexrk_module.f90 sourcefile~atl_global_time_integration_module.f90->sourcefile~atl_ssprk2_module.f90 sourcefile~atl_global_time_integration_module.f90->sourcefile~atl_fwdeuler_module.f90 sourcefile~atl_global_time_integration_module.f90->sourcefile~atl_predcor_cerk4_module.f90 sourcefile~atl_global_time_integration_module.f90->sourcefile~atl_rk4_module.f90 sourcefile~atl_global_time_integration_module.f90->sourcefile~atl_rktaylor_module.f90 sourcefile~atl_program_module.f90 atl_program_module.f90 sourcefile~atl_program_module.f90->sourcefile~atl_global_time_integration_module.f90 sourcefile~atl_container_module.f90 atl_container_module.f90 sourcefile~atl_container_module.f90->sourcefile~atl_global_time_integration_module.f90

Contents


Source Code

! Copyright (c) 2012, 2014 Jens Zudrop <j.zudrop@grs-sim.de>
! Copyright (c) 2013 Melven Zoellner <yameta@freenet.de>
! Copyright (c) 2013-2016, 2018 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2013, 2015-2018 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2014 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2015-2016 Nikhil Anand <nikhil.anand@uni-siegen.de>
! Copyright (c) 2016 Parid Ndreka
! Copyright (c) 2016 Tobias Girresser <tobias.girresser@student.uni-siegen.de>
! Copyright (c) 2017 Daniel PetrĂ³ <daniel.petro@student.uni-siegen.de>
! Copyright (c) 2019 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.
! **************************************************************************** !

! Copyright (c) 2014,2016-2017 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2014 Harald Klimach <harald.klimach@uni-siegen.de>
!
! Parts of this file were written by Peter Vitt and Harald Klimach for
! University of Siegen.
!
! 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.
! **************************************************************************** !
!
! Return the position of a given ansatz function combination in the
! linearized list of modal coefficients for Q-Tensor product polynomials.
! You must provide
! * Ansatzfunction index in x direction. Index starts with 1.
! * Ansatzfunction index in y direction. Index starts with 1.
! * Ansatzfunction index in z direction. Index starts with 1.
! * The maximal polynomial degree per spatial direction.
! * The variable to store the position of the modal coefficient in the list of
!   modal coefficients in.


! Return the position of a given ansatz function combination in the
! linearized list of modal coefficients for Q-Tensor product polynomials.
! You must provide
! * Ansatzfunction index in x direction. Index starts with 1.
! * Ansatzfunction index in y direction. Index starts with 1.
! * The maximal polynomial degree per spatial direction.
! * The variable to store the position of the modal coefficient in the list of
!   modal coefficients in.


! Return the position of a given ansatz function combination in the
! linearized list of modal coefficients for Q-Tensor product polynomials.
! You must provide
! * Ansatzfunction index in x direction. Index starts with 1.
! * The variable to store the position of the modal coefficient in the list of
!   modal coefficients in.


! Return the position of a given ansatz function combination in the
! linearized list of modal coefficients for P-Tensor product polynomials.
! You must provide
! * Ansatzfunction index in x direction. Index starts with 1.
! * Ansatzfunction index in y direction. Index starts with 1.
! * Ansatzfunction index in z direction. Index starts with 1.
! * The maximal polynomial degree per spatial direction.
! * The variable to store the position of the modal coefficient in the list of
!   modal coefficients in.


! Return the position of a given ansatz function combination in the
! linearized list of modal coefficients for P-Tensor product polynomials.
! You must provide
! * Ansatzfunction index in x direction. Index starts with 1.
! * Ansatzfunction index in y direction. Index starts with 1.
! * The maximal polynomial degree per spatial direction.
! * The variable to store the position of the modal coefficient in the list of
!   modal coefficients in.


! Return the position of a given ansatz function combination in the
! linearized list of modal coefficients for P-Tensor product polynomials.
! You must provide
! * Ansatzfunction index in x direction. Index starts with 1.
! * The variable to store the position of the modal coefficient in the list of
!   modal coefficients in.


! Return the number of degrees of freedom for Q polynomial space
! Your must provide:
! * The maximal polynomial degree per spatial direction
! * The variable to store the number of degrees of freedom for a Q tensor
!   product polynomial


! Return the number of degrees of freedom for broken polynomial space
! Your must provide:
! * The maximal polynomial degree per spatial direction (for P Tensor product
!   polynomials this assumed to be the same for each spatial direction).
! * The variable to store the number of degrees of freedom for a P tensor
!   product polynomial


! Return the number of degrees of freedom for Q polynomial space
! You must provide:
! * The maximal polynomial degree per spatial direction
! * The variable to store the number of degrees of freedom for a Q tensor
!   product polynomial


! Return the number of degrees of freedom for broken polynomial space
! You must provide:
! * The maximal polynomial degree per spatial direction (for P Tensor product
!   polynomials this assumed to be the same for each spatial direction).
! * A variable to store the number of degrees of freedom for a P tensor product
!   polynomial


! Return the number of degrees of freedom for Q polynomial space
! You must provide:
! * The maximal polynomial degree per spatial direction
! * The variable to store the number of degrees of freedom for a Q tensor
!   product polynomial


! Return the number of degrees of freedom for broken polynomial space
! You must provide:
! * The maximal polynomial degree per spatial direction (for P Tensor product
!   polynomials this assumed to be the same for each spatial direction).
! * The variable to store the number of degrees of freedom for a P tensor
!   product polynomial

! The x, y and z ansatz degrees are turned into the degrees of the next
! ansatz function in the layered P list
! You must provide:
! * Ansatz function index in x direction. First ansatz function has index 1.
! * Ansatz function index in y direction. First ansatz function has index 1.
! * Ansatz function index in z direction. First ansatz function has index 1.

! The x and y ansatz degrees are turned into the degrees of the next
! ansatz function in the layered P list
! You must provide:
! * Ansatz function index in x direction. First ansatz function has index 1.
! * Ansatz function index in y direction. First ansatz function has index 1.

! The x ansatz degree is turned into the degree of the next
! ansatz function in the layered P list
! You must provide:
! * Ansatz function index in x direction. First ansatz function has index 1.

! The x, y and z ansatz degrees are turned into the degrees of the next
! ansatz function in the linearized Q tensor
! You must provide:
! * Ansatz function index in x direction. First ansatz function has index 1.
! * Ansatz function index in y direction. First ansatz function has index 1.
! * Ansatz function index in z direction. First ansatz function has index 1.
! * Maximal polynomial degree

! The x and y ansatz degrees are turned into the degrees of the next
! ansatz function in the linearized Q tensor
! You must provide:
! * Ansatz function index in x direction. First ansatz function has index 1.
! * Ansatz function index in y direction. First ansatz function has index 1.
! * Maximal polynomial degree

! The x ansatz degree is turned into the degree of the next
! ansatz function in the linearized Q tensor
! You must provide:
! * Ansatz function index in x direction. First ansatz function has index 1.
!> author: Jens Zudrop
!! Collections of operations and datatypes related to multilevel simulations for
!! the MODG scheme.
!!
!! Due to a compiler problem in Intel 15, we had to replace array syntax by do
!! loops below, and use OMP DO instead of WORKSHARE constructs.
module atl_modg_2d_multilevel_module
  use mpi
  use env_module,                only: rk

  use tem_topology_module,       only: tem_childNumber
  use tem_faceData_module,       only: tem_invFace_map
  use tem_element_module,        only: eT_fluid,            &
    &                                  eT_ghostFromCoarser, &
    &                                  eT_ghostFromFiner
  use tem_timer_module,         only: tem_startTimer, &
    &                                 tem_stopTimer

  use atl_scheme_module,         only: atl_scheme_type
  use atl_modg_2d_scheme_module, only: atl_modg_2d_scheme_type
  use atl_cube_elem_module,      only: atl_cube_elem_type
  use atl_facedata_module,       only: atl_facedata_type
  use atl_kerneldata_module,     only: atl_statedata_type
  use atl_timer_module,          only: atl_elemTimers
  use ply_modg_basis_module,     only: ply_modg_basis_type

  implicit none
  private

  public :: atl_modg_2d_coarseToFineFace, atl_modg_2d_fineToCoarseFace
  public :: atl_modg_coarseToFineElem_2d, atl_modg_fineToCoarseElem_2d


contains


  !> summary: Interpolate modal face representation from coarse to next finer faces (level
  !! difference between coarser and finer faces has to be 1).
  !!
  !! Interpolates functions defined on faces of the current level to faces
  !! of the next finer level. \n
  !! \n
  !!      faces on fine                            face on current           \n
  !!          level                                     level                \n
  !! ------------------------                 ------------------------       \n
  !! |          |           |                 |                      |       \n
  !! |    3     |     4     |                 |                      |       \n
  !! |          |           |                 |                      |       \n
  !! ------------------------   <<----------  |          5           |       \n
  !! |          |           |                 |                      |       \n
  !! |    1     |     2     |                 |                      |       \n
  !! |          |           |                 |                      |       \n
  !! ------------------------                 ------------------------       \n
  !! \n
  !! This is accomplished with lower complexity (with respect to the polynomial
  !! degree) by a dimension by dimension approach: \n
  !! \n
  !!      faces on fine                            face on current           \n
  !!          level                                     level                \n
  !! ------------------------                 ------------------------       \n
  !! |          |           |                 |                      |       \n
  !! |    3     |     4     |                 |                      |       \n
  !! |          |           |                 |                      |       \n
  !! ------------------------   <<----------  |          5           |       \n
  !! |          |           |                 |                      |       \n
  !! |    1     |     2     |                 |                      |       \n
  !! |          |           |                 |                      |       \n
  !! ------------------------                 ------------------------       \n
  !!                 \                               /                       \n
  !!                  \                             /                        \n
  !!                   \                           /                         \n
  !!                    \                         /                          \n
  !!                     ------------------------                            \n
  !!                     |                      |                            \n
  !!                     |          b           |                            \n
  !!                     |                      |                            \n
  !!                     ------------------------                            \n
  !!                     |                      |                            \n
  !!                     |          a           |                            \n
  !!                     |                      |                            \n
  !!                     ------------------------                            \n
  !!
  subroutine atl_modg_2d_coarseToFineFace( minLevel, maxLevel, currentLevel, &
                                  & mesh, facedata, scheme , nScalars)
    ! --------------------------------------------------------------------------
    !> The minumum level of the mesh.
    integer, intent(in) :: minLevel
    !> The maximum level of the mesh.
    integer, intent(in) :: maxLevel
    !> The current level (i.e. the coarse level).
    integer, intent(in) :: currentLevel
    !> The mesh representation.
    type(atl_cube_elem_type), intent(in) :: mesh(minLevel:maxLevel)
    !> The face representations (finer faces are interpolated from coarser ones).
    type(atl_facedata_type), intent(inout) :: facedata(minLevel:maxLevel)
    !> The schemes on the different levels.
    type(atl_scheme_type), intent(in) :: scheme(minLevel:maxLevel)
    !> The number of scalar variables in your equation system.
    integer, intent(in) :: nScalars
    ! --------------------------------------------------------------------------
    integer :: iDir, iAlign, iFace, invAlign
    integer :: iChild
    integer :: i, iDof, iVar
    integer :: nDofsCoarse, nDofsFine, nDofsInter
    integer :: nFluids
    integer :: elemPos, childPos(4)
    integer :: totpos
    integer :: coff, yshift
    integer :: nFaces
    real(kind=rk), allocatable :: faceDat(:,:)
    real(kind=rk), allocatable :: aFace(:,:)
    real(kind=rk), allocatable :: firstFace(:,:), secondFace(:,:)
    ! --------------------------------------------------------------------------

    ! The number of degrees of freedom on a face for the coarse, fine and
    ! intermediate level
    nDofsCoarse = scheme(currentLevel)%modg_2d%maxPolyDegree+1
    nDofsFine = scheme(currentLevel+1)%modg_2d%maxPolyDegree+1
    nDofsInter = (scheme(currentLevel)%modg_2d%maxPolyDegree+1)*1
    nfluids = mesh(currentLevel)%descriptor%elem%nElems(eT_fluid)

    ! Create the intermediate and resulting arrays
    allocate( faceDat(nDofsCoarse,nScalars), &
            & aFace(nDofsInter,nScalars), &
            & firstFace(nDofsFine,nScalars), secondFace(nDofsFine,nScalars) )

    ! Iterate over all the from finer faces of the current level and project
    ! the polynomials downwards to the next finer level.
    do iDir = 1,2
      do iAlign = 1,2
        nFaces = size( mesh(currentLevel)%faces%faces(iDir)           &
          &                                    %fromFinerFace(iAlign) &
          &                                    %elemPos               )
        do iFace = 1,nFaces

          ! The face we have to interpolate. If the face is refined from its
          ! left element (i.e. the right face of the element element is refined,
          !               -> iAlign == 2), we have to work on the left face of
          ! the right element.
          invAlign = tem_invFace_map(iAlign)

          ! Get the face representation (i.e. face 5)
          elemPos = mesh(currentLevel)%faces%faces(iDir)           &
            &                               %fromFinerFace(iAlign) &
            &                               %elemPosOp(iFace)
          if (elempos <= nFluids) then
            totpos = mesh(currentLevel)%descriptor%pntTID(elempos)
            call tem_startTimer( me          = atl_elemTimers, &
              &                  timerHandle = totPos          )
          end if
          do iChild=1,4
            childPos(iChild) = mesh(currentLevel)%faces%faces(iDir)      &
              &                                  %fromFinerFace(iAlign)  &
              &                                  %childPosOp(iChild,iFace)
          end do

          do i=1,nDofsCoarse*nScalars
            iDof = mod(i-1,nDofsCoarse)+1
            iVar = (i-1)/nDofsCoarse + 1
            faceDat(iDof,iVar) = facedata(currentLevel)%faceRep(iDir) &
              &                  %dat(elemPos,iDof,iVar,invAlign)
          end do

          ! 2 Slices of refinement in y
          do yshift=1,2
            do i=1,nDofsInter*nScalars
              iDof = mod(i-1,nDofsInter)+1
              iVar = (i-1)/nDofsInter + 1
              aFace(iDof,iVar) = 0.0_rk
            end do

            do i=1,nDofsFine*nScalars
              iDof = mod(i-1,nDofsFine)+1
              iVar = (i-1)/nDofsFine + 1
              firstFace(iDof,iVar)  = 0.0_rk
              secondFace(iDof,iVar) = 0.0_rk
            end do

            ! Project face 5 to intermediate face:
            call modg_2d_semiRefineFace(                         &
              & modalRepFace  = faceDat,                         &
              & modg_basis    = scheme(currentLevel)%modg_basis, &
              & schemeCoarse  = scheme(currentLevel)%modg_2d,    &
              & schemeFine    = scheme(currentLevel+1)%modg_2d,  &
              & refineDir     = 2,                               &
              & fineFaceShift = yshift,                          &
              & modalRefined  = aFace                            )

            ! Project intermediate face a to fine faces:

            ! ... project a to 1
            call modg_2d_semiRefineFace(                         &
              & modalRepFace  = aFace,                           &
              & modg_basis    = scheme(currentLevel)%modg_basis, &
              & schemeCoarse  = scheme(currentLevel)%modg_2d,    &
              & schemeFine    = scheme(currentLevel+1)%modg_2d,  &
              & refineDir     = 1,                               &
              & fineFaceShift = 1,                               &
              & modalRefined  = firstFace                        )

            ! ... project a to 2
            call modg_2d_semiRefineFace(                         &
              & modalRepFace  = aFace,                           &
              & modg_basis    = scheme(currentLevel)%modg_basis, &
              & schemeCoarse  = scheme(currentLevel)%modg_2d,    &
              & schemeFine    = scheme(currentLevel+1)%modg_2d,  &
              & refineDir     = 1,                               &
              & fineFaceShift = 2,                               &
              & modalRefined  = secondFace                       )

            ! Assign the interpolated data to the child faces.
            coff = (yshift-1)*2
            do i=1,nDofsFine*nScalars
              iDof = mod(i-1,nDofsFine)+1
              iVar = (i-1)/nDofsFine + 1
              facedata(currentLevel+1)%faceRep(iDir)           &
                &                     %dat(childPos(1+coff),   &
                &                          iDof,iVar,invAlign) &
                &  = firstFace(iDof,iVar)
              facedata(currentLevel+1)%faceRep(iDir)           &
                &                     %dat(childPos(2+coff),   &
                &                          iDof,iVar,invAlign) &
                &  = secondFace(iDof,iVar)
            end do
          end do

          if (elemPos <= nFluids) then
            call tem_stopTimer( me          = atl_elemTimers, &
              &                 timerHandle = totPos          )
          end if

        end do
      end do
    end do

  end subroutine atl_modg_2d_coarseToFineFace



  !> summary: Project modal representation on a face to a semi-refined face (i.e. on a face
  !! that is refined in one of the spatial directions.
  !!
  !! Project modal representation on a face to a semi-refined face (i.e. on a face
  !! that is refined in one of the spatial directions. The result is a modal representation
  !! on the semi-refined element. \n
  !! \n
  !! The function is executing one of the following projections: \n
  !!
  !!       face on current                      semi-refined                            \n
  !!            level                               face                                \n
  !!
  !!  ------------------------             ------------------------                     \n
  !!  |                      |             |          |           |                     \n
  !!  |                      |  refineDir  |          |           |                     \n
  !!  |                      |      == 1   | fineFace | fineFace  |                     \n
  !!  |                      | -------->>  |  Shift   |  Shift    |                     \n
  !!  |                      |             |   == 1   |   == 2    |                     \n
  !!  |                      |             |          |           |                     \n
  !!  |                      |             |          |           |                     \n
  !!  ------------------------             ------------------------                     \n
  !!
  !!  or:
  !!
  !!  ------------------------             ------------------------                     \n
  !!  |                      |             |                      |                     \n
  !!  |                      |  refineDir  |     fineFaceShift    |                     \n
  !!  |                      |      == 2   |        == 2          |                     \n
  !!  |                      | -------->>  ------------------------                     \n
  !!  |                      |             |                      |                     \n
  !!  |                      |             |     fineFaceShift    |                     \n
  !!  |                      |             |        == 1          |                     \n
  !!  ------------------------             ------------------------                     \n
  !!
  subroutine modg_2d_semiRefineFace( modalRepFace, modg_basis, schemeCoarse, &
    &                                schemeFine, refineDir, fineFaceShift,   &
    &                                modalRefined                            )
    ! --------------------------------------------------------------------------
    !> Modal representation of a function on the non-refined face.
    !! Dimensions are: (modg%maxPolyDegree+1)^2 for the first dimension
    !! and nScalars for the second dimension.
    real(kind=rk), intent(in) :: modalRepFace(:,:)
    !> Informations about the polynomial basis of the MODG scheme.
    type(ply_modg_basis_type) :: modg_basis
    !> The parameters of your MODG scheme on the coarse level.
    type(atl_modg_2d_scheme_type), intent(in) :: schemeCoarse
    !> The parameters of your MODG scheme on the fine level.
    type(atl_modg_2d_scheme_type), intent(in) :: schemeFine
    !> The direction of the semi-refinement. Either 1 or 2. Have a look at the
    !! function description.
    integer, intent(in) :: refineDir
    !> The semi-refined element you want to obtain.
    integer, intent(in) :: fineFaceShift
    !> The modal representation of modalRepFace restricted to the semi-refined
    !! element.
    real(kind=rk), intent(inout) :: modalRefined(:,:)
    ! --------------------------------------------------------------------------
    integer :: iFunc, iCoarseFunc, comFunc
    integer :: funcPos, coarsePos
    real(kind=rk) :: fineSqNorm
    ! --------------------------------------------------------------------------

    ! Loop over all the ansatz functions of the semi-refined element and calculate
    ! all the modal coefficients for it.

    ! We treat the refinement directions separately to have the if condition
    ! outside of the loop.
    if( refineDir .eq. 1 ) then

      ! ... loop over the ansatz function that are common for the fine and current element
      comFunc = 1
      ! ... loop over the ansatz function of the semi-refined element.
      do iFunc = 1, schemeFine%maxPolyDegree+1
  funcpos = ifunc                                      &
    &      + ( ( comfunc-1)                             &
    &      + (1-1)*(schemefine%maxpolydegree+1))*(schemefine%maxpolydegree+1)
        fineSqNorm = 2.0_rk /(2.0_rk * iFunc - 1.0_rk)
        ! ... loop over the ansatz functions of the current (non-refined) element
        do iCoarseFunc = 1, schemeCoarse%maxPolyDegree+1
          coarsePos = 1                                      &
            & + (iCoarseFunc-1)                              &
            & + (comFunc-1) * (schemeCoarse%maxPolyDegree + 1)
          ! Weight and add it to the semi refined face representation
          modalRefined(funcPos,:) = modalRefined(funcPos,:)                &
            & + modg_basis%refineBaseCoeff                                 &
            &             %anz_anzShift(iFunc, iCoarseFunc, fineFaceShift) &
            & * modalRepFace(coarsePos,:) / ( fineSqNorm )
        end do
      end do

    else

      ! ... loop over the ansatz function that are common for the fine and current element
      do comFunc = 1, schemeCoarse%maxPolyDegree+1
        ! ... loop over the ansatz function of the semi-refined element.
        iFunc = 1
        funcPos = 1 + (comFunc-1) + (iFunc-1)*(schemeCoarse%maxPolyDegree+1)
        fineSqNorm = 2.0_rk /(2.0_rk * iFunc - 1.0_rk)
        ! ... loop over the ansatz functions of the current (non-refined) element
        iCoarseFunc = 1
  coarsepos = comfunc                                      &
    &      + ( ( icoarsefunc-1)                             &
    &      + (1-1)*(schemecoarse%maxpolydegree+1))*(schemecoarse%maxpolydegree+1)
        modalRefined(funcPos,:) = modalRefined(funcPos,:)                &
          & + modg_basis%refineBaseCoeff                                 &
          &             %anz_anzShift(iFunc, iCoarseFunc, fineFaceShift) &
          & * modalRepFace(coarsePos,:) / ( fineSqNorm )
      end do

    end if


  end subroutine modg_2d_semiRefineFace




  !> summary: Interpolate modal face representation from next finer faces to coarse level (level
  !! difference between coarser and finer faces has to be 1).
  !!
  !! Interpolates functions defined on finer faces to faces of the current level.\n
  !! \n
  !!      faces on fine                            face on current
  !!          level                                     level
  !! ------------------------                 ------------------------
  !! |          |           |                 |                      |
  !! |    3     |     4     |                 |                      |
  !! |          |           |                 |                      |
  !! ------------------------   ---------->>  |          5           |
  !! |          |           |                 |                      |
  !! |    1     |     2     |                 |                      |
  !! |          |           |                 |                      |
  !! ------------------------                 ------------------------
  !! \n
  !! This is accomplished with lower complexity (with respect to the polynomial
  !! degree) by a dimension by dimension approach: \n
  !! \n
  !!      faces on fine                            face on current           \n
  !!          level                                     level                \n
  !! ------------------------                 ------------------------       \n
  !! |          |           |                 |                      |       \n
  !! |    3     |     4     |                 |                      |       \n
  !! |          |           |                 |                      |       \n
  !! ------------------------   ---------->>  |          5           |       \n
  !! |          |           |                 |                      |       \n
  !! |    1     |     2     |                 |                      |       \n
  !! |          |           |                 |                      |       \n
  !! ------------------------                 ------------------------       \n
  !!                 \                               /                       \n
  !!                  \                             /                        \n
  !!                   \                           /                         \n
  !!                    \                         /                          \n
  !!                     ------------------------                            \n
  !!                     |                      |                            \n
  !!                     |          b           |                            \n
  !!                     |                      |                            \n
  !!                     ------------------------                            \n
  !!                     |                      |                            \n
  !!                     |          a           |                            \n
  !!                     |                      |                            \n
  !!                     ------------------------                            \n
  !!
  subroutine atl_modg_2d_fineToCoarseFace( minLevel, maxLevel, currentLevel, &
                                  & mesh, facedata, scheme, nScalars)
    ! --------------------------------------------------------------------------
    !> The minumum level of the mesh.
    integer, intent(in) :: minLevel
    !> The maximum level of the mesh.
    integer, intent(in) :: maxLevel
    !> The current level (i.e. the coarse level).
    integer, intent(in) :: currentLevel
    !> The mesh representation.
    type(atl_cube_elem_type), intent(in) :: mesh(minLevel:maxLevel)
    !> The face representations (finer faces are interpolated from coarser ones).
    type(atl_facedata_type), intent(inout) :: facedata(minLevel:maxLevel)
    !> The schemes on the different levels.
    type(atl_scheme_type), intent(in) :: scheme(minLevel:maxLevel)
    !> The number of scalar variables in your equation system.
    integer, intent(in) :: nScalars
    ! --------------------------------------------------------------------------
    integer :: iDir, iAlign, iFace, invAlign
    integer :: i, iVar, iDof, yshift
    integer :: cp1, cp2
    integer :: nDofsCoarse, nDofsFine, nDofsInter
    integer :: nFluids
    integer :: elemPos, childPos(4)
    integer :: totpos
    real(kind=rk), allocatable :: faceDat(:,:), aFace(:,:)
    real(kind=rk), allocatable :: fineDat(:,:)
    ! --------------------------------------------------------------------------

    nfluids = mesh(currentLevel)%descriptor%elem%nElems(eT_fluid)

    ! The number of degrees of freedom on a face
    nDofsCoarse = scheme(currentLevel)%modg_2d%maxPolyDegree+1
    nDofsFine = scheme(currentLevel+1)%modg_2d%maxPolyDegree+1
    nDofsInter = (scheme(currentLevel)%modg_2d%maxPolyDegree+1)*1

    ! Create the intermediate and resulting arrays
    allocate( faceDat(nDofsCoarse,nScalars), &
            & aFace(nDofsInter,nScalars),    &
            & fineDat(nDofsFine,nScalars)    )

    ! Iterate over all the faces and project from the faces of the finer level
    ! to the current level.
    do iDir = 1,2
      do iAlign = 1,2
        do iFace = 1, size(mesh(currentLevel)%faces%faces(iDir)%fromFinerFace(iAlign)%elemPos)
          ! get elempos
          elemPos = mesh(currentLevel)%faces%faces(iDir)%fromFinerFace(iAlign)%elemPosOp(iFace)

          if (elempos <= nFluids) then
            totpos = mesh(currentLevel)%descriptor%pntTID(elempos)
            ! start element wise timer for LB weights
            call tem_startTimer( me          = atl_elemTimers, &
              &                  timerHandle = totPos          )
          end if


          ! The face we have to interpolate. If the face is refined from its left element
          ! (i.e. the right face of the element element is refined, -> iAlign == 2), we have
          ! to work on the left face of the right element.
          invAlign = tem_invFace_map(iAlign)

          ! Get the data of the faces 1,2,3 and 4
          ! ... face 1
          do i=1,4
            childPos(i) = mesh(currentLevel)%faces%faces(iDir)           &
              &                                   %fromFinerFace(iAlign) &
              &                                   %childPosOp(i,iFace)
          end do

          do i=1,nDofsCoarse*nScalars
            iDof = mod(i-1,nDofsCoarse)+1
            iVar = (i-1)/nDofsCoarse + 1
            faceDat(iDof,iVar) = 0.0_rk
          end do

          do yshift=0,1

            do i=1,nDofsInter*nScalars
              iDof = mod(i-1,nDofsInter)+1
              iVar = (i-1)/nDofsInter + 1
              aFace(iDof,iVar) = 0.0_rk
            end do

            cp1 = childpos(1 + yshift*2)
            cp2 = childpos(2 + yshift*2)
            do i=1,nDofsFine*nScalars
              iDof = mod(i-1,nDofsFine)+1
              iVar = (i-1)/nDofsFine + 1
              fineDat(iDof,iVar) = facedata(currentLevel+1)%faceFlux(iDir) &
                &                   %dat(cp1,iDof,iVar,invAlign)
            end do

            ! ... project 1 to a
            call modg_2d_semiCoarseFace(                         &
              & modalRepFace  = fineDat,                         &
              & modg_basis    = scheme(currentLevel)%modg_basis, &
              & schemeCoarse  = scheme(currentLevel)%modg_2d,    &
              & schemeFine    = scheme(currentLevel+1)%modg_2d,  &
              & coarseDir     = 1,                               &
              & fineFaceShift = 1,                               &
              & modalCoarsed  = aFace                            )

            do i=1,nDofsFine*nScalars
              iDof = mod(i-1,nDofsFine)+1
              iVar = (i-1)/nDofsFine + 1
              fineDat(iDof,iVar) = facedata(currentLevel+1)%faceFlux(iDir) &
                &                   %dat(cp2,iDof,iVar,invAlign)
            end do

            ! ... project 2 to a
            call modg_2d_semiCoarseFace(                         &
              & modalRepFace  = fineDat,                         &
              & modg_basis    = scheme(currentLevel)%modg_basis, &
              & schemeCoarse  = scheme(currentLevel)%modg_2d,    &
              & schemeFine    = scheme(currentLevel+1)%modg_2d,  &
              & coarseDir     = 1,                               &
              & fineFaceShift = 2,                               &
              & modalCoarsed  = aFace                            )

            ! ... project a to coarse
            call modg_2d_semiCoarseFace(                         &
              & modalRepFace  = aFace,                           &
              & modg_basis    = scheme(currentLevel)%modg_basis, &
              & schemeCoarse  = scheme(currentLevel)%modg_2d,    &
              & schemeFine    = scheme(currentLevel+1)%modg_2d,  &
              & coarseDir     = 2,                               &
              & fineFaceShift = 1+yshift,                        &
              & modalCoarsed  = faceDat                          )
          end do

          do i=1,nDofsCoarse*nScalars
            iDof = mod(i-1,nDofsCoarse)+1
            iVar = (i-1)/nDofsCoarse + 1
            facedata(currentLevel)%faceFlux(iDir)%dat(elemPos,iDof,iVar,invAlign) = faceDat(iDof,iVar)
          end do

          if (elempos <= nFluids) then
            ! start element wise timer for LB weights
            call tem_stopTimer( me          = atl_elemTimers, &
              &                 timerHandle = totPos          )
          end if

        end do
      end do
    end do

  end subroutine atl_modg_2d_fineToCoarseFace



  !> summary: Project modal representation of a semi-refined face (i.e. on a face
  !! that is refined in one of the spatial directions) to a coarse element.
  !!
  !! Project modal representation of a semi-refined face (i.e. on a face
  !! that is refined in one of the spatial directions) to a coarser element.
  !! The result is a modal representation on the coarse element that is representing
  !! the polynomial function of the refined element on the coarse element when restricted
  !! to the refined element again. However the modal representation is done in terms
  !! of the ansatz functions of the coarse element. \n
  !! \n
  !! The function is executing one of the following projections: \n
  !!
  !!       face on current                      semi-refined                            \n
  !!            level                               face                                \n
  !!                                                                                    \n
  !!  ------------------------             ------------------------                     \n
  !!  |                      |             |          |           |                     \n
  !!  |                      |  coarseDir  |          |           |                     \n
  !!  |                      |      == 1   | fineFace | fineFace  |                     \n
  !!  |                      | <<--------  |  Shift   |  Shift    |                     \n
  !!  |                      |             |   == 1   |   == 2    |                     \n
  !!  |                      |             |          |           |                     \n
  !!  |                      |             |          |           |                     \n
  !!  ------------------------             ------------------------                     \n
  !!
  !!  or:
  !!
  !!  ------------------------             ------------------------                     \n
  !!  |                      |             |                      |                     \n
  !!  |                      |  coarseDir  |     fineFaceShift    |                     \n
  !!  |                      |      == 2   |        == 2          |                     \n
  !!  |                      | <<--------  ------------------------                     \n
  !!  |                      |             |                      |                     \n
  !!  |                      |             |     fineFaceShift    |                     \n
  !!  |                      |             |        == 1          |                     \n
  !!  ------------------------             ------------------------                     \n
  !!
  subroutine modg_2d_semiCoarseFace( modalRepFace, modg_basis, schemeCoarse, &
    &                                schemeFine, coarseDir, fineFaceShift,   &
    &                                modalCoarsed                            )
    ! --------------------------------------------------------------------------
    !> Modal representation of a function on one of refined face. Which fine face
    !! is determined by fineFaceShift
    !! Dimensions are: (modg%maxPolyDegree+1)^2 for the first dimension
    !! and nScalars for the second dimension.
    real(kind=rk), intent(in) :: modalRepFace(:,:)
    !> Informations about the polynomial basis of the MODG scheme.
    type(ply_modg_basis_type) :: modg_basis
    !> The parameters of your MODG scheme on the coarse level.
    type(atl_modg_2d_scheme_type), intent(in) :: schemeCoarse
    !> The parameters of your MODG scheme on the fint level.
    type(atl_modg_2d_scheme_type), intent(in) :: schemeFine
    !> The direction of the semi-coarsening. Either 1 or 2. Have a look at the
    !! function description.
    integer, intent(in) :: coarseDir
    !> The semi-refined element you want to obtain.
    integer, intent(in) :: fineFaceShift
    !> The modal representation of modalRepFace on the coarser element, restricted
    !! to the given fine element.
    real(kind=rk), intent(inout) :: modalCoarsed(:,:)
    ! --------------------------------------------------------------------------
    integer :: comFunc, iFunc, iCoarseFunc
    integer :: funcPos, coarsePos
    real(kind=rk) :: coarseSqNorm, jacobiDetFineToCoarse
    ! --------------------------------------------------------------------------


    ! Ratio of determinants of the Jacobians of the mapping from fine to reference and
    ! coarse to reference (one-dimensional)
    jacobiDetFineToCoarse = 0.5_rk

    ! Check for the spatial direction of the coarsening.
    if(coarseDir.eq.1)  then

      ! ... loop over the ansatz function that are common for the fine and current element
      comFunc = 1
      ! ... loop over the ansatz functions of the current (non-refined) element
      do iCoarseFunc = 1, schemeCoarse%maxPolyDegree+1
        coarsePos = 1 &
          & + (iCoarseFunc-1) &
          & + (comFunc-1) * (schemeCoarse%maxPolyDegree + 1)
        coarseSqNorm = 2.0_rk / ( 2.0_rk * iCoarseFunc - 1.0_rk)
        ! ... loop over the ansatz function of the semi-refined element.
        do iFunc = 1, schemeFine%maxPolyDegree+1
  funcpos = ifunc                                      &
    &      + ( ( comfunc-1)                             &
    &      + (1-1)*(schemefine%maxpolydegree+1))*(schemefine%maxpolydegree+1)
          modalCoarsed(coarsePos,:) = modalCoarsed(coarsePos,:)            &
            & + modg_basis%refineBaseCoeff                                 &
            &             %anz_anzShift(iFunc, iCoarseFunc, fineFaceShift) &
            & * modalRepFace(funcPos,:)                                    &
            & * jacobiDetFineToCoarse / ( coarseSqNorm )
        end do
      end do

    else

      ! ... loop over the ansatz function that are common for the fine and current element
      do comFunc = 1, schemeCoarse%maxPolyDegree+1
        ! ... loop over the ansatz functions of the current (non-refined) element
        iCoarseFunc = 1
  coarsepos = comfunc                                      &
    &      + ( ( icoarsefunc-1)                             &
    &      + (1-1)*(schemecoarse%maxpolydegree+1))*(schemecoarse%maxpolydegree+1)
        coarseSqNorm = 2.0_rk / ( 2.0_rk * iCoarseFunc - 1.0_rk)
        ! ... loop over the ansatz function of the semi-refined element.
        iFunc = 1
        funcPos = 1 + (comFunc-1) + (iFunc-1)*(schemeCoarse%maxPolyDegree+1)
        modalCoarsed(coarsePos,:) = modalCoarsed(coarsePos,:)                  &
          & + modg_basis%refineBaseCoeff                                       &
          &             %anz_anzShift(iFunc,iCoarseFunc,fineFaceShift)         &
          & * modalRepFace(funcPos,:) * jacobiDetFineToCoarse / ( coarseSqNorm )
      end do

    end if

  end subroutine modg_2d_semiCoarseFace

  !> Project data from 8 smaller elements to its parent element in terms
  !! of L2 projections.
  subroutine atl_modg_fineToCoarseElem_2d( minLevel, maxLevel, currentLevel, &
    &                                      iDir, mesh, state_stab, scheme,   &
    &                                      nScalars                          )
    ! --------------------------------------------------------------------------
    !> The minumum level of the mesh.
    integer, intent(in) :: minLevel
    !> The maximum level of the mesh.
    integer, intent(in) :: maxLevel
    !> The current level (i.e. the coarse level).
    integer, intent(in) :: currentLevel
    !> The direction to interpolate.
    integer, intent(in) :: iDir
    !> The mesh representation.
    type(atl_cube_elem_type), intent(in) :: mesh(minLevel:maxLevel)
    !> The face representations (finer faces are interpolated from coarser ones).
    type(atl_statedata_type), intent(inout) :: state_stab(minLevel:maxLevel,1:3)
    !> The schemes on the different levels.
    type(atl_scheme_type), intent(in) :: scheme(minLevel:maxLevel)
    !> The number of scalar variables in your equation system.
    integer, intent(in) :: nScalars
    ! --------------------------------------------------------------------------
    integer :: iElem, iChild
    integer :: iRefineX, iRefineY
    integer :: nFluids
    integer :: nDofsCoarse, nDofsFine, nDofsInter_firstCoarse, &
      &        nDofsInter_secondCoarse
    integer :: elemPos, childPos
    integer :: totpos
    real(kind=rk) :: starttime, stoptime
    real(kind=rk), allocatable :: faceDat(:,:)
    real(kind=rk), allocatable :: childFace(:,:,:)
    real(kind=rk), allocatable :: firstCoarse(:,:,:), secondCoarse(:,:,:)
    ! --------------------------------------------------------------------------

    nfluids = mesh(currentLevel)%descriptor%elem%nElems(eT_fluid)

    ! The number of degrees of freedom on a face
    nDofsCoarse = scheme(currentLevel)%nDofs
    nDofsFine = scheme(currentLevel+1)%nDofs
    nDofsInter_firstCoarse =  (scheme(currentLevel+1)%modg_2d%maxPolyDegree+1) &
                           & *(scheme(currentLevel+1)%modg_2d%maxPolyDegree+1)
    nDofsInter_secondCoarse =  (scheme(currentLevel+1)%modg_2d%maxPolyDegree+1) &
                            & *(scheme(currentLevel)%modg_2d%maxPolyDegree+1)

    ! Create the intermediate and resulting arrays
    allocate( childFace(nDofsFine,nScalars,8) )
    allocate( firstCoarse(nDofsInter_firstCoarse,nScalars,4) )
    allocate( secondCoarse(nDofsInter_firstCoarse,nScalars,2) )
    allocate( faceDat(nDofsCoarse,nScalars) )

    ! Iterate over all the from finer elements and project from the elements of the finer level
    ! to the current level.
    do iElem = 1, mesh(currentLevel)%faces_stab              &
      &                             %dimByDimDesc(iDir)      &
      &                             %elem                    &
      &                             %nElems(eT_ghostFromFiner)
      elemPos = iElem + mesh(currentLevel)%faces_stab                &
        &                                 %dimByDimDesc(iDir)        &
        &                                 %elem                      &
        &                                 %nElems(eT_fluid)          &
        &             + mesh(currentLevel)%faces_stab                &
        &                                 %dimByDimDesc(iDir)        &
        &                                 %elem                      &
        &                                 %nElems(eT_ghostFromCoarser)

      starttime = MPI_Wtime()

      do iChild = 1,8
        childPos = mesh(currentLevel)%faces_stab%dimByDimDesc(iDir)%depFromFiner(iElem)%elem%val(iChild)
        childFace(:,:,iChild) = state_stab(currentLevel+1,iDir)%state(childPos,:,:)
      end do

      ! First coarsening step in z direction (just average in this direction)
      firstCoarse(:,:,:) = 0.0_rk
      do iRefineY = 1,2
        do iRefineX = 1,2
          firstCoarse(:,:,(iRefineY-1)*2+iRefineX) = &
               & ( childFace(:,:,(iRefineY-1)*2+iRefineX) &
               & + childFace(:,:,4+(iRefineY-1)*2+iRefineX) )/2.0_rk
        end do
      end do

      ! Second coarsening step in y direction
      secondCoarse(:,:,:) = 0.0_rk
      do iRefineY = 1,2
        do iRefineX = 1,2
          call modg_semiCoarseElem_2d(                                  &
            & modalRepFace  = firstCoarse(:,:,(iRefineY-1)*2+iRefineX), &
            & modg_basis    = scheme(currentLevel)%modg_basis,          &
            & schemeCoarse  = scheme(currentLevel)%modg_2d,             &
            & schemeFine    = scheme(currentLevel+1)%modg_2d,           &
            & coarseDir     = 2,                                        &
            & fineElemShift = iRefineY,                                 &
            & modalCoarsed  = secondCoarse(:,:,iRefineX)                )
        end do
      end do

      ! Third coarsening step in x direction
      do iChild=lbound(faceDat,2),ubound(faceDat,2)
        faceDat(:,iChild) = 0.0_rk
      end do
      do iRefineX = 1,2
        call modg_semiCoarseElem_2d(                         &
          & modalRepFace  = secondCoarse(:,:,iRefineX),      &
          & modg_basis    = scheme(currentLevel)%modg_basis, &
          & schemeCoarse  = scheme(currentLevel)%modg_2d,    &
          & schemeFine    = scheme(currentLevel+1)%modg_2d,  &
          & coarseDir     = 1,                               &
          & fineElemShift = iRefineX,                        &
          & modalCoarsed  = faceDat(:,:)                     )
      end do

      ! Assign the data (same for both sides of the face)
      state_stab(currentLevel,iDir)%state(elemPos,:,:) = faceDat(:,:)
      stoptime = MPI_Wtime()

      do iChild = 1,4
        childPos = mesh(currentLevel)%faces_stab%dimByDimDesc(iDir)  &
          &                                     %depFromFiner(iElem) &
          &                                     %elem%val(iChild)

        if ( childpos <= nFluids ) then
          totpos = mesh(currentLevel)%descriptor%pntTID(childpos)
          atl_elemTimers%duration%val(totPos)        &
            &  = atl_elemTimers%duration%val(totPos) &
            &    + 0.25_rk*(stoptime-starttime)
        end if
      end do

    end do

  end subroutine atl_modg_fineToCoarseElem_2d


  !> Subroutine to semi-coarsen an element with modal polynomial representation
  !! to its semi-parent.
  subroutine modg_semiCoarseElem_2d( modalRepFace, modg_basis, schemeCoarse, &
    &                                schemeFine, coarseDir, fineElemShift,   &
    &                                modalCoarsed                            )
    ! --------------------------------------------------------------------------
    !> Modal representation of a function on one of refined element.
    !! Which fine element is determined by fineFaceShift
    !! Dimensions are: (modg%maxPolyDegree+1)^3 for the first dimension
    !! and nScalars for the second dimension.
    real(kind=rk), intent(in) :: modalRepFace(:,:)
    !> Informations about the polynomial basis of the MODG scheme.
    type(ply_modg_basis_type) :: modg_basis
    !> The parameters of your MODG scheme on the coarse level.
    type(atl_modg_2d_scheme_type), intent(in) :: schemeCoarse
    !> The parameters of your MODG scheme on the fint level.
    type(atl_modg_2d_scheme_type), intent(in) :: schemeFine
    !> The direction of the semi-coarsening. Either 1 or 2 or 3. Have a look at the
    !! function description.
    integer, intent(in) :: coarseDir
    !> The semi-refined element you want to obtain.
    integer, intent(in) :: fineElemShift
    !> The modal representation of modalRepFace on the coarser element, restricted
    !! to the given fine element.
    real(kind=rk), intent(inout) :: modalCoarsed(:,:)
    ! --------------------------------------------------------------------------
    integer :: iFunc, iDegX, iDegY
    integer :: funcPos, dof
    real(kind=rk) :: coarseSqNorm, jacobiDetFineToCoarse
    integer :: mpd1, mpd1_square
    ! --------------------------------------------------------------------------

    ! Ratio of determinants of the Jacobians of the mapping from fine to reference and
    ! coarse to reference (one-dimensional)
    jacobiDetFineToCoarse = 0.5_rk

    ! Check for the spatial direction of the coarsening.
    if(coarseDir.eq.1)  then

      mpd1 = schemeCoarse%maxPolyDegree+1
      mpd1_square = (schemeFine%maxPolyDegree+1)*mpd1

      do dof = 1, mpd1_square
        iDegX = mod(dof-1,mpd1)+1
        iDegY = (dof-1)/mpd1+1
        coarseSqNorm = 2.0_rk / ( 2.0_rk * iDegX - 1.0_rk)
        ! ... loop over the ansatz function of the semi-refined element.
        do iFunc = 1, schemeFine%maxPolyDegree+1
  funcpos = ifunc                                      &
    &      + ( ( idegy-1)                             &
    &      + (1-1)*(schemefine%maxpolydegree+1))*(schemefine%maxpolydegree+1)
          modalCoarsed(dof,:) = modalCoarsed(dof,:)                  &
            & + modg_basis%refineBaseCoeff                           &
            &             %anz_anzShift(iFunc, iDegX, fineElemShift) &
            & * modalRepFace(funcPos,:)                              &
            & * jacobiDetFineToCoarse / ( coarseSqNorm )
        end do
      end do

    else

      mpd1 = schemeCoarse%maxPolyDegree+1
      mpd1_square = (schemeCoarse%maxPolyDegree+1)*mpd1

      do dof = 1, mpd1_square
        iDegX = mod(dof-1,mpd1)+1
        iDegY = (dof-1)/mpd1+1
        coarseSqNorm = 2.0_rk / ( 2.0_rk * iDegY - 1.0_rk)
        ! ... loop over the ansatz function of the semi-refined element.
        do iFunc = 1, schemeFine%maxPolyDegree+1
  funcpos = idegx                                      &
    &      + ( ( ifunc-1)                             &
    &      + (1-1)*(schemefine%maxpolydegree+1))*(schemefine%maxpolydegree+1)
          modalCoarsed(dof,:) = modalCoarsed(dof,:)                  &
            & + modg_basis%refineBaseCoeff                           &
            &             %anz_anzShift(iFunc, iDegY, fineElemShift) &
            & * modalRepFace(funcPos,:)                              &
            & * jacobiDetFineToCoarse / ( coarseSqNorm )
        end do
      end do

    end if

  end subroutine modg_semiCoarseElem_2d


  !> Project coarse parent element to its 8 finer child elements
  !! by a simple L2 projection.
  subroutine atl_modg_coarseToFineElem_2d( minLevel, maxLevel, currentLevel, &
    &                                      iDir, mesh, state_stab, scheme,   &
    &                                      nScalars                          )
    ! --------------------------------------------------------------------------
    !> The minumum level of the mesh.
    integer, intent(in) :: minLevel
    !> The maximum level of the mesh.
    integer, intent(in) :: maxLevel
    !> The current level (i.e. the coarse level).
    integer, intent(in) :: currentLevel
    !> The direction to project
    integer, intent(in) :: iDir
    !> The mesh representation.
    type(atl_cube_elem_type), intent(in) :: mesh(minLevel:maxLevel)
    !> The face representations (finer faces are interpolated from coarser ones).
    type(atl_statedata_type), intent(inout) :: state_stab(minLevel:maxLevel,1:3)
    !> The schemes on the different levels.
    type(atl_scheme_type), intent(in) :: scheme(minLevel:maxLevel)
    !> The number of scalar variables in your equation system.
    integer, intent(in) :: nScalars
    ! --------------------------------------------------------------------------
    integer :: iElem, childPos, childNum, parentPos, ii
    integer :: totpos
    integer :: iRefineX, iRefineY, iRefineZ
    integer :: nDofsCoarse, nDofsFine
    integer :: nDofsInter_firstRefine, nDofsInter_secondRefine
    integer :: nFluids
    real(kind=rk), allocatable :: faceDat(:,:)
    real(kind=rk), allocatable :: firstRefine(:,:), secondRefine(:,:)
    ! --------------------------------------------------------------------------

    ! The number of degrees of freedom on an element for the coarse, fine and intermediate level
    nDofsCoarse = scheme(currentLevel-1)%nDofs
    nDofsFine = scheme(currentLevel)%nDofs
    nDofsInter_firstRefine = &
                         &   (scheme(currentLevel)%modg_2d%maxPolyDegree+1) &
                         & * (scheme(currentLevel-1)%modg_2d%maxPolyDegree+1)
    nDofsInter_secondRefine = &
                         &   (scheme(currentLevel)%modg_2d%maxPolyDegree+1) &
                         & * (scheme(currentLevel)%modg_2d%maxPolyDegree+1)

    nfluids = mesh(currentLevel)%descriptor%elem%nElems(eT_fluid)

    ! Create the intermediate and resulting arrays
    allocate( faceDat(nDofsCoarse,nScalars) )
    allocate( firstRefine(nDofsInter_firstRefine,nScalars) )
    allocate( secondRefine(nDofsInter_secondRefine,nScalars) )

    do iElem = 1, mesh(currentLevel)%faces_stab                &
      &                             %dimByDimDesc(iDir)        &
      &                             %elem                      &
      &                             %nElems(eT_ghostFromCoarser)

      ! Get position and number of the fine (child) element
      childPos = iElem + mesh(currentLevel)%faces_stab         &
        &                                  %dimByDimDesc(iDir) &
        &                                  %elem               &
        &                                  %nElems(eT_fluid)
      childNum = tem_childNumber(mesh(currentLevel)%faces_stab         &
        &                                          %dimByDimDesc(iDir) &
        &                                          %total(childPos)    )

      ! For 2D the if childpos would return childnum > 4 then it should
      ! just project the corresponding child from 1-4
      childNum = mod(childNum-1, 4) + 1

      ! Get position of the parent element
      parentPos = mesh(currentLevel)%faces_stab%dimByDimDesc(iDir)%depFromCoarser(iElem)%elem%val(1)

      if (parentPos <= nfluids) then
        totpos = mesh(currentLevel)%descriptor%pntTID(parentpos)
        ! start timer for LB of parent element
        call tem_startTimer( me          = atl_elemTimers, &
          &                  timerHandle = totPos          )
      end if

      ! Get state of the coarser element

      faceDat(:,:) = state_stab(currentLevel-1,iDir)%state(parentPos,:,:)

      ! Get the fine element shift from the child number
      iRefineZ = (childNum-1)/4 + 1
      iRefineY = (childNum-1-(iRefineZ-1)*4)/2+1
      iRefineX = mod(childNum-1,2)+1

      ! Do the first refinement step in x direction
      do ii = lbound(firstRefine,2),ubound(firstRefine,2)
        firstRefine(:, ii) = 0.0_rk
      end do
      call modg_semiRefineElem_2d(                         &
        & modalRepFace  = faceDat,                         &
        & modg_basis    = scheme(currentLevel)%modg_basis, &
        & schemeCoarse  = scheme(currentLevel-1)%modg_2d,  &
        & schemeFine    = scheme(currentLevel)%modg_2d,    &
        & refineDir     = 1,                               &
        & fineElemShift = iRefineX,                        &
        & modalRefined  = firstRefine(:,:)                 )

      ! Do the second refinement step in y direction
      do ii=lbound(secondRefine,2),ubound(secondRefine,2)
        secondRefine(:,ii) = 0.0_rk
      end do
      call modg_semiRefineElem_2d(                         &
        & modalRepFace  = firstRefine(:,:),                &
        & modg_basis    = scheme(currentLevel)%modg_basis, &
        & schemeCoarse  = scheme(currentLevel-1)%modg_2d,  &
        & schemeFine    = scheme(currentLevel)%modg_2d,    &
        & refineDir     = 2,                               &
        & fineElemShift = iRefineY,                        &
        & modalRefined  = secondRefine(:,:)                )

      ! Assign element state for the fine (child) element
      state_stab(currentLevel,iDir)%state(childPos,:,:) = secondRefine(:,:)

      if (parentPos <= nFluids) then
        call tem_stopTimer( me          = atl_elemTimers, &
          &                 timerHandle = totPos          )
      end if

    end do

  end subroutine atl_modg_coarseToFineElem_2d


  !> Subroutine to semi-refine an element with modal polynomial representation
  !! into its semi-children.
  subroutine modg_semiRefineElem_2d( modalRepFace, modg_basis, schemeCoarse, &
    &                                schemeFine, refineDir, fineElemShift,   &
    &                                modalRefined                            )
    ! --------------------------------------------------------------------------
    !> Modal representation of a function on the non-refined face.
    !! Dimensions are: (modg%maxPolyDegree+1)^2 for the first dimension
    !! and nScalars for the second dimension.
    real(kind=rk), intent(in) :: modalRepFace(:,:)
    !> Informations about the polynomial basis of the MODG scheme.
    type(ply_modg_basis_type) :: modg_basis
    !> The parameters of your MODG scheme on the coarse level.
    type(atl_modg_2d_scheme_type), intent(in) :: schemeCoarse
    !> The parameters of your MODG scheme on the fine level.
    type(atl_modg_2d_scheme_type), intent(in) :: schemeFine
    !> The direction of the semi-refinement. Either 1 or 2. Have a look at the
    !! function description.
    integer, intent(in) :: refineDir
    !> The semi-refined element you want to obtain.
    integer, intent(in) :: fineElemShift
    !> The modal representation of modalRepFace restricted to the semi-refined
    !! element.
    real(kind=rk), intent(inout) :: modalRefined(:,:)
    ! --------------------------------------------------------------------------
    integer :: iDegX, iDegY, iCoarseFunc
    integer :: dof, coarsePos
    real(kind=rk) :: fineSqNorm
    integer :: mpd1, mpd1_square
    ! --------------------------------------------------------------------------

    ! We treat the refinement directions separately to have the if condition
    ! outside of the loop.
    if( refineDir .eq. 1 ) then

      mpd1 = schemeCoarse%maxPolyDegree+1
      mpd1_square = mpd1 * (schemeFine%maxPolyDegree+1)

      ! loop over all degrees of freedoms on the (semi-)refined element
      do dof = 1, mpd1_square
        iDegX = mod(dof-1,mpd1)+1
        iDegY = (dof-1)/mpd1+1
        fineSqNorm = 2.0_rk /(2.0_rk * iDegX - 1.0_rk)
        ! ... loop over the ansatz functions of the current (non-refined) element
        do iCoarseFunc = 1, schemeCoarse%maxPolyDegree+1
  coarsepos = icoarsefunc                                      &
    &      + ( ( idegy-1)                             &
    &      + (1-1)*(schemecoarse%maxpolydegree+1))*(schemecoarse%maxpolydegree+1)
          ! Weight and add it to the semi refined face representation
          modalRefined(dof,:) = modalRefined(dof,:)                        &
            & + modg_basis%refineBaseCoeff                                 &
            &             %anz_anzShift(iDegX, iCoarseFunc, fineElemShift) &
            & * modalRepFace(coarsePos,:) / ( fineSqNorm )
        end do
      end do

    elseif( refineDir .eq. 2 ) then

      mpd1 = schemeFine%maxPolyDegree+1
      mpd1_square = mpd1 * (schemeFine%maxPolyDegree+1)

      ! loop over all degrees of freedoms on the (semi-)refined element
      do dof = 1, mpd1_square
        iDegX = mod(dof-1,mpd1)+1
        iDegY = (dof-1)/mpd1+1
        fineSqNorm = 2.0_rk /(2.0_rk * iDegY - 1.0_rk)
        ! ... loop over the ansatz functions of the current (non-refined) element
        do iCoarseFunc = 1, schemeCoarse%maxPolyDegree+1
  coarsepos = idegx                                      &
    &      + ( ( icoarsefunc-1)                             &
    &      + (1-1)*(schemecoarse%maxpolydegree+1))*(schemecoarse%maxpolydegree+1)
          ! Weight and add it to the semi refined face representation
          modalRefined(dof,:) = modalRefined(dof,:)                        &
            & + modg_basis%refineBaseCoeff                                 &
            &             %anz_anzShift(iDegY, iCoarseFunc, fineElemShift) &
            & * modalRepFace(coarsePos,:) / ( fineSqNorm )
        end do
      end do

    end if

  end subroutine modg_semiRefineElem_2d

end module atl_modg_2d_multilevel_module