ply_poly_project_module.f90 Source File


This file depends on

sourcefile~~ply_poly_project_module.f90~~EfferentGraph sourcefile~ply_poly_project_module.f90 ply_poly_project_module.f90 sourcefile~ply_dof_module.f90 ply_dof_module.f90 sourcefile~ply_poly_project_module.f90->sourcefile~ply_dof_module.f90 sourcefile~ply_nodes_header_module.f90 ply_nodes_header_module.f90 sourcefile~ply_poly_project_module.f90->sourcefile~ply_nodes_header_module.f90 sourcefile~ply_nodes_module.f90 ply_nodes_module.f90 sourcefile~ply_poly_project_module.f90->sourcefile~ply_nodes_module.f90 sourcefile~ply_l2p_module.f90 ply_l2p_module.f90 sourcefile~ply_poly_project_module.f90->sourcefile~ply_l2p_module.f90 sourcefile~ply_dynarray_project_module.f90 ply_dynArray_project_module.f90 sourcefile~ply_poly_project_module.f90->sourcefile~ply_dynarray_project_module.f90 sourcefile~ply_legfpt_module.f90 ply_legFpt_module.f90 sourcefile~ply_poly_project_module.f90->sourcefile~ply_legfpt_module.f90 sourcefile~ply_legfpt_3d_module.f90 ply_legFpt_3D_module.f90 sourcefile~ply_poly_project_module.f90->sourcefile~ply_legfpt_3d_module.f90 sourcefile~ply_prj_header_module.f90 ply_prj_header_module.f90 sourcefile~ply_poly_project_module.f90->sourcefile~ply_prj_header_module.f90 sourcefile~ply_legfpt_2d_module.f90 ply_legFpt_2D_module.f90 sourcefile~ply_poly_project_module.f90->sourcefile~ply_legfpt_2d_module.f90 sourcefile~ply_fxt_module.f90 ply_fxt_module.f90 sourcefile~ply_poly_project_module.f90->sourcefile~ply_fxt_module.f90 sourcefile~ply_nodes_module.f90->sourcefile~ply_nodes_header_module.f90 sourcefile~ply_nodeset_module.f90 ply_nodeset_module.f90 sourcefile~ply_nodes_module.f90->sourcefile~ply_nodeset_module.f90 sourcefile~ply_l2p_module.f90->sourcefile~ply_nodeset_module.f90 sourcefile~ply_space_integration_module.f90 ply_space_integration_module.f90 sourcefile~ply_l2p_module.f90->sourcefile~ply_space_integration_module.f90 sourcefile~ply_modg_basis_module.f90 ply_modg_basis_module.f90 sourcefile~ply_l2p_module.f90->sourcefile~ply_modg_basis_module.f90 sourcefile~ply_l2p_header_module.f90 ply_l2p_header_module.f90 sourcefile~ply_l2p_module.f90->sourcefile~ply_l2p_header_module.f90 sourcefile~ply_lagrange_module.f90 ply_lagrange_module.f90 sourcefile~ply_l2p_module.f90->sourcefile~ply_lagrange_module.f90 sourcefile~ply_dynarray_project_module.f90->sourcefile~ply_prj_header_module.f90 sourcefile~ply_polybaseexc_module.f90 ply_polyBaseExc_module.f90 sourcefile~ply_legfpt_module.f90->sourcefile~ply_polybaseexc_module.f90 sourcefile~ply_fpt_header_module.f90 ply_fpt_header_module.f90 sourcefile~ply_legfpt_module.f90->sourcefile~ply_fpt_header_module.f90 sourcefile~ply_legfpt_3d_module.f90->sourcefile~ply_legfpt_module.f90 sourcefile~ply_prj_header_module.f90->sourcefile~ply_fpt_header_module.f90 sourcefile~ply_fxt_header_module.f90 ply_fxt_header_module.f90 sourcefile~ply_prj_header_module.f90->sourcefile~ply_fxt_header_module.f90 sourcefile~ply_prj_header_module.f90->sourcefile~ply_l2p_header_module.f90 sourcefile~ply_legfpt_2d_module.f90->sourcefile~ply_legfpt_module.f90 sourcefile~ply_fxt_module.f90->sourcefile~ply_fxt_header_module.f90 sourcefile~ply_polybaseexc_module.f90->sourcefile~ply_fpt_header_module.f90 sourcefile~ply_fpt_header_module.f90->sourcefile~ply_nodes_header_module.f90 sourcefile~ply_fxt_header_module.f90->sourcefile~ply_nodes_header_module.f90 sourcefile~ply_modg_basis_module.f90->sourcefile~ply_dof_module.f90 sourcefile~ply_modg_basis_module.f90->sourcefile~ply_space_integration_module.f90 sourcefile~ply_l2p_header_module.f90->sourcefile~ply_nodes_header_module.f90 sourcefile~ply_lagrange_module.f90->sourcefile~ply_nodeset_module.f90

Files dependent on this one

sourcefile~~ply_poly_project_module.f90~~AfferentGraph sourcefile~ply_poly_project_module.f90 ply_poly_project_module.f90 sourcefile~ply_oversample_module.f90 ply_oversample_module.f90 sourcefile~ply_oversample_module.f90->sourcefile~ply_poly_project_module.f90 sourcefile~sdr_proto2treelm_module.f90 sdr_proto2treelm_module.f90 sourcefile~sdr_proto2treelm_module.f90->sourcefile~ply_poly_project_module.f90 sourcefile~sdr_proto2treelm_module.f90->sourcefile~ply_oversample_module.f90 sourcefile~seeder.f90 seeder.f90 sourcefile~seeder.f90->sourcefile~sdr_proto2treelm_module.f90

Contents


Source Code

! Copyright (c) 2013-2014 Jens Zudrop <j.zudrop@grs-sim.de>
! Copyright (c) 2013-2014, 2016-2017, 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2013-2015 Nikhil Anand <nikhil.anand@uni-siegen.de>
! Copyright (c) 2013-2018, 2020 Harald Klimach <harald.klimach@uni-siegen.de.de>
! Copyright (c) 2013-2014 Verena Krupp <v.krupp@grs-sim.de>
! Copyright (c) 2015 Kay Langhammer <kay.langhammer@student.uni-siegen.de>
! Copyright (c) 2016-2017, 2019 Neda Ebrahimi Pour <neda.epour@uni-siegen.de>
! Copyright (c) 2016 Tobias Girresser <tobias.girresser@student.uni-siegen.de>
! Copyright (c) 2016 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
!
! Parts of this file were written by Jens Zudrop for German Research School
! for Simulation Sciences GmbH.
!
! Parts of this file were written by Harald Klimach, Verena Krupp, Peter Vitt,
! Tobias Girresser, Nikhil Anand, Kay Langhammer, Kannan Masilamani and Neda
! Ebrahimi Pour 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.
! **************************************************************************** !
! 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.
!> The polynomial projection describes the change from modal to nodal
!! representation and vice-versa.
!!
!! The method to be used for this transformation is configured in a table that
!! contains the `kind` of transformation and then further options to the
!! transformation.
!! There are three kinds of transformation available:
!!
!! * 'fpt' Fast polynomial transformation, only available if linked against
!!   FFTW.
!! * 'l2p' Standard projection, evaluate each polynomial for each point.
!! * 'fxt' Fast transformation based on multipole method, enabled by the
!!   (included) FXTPACK library.
!!
!! By default the nodal representation uses as many nodes as there are modes
!! in the modal representation. This can be changed with the `factor` setting.
!! There will `factor` times number of modes nodes used in the transformation.
!! This oversampling allows for dealiasing of the nodal representation.
!! A different set of points will be used depending on the kind of
!! transformation to be used:
!!
!! * FPT uses Chebyshev integration nodes
!! * L2P may make use of either Chebyshev or Gauss-Legendre nodes
!! * FXT Gauss-Legendre integration nodes
!!
!! If the Chebyshev nodes are used for the nodal representation, the interval
!! boundaries can be included by setting `lobattoPoints=true`. In this case
!! there will also be at least one point more be used in the nodal
!! representation than there are modes in the Legendre series.
!!
!! If FPT is configured but the executable is not linked against the FFTW, there
!! will be a warning, and the simulation uses the L2P method instead.
!! The L2P method is also the default when no kind is provided at all.
!!
!! For further options of the individual projection kinds, see their respective
!! descriptions:
!!
!! * L2P: [[ply_l2p_header_module]]
!! * FPT: [[ply_fpt_header_module]]
!! * FXT: [[ply_fxt_header_module]]
!!
!! A minimalistic configuration for the projection method looks as follows:
!!
!!```lua
!!  projection = {
!!    kind = 'l2p',
!!    factor = 1.5
!!  }
!!```
!!
!! The name of the table (here `projection`) is arbitrary and defined by the
!! application.
!! An application may load multiple projection definitions with differing names.
!!
module ply_poly_project_module
  use env_module,                only: rk, labelLen

  use tem_aux_module,            only: tem_abort
  use tem_logging_module,        only: logUnit
  use tem_tools_module,          only: tem_horizontalSpacer
  use ply_dof_module,            only: Q_space, &
   &                                   P_space
  use ply_prj_header_module
  use ply_dynArray_project_module, only: dyn_ProjectionArray_type, &
    &                                    ply_fill_dynProjectArray, &
    &                                    ply_prj_init_type

  use ply_LegFpt_module,           only: ply_legFpt_type, &
    &                                    ply_init_legFpt, &
    &                                    assignment(=)
  use ply_legFpt_2D_module,        only: ply_pntToLeg_2D, &
    &                                    ply_legToPnt_2D
  use ply_legFpt_3D_module,        only: ply_pntToLeg_3D, &
    &                                    ply_legToPnt_3D
  use ply_l2p_module,              only: ply_l2p_type,     &
    &                                    ply_init_l2p,     &
    &                                    assignment(=),    &
    &                                    ply_l2p_trafo_3d, &
    &                                    ply_l2p_trafo_2d, &
    &                                    ply_l2p_trafo_1d
  use ply_nodes_header_module,     only: ply_nodes_header_type, &
    &                                    assignment(=)
  use ply_nodes_module,            only: ply_nodes_create, &
    &                                    ply_facenodes_type
  use ply_fxt_module,              only: ply_fxt_type,   &
    &                                    ply_init_fxt,   &
    &                                    ply_fxt_m2n_1D, &
    &                                    ply_fxt_m2n_3D, &
    &                                    ply_fxt_m2n_2D, &
    &                                    ply_fxt_n2m_1D, &
    &                                    ply_fxt_n2m_3D, &
    &                                    ply_fxt_n2m_2D

  implicit none

  private

  !> Additional data, required for the projection.
  type ply_prj_body_type
    !> The fast polynomial transformation which will be used in case
    !! of nonlinear equations. It is used if fpt is choses as projection
    !! method in the lua file
    type(ply_legFpt_type)  :: fpt
    !> The Legendre Polynomial type for the Fast Orthogonal Function
    !! Transform via fxtpack. It is used if 'fxt' is chosen as projection
    !! method in the lua file
    type(ply_fxt_type) :: fxt
    !> Projection method which cam be used for transfoamtion from modal to
    !! nodal space and vice versa. It is used if 'l2p' is chosen as projection
    !! method in the lua file
    type(ply_l2p_type) :: l2p
    !> Volume quadrature points in the reference element
    real(kind=rk), allocatable :: nodes(:,:)
    !> Facial quadrature nodes (reference element) for all three spatial
    !! direction and left and right face.
    !! These points are necessary to transfer boundary conditions given
    !! in physical space to modal space by projection (l2p or fpt)
    type(ply_faceNodes_type), allocatable :: faces(:,:)
    !> quadrature points including oversampling factor
    integer                    :: nQuadPoints
    !> degree of freedom of the scheme depending on maxPolyDegree
    integer                    :: ndofs
    ! the oversamp_dofs are the degrees of freedom for this 'oversampling
    ! degree and equal to (s*(m+1))**d
    integer                    :: oversamp_dofs
    ! minimal number of dofs, used to enable the flexibilty to
    ! set the oversampling factor < 1
    integer                    :: min_dofs

  end type ply_prj_body_type


  !> Projection definition.
  type ply_poly_project_type
    !> Polynomial basis type.
    !!
    !! 3D Monomials have the form x^i * y^j * z^k
    !! - Q_space: quadratic polynomial space (i,j,k) <= maxPolyDegree
    !! - P_space: polynomial space i+j+k <= maxPolyDegree
    integer                       :: basisType

    !> Kind of projection. Currently available:
    !! - 'l2p', L2-Projection
    !! - 'fpt', Fast Polynomial Transformation. Requires the FFTW.
    !! - 'fxt', Fast Polynomial Transformation. uses FXTPACK
    character(len=labelLen)       :: kind

    !> The maximal polynomial degree per spatial direction.
    integer                       :: maxPolyDegree

    !> Using oversampling, the modal space need to be extended according
    ! to the oversampling factor, thus the oversampling degree is (s*m+1)-1
    integer                       :: oversamp_degree

    ! minimal number of dofs, used to enable the flexibilty to
    ! set the oversampling factor < 1
    integer                       :: min_degree

    !> quadrature points including oversampling factor
    ! per spatial direction
    integer                       :: nQuadPointsPerDir

    !> Logical to indicate whether Chebyshev-Lobatto points or simple
    !! Chebyshev points are used
    logical                       :: lobattoPoints = .false.

    !> projection header consits of general information like which kind
    !! of projection is used

    !> In the body datatype, there is for each dimension the main data
    !! for the projection method stored
    type(ply_prj_body_type)   :: body_1d
    type(ply_prj_body_type)   :: body_2d
    type(ply_prj_body_type)   :: body_3d

  end type ply_poly_project_type


  interface assignment(=)
    module procedure Copy_poly_project
    module procedure Copy_poly_project_body
  end interface

  interface ply_poly_project_m2n
    module procedure ply_poly_project_m2n_multiVar
  end interface

  interface ply_poly_project_n2m
    module procedure ply_poly_project_n2m_multiVar
  end interface

  public :: assignment(=)
  public :: ply_fill_project_list
  public :: ply_poly_project_fillbody
  public :: ply_poly_project_m2n
  public :: ply_poly_project_n2m
  public :: ply_poly_project_type
  public :: ply_faceNodes_type
  public :: ply_get_quadpoints_faces
  public :: ply_get_quadpoints_faces_2d
  public :: ply_get_quadpoints_faces_1d
  public :: ply_prj_body_type


contains


  ! ------------------------------------------------------------------------ !
  subroutine Copy_poly_project(left,right)
    ! -------------------------------------------------------------------- !
    !> fpt to copy to
    type(ply_poly_project_type), intent(out) :: left
    !> fpt to copy from
    type(ply_poly_project_type), intent(in) :: right
    ! -------------------------------------------------------------------- !
    ! -------------------------------------------------------------------- !

    left%body_1d = right%body_1d
    left%body_2d = right%body_2d
    left%body_3d = right%body_3d

    left%kind = right%kind
    left%maxPolyDegree = right%maxPolyDegree
    left%basisType = right%basisType
    left%oversamp_degree = right%oversamp_degree
    left%min_degree = right%min_degree
    left%nquadpointsPerDir = right%nquadpointsPerDir
    left%lobattoPoints = right%lobattoPoints

  end subroutine copy_poly_project
  ! ------------------------------------------------------------------------ !


  ! ------------------------------------------------------------------------ !
  subroutine Copy_poly_project_body(left,right)
    ! -------------------------------------------------------------------- !
    ! fpt to copy to
    type(ply_prj_body_type), intent(out) :: left
    ! fpt to copy from
    type(ply_prj_body_type), intent(in) :: right
    ! -------------------------------------------------------------------- !
    ! -------------------------------------------------------------------- !

    left%fpt = right%fpt
    left%l2p = right%l2p
    left%fxt = right%fxt
    left%nodes = right%nodes
    left%faces = right%faces
    left%nquadpoints = right%nquadpoints
    left%ndofs = right%ndofs
    left%oversamp_dofs = right%oversamp_dofs

  end subroutine copy_poly_project_body
  ! ------------------------------------------------------------------------ !


  ! ------------------------------------------------------------------------ !
  !> Fill ups the bodys accroding to the DA.
  subroutine ply_fill_project_list( proj_list, dyn_projectArray, scheme_dim )
    ! -------------------------------------------------------------------- !
    type(ply_poly_project_type), intent(out), allocatable :: proj_list(:)
    type(dyn_ProjectionArray_type), intent(in) :: dyn_projectArray
    integer, intent(in) :: scheme_dim
    ! -------------------------------------------------------------------- !
    integer :: ipos
    ! -------------------------------------------------------------------- !
    call tem_horizontalSpacer(fUnit=logUnit(2))
    write(logUnit(2),*) 'Loading list of projection methods ... '

    ! allocate the poly_proj_list, the maximum value of all positions will give
    ! the number of elements in the DA
    allocate (proj_list(dyn_projectArray%nVals) )
    write(logUnit(5),*) 'the number of elements in proj_pos is =', &
      &                 dyn_projectArray%nVals
    do ipos=1, dyn_projectArray%nVals
      write(logUnit(5),*) 'for pos=', ipos, 'dyn array is',         &
        &                 dyn_projectArray%val(ipos)%basisType,     &
        &                 dyn_projectArray%val(ipos)%maxpolydegree, &
        &                 dyn_projectArray%val(ipos)%header%kind

      call ply_poly_project_fillbody(                 &
        &    me         = proj_list(ipos),            &
        &    proj_init  = dyn_projectArray%val(ipos), &
        &    scheme_dim = scheme_dim                  )

      write(logUnit(5),*) ' for position', ipos, &
        &                 'of projection list, projection type is'
      write(logUnit(5),*) '  kind = ', proj_list(ipos)%kind
      write(logUnit(5),*) '  maxPolyDegree =', proj_list(ipos)%maxPolyDegree
      write(logUnit(5),*) '  basisType = ', proj_list(ipos)%basisType
      write(logUnit(5),*) '  oversamp_degree = ', &
        &                 proj_list(ipos)%oversamp_degree
      write(logUnit(5),*) '  min_degree = ', proj_list(ipos)%min_degree
      write(logUnit(5),*) '  nquadpointsPerDir = ', &
        &                 proj_list(ipos)%nquadpointsPerDir
      write(logUnit(5),*) '  lobattoPoints = ', proj_list(ipos)%lobattoPoints

    end do

  end subroutine ply_fill_project_list
  ! ------------------------------------------------------------------------ !


  ! ------------------------------------------------------------------------ !
  !> Fill the body of the projection with all required data,
  !! ply_poly_project_define has to be used beforehand to set necessary header
  !! information.
  subroutine ply_poly_project_fillbody(me, proj_init, scheme_dim)
    ! -------------------------------------------------------------------- !
    type(ply_poly_project_type), intent(inout)  :: me
    type(ply_prj_init_type), intent (in)    :: proj_init
    integer, intent(in) :: scheme_dim
    ! -------------------------------------------------------------------- !
    ! the oversampling order, need to get the number of modes in when
    ! oversampling is used
    integer :: oversampling_order
    ! the number of volume quadrature points per spatial direction
    integer :: numQuadPointsPerDir
    real(kind=rk) :: log_order, rem_log
    real(kind=rk) :: over_factor
    integer :: lb_log
    type(ply_nodes_header_type) :: nodes_header
    ! -------------------------------------------------------------------- !
    ! set the kind in the final projection type
    me%kind = trim(proj_init%header%kind)
    me%maxPolyDegree = proj_init%maxPolyDegree
    me%basisType = proj_init%basisType

    select case(trim(proj_init%header%kind))
    case('fpt')
      over_factor = proj_init%header%fpt_header%factor
      nodes_header = proj_init%header%fpt_header%nodes_header
      me%lobattopoints = nodes_header%lobattopoints
    case('l2p')
      over_factor = proj_init%header%l2p_header%factor
      nodes_header = proj_init%header%l2p_header%nodes_header
      me%lobattopoints = nodes_header%lobattopoints
    case('fxt')
      over_factor = proj_init%header%fxt_header%factor
      nodes_header = proj_init%header%fxt_header%nodes_header
      me%lobattopoints = .false.
    end select

    ! Find the oversampling order
    oversampling_order = ceiling( over_factor*(me%maxpolyDegree+1) )
    if (trim(proj_init%header%kind) == 'fpt') then
      ! Ensure an appropriate order in the oversampled polynomial
      ! representation.
      if ( proj_init%header%fpt_header%adapt_factor_pow2               &
        &  .and. (iand(oversampling_order, oversampling_order-1) /= 0) &
        &  ) then

        write(logUnit(1),*) '*** NOTE: oversampling order is increased to' &
          &                 //' next power of 2! ***'
        write(logUnit(2),*) '          original oversampling order would' &
          &                 // ' have been: ', oversampling_order
        ! Oversampling_order is not a power of 2, find the next power of 2
        ! and use that one instead...
        log_order = log(real(oversampling_order, kind=rk))/log(2.0_rk)
        lb_log = max(floor(log_order), 0)
        rem_log = log_order - lb_log
        if (rem_log > epsilon(log_order)*lb_log) then
          oversampling_order = 2**(lb_log+1)
        else
          oversampling_order = 2**lb_log
        end if

      end if
    end if

    if (me%lobattoPoints) then
      ! If lobatto Points are to be used, use at least one point more
      ! in the nodal representation.
      oversampling_order = max(oversampling_order, me%maxPolyDegree+2)
    end if

    if (trim(proj_init%header%kind) == 'fxt') then
      ! Ensure an even order for the oversampled polynomial representation.
      ! There is a bug in FXTPACK resulting faulty conversions for odd
      ! numbers of evaluation points.
      ! To avoid this, we increase the oversampling order by 1, if it is odd.
      oversampling_order = oversampling_order + mod(oversampling_order,2)
    end if

    write(logUnit(1),*) 'Using an oversampled order of: ', oversampling_order
    write(logUnit(1),*) 'maxpolydegree is: ', me%maxPolyDegree
    write(logUnit(2),*) '(Actual oversampling factor: ',                   &
      &                 real(oversampling_order)/real(me%maxpolydegree+1), &
      &                 ')'

    numQuadPointsPerDir = oversampling_order

    me%basisType = proj_init%basisType
    me%maxPolyDegree = proj_init%maxPolyDegree
    me%nquadpointsPerDir = numQuadPointsPerDir
    me%oversamp_degree = oversampling_order-1
    me%min_degree = min(me%maxPolyDegree, me%oversamp_degree)

    ! number of dof depending on q_space or p_space
    if (me%basisType == Q_space) then
       me%body_3d%ndofs = (me%maxPolyDegree+1)**3
       me%body_2d%ndofs = (me%maxPolyDegree+1)**2
       me%body_1d%ndofs = me%maxPolyDegree+1
       me%body_3d%min_dofs = (me%min_degree+1)**3
       me%body_2d%min_dofs = (me%min_degree+1)**2
       me%body_1d%min_dofs = me%min_degree+1
    else !p_space
  me%body_3d%ndofs = (((me%maxpolydegree) + 1) &
    &   * ((me%maxpolydegree) + 2) &
    &   * ((me%maxpolydegree) + 3)) &
    & / 6
  me%body_2d%ndofs = ((me%maxpolydegree)+1)*((me%maxpolydegree)+2)/2
  me%body_1d%ndofs = ((me%maxpolydegree)+1)
  me%body_3d%min_dofs = (((me%min_degree) + 1) &
    &   * ((me%min_degree) + 2) &
    &   * ((me%min_degree) + 3)) &
    & / 6
  me%body_2d%min_dofs = ((me%min_degree)+1)*((me%min_degree)+2)/2
  me%body_1d%min_dofs = ((me%min_degree)+1)
    end if

    me%body_3d%nquadpoints = numQuadPointsPerDir**3
    me%body_2d%nquadpoints = numQuadPointsPerDir**2
    me%body_1d%nquadpoints = numQuadPointsPerDir

    me%body_3d%oversamp_dofs = (oversampling_order)**3
    me%body_2d%oversamp_dofs = (oversampling_order)**2
    me%body_1d%oversamp_dofs = oversampling_order

    select case (trim(proj_init%header%kind))
    case('fpt')
      ! Fill fpt datatype
      call ply_init_legfpt(                                 &
        &    maxPolyDegree    = me%oversamp_degree,         &
        &    nIndeps          = 1,                          &
        &    fpt              = me%body_1d%fpt,             &
        &    header           = proj_init%header%fpt_header )

      if (scheme_dim >= 2) then
        call ply_init_legfpt(                                 &
          &    maxPolyDegree    = me%oversamp_degree,         &
          &    nIndeps          = me%oversamp_degree+1,       &
          &    fpt              = me%body_2d%fpt,             &
          &    header           = proj_init%header%fpt_header )
      end if

      if (scheme_dim >= 3) then
        call ply_init_legfpt(                                 &
          &    maxPolyDegree    = me%oversamp_degree,         &
          &    nIndeps          = (me%oversamp_degree+1)**2,  &
          &    fpt              = me%body_3D%fpt,             &
          &    header           = proj_init%header%fpt_header )
      end if

    case('l2p')
      ! Fill the L2 projection datatype
      if (scheme_dim >= 3) then
        call ply_init_l2p( l2p    = me%body_3d%l2p,              &
          &                header = proj_init%header%l2p_header, &
          &                degree = me%oversamp_degree           )
      end if

      if (scheme_dim >= 2) then
        call ply_init_l2p( l2p    = me%body_2d%l2p,              &
          &                header = proj_init%header%l2p_header, &
          &                degree = me%oversamp_degree           )
      end if

      call ply_init_l2p( l2p    = me%body_1d%l2p,              &
        &                header = proj_init%header%l2p_header, &
        &                degree = me%oversamp_degree           )

    case ('fxt')
      ! Fill the fxt Legendre Polynomial datatype
      if (scheme_dim >= 3) then
        call ply_init_fxt( fxt    = me%body_3d%fxt,              &
          &                header = proj_init%header%fxt_header, &
          &                degree = me%oversamp_degree           )
      end if

      if (scheme_dim >= 2) then
        call ply_init_fxt( fxt    = me%body_2d%fxt,              &
          &                header = proj_init%header%fxt_header, &
          &                degree = me%oversamp_degree           )
      end if
        call ply_init_fxt( fxt    = me%body_1d%fxt,              &
          &                header = proj_init%header%fxt_header, &
          &                degree = me%oversamp_degree           )

    case default
      write(logUnit(1),*) 'ERROR in initializing projection:'
      write(logUnit(1),*) 'Unknown projection method <', &
        &                 trim(proj_init%header%kind), &
        &                 '>'
      write(logUnit(1),*) 'Stopping....'
      call tem_abort()

    end select

    !> Initialization/Create of the volume quadrature nodes and the
    !! quadrature points on the face
    call ply_nodes_create(                           &
      &    me                = nodes_header,         &
      &    nodes             = me%body_1d%nodes,     &
      &    faces             = me%body_1d%faces,     &
      &    nQuadPointsPerDir = me%nQuadPointsPerDir, &
      &    nDims             = 1                     )

    if (scheme_dim >= 2) then
      call ply_nodes_create(                           &
        &    me                = nodes_header,         &
        &    nodes             = me%body_2d%nodes,     &
        &    faces             = me%body_2d%faces,     &
        &    nQuadPointsPerDir = me%nQuadPointsPerDir, &
        &    nDims             = 2                     )
    end if

    if (scheme_dim >= 3) then
      call ply_nodes_create(                           &
        &    me                = nodes_header,         &
        &    nodes             = me%body_3d%nodes,     &
        &    faces             = me%body_3d%faces,     &
        &    nQuadPointsPerDir = me%nQuadPointsPerDir, &
        &    nDims             = 3                     )
    end if

  end subroutine ply_poly_project_fillbody
  ! ------------------------------------------------------------------------ !


  ! ***************** MODAL to NODAL ***************************************** !
  ! ------------------------------------------------------------------------ !
  !> Convert nDoF modes to nodal values.
  subroutine ply_poly_project_m2n_multiVar(me, dim, nVars, modal_data, &
    &                                      nodal_data)
    ! -------------------------------------------------------------------- !
    type(ply_poly_project_type), intent(inout) :: me
    integer, intent(in) :: dim
    !> The number of variables to project. If a variable consists of more than
    !! one component, the number of components has to be passed. If there are
    !! more than one variable, the sum of all components has to be passed (e.g.
    !! 6 when there are two three-dimensional vectors).
    integer, intent(in) :: nVars
    real(kind=rk), intent(inout) :: modal_data(:,:)
    real(kind=rk), intent(inout) :: nodal_data(:,:)
    ! -------------------------------------------------------------------- !
    integer :: iVar
    ! -------------------------------------------------------------------- !

    select case(trim(me%kind))
    case ('l2p')
      ! for the projection modal to nodal, we do not need to distingusih
      ! between the spaces since the modal values results from the computation
      ! and the projection is on evluation of the nodes at the points with
      ! additional summation
      select case(dim)
      case (1)
        do iVar = 1, nVars
          call ply_l2p_trafo_1D( trafo = me%body_1D%l2p%leg2node, &
            &                    projected = nodal_data(:,iVar),  &
            &                    original  = modal_data(:,iVar)   )
        end do
      case (2)
        do iVar = 1, nVars
          call ply_l2p_trafo_2D( trafo = me%body_2D%l2p%leg2node, &
            &                    projected = nodal_data(:,iVar),  &
            &                    original  = modal_data(:,iVar)   )
        end do
      case (3)
        do iVar = 1, nVars
          call ply_l2p_trafo_3D( trafo = me%body_3D%l2p%leg2node, &
            &                    projected = nodal_data(:,iVar),  &
            &                    original  = modal_data(:,iVar)   )
        end do
      end select

    case ('fpt')

      select case (dim)
      case (3)
        call ply_LegToPnt_3D( fpt       = me%body_3d%fpt, &
          &                   pntVal    = nodal_data,     &
          &                   legCoeffs = modal_data,     &
          &                   nVars     = nVars           )
      case (2)
        call ply_LegToPnt_2D( fpt       = me%body_2d%fpt, &
          &                   pntVal    = nodal_data,     &
          &                   legCoeffs = modal_data,     &
          &                   nVars     = nVars           )
      case (1)
        do iVar = 1,nVars
          call me%body_1d%fpt%LegToPnt(          &
            &    pntVal    = nodal_data(:,iVar), &
            &    legCoeffs = modal_data(:,iVar), &
            &    nIndeps   = 1                   )
        end do
      end select

    case ('fxt')
      select case (dim)
      case (3)
        do iVar = 1,nVars
          call ply_fxt_m2n_3D( fxt = me%body_3d%fxt,            &
            &               modal_data = modal_data(:,iVar),    &
            &               nodal_data = nodal_data(:,iVar),    &
            &              oversamp_degree = me%oversamp_degree )
        end do
      case (2)
        do iVar = 1,nVars
          call ply_fxt_m2n_2D( fxt = me%body_2d%fxt,            &
            &               modal_data = modal_data(:,iVar),    &
            &               nodal_data = nodal_data(:,iVar),    &
            &              oversamp_degree = me%oversamp_degree )
        end do

      case (1)
        do iVar = 1,nVars
          call ply_fxt_m2n_1D( fxt = me%body_1d%fxt,         &
            &               modal_data = modal_data(:,iVar), &
            &               nodal_data = nodal_data(:,iVar)  )
        end do
      end select
    end select

  end subroutine ply_poly_project_m2n_multivar
  ! ------------------------------------------------------------------------ !



  ! **************** NODAL to MODAL ****************************************** !
  ! ------------------------------------------------------------------------ !
  !> Convert nodal values to nDoFs modes.
  subroutine ply_poly_project_n2m_multiVar(me, dim, nVars, nodal_data, &
    &                                      modal_data)
    ! -------------------------------------------------------------------- !
    type(ply_poly_project_type), intent(inout) :: me
    integer, intent(in) :: dim
    integer, intent(in) :: nVars
    real(kind=rk), intent(inout) :: nodal_data(:,:)
    real(kind=rk), intent(inout) :: modal_data(:,:)
    ! -------------------------------------------------------------------- !
    integer :: iVar
    ! -------------------------------------------------------------------- !

    select case(trim(me%kind))
    case ('l2p')
      select case (dim)
      case (1)
        do iVar = 1, nVars
          call ply_l2p_trafo_1D( trafo = me%body_1D%l2p%node2leg, &
            &                    projected = modal_data(:,iVar),  &
            &                    original  = nodal_data(:,iVar)   )
        end do
      case (2)
        do iVar = 1, nVars
          call ply_l2p_trafo_2D( trafo = me%body_2D%l2p%node2leg, &
            &                    projected = modal_data(:,iVar),  &
            &                    original  = nodal_data(:,iVar)   )
        end do
      case (3)
        do iVar = 1, nVars
          call ply_l2p_trafo_3D( trafo = me%body_3D%l2p%node2leg, &
            &                    projected = modal_data(:,iVar),  &
            &                    original  = nodal_data(:,iVar)   )
        end do
      end select

     case ('fpt')
      !projection via fpt
      select case (dim)
      case (3)
        call ply_pntToLeg_3D( fpt       = me%body_3d%fpt, &
          &                   nVars     = nVars,          &
          &                   pntVal    = nodal_data,     &
          &                   legCoeffs = modal_data      )
      case (2)
        call ply_pntToLeg_2D( fpt       = me%body_2d%fpt, &
          &                   nVars     = nVars,          &
          &                   pntVal    = nodal_data,     &
          &                   legCoeffs = modal_data      )
      case (1)
        do iVar = 1,nVars
          call me%body_1d%fpt%pntToLeg(          &
            &    nIndeps   = 1,                  &
            &    pntVal    = nodal_data(:,iVar), &
            &    legCoeffs = modal_data(:,iVar)  )
        end do
      end select

    case ('fxt')
      select case (dim)
      case (3)
        do iVar = 1, nVars
          call ply_fxt_n2m_3D(                          &
            &    fxt              = me%body_3d%fxt,     &
            &    nodal_data       = nodal_data(:,iVar), &
            &    modal_data       = modal_data(:,iVar), &
            &    oversamp_degree  = me%oversamp_degree  )
        end do

      case (2)
        do iVar = 1, nVars
          call ply_fxt_n2m_2D(                          &
            &    fxt              = me%body_2d%fxt,     &
            &    nodal_data       = nodal_data(:,iVar), &
            &    modal_data       = modal_data(:,iVar), &
            &    oversamp_degree  = me%oversamp_degree  )
        end do

      case (1)
        do iVar = 1, nVars
          call ply_fxt_n2m_1D(                          &
            &    fxt              = me%body_1d%fxt,     &
            &    nodal_data       = nodal_data(:,iVar), &
            &    modal_data       = modal_data(:,iVar)  )
        end do
      end select

    case default
       write(logUnit(1),*) 'ERROR in projection nodal to modal'
    end select

  end subroutine ply_poly_project_n2m_multivar
  ! ------------------------------------------------------------------------ !


  ! ------------------------------------------------------------------------ !
  !> function to provide the coordinates from the quadrature points on the faces
  ! idir and ialign are inputs which identify which face is needed
  ! faces is allocated as face(dir,align)
  subroutine ply_get_quadpoints_faces(poly_proj, idir, ialign, points)
    ! -------------------------------------------------------------------- !
    type(ply_poly_project_type), intent(in) :: poly_proj
    integer, intent(in)                     :: idir
    integer, intent(in)                     :: ialign
    real(kind=rk), allocatable, intent(inout) :: points (:,:)
    ! -------------------------------------------------------------------- !
    ! -------------------------------------------------------------------- !

     allocate (points(poly_proj%body_3d%faces(idir,iAlign)%nquadpoints,3))
     points = poly_proj%body_3d%faces(idir,iAlign)%points

  end subroutine ply_get_quadpoints_faces
  ! ------------------------------------------------------------------------ !


  ! ------------------------------------------------------------------------ !
  subroutine ply_get_quadpoints_faces_2d(poly_proj, idir, ialign, points)
    ! -------------------------------------------------------------------- !
    type(ply_poly_project_type), intent(in) :: poly_proj
    integer, intent(in)                     :: idir
    integer, intent(in)                     :: ialign
    real(kind=rk), allocatable, intent(out) :: points (:,:)
    ! -------------------------------------------------------------------- !
    ! -------------------------------------------------------------------- !

     allocate (points(poly_proj%body_2d%faces(idir,iAlign)%nquadpoints,3))
     points = poly_proj%body_2d%faces(idir,iAlign)%points

  end subroutine ply_get_quadpoints_faces_2d
  ! ------------------------------------------------------------------------ !


  ! ------------------------------------------------------------------------ !
  subroutine ply_get_quadpoints_faces_1d(poly_proj, idir, ialign, points)
    ! -------------------------------------------------------------------- !
    type(ply_poly_project_type), intent(in) :: poly_proj
    integer, intent(in)                     :: idir
    integer, intent(in)                     :: ialign
    real(kind=rk), allocatable, intent(out) :: points (:,:)
    ! -------------------------------------------------------------------- !
    ! -------------------------------------------------------------------- !
     allocate (points(poly_proj%body_1d%faces(idir,iAlign)%nquadpoints,3))
     points = poly_proj%body_1d%faces(idir,iAlign)%points
  end subroutine ply_get_quadpoints_faces_1d
  ! ------------------------------------------------------------------------ !

end module ply_poly_project_module