! Copyright (c) 2017 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2017 Daniel Fleischer <daniel.fleischer@student.uni-siegen.de>
! Copyright (c) 2017 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Parts of this file were written by Harald Klimach, Daniel Fleischer and
! Peter Vitt 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.
! **************************************************************************** !
!! This module is used for the projection of QLegendre polynomials from parent
!! to child elements.

module ply_poly_transformation_module
  use env_module,                   only: rk
  use treelmesh_module,             only: treelmesh_type

  implicit none

  private

  ! ************************************************************************ !
  type ply_subsample_type
    !> Is subsampling active
    logical :: isActive = .false.

    !> The current sampling lvl.
    integer :: sampling_lvl

    !> Maximal Level down to which subsampling should be done.
    integer :: caplevel = 20

    !> Minimal subsampling depth:
    integer :: minsub = 0

    !> Maximal subsampling depth:
    integer :: maxsub = 0

    !> Maximum allowed oscillation of the solution.
    !! For adaptive subsampling only.
    real(kind=rk) :: eps_osci

    !> Factor for the reduction of the degrees of freedom in one subsampling
    !! step (per spatial direction).
    real(kind=rk) :: dofReducFactor

    !> Indicator for limitation of total memory consumption
    logical :: adaptiveDofReduction

    !> Absolute upper bound level to refine to.
    integer :: AbsUpperBoundLevel
  end type ply_subsample_type
  ! ************************************************************************ !

  ! ************************************************************************ !
  type ply_array_type
    real(kind=rk),allocatable :: dat(:)
  end type ply_array_type
  ! ************************************************************************ !

  public :: ply_Poly_Transformation
  public :: ply_subsample_type
  public :: ply_array_type

contains

  ! ************************************************************************ !
  !> Projection of polynomial data from parent elements to child elements.
  !! The projection is done by a direct transformation of the modal
  !! coeffiecients to another coordinate system with z=ax+b.
  subroutine ply_Poly_Transformation( subsamp, dofReduction, mesh, meshData,   &
    &                                 varDofs, varComps, ndims, refine_tree,   &
    &                                 new_refine_tree, newMeshData, newVarDofs )
    ! -------------------------------------------------------------------- !
    !> Parameters for the subsampling
    type(ply_subsample_type), intent(in) :: subsamp

    !> Factor for reduction of degrees of freedom.
    real(kind=rk), intent(in) :: dofReduction(:)

    !> The mesh related to meshData.
    type(treelmesh_type), intent(in) :: mesh

    !> The data for subsampling.
    type(ply_array_type), intent(in) :: meshData(:)

    !> The number of degrees of freedom for every variable.
    integer, intent(in) :: varDofs(:)

    !> The number of components for every variable.
    integer, intent(in) :: varComps(:)

    !> Number of dimensions in the polynomial representation.
    integer, intent(in) :: ndims

    !> Logical array that marks elements for refinement
    !! of the previous sampling level.
    logical, intent(in) :: refine_tree(:)

    !> Logical array that marks elements for refinement.
    logical, intent(in) :: new_refine_tree(:)

    !> The subsampled data for new_refine_tree.
    type(ply_array_type), allocatable, intent(out) :: newMeshData(:)

    !> The number of dofs for the subsampled data.
    integer, allocatable, intent(out) :: newVarDofs(:)
    ! -------------------------------------------------------------------- !
    real(kind=rk), allocatable :: workData(:)
    real(kind=rk), allocatable :: newWorkData(:)
    integer :: nVars, nDofs, nComponents, nChildDofs
    integer :: iVar
    ! -------------------------------------------------------------------- !

    nVars = size(varDofs)
    allocate(newVarDofs(nVars))
    allocate(newMeshData(nVars))

    varLoop: do iVar = 1, nVars
      nComponents = varComps(iVar)
      nDofs = vardofs(iVar)

      allocate(workData(size(meshData(iVar)%dat)))
      workData = meshData(iVar)%dat

      if(subsamp%sampling_lvl .eq. subsamp%maxsub) then
        nChildDofs = 1
      else
        nChildDofs = (ceiling(nint(nDofs**(1.0_rk/real(ndims, kind=rk))) &
          &                        * dofReduction(iVar)))**ndims

        if (nChildDofs > nDofs) then
          nChildDofs = nDofs
        elseif (nChildDofs<1) then
          nChildDofs = 1
        end if
      end if

      call ply_subsampleData( mesh             = mesh,             &
        &                     meshData         = workData,         &
        &                     nDofs            = nDofs,            &
        &                     nChildDofs       = nChildDofs,       &
        &                     nComponents      = nComponents,      &
        &                     refine_tree      = refine_tree,      &
        &                     new_refine_tree  = new_refine_tree,  &
        &                     ndims            = ndims,            &
        &                     subsamp          = subsamp,          &
        &                     newMeshData      = newWorkData       )

      allocate(newMeshData(iVar)%dat(size(newWorkData)))
      newMeshData(iVar)%dat = newWorkData

      newVarDofs(iVar) = nChildDofs

      deallocate(workData)

    end do varLoop

  end subroutine ply_Poly_Transformation
  ! ************************************************************************ !

  ! ************************************************************************ !
  subroutine ply_subsampleData( mesh, meshData, nDofs, nChildDofs,         &
    &                           nComponents, refine_tree, new_refine_tree, &
    &                           nDims, subsamp, newMeshData                )
    ! -------------------------------------------------------------------- !
    !> The mesh for the data.
    type(treelmesh_type), intent(in) :: mesh

    !> The data to subsample
    real(kind=rk), intent(in) :: meshData(:)

    !> The number of degrees of freedom.
    integer, intent(in) :: nDofs

    !> The number of degrees of freedom for the child elements.
    integer, intent(in) :: nChildDofs

    !> Number of Components.
    integer, intent(in) :: nComponents

    !> Logical array that marks all elements for refinement for the previous
    !! sampling level.
    logical, intent(in) :: refine_tree(:)

    !> Logical array that marks all elements for refinement for the current
    !! sampling level.
    logical, intent(in) :: new_refine_tree(:)

    !> The number of dimensions in the polynomial representation.
    integer, intent(in) :: nDims

    !> Parameters for subsampling.
    type(ply_subsample_type), intent(in) :: subsamp

    !> The subsampled Data.
    real(kind=rk), allocatable, intent(out) :: newMeshData(:)
    ! -------------------------------------------------------------------- !
    integer :: nChilds, nElems, nElemsToRefine, nElemsNotToRefine
    integer :: nParentElems, childPos, lowElemIndex, upElemIndex
    integer :: iParentElem, iChild, iElem, upChildIndex, lowChildIndex
    integer :: oneDof, noChilds, max_modes
    real(kind=rk), allocatable :: transform_matrix(:,:)
    real(kind=rk), allocatable :: childData(:)
    ! -------------------------------------------------------------------- !
    nChilds = 2**nDims

    max_modes = nint(real(nDofs, kind=rk)**(1.0_rk/real(nDims, kind=rk)))

    ! Get the transformation matrix for the maximal polynomial degree.
    allocate(transform_matrix(max_modes, max_modes))
    call ply_transform_matrix( max_modes = max_modes,       &
      &                        v         = transform_matrix )

    nElems = mesh%nElems
    nElemsToRefine = count(new_refine_tree)
    nElemsNotToRefine = nElems - nElemsToRefine
    nParentElems = size(refine_tree)

    ! Now, we set the correct data for the newMeshData.
    allocate(newMeshData((nElemsToRefine * nChilds * nChildDofs &
      &                  + nElemsNotToRefine) * nComponents))

    newMeshData = 0.0_rk
    upChildIndex = 0
    upElemIndex = 0
    childPos = 0

    if (subsamp%sampling_lvl > 1) then

      elementLoop: do iParentElem=1,nParentElems
        ! Check if the parent cell was already refined...
        if (refine_tree(iParentElem)) then
          ! Parent cell was already refined so it contains data with nDofs
          childLoop: do iChild=1,nChilds
            childPos = childPos + 1
            ! Check if the child elems will be refined...
            if (new_refine_tree(childPos)) then

              ! Child cell will be refined.
              !
              ! Need to project current elem data to new childs
              ! with reduced dofs.
              ! Create lower and upper indices for all data of
              ! iElem in meshData.
              lowElemIndex = upElemIndex + 1
              upElemIndex = (lowElemIndex-1) + nDofs * nComponents

              ! Project these dofs from the coarse element to the
              ! finer elements.
              call ply_projDataToChild(                                  &
                & parentData       = meshData(lowElemIndex:upElemIndex), &
                & nParentDofs      = nDofs,                              &
                & nChildDofs       = nChildDofs,                         &
                & nComponents      = nComponents,                        &
                & nDimensions      = nDims,                              &
                & nChilds          = nChilds,                            &
                & transform_matrix = transform_matrix,                   &
                & childData        = childData                           )

              ! Set the data correctly in newMeshData.
              lowChildIndex = upChildIndex + 1
              upChildIndex = (lowChildIndex-1) + &
                &            nChilds * nChildDofs * nComponents

              newMeshData(lowChildIndex:upChildIndex) = childData
              deallocate(childData)
            else
              ! Child cell won't be refined.
              !
              ! Need projecton from current dofs to 1 dof.
              ! Create lower and upper indices for all data of
              ! iElem in meshData.
              lowElemIndex = upElemIndex + 1
              upElemIndex = (lowElemIndex-1) + nDofs * nComponents

              ! Projection from nDofs to oneDof (integral mean valuea).
              oneDof = 1
              noChilds = 1
              call ply_projDataToChild(                                  &
                & parentData       = meshData(lowElemIndex:upElemIndex), &
                & nParentDofs      = oneDof,                             &
                & nChildDofs       = oneDof,                             &
                & nComponents      = nComponents,                        &
                & nDimensions      = nDims,                              &
                & nChilds          = noChilds,                           &
                & transform_matrix = transform_matrix,                   &
                & childData        = childData                           )

              ! Iterate over all childDofs and set the data corectly
              ! in newMeshData.
              lowChildIndex = upChildIndex + 1
              upChildIndex = (lowChildIndex-1) + nComponents

              newMeshData(lowChildIndex:upChildIndex) = childData
              deallocate(childData)
            end if
          end do childLoop
        else
          ! Parent cell wasn't refined so it contains data with only one dof.
          ! Simple copying.
          allocate(childData(nComponents))

          lowElemIndex = upElemIndex + 1
          upElemIndex = (lowElemIndex-1) + nComponents

          childData = meshData(lowElemIndex:upElemIndex)

          lowChildIndex = upChildIndex + 1
          upChildIndex = (lowChildIndex-1) + nComponents

          newMeshData(lowChildIndex:upChildIndex) = childData
          childpos = childpos + 1
          deallocate(childData)
        end if
      end do elementLoop
    else
      elemLoop: do iElem=1,nElems
        if (new_refine_tree(iElem)) then
          ! Create lower and upper indices for all data of iElem in meshData.
          lowElemIndex = upElemIndex + 1
          upElemIndex = (lowElemIndex-1) + nDofs * nComponents

          ! Project these dofs from the coarse element to the
          ! finer elements.
          call ply_projDataToChild(                                  &
            & parentData       = meshData(lowElemIndex:upElemIndex), &
            & nParentDofs      = nDofs,                              &
            & nChildDofs       = nChildDofs,                         &
            & nComponents      = nComponents,                        &
            & nDimensions      = nDims,                              &
            & nChilds          = nChilds,                            &
            & transform_matrix = transform_matrix,                   &
            & childData        = childData                           )

          ! Iterate over all childDofs and set the data corectly in newMeshData
          lowChildIndex = upChildIndex + 1
          upChildIndex = (lowChildIndex-1) + nChilds * nChildDofs * nComponents

          newMeshData(lowChildIndex:upChildIndex) = childData
          deallocate(childData)
        else
          ! Create lower and upper indices for all data of iElem in meshData.
          lowElemIndex = upElemIndex + 1
          upElemIndex = (lowElemIndex-1) + nDofs * nComponents

          ! Projection from nDofs to oneDof (integral mean value).
          oneDof = 1
          noChilds = 1
          call ply_projDataToChild(                                  &
            & parentData       = meshData(lowElemIndex:upElemIndex), &
            & nParentDofs      = oneDof,                             &
            & nChildDofs       = oneDof,                             &
            & nComponents      = nComponents,                        &
            & nDimensions      = nDims,                              &
            & nChilds          = noChilds,                           &
            & transform_matrix = transform_matrix,                   &
            & childData        = childData                           )

          ! Iterate over all childDofs and set the data corectly in newMeshData
          lowChildIndex = upChildIndex + 1
          upChildIndex = (lowChildIndex-1) + nComponents

          newMeshData(lowChildIndex:upChildIndex) = childData
          deallocate(childData)
        end if
      end do elemLoop
    end if

    deallocate(transform_matrix)

  end subroutine ply_subsampleData
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Subroutine to project element data from a parent cell to its children.
  subroutine ply_projDataToChild( parentData, nParentDofs, nChildDofs, &
    &                             nComponents, nDimensions, nChilds,   &
    &                             transform_matrix, childData          )
    ! -------------------------------------------------------------------- !
    !> The polynomial data for a single parent element.
    real(kind=rk), intent(in) :: parentData(:)

    !> The number of dofs of the parent element.
    integer, intent(in) :: nParentDofs

    !> The total number of dofs for the child cells.
    integer, intent(in) :: nChildDofs

    !> The number of componentns for the given variable.
    integer, intent(in) :: nComponents

    !> The number of dimensions.
    integer, intent(in) :: nDimensions

    !> The number of child elements.
    integer, intent(in) :: nChilds

    !> The transformation matrix for the linear coordinate transformation.
    real(kind=rk), intent(in) :: transform_matrix(:,:)

    !> The new data representation for all child cell of the parent cell.
    real(kind=rk), allocatable, intent(out) :: childData(:)
    ! -------------------------------------------------------------------- !
    integer :: iDimension, iIndep, iComponent, iMode, jMode
    integer :: nSubElems, nIndeps, lMode, kMode, iChildElem_prev
    integer :: nChildElems_prev, nChildElems_cur
    integer :: child_dofPos, childElem
    integer :: parent_modes, child_modes
    integer :: lower_bound, upper_bound, stride
    real(kind=rk), allocatable :: temp_data(:)
    real(kind=rk), allocatable :: childData_prev(:)
    ! -------------------------------------------------------------------- !
    parent_modes = nint(real(nParentDofs,kind=rk)        &
      &                 **(1/real(nDimensions,kind=rk)))

    child_modes = nint(real(nChildDofs,kind=rk)         &
      &                **(1/real(nDimensions,kind=rk)))

    allocate(temp_Data(parent_modes))

    do iDimension = 1,nDimensions

      if (nChilds > 1) then
        nSubElems = 2
        nChildElems_cur = 2**(iDimension)
        nChildElems_prev = 2**(iDimension-1)
      else
        nChildElems_cur = nChilds
        nChildElems_prev = nChilds
        nSubElems = nChilds
      end if

      stride = child_modes**(iDimension-1) * nComponents

      if (iDimension .eq. 1) then

        nIndeps = parent_modes**(nDimensions-1)

        ! Allocate memory for two childs in x-direction
        if (iDimension .eq. nDimensions) then
          allocate(childData(nChildElems_cur * child_modes**nDimensions &
            &                * nComponents))
        else
          allocate(childData(nChildElems_cur                   &
            &                  * parent_modes**(nDimensions-1) &
            &                  * child_modes**(nDimensions-2)  &
            &                  * nComponents                   ))
        end if
        childData = 0.0_rk

        call splitx()

      elseif (iDimension .eq. 2) then

        allocate(childData_prev(size(childData)))
        childData_prev = childData
        deallocate(childData)

        ! Allocate memory for the four childs in y-direction
        if (iDimension .eq. nDimensions) then
          allocate(childData(nChildElems_cur * child_modes**nDimensions &
            &                * nComponents) )
        else
          allocate(childData(nChildElems_cur                                &
            &                  * parent_modes**(nDimensions-2)              &
            &                  * child_modes**(nDimensions-1) * nComponents))
        end if
        childData = 0.0_rk

        call splity()

      else

        deallocate(childData_prev)
        allocate(childData_prev(size(childData)))
        childData_prev = childData

        deallocate(childData)

        ! Allocate memory for the eight childs in z-direction
        allocate(childData(nChildElems_cur * child_modes**3 &
          &                 * nComponents))
        childData = 0.0_rk

        call splitz()

      end if
    end do

    deallocate(childData_prev)
    deallocate(temp_Data)

  contains

    ! NOTE: these internal subroutines were introduced to help the compiler
    !       and avoid internal compiler errors when compiling for AVX
    !       instructions.

    ! Perform the projection in x direction onto two child elements
    subroutine splitx()
      ! iSubElem = 1
      do iComponent = 1, nComponents
        do iIndep = 1, nIndeps

          lower_bound = iComponent + (iIndep - 1) * parent_modes * nComponents
          upper_bound = iIndep * parent_modes * nComponents

          temp_data = parentData(lower_bound:upper_bound:stride)

          do iMode = 1, child_modes
              child_dofPos = iComponent                                   &
                &              + (iMode-1) * nComponents                  &
                &              + (iIndep - 1) * child_modes * nComponents

            do jMode = iMode, parent_modes

              childData(child_dofpos) = childData(child_dofpos)         &
                &                         + temp_Data(jMode)            &
                &                         * transform_matrix(jMode,iMode)

            end do
          end do
        end do
      end do

      if (nSubElems > 1) then
      ! iSubElem = 2
        do iComponent = 1, nComponents
          do iIndep = 1, nIndeps

            lower_bound = iComponent &
              &             + (iIndep - 1) * parent_modes * nComponents
            upper_bound = iIndep * parent_modes * nComponents

            temp_data = parentData(lower_bound:upper_bound:stride)

            do iMode = 1, child_modes
                child_dofPos = iComponent                                   &
                  &              + (iMode-1) * nComponents                  &
                  &              + (iIndep - 1) * child_modes * nComponents &
                  &              + nComponents                              &
                  &                * child_modes**(nDimensions-2)           &
                  &                * parent_modes**(nDimensions-1)

              do jMode = iMode, parent_modes

                childData(child_dofpos) = childData(child_dofpos)           &
                  &                         + temp_Data(jMode)              &
                  &                           * transform_matrix(iMode,jMode)

              end do
            end do
          end do
        end do
      end if
    end subroutine splitx

    ! Perform the projection in y direction onto the two child elements
    ! for each of the  two parents
    ! (the projected elements split in x direction)
    subroutine splity()
      do iChildElem_prev = 1, nChildElems_prev
        ! iSubElem = 1
        childElem = iChildElem_prev
        do iComponent = 1, nComponents
          do kMode = 1, parent_modes
            do lMode = 1, child_modes

              lower_bound = iComponent                                 &
                & + (lMode-1) * nComponents                            &
                & + (kMode-1) * child_modes*parent_modes * nComponents &
                & + (iChildElem_prev-1) * nComponents                  &
                &   * child_modes**(nDimensions-2)                     &
                &   * parent_modes**(nDimensions-1)

              upper_bound = kMode * parent_modes * child_modes * nComponents &
                &             + (iChildElem_prev-1) * nComponents            &
                &               * child_modes**(nDimensions-2)               &
                &               * parent_modes**(nDimensions-1)

              temp_data = childData_prev(lower_bound:upper_bound:stride)

              do iMode = 1, child_modes
                  child_dofPos = iComponent                                &
                    &              + (iMode-1) * child_modes * nComponents &
                    &              + (lMode-1) * nComponents               &
                    &              + (kMode-1) * nComponents               &
                    &                * child_modes**2                      &
                    &              + (childElem-1) * nComponents           &
                    &                * child_modes**(nDimensions-1)        &
                    &                * parent_modes**(nDimensions-2)

                do jMode = iMode, parent_modes

                  childData(child_dofpos) = childData(child_dofpos)         &
                    & + temp_Data(jMode) * transform_matrix(jMode,iMode)

                end do
              end do
            end do
          end do
        end do

        if (nSubElems > 1) then
          ! iSubElem = 2
          childElem = iChildElem_prev + 2
          do iComponent = 1, nComponents
            do kMode = 1, parent_modes
              do lMode = 1, child_modes

                lower_bound = iComponent                                 &
                  & + (lMode-1) * nComponents                            &
                  & + (kMode-1) * child_modes*parent_modes * nComponents &
                  & + (iChildElem_prev-1) * nComponents                  &
                  &   * child_modes**(nDimensions-2)                     &
                  &   * parent_modes**(nDimensions-1)

                upper_bound = kMode * parent_modes * child_modes &
                  &   * nComponents                              &
                  & + (iChildElem_prev-1) * nComponents          &
                  &   * child_modes**(nDimensions-2)             &
                  &   * parent_modes**(nDimensions-1)

                temp_data = childData_prev(lower_bound:upper_bound:stride)

                do iMode = 1, child_modes
                    child_dofPos = iComponent                      &
                      & + (iMode-1) * child_modes * nComponents    &
                      & + (lMode-1) * nComponents                  &
                      & + (kMode-1) * nComponents * child_modes**2 &
                      & + (childElem-1) * nComponents              &
                      &   * child_modes**(nDimensions-1)           &
                      &   * parent_modes**(nDimensions-2)

                  do jMode = iMode, parent_modes

                    childData(child_dofpos) = childData(child_dofpos) &
                      & + temp_Data(jMode) * transform_matrix(iMode,jMode)

                  end do
                end do
              end do
            end do
          end do
        end if
      end do
    end subroutine splity

    ! Perform the projection in z direction onto the two child elements
    ! for each of the four parents
    ! (the projected elements split in x and y direction)
    subroutine splitz()
      do iChildElem_prev = 1, nChildElems_prev
        ! iSubElem = 1
        childElem = iChildElem_prev
        do iComponent = 1, nComponents
          iIndep = 0
          do kMode = 1,child_modes
            do lMode = 1,child_modes
              iIndep = iIndep + 1

              lower_bound = iComponent                &
                & + (iIndep - 1) * nComponents        &
                & + (iChildElem_prev-1) * nComponents &
                &   * child_modes**(nDimensions-1)    &
                &   * parent_modes**(nDimensions-2)

              upper_bound = parent_modes * child_modes**(nDimensions-1) &
                &   * nComponents                                       &
                & - (nComponents - iComponent)                          &
                & + (iChildElem_prev-1) * nComponents                   &
                &   * child_modes**(nDimensions-1)                      &
                &   * parent_modes**(nDimensions-2)

              temp_data(:) = childData_prev &
                &            (lower_bound:upper_bound:stride)

              do iMode = 1, child_modes
                  child_dofPos = iComponent                      &
                    & + (iMode-1) * child_modes**2 * nComponents &
                    & + (lMode - 1) * nComponents                &
                    & + (kMode - 1) * child_modes * nComponents  &
                    & + (childElem - 1) * nComponents            &
                    &   * child_modes**nDimensions

                do jMode = iMode, parent_modes

                  childData(child_dofpos) = childData(child_dofpos) &
                    & + temp_Data(jMode) * transform_matrix(jMode,iMode)

                end do
              end do
            end do
          end do
        end do

        if (nSubElems > 1) then
          ! iSubElem = 2
          childElem = iChildElem_prev + 4
          do iComponent = 1, nComponents
            iIndep = 0
            do kMode = 1,child_modes
              do lMode = 1,child_modes
                iIndep = iIndep + 1

                lower_bound = iComponent                &
                  & + (iIndep - 1) * nComponents        &
                  & + (iChildElem_prev-1) * nComponents &
                  &   * child_modes**(nDimensions-1)    &
                  &   * parent_modes**(nDimensions-2)

                upper_bound = parent_modes * child_modes**(nDimensions-1) &
                  &   * nComponents                                       &
                  & - (nComponents - iComponent)                          &
                  & + (iChildElem_prev-1) * nComponents                   &
                  &   * child_modes**(nDimensions-1)                      &
                  &   * parent_modes**(nDimensions-2)

                temp_data(:) = childData_prev(lower_bound:upper_bound:stride)

                do iMode = 1, child_modes
                    child_dofPos = iComponent                      &
                      & + (iMode-1) * child_modes**2 * nComponents &
                      & + (lMode - 1) * nComponents                &
                      & + (kMode - 1) * child_modes * nComponents  &
                      & + (childElem - 1) * nComponents            &
                      &   * child_modes**nDimensions

                  do jMode = iMode, parent_modes

                    childData(child_dofpos) = childData(child_dofpos) &
                      & + temp_Data(jMode) * transform_matrix(iMode,jMode)

                  end do
                end do
              end do
            end do
          end do
        end if
      end do
    end subroutine splitz

  end subroutine ply_projDataToChild
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Compute the transformation matrix for a projection to the left and right
  !! half-interval of Legendre polynomials for the given maximal number of
  !! modes.
  !!
  !! Note: The transformation matrices to each subinterval are triangular, and
  !!       the diagonal entries are the same. To save memory both matrices are
  !!       stored in a single 2 dimensional array of size
  !!       (max_modes, max_modes).
  !!
  !! This matrix only needs to be computed once for a sufficiently high order,
  !! as submatices out of it can by used to perform the transformation for
  !! any lower polynomial degree.
  subroutine ply_transform_matrix(max_modes, v)
    ! -------------------------------------------------------------------- !
    !> The maximal number of modes to compute the transformation for.
    !!
    !! The resulting matrix v will be max_modes x max_modes large and can
    !! be used for the transformation of all polynomials with up to this
    !! many modes.
    integer, intent(in) :: max_modes

    !> The transformation matrix.
    !!
    !! Upper triangular matrix is created for shifting and lower triangular
    !! for (-1) * shifting.
    !! For the right interval we interpret the first index as row index
    !! and the second as column. For the left interval this is reverted and
    !! we interpret the first index as columns of the matrix.
    real(kind=rk), allocatable, intent(out) :: v(:,:)
    ! -------------------------------------------------------------------- !
    integer :: m, orig
    real(kind=rk) :: shifting, scaling
    ! -------------------------------------------------------------------- !

    ! transformation matrix looks like this:
    ! [1.0  --  --     shift=0.5   ]
    ! | |  0.5  --                 |
    ! | |   |  0.25             ...|
    ! |             0.125          |
    ! | shift=-0.5        0.0625   |
    ! [    :                    ...]

    allocate(v(max_modes,max_modes))
    v = 0.0_rk

    scaling = 0.5_rk
    shifting = 0.5_rk

    ! Set the first two entries of v manually.
    v(1,1) = 1.0_rk
    if (max_modes > 1) then
      v(1,2) = shifting
      v(2,2) = scaling

      if (max_modes > 2) then
        do orig = 3,max_modes
          v(1,orig) = ply_beta(orig-1) * v(1,orig-2)                    &
            &       + ply_alpha(orig-1) * shifting * v(1,orig-1)        &
            &       - scaling * ply_alpha_beta(2,orig-1) * v(2, orig-1)
          do m = 2,orig
            if (m < max_modes) then
              v(m,orig) = ply_beta(orig-1) * v(m,orig-2)                     &
                &       + ply_alpha(orig-1) * shifting * v(m,orig-1)         &
                &       - scaling * ply_alpha_beta(m+1,orig-1)               &
                &                 * v(m+1, orig-1)                           &
                &       + scaling * ply_alpha_frac(m-1,orig-1) * v(m-1,orig-1)
            else
              ! Need to skip one summand for v(max_modes,max_modes).
              v(m,orig) = scaling * ply_alpha_frac(m-1,orig-1) * v(m-1,orig-1)
            end if
          end do
        end do
      end if

      ! Due to the symmetry of the problem (the left subinterval has just
      ! the shifting with a changed sign), we can fill the other half of
      ! the matrix by copying the already computed values accordingly with
      ! a change in sign, as needed (alternatingly).
      do m = 1 , max_modes
        do orig = 1, m-1
          if (mod((m+orig),2) /= 0) then
            v(m, orig) = - v(orig, m)
          else
            v(m, orig) = v(orig, m)
          end if
        end do
      end do
    end if

  end subroutine ply_transform_matrix
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Coefficients from the recursive formulation of legendre polynomials.
  !! L_n = alpha * x * L_n-1 + beta * L_n-2
  function ply_alpha(mode) result(alpha)
    ! -------------------------------------------------------------------- !
    !> The current mode in the polynomial representation.
    integer, intent(in) :: mode

    !> Alpha coefficient from the recursive formulation of legendre
    !! polynomials.
    real(kind=rk) :: alpha
    ! -------------------------------------------------------------------- !
    if (mode > 0) then
      alpha = real((2 * mode - 1), kind=rk) / real(mode, kind=rk)
    else
      alpha = 0.0
    end if
  end function ply_alpha
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Coefficients from the recursive formulation of legendre polynomials.
  !! L_n = alpha * x * L_n-1 + beta * L_n-2
  !!
  function ply_beta(mode) result(beta)
    ! -------------------------------------------------------------------- !
    !> The current mode in the polynomial representation.
    integer, intent(in) :: mode

    !> Beta coefficient from the recursive formulation of legendre
    !! polynomials.
    real(kind=rk) :: beta
    ! -------------------------------------------------------------------- !
    if (mode > 0) then
      beta = real((1 - mode), kind=rk) / real(mode, kind=rk)
    else
      beta = 0.0
    end if
  end function ply_beta
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Quotient of two alpha values.
  function ply_alpha_frac(denominator, numerator) result(alpha_frac)
    ! -------------------------------------------------------------------- !
    !> Numerator
    integer, intent(in) :: numerator

    !> Denominator
    integer, intent(in) :: denominator

    !> The quotient of two alpha values.
    real(kind=rk) :: alpha_frac
    ! -------------------------------------------------------------------- !
    if ( denominator > 0 .and. numerator > 0 ) then
      if ( denominator == numerator ) then
        alpha_frac = 1.0_rk
      else
        alpha_frac = real((2*numerator-1)*denominator, kind=rk)  &
          &            / real((2*denominator-1)*numerator, kind=rk)
      end if
    else
      alpha_frac = 0.0_rk
    end if
  end function ply_alpha_frac
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Prodcut of alpha(numerator) * beta(denominator) / alpha(denominator)
  function ply_alpha_beta(denominator, numerator) result(alpha_beta)
    ! -------------------------------------------------------------------- !
    !> Numerator
    integer, intent(in) :: numerator

    !> Denominator
    integer, intent(in) :: denominator

    !> The product of alpha(n) * beta(d) / alpha(d)
    real(kind=rk) :: alpha_beta
    ! -------------------------------------------------------------------- !
    if (numerator > 0) then
      if ( denominator == numerator ) then
        alpha_beta = real((1-denominator), kind=rk)/ real(denominator, kind=rk)
      else
        alpha_beta = real((2*numerator-1)*(1-denominator), kind=rk) &
        &            / real(numerator*(2*denominator-1), kind=rk)
      end if
    else
      alpha_beta = 0.0_rk
    end if
  end function ply_alpha_beta
  ! ************************************************************************ !

end module ply_poly_transformation_module
