scale_element_base.F90

Go to the documentation of this file.

!> module FElib / Element / Base
!!
!! @par Description
!!           A base module for finite element
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_element_base
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision  
  
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
 
  !- Base
 
  !> Derived type representing an arbitrary finite element
  type, public :: elementbase
    integer :: np                       !< Number of nodes within an element
    integer :: nfaces                   !< Number of faces
    integer :: nfptot                   !< Total number of nodes on faces
    integer :: nv                       !< Number of vertices with an element
    logical, private :: lumpedmatflag   !< Flag whether the lumped mass matrix is used
 
    real(rp), allocatable :: v(:,:)             !< The Vandermonde matrix (V) whose size is Np x Np 
    real(rp), allocatable :: invv(:,:)          !< Inversion of the Vandermonde matrix (V^-1) whose size is Np x Np
    real(rp), allocatable :: m(:,:)             !< Mass matrix (M) whose size is Np x Np
    real(rp), allocatable :: invm(:,:)          !< Inversion of the mass matrix (M^-1) whose size Np x NP
    real(rp), allocatable :: lift(:,:)          !< Lifting matrix with element boundary integrals whose size Np x NfpTot
    real(rp), allocatable :: intweight_lgl(:)   !< Weights of gaussian quadrature with the LGL nodes
  contains
    procedure :: islumpedmatrix => elementbase_islumpedmatrix
  end type elementbase
  public :: elementbase_construct_massmat
  public :: elementbase_construct_stiffmat
  public :: elementbase_construct_liftmat
 
  !- 1D
 
  !> Derived type representing a 1D reference element
  type, public, extends(elementbase) :: elementbase1d
    integer :: polyorder                 !< Polynomial order
    integer :: nfp                       !< Number of nodes on an element face
    integer, allocatable :: fmask(:,:)   !< Array saving indices to extract nodal values on the faces
 
    real(rp), allocatable :: x1(:)       !< Array saving x1-coordinate of nodes within the reference element
 
    real(rp), allocatable :: dx1(:,:)    !< Elementwise differential matrix for the x1-coordinate direction (Dx1 = M^-1 Sx1)
 
    real(rp), allocatable :: sx1(:,:)    !< Elementwise stiffness matrix for the x1-coordinate direction
  end type elementbase1d
  
  public :: elementbase1d_init
  public :: elementbase1d_final
 
  !- 2D
 
  !> Derived type representing a 2D reference element
  type, public, extends(elementbase) :: elementbase2d
    integer :: polyorder               !< Polynomial order
    integer :: nfp                     !< Number of nodes on an element face
    integer, allocatable :: fmask(:,:) !< Array saving indices to extract nodal values on the faces
 
    real(rp), allocatable :: x1(:)     !< Array saving x1-coordinate of nodes within the reference element
    real(rp), allocatable :: x2(:)     !< Array saving x2-coordinate of nodes within the reference element
 
    real(rp), allocatable :: dx1(:,:)  !< Elementwise differential matrix for the x1-coordinate direction (Dx1 = M^-1 Sx1)
    real(rp), allocatable :: dx2(:,:)  !< Elementwise differential matrix for the x2-coordinate direction (Dx2 = M^-1 Sx2)
 
    real(rp), allocatable :: sx1(:,:)  !< Elementwise stiffness matrix for the x1-coordinate direction
    real(rp), allocatable :: sx2(:,:)  !< Elementwise stiffness matrix for the x2-coordinate direction
  contains
    procedure :: generate_l2projmat => elementbase2d_gen_l2projmat
    procedure :: generate_interpmat => elementbase2d_gen_interpmat
    procedure :: generate_modaltruncationmat => elementbase2d_gen_modaltruncationmat  
  end type elementbase2d
 
  public :: elementbase2d_init
  public :: elementbase2d_final
 
  interface
    function elementbase2d_genintgausslegendreintrpmat(this, IntrpPolyOrder, &
      intw_intrp, x_intrp, y_intrp ) result(IntrpMat)
  
      import elementbase2d
      import rp
      class(elementbase2d), intent(in) :: this
      integer, intent(in) :: intrppolyorder
      real(rp), intent(out), optional :: intw_intrp(intrppolyorder**2)
      real(rp), intent(out), optional :: x_intrp(intrppolyorder**2)
      real(rp), intent(out), optional :: y_intrp(intrppolyorder**2)
      real(rp) :: intrpmat(intrppolyorder**2,this%np)
    end function elementbase2d_genintgausslegendreintrpmat
  end interface
 
  !- 3D
 
  !> Derived type representing a 3D reference element
  type, public, extends(elementbase) :: elementbase3d    
    integer :: polyorder_h               !< Polynomial order with the horizontal direction
    integer :: nnode_h1d                 !< Number of nodes along the horizontal coordinate
    integer :: nfaces_h                  !< Number of nodes on an horizontal face of the reference element
    integer :: nfp_h                     !< Number of horizontal faces of the reference element
    integer, allocatable :: fmask_h(:,:) !< Array saving indices to extract nodal values on the horizontal faces
 
    integer :: polyorder_v               !< Polynomial order with the vertical direction
    integer :: nnode_v                   !< Number of nodes along the vertical coordinate
    integer :: nfaces_v                  !< Number of nodes on an vertical face of the reference element
    integer :: nfp_v                     !< Number of vertical faces of the reference element
    integer, allocatable :: fmask_v(:,:) !< Array saving indices to extract nodal values on the vertical faces
 
    integer, allocatable :: colmask(:,:)        !< Array saving indices to extract nodal values on the vertical columns
    integer, allocatable :: hslice(:,:)         !< Array saving indices to extract nodal values on the horizontal plane
    integer, allocatable :: indexh2dto3d(:)     !< Array saving indices to expand 2D horizontal nodal values into the 3D nodal values
    integer, allocatable :: indexh2dto3d_bnd(:) !< Array saving indices to expand 2D horizontal nodal values into the 3D nodal values on element faces
    integer, allocatable :: indexz1dto3d(:)     !< Array saving indices to expand 1D vertical nodal values into the 3D nodal values
 
    real(rp), allocatable :: x1(:) !< Array saving x1-coordinate of nodes within the reference element
    real(rp), allocatable :: x2(:) !< Array saving x2-coordinate of nodes within the reference element
    real(rp), allocatable :: x3(:) !< Array saving x3-coordinate of nodes within the reference element
    
    real(rp), allocatable :: dx1(:,:) !< Elementwise differential matrix for the x1-coordinate direction (Dx1 = M^-1 Sx1)
    real(rp), allocatable :: dx2(:,:) !< Elementwise differential matrix for the x2-coordinate direction (Dx2 = M^-1 Sx2)
    real(rp), allocatable :: dx3(:,:) !< Elementwise differential matrix for the x3-coordinate direction (Dx3 = M^-1 Sx3)
  
    real(rp), allocatable :: sx1(:,:) !< Elementwise stiffness matrix for the x1-coordinate direction
    real(rp), allocatable :: sx2(:,:) !< Elementwise stiffness matrix for the x2-coordinate direction    
    real(rp), allocatable :: sx3(:,:) !< Elementwise stiffness matrix for the x3-coordinate direction
  contains
    procedure :: generate_l2projmat => elementbase3d_gen_l2projmat
    procedure :: generate_interpmat => elementbase3d_gen_interpmat
    procedure :: generate_modaltruncationmat => elementbase3d_gen_modaltruncationmat
  end type elementbase3d
  
  public :: elementbase3d_init
  public :: elementbase3d_final
 
  !-----------------------------------------------------------------------------
  !
  !++ Private type & procedure
  ! 
 
  private :: elementbase_init
  private :: elementbase_final
 
  private :: elementbase2d_gen_nodaltransfermat
  private :: elementbase3d_gen_nodaltransfermat
 
contains
  !-- Base Element ------------------------------------------------------------------------------
 
  !> Initialize a base object to manage a reference element
!OCL SERIAL
  subroutine elementbase_init( elem, lumpedmat_flag )
    implicit none
    class(elementbase), intent(inout) :: elem
    logical, intent(in) :: lumpedmat_flag      !< Flag whether mass lumping is considered
    !-----------------------------------------------------------------------------
 
    allocate( elem%M(elem%Np, elem%Np) )
    allocate( elem%invM(elem%Np, elem%Np) )
    allocate( elem%V(elem%Np, elem%Np) )
    allocate( elem%invV(elem%Np, elem%Np) )
    allocate( elem%Lift(elem%Np, elem%NfpTot) )    
    
    allocate( elem%IntWeight_lgl(elem%Np) )
    !$acc enter data create( elem%M, elem%invM, elem%V, elem%invV, elem%Lift, elem%IntWeight_lgl )
    
    elem%LumpedMatFlag = lumpedmat_flag
    return
  end subroutine elementbase_init
 
  !> Finalize a base object to manage a reference element
!OCL SERIAL
  subroutine elementbase_final( elem )
    implicit none
    class(elementbase), intent(inout) :: elem
    !-----------------------------------------------------------------------------
    if ( allocated(elem%M) ) then
      !$acc exit data delete( elem%M, elem%invM, elem%V, elem%invV, elem%Lift, elem%IntWeight_lgl )
      deallocate( elem%M )
      deallocate( elem%invM )
      deallocate( elem%V )
      deallocate( elem%invV )
      deallocate( elem%Lift )
      deallocate( elem%IntWeight_lgl )
    end if
 
    return
  end subroutine elementbase_final
 
  !> Get a flag whether the lumped mass matrix is used
!OCL SERIAL
  function elementbase_islumpedmatrix( elem ) result(lumpedmat_flag)
    implicit none
    class(elementbase), intent(in) :: elem
    logical :: lumpedmat_flag
    !---------------------------------------------
 
    lumpedmat_flag = elem%LumpedMatFlag
    return
  end function elementbase_islumpedmatrix
 
 
  !> Construct mass matrix
  !!  M^-1 = V V^T
  !!  M = ( M^-1 )^-1
!OCL SERIAL
  subroutine elementbase_construct_massmat( V, Np, &
    MassMat, invMassMat )
    use scale_linalgebra, only: linalgebra_inv
    implicit none
    integer, intent(in) :: np
    real(rp), intent(in) :: v(np,np)
    real(rp), intent(out) :: massmat(np,np)
    real(rp), intent(out), optional :: invmassmat(np,np)
 
    real(rp) :: tmpmat(np,np)
    real(rp) :: invm(np,np)
    !------------------------------------
 
    tmpmat(:,:) = transpose(v)
    invm(:,:) = matmul( v, tmpmat )
    massmat(:,:) = linalgebra_inv( invm )
 
    if ( present(invmassmat) ) invmassmat(:,:) = invm(:,:)
    return
  end subroutine elementbase_construct_massmat
 
  !> Construct stiffness matrix
  !!  StiffMat_i = M^-1 ( M D_xi )^T
!OCL SERIAL
  subroutine elementbase_construct_stiffmat( MassMat, invMassMat, DMat, Np, &
      StiffMat )
    implicit none
    integer, intent(in) :: np
    real(rp), intent(in) :: massmat(np,np)
    real(rp), intent(in) :: invmassmat(np,np)
    real(rp), intent(in) :: dmat(np,np)
    real(rp), intent(out) :: stiffmat(np,np)
 
    real(rp) :: tmpmat1(np,np)
    real(rp) :: tmpmat2(np,np)
    !------------------------------------
 
    tmpmat1(:,:) = matmul( massmat, dmat )
    tmpmat2(:,:) = transpose( tmpmat1 )
    stiffmat(:,:) = matmul( invmassmat, tmpmat2 )
 
    return
  end subroutine elementbase_construct_stiffmat
 
  !> Construct stiffness matrix
  !!  StiffMat_i = M^-1 ( M D_xi )^T
!OCL SERIAL
  subroutine elementbase_construct_liftmat( invM, EMat, Np, NfpTot, &
      LiftMat )
    implicit none
    integer, intent(in) :: np
    integer, intent(in) :: nfptot
    real(rp), intent(in) :: invm(np,np)
    real(rp), intent(in) :: emat(np,nfptot)
    real(rp), intent(out) :: liftmat(np,nfptot)
    !------------------------------------
 
    liftmat(:,:) = matmul( invm, emat )
    return
  end subroutine elementbase_construct_liftmat 
 
  !-- 1D Element ------------------------------------------------------------------------------
 
  !> Initialize an object to manage a 1D reference element
  !!
  !! @param elem Object of finite element
  !! @param elem Flag whether mass lumping is considered
!OCL SERIAL
  subroutine elementbase1d_init( elem, lumpedmat_flag )
    implicit none
    class(elementbase1d), intent(inout) :: elem
    logical, intent(in) :: lumpedmat_flag    
    !-----------------------------------------------------------------------------
 
    call elementbase_init( elem, lumpedmat_flag )
 
    !$acc enter data create( elem )
    !$acc update device( elem%Np, elem%Nfaces, elem%NfpTot, elem%Nv )
 
    allocate( elem%x1(elem%Np) )  
    allocate( elem%Fmask(elem%Nfp, elem%Nfaces) )
    
    allocate( elem%Dx1(elem%Np, elem%Np) )
    allocate( elem%Sx1(elem%Np, elem%Np) )
    !$acc enter data create( elem%x1, elem%Dx1, elem%Sx1, elem%Fmask )
    return
  end subroutine elementbase1d_init
 
!> Finalize an object to manage a 1D reference element
!OCL SERIAL
  subroutine elementbase1d_final( elem )
    implicit none
    class(elementbase1d), intent(inout) :: elem
    !-----------------------------------------------------------------------------
 
    if ( allocated( elem%x1 ) )  then
      !$acc exit data delete( elem%x1, elem%Dx1, elem%Sx1, elem%Fmask )
      !$acc exit data delete( elem )
      deallocate( elem%x1 )
      deallocate( elem%Dx1 )
      deallocate( elem%Sx1 )
      deallocate( elem%Fmask )
    end if
 
    call elementbase_final( elem )
 
    return
  end subroutine elementbase1d_final
 
  !-- 2D Element ------------------------------------------------------------------------------
 
  !> Initialize an object to manage a 2D reference element
  !!
  !! @param elem Object of finite element
  !! @param elem Flag whether mass lumping is considered
!OCL SERIAL
  subroutine elementbase2d_init( elem, lumpedmat_flag )
    implicit none
    class(elementbase2d), intent(inout) :: elem
    logical, intent(in) :: lumpedmat_flag    
    !-----------------------------------------------------------------------------
 
    call elementbase_init( elem, lumpedmat_flag )
 
    !$acc enter data create( elem )
    !$acc update device( elem%Np, elem%Nfaces, elem%NfpTot, elem%Nv )
 
    allocate( elem%x1(elem%Np), elem%x2(elem%Np) )
    allocate( elem%Fmask(elem%Nfp, elem%Nfaces) )
 
    allocate( elem%Dx1(elem%Np, elem%Np), elem%Dx2(elem%Np, elem%Np) )
    allocate( elem%Sx1(elem%Np, elem%Np), elem%Sx2(elem%Np, elem%Np) )
    !$acc enter data create( elem%x1, elem%x2, elem%Dx1, elem%Dx2, elem%Sx1, elem%Sx2, elem%Fmask )
 
    return
  end subroutine elementbase2d_init
 
  !> Finalize an object to manage a 2D reference element
!OCL SERIAL
  subroutine elementbase2d_final( elem )
    implicit none
    class(elementbase2d), intent(inout) :: elem
    !-----------------------------------------------------------------------------
 
    if ( allocated( elem%x1 ) )  then
      !$acc exit data delete( elem%x1, elem%x2, elem%Dx1, elem%Dx2, elem%Sx1, elem%Sx2, elem%Fmask )
      !$acc exit data delete( elem )
      deallocate( elem%x1, elem%x2 )
      deallocate( elem%Fmask )
 
      deallocate( elem%Dx1, elem%Dx2 )
      deallocate( elem%Sx1, elem%Sx2 )
    end if
 
    call elementbase_final( elem )
 
    return
  end subroutine elementbase2d_final
 
  !> Generate a projection matrix for L2 projection.
  !! This matrix maps nodal values on elem_in to nodal values on elem. It is intended for p-restriction, i.e. elem order <= elem_in order.
  !!
  !! With an orthogonal Legendre modal basis, this corresponds to L2 projection onto the polynomial space represented by elem.
  !!
!OCL SERIAL
  subroutine elementbase2d_gen_l2projmat( elem, &
      elem_in,  &
      L2ProjMat )
    implicit none
    class(elementbase2d), intent(in) :: elem
    class(elementbase2d), intent(in) :: elem_in
    real(rp), intent(out) :: l2projmat(elem%np,elem_in%np)
    !---------------------------------------------
 
    call elementbase2d_gen_nodaltransfermat( elem, elem_in, &
        elem%PolyOrder,                                     &
        l2projmat )
    return
  end subroutine elementbase2d_gen_l2projmat
 
  !> Generate an interpolation matrix by p-prolongation.
  !! This matrix maps nodal values on elem_in to nodal values on elem. It is intended for p-prolongation, i.e. elem order >= elem_in order.
  !!
  !! In Legendre modal space, higher modes not present in elem_in are set to zero.
!OCL SERIAL
  subroutine elementbase2d_gen_interpmat( elem, &
      elem_in,  &
      InterpMat )
    implicit none
    class(elementbase2d), intent(in) :: elem
    class(elementbase2d), intent(in) :: elem_in
    real(rp), intent(out) :: interpmat(elem%np,elem_in%np)
    !---------------------------------------------
 
    call elementbase2d_gen_nodaltransfermat( elem, elem_in, &
        elem_in%PolyOrder,                                  &
        interpmat )
    return
  end subroutine elementbase2d_gen_interpmat  
  
  !> Generate a nodal-to-nodal transfer matrix to remove several high modes
!OCL SERIAL
  subroutine elementbase2d_gen_modaltruncationmat( elem, &
      polyOrder_tr, &
      TruncateMat )
    implicit none
    class(elementbase2d), intent(in) :: elem
    integer, intent(in) :: polyorder_tr
    real(rp), intent(out) :: truncatemat(elem%np,elem%np)
    !---------------------------------------------
 
    call elementbase2d_gen_nodaltransfermat( elem, elem, &
      polyorder_tr,                                      &
      truncatemat )
    return
  end subroutine elementbase2d_gen_modaltruncationmat
 
  !-- 3D Element ------------------------------------------------------------------------------
  
!> Initialize an object to manage a 3D reference element
!!
!! @param elem Object of finite element
!! @param elem Flag whether mass lumping is considered
!OCL SERIAL
  subroutine elementbase3d_init( elem, lumpedmat_flag )
    implicit none
    class(elementbase3d), intent(inout) :: elem
    logical, intent(in) :: lumpedmat_flag
    !-----------------------------------------------------------------------------
 
    call elementbase_init( elem, lumpedmat_flag )
 
    !$acc enter data create( elem )
    !$acc update device( elem%Np, elem%Nfaces, elem%NfpTot, elem%Nv )
 
    allocate( elem%x1(elem%Np), elem%x2(elem%Np), elem%x3(elem%Np) )
    allocate( elem%Dx1(elem%Np, elem%Np), elem%Dx2(elem%Np, elem%Np), elem%Dx3(elem%Np, elem%Np) )
    allocate( elem%Sx1(elem%Np, elem%Np), elem%Sx2(elem%Np, elem%Np), elem%Sx3(elem%Np, elem%Np) )
    allocate( elem%Fmask_h(elem%Nfp_h, elem%Nfaces_h), elem%Fmask_v(elem%Nfp_v, elem%Nfaces_v) )
    allocate( elem%Colmask(elem%Nnode_v,elem%Nfp_v))
    allocate( elem%Hslice(elem%Nfp_v,elem%Nnode_v) )
    allocate( elem%IndexH2Dto3D(elem%Np) )
    allocate( elem%IndexH2Dto3D_bnd(elem%NfpTot) )
    allocate( elem%IndexZ1Dto3D(elem%Np) )
    !$acc enter data create( elem%x1, elem%x2, elem%x3,                &
    !$acc  elem%Dx1, elem%Dx2, elem%Dx3, elem%Sx1, elem%Sx2, elem%Sx3, &
    !$acc  elem%Fmask_h, elem%Fmask_v, elem%Colmask, elem%Hslice,      &
    !$acc  elem%IndexH2Dto3D, elem%IndexH2Dto3D_bnd, elem%IndexZ1Dto3D )
    return
  end subroutine elementbase3d_init
 
!> Finalize an object to manage a 3D reference element
!OCL SERIAL
  subroutine elementbase3d_final( elem )
    implicit none
    class(elementbase3d), intent(inout) :: elem
    !-----------------------------------------------------------------------------
 
    if ( allocated( elem%x1 ) )  then
      !$acc exit data delete( elem%x1, elem%x2, elem%x3,                 &
      !$acc  elem%Dx1, elem%Dx2, elem%Dx3, elem%Sx1, elem%Sx2, elem%Sx3, &
      !$acc  elem%Fmask_h, elem%Fmask_v, elem%Colmask, elem%Hslice,      &
      !$acc  elem%IndexH2Dto3D, elem%IndexH2Dto3D_bnd, elem%IndexZ1Dto3D )
      !$acc exit data delete( elem )
      deallocate( elem%x1, elem%x2, elem%x3 )
      deallocate( elem%Dx1, elem%Dx2, elem%Dx3 )
      deallocate( elem%Sx1, elem%Sx2, elem%Sx3 )
      deallocate( elem%Fmask_h, elem%Fmask_v )
      deallocate( elem%Colmask, elem%Hslice )
      deallocate( elem%IndexH2Dto3D, elem%IndexH2Dto3D_bnd )
      deallocate( elem%IndexZ1Dto3D )
    end if
 
    call elementbase_final( elem )    
 
    return
  end subroutine elementbase3d_final
 
  !> Generate a projection matrix for L2 projection.
  !! This matrix maps nodal values on elem_in to nodal values on elem. It is intended for p-restriction, i.e. elem order <= elem_in order.
  !!
  !! With an orthogonal Legendre modal basis, this corresponds to L2 projection onto the polynomial space represented by elem.
  !!
!OCL SERIAL
  subroutine elementbase3d_gen_l2projmat( elem, &
      elem_in,  &
      L2ProjMat )
    implicit none
    class(elementbase3d), intent(in) :: elem
    class(elementbase3d), intent(in) :: elem_in
    real(rp), intent(out) :: l2projmat(elem%np,elem_in%np)
    !---------------------------------------------
 
    call elementbase3d_gen_nodaltransfermat( elem, elem_in, &
        elem%PolyOrder_h, elem%PolyOrder_v,                 &
        l2projmat )
    return
  end subroutine elementbase3d_gen_l2projmat
 
  !> Generate an interpolation matrix by p-prolongation.
  !! This matrix maps nodal values on elem_in to nodal values on elem. It is intended for p-prolongation, i.e. elem order >= elem_in order.
  !!
  !! In Legendre modal space, higher modes not present in elem_in are set to zero.
!OCL SERIAL
  subroutine elementbase3d_gen_interpmat( elem, &
      elem_in,  &
      InterpMat )
    implicit none
 
    class(elementbase3d), intent(in) :: elem
    class(elementbase3d), intent(in) :: elem_in
    real(rp), intent(out) :: interpmat(elem%np,elem_in%np)
    !---------------------------------------------
 
    call elementbase3d_gen_nodaltransfermat( elem, elem_in, &
        elem_in%PolyOrder_h, elem_in%PolyOrder_v,           &
        interpmat )
    return
  end subroutine elementbase3d_gen_interpmat  
  
  !> Generate a nodal-to-nodal transfer matrix to remove several high modes
!OCL SERIAL
  subroutine elementbase3d_gen_modaltruncationmat( elem, &
      polyOrder_h_tr, polyOrder_v_tr, &
      TruncateMat )
    implicit none
 
    class(elementbase3d), intent(in) :: elem
    integer, intent(in) :: polyorder_h_tr
    integer, intent(in) :: polyorder_v_tr
    real(rp), intent(out) :: truncatemat(elem%np,elem%np)
    !---------------------------------------------
 
    call elementbase3d_gen_nodaltransfermat( elem, elem, &
      polyorder_h_tr, polyorder_v_tr,                    &
      truncatemat )
    return
  end subroutine elementbase3d_gen_modaltruncationmat
 
!--- private procedures ------------------------------------------------------------------------------
 
 
  !> Generate a nodal-to-nodal transfer matrix for 2D element.
  !!   Mat = V_elem  * S * invV_elem_in
  !! This matrix maps nodal values on elem_in to nodal values on elem.
  !! In modal space, common modes are copied and the other modes are zero.
  !!
  !! If elem has lower order than elem_in, this corresponds to modal truncation, i.e. L2 projection under an orthogonal Legendre modal basis.
  !! If elem has higher order than elem_in, this corresponds to p-prolongation by modal zero-padding.
!OCL SERIAL
  subroutine elementbase2d_gen_nodaltransfermat( elem, elem_in, &
    pmax, &
    TransferMat )
    implicit none
    class(elementbase2d), intent(in) :: elem                  !< Object to manage a reference element
    class(elementbase2d), intent(in) :: elem_in               !< Object to manage a reference element
    integer, intent(in) :: pmax                               !< Maximum polynomial order of the common modes between elem and elem_in
    real(rp), intent(out) :: transfermat(elem%np,elem_in%np)  !< Nodal-to-nodal transfer matrix from elem_in to elem
 
    integer :: p1, p2
    integer :: p_out, p_in
 
    real(rp) :: invv_in(elem%np,elem_in%np)
    !---------------------------------------------
 
    invv_in(:,:) = 0.0_rp
    do p2=1, pmax+1
    do p1=1, pmax+1
      p_out = p1 + (p2-1)*(elem%PolyOrder    + 1)
      p_in  = p1 + (p2-1)*(elem_in%PolyOrder + 1)
      invv_in(p_out,:) = elem_in%invV(p_in,:)
    end do
    end do
 
    transfermat(:,:) = matmul(elem%V, invv_in)
    return
  end subroutine elementbase2d_gen_nodaltransfermat
  
  !> Generate a nodal-to-nodal transfer matrix for 3D element.
  !!   Mat = V_elem  * S * invV_elem_in
  !! This matrix maps nodal values on elem_in to nodal values on elem.
  !! In modal space, common modes are copied and the other modes are zero.
  !!
  !! If elem has lower order than elem_in, this corresponds to modal truncation, i.e. L2 projection under an orthogonal Legendre modal basis.
  !! If elem has higher order than elem_in, this corresponds to p-prolongation by modal zero-padding.
!OCL SERIAL
  subroutine elementbase3d_gen_nodaltransfermat( elem, elem_in, &
    pmax_h, pmax_v, &
    TransferMat )
    implicit none
    class(elementbase3d), intent(in) :: elem                  !< Object to manage a reference element
    class(elementbase3d), intent(in) :: elem_in               !< Object to manage a reference element
    integer, intent(in) :: pmax_h                             !< Maximum polynomial order of the common modes between elem and elem_in
    integer, intent(in) :: pmax_v                             !< Maximum polynomial order of the common modes between elem and elem_in
    real(rp), intent(out) :: transfermat(elem%np,elem_in%np)  !< Nodal-to-nodal transfer matrix from elem_in to elem
 
    integer :: p1, p2, p3
    integer :: p_out, p_in
 
    real(rp) :: invv_in(elem%np,elem_in%np)
    !---------------------------------------------
 
    invv_in(:,:) = 0.0_rp
    do p3=1, pmax_v+1
    do p2=1, pmax_h+1
    do p1=1, pmax_h+1
      p_out = p1 + (p2-1)*(elem%PolyOrder_h    + 1) + (p3-1)*(elem%PolyOrder_h    + 1)**2
      p_in  = p1 + (p2-1)*(elem_in%PolyOrder_h + 1) + (p3-1)*(elem_in%PolyOrder_h + 1)**2
      invv_in(p_out,:) = elem_in%invV(p_in,:)
    end do
    end do
    end do
 
    transfermat(:,:) = matmul(elem%V, invv_in)
    return
  end subroutine elementbase3d_gen_nodaltransfermat
 
end module scale_element_base