scale_element_operation_tensorprod3D_interp.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Element / Interpolation with 3D tensor product elements
!!
!! @par Description
!!           A module for providing interpolation operations assuming a 3D tensor product element with (p+1)^3 DOF
!!
!! @author Yuta Kawai, Xuanzhengbo Ren, and Team SCALE
!!
!<
#include "scaleFElib.h"
module scale_element_operation_tensorprod3d_interp
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_prc, only: prc_abort 
 
  use scale_element_base, only: elementbase3d
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
 
  !> Derived type for interpolation operation with 3D tensor product element
  type, public :: elementoptrtensorprod3dinterp
    class(elementbase3d), pointer :: elem  
    integer :: polyorder_h_intrp
    integer :: nnodeh1d_gl
    integer :: np3d_hglvgll
    integer :: np3d
    real(rp), allocatable :: intrpmat1d_gll2gl(:,:)
    real(rp), allocatable :: intrpmat1d_gll2gl_tr(:,:)
 
    real(rp),  allocatable :: v1d_gl(:,:)
    real(rp),  allocatable :: inv_v1d_gl(:,:)
    real(rp), allocatable :: intrpmat1d_gl2gll(:,:)
    real(rp), allocatable :: intrpmat1d_gl2gll_tr(:,:)
 
  contains
    procedure :: init => elementoptrtensorprod3dinterp_init
    procedure :: final => elementoptrtensorprod3dinterp_final
    procedure :: do_gll2gl => elementoptrtensorprod3dinterp_gll2gl
    procedure :: do_gl2gll => elementoptrtensorprod3dinterp_gl2gll
  end type elementoptrtensorprod3dinterp
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  private :: intrp_core
  
contains
!OCL SERIAL
  subroutine elementoptrtensorprod3dinterp_init( this, elem, PolyOrder_h_ip )
    use scale_linalgebra, only: linalgebra_inv
    use scale_polynomial, only: &
      polynomial_genlegendrepoly,    &
      polynomial_gengausslegendrept, &
      polynomial_gengausslegendreptintweight
    use scale_element_line, only: lineelement
    use scale_prc
    implicit none
    class(elementoptrtensorprod3dinterp), intent(inout) :: this
    type(elementbase3d), intent(in), target :: elem
    integer, intent(in) :: Polyorder_h_ip
 
    real(RP) :: r_int1D_i(Polyorder_h_ip+1)
    real(RP) :: r_int1Dw_i(Polyorder_h_ip+1)
    real(RP) :: P_int1D_ori_h(Polyorder_h_ip+1,elem%PolyOrder_h+1)
    real(RP) :: Vint1D(Polyorder_h_ip+1,elem%PolyOrder_h+1)
 
    real(RP) :: P1D_gl_h(Polyorder_h_ip+1,Polyorder_h_ip+1)
    real(RP) :: P1D_gll_h(elem%Nnode_h1D,Polyorder_h_ip+1)
    real(RP) :: Vint1D_gl(elem%Nnode_h1D,Polyorder_h_ip+1)
 
    integer :: p1, p1_
    integer :: n_, l_
 
    type(lineelement) :: elem1D
    !----------------------------------------
 
    this%elem => elem
    this%PolyOrder_h_intrp = polyorder_h_ip
    this%NnodeH1D_GL = polyorder_h_ip + 1
    this%Np3D_hGLvGLL = this%NnodeH1D_GL**2 * elem%Nnode_v
    this%Np3D = elem%Np
 
    allocate( this%IntrpMat1D_gll2gl(this%NnodeH1D_GL,elem%Nnode_h1D) )
    allocate( this%IntrpMat1D_gll2gl_tr(elem%Nnode_h1D,this%NnodeH1D_GL) )
 
    !-
    call elem1d%Init( elem%PolyOrder_h, .false. )
    r_int1d_i(:) = polynomial_gengausslegendrept( this%NnodeH1D_GL )
    r_int1dw_i(:) = polynomial_gengausslegendreptintweight(this%NnodeH1D_GL )
    p_int1d_ori_h(:,:) = polynomial_genlegendrepoly( elem%PolyOrder_h, r_int1d_i )
 
    do p1_=1,this%NnodeH1D_GL
      n_ = p1_
      do p1=1, elem%Nnode_h1D
        l_ = p1
        vint1d(n_,l_) = p_int1d_ori_h(p1_,p1) * sqrt( real(p1-1,kind=rp) + 0.5_rp )
      end do
    end do
    this%IntrpMat1D_gll2gl(:,:) = matmul( vint1d, elem1d%invV )
    this%IntrpMat1D_gll2gl_tr(:,:) = transpose(this%IntrpMat1D_gll2gl)
 
    !-
    allocate( this%V1D_GL(this%NnodeH1D_GL,polyorder_h_ip+1) )
    allocate( this%inv_V1D_GL(polyorder_h_ip+1,this%NnodeH1D_GL) )
    allocate( this%IntrpMat1D_gl2gll(elem%Nnode_h1D,this%NnodeH1D_GL) )
    allocate( this%IntrpMat1D_gl2gll_tr(this%NnodeH1D_GL,elem%Nnode_h1D) )
 
    p1d_gl_h(:,:) = polynomial_genlegendrepoly( polyorder_h_ip, r_int1d_i )
    p1d_gll_h(:,:) = polynomial_genlegendrepoly( polyorder_h_ip, elem1d%x1 )
 
    do p1_=1, this%NnodeH1D_GL
      n_ = p1_
      do p1=1, polyorder_h_ip+1
        l_ = p1
        this%V1D_GL(n_,l_) = p1d_gl_h(p1_,p1) * sqrt( real(p1-1,kind=rp) + 0.5_rp )
      end do
    end do
    this%inv_V1D_GL(:,:) = linalgebra_inv( this%V1D_GL(:,:) )
 
    vint1d_gl(:,:) = 0.0_rp
    do p1_=1,this%elem%Nnode_h1D
      n_ = p1_
      do p1=1, this%NnodeH1D_GL
        l_ = p1
        vint1d_gl(n_,l_) = p1d_gll_h(p1_,p1) * sqrt( real(p1-1,kind=rp) + 0.5_rp )
      end do
    end do
    this%IntrpMat1D_gl2gll(:,:) = matmul( vint1d_gl, this%inv_V1D_GL )
    this%IntrpMat1D_gl2gll_tr(:,:) = transpose(this%IntrpMat1D_gl2gll(:,:))
 
    !-
    call elem1d%Final()
    return
  end subroutine elementoptrtensorprod3dinterp_init
 
!OCL SERIAL
  subroutine elementoptrtensorprod3dinterp_final( this )
    implicit none
    class(elementoptrtensorprod3dinterp), intent(inout) :: this
    !----------------------------------------
    deallocate( this%IntrpMat1D_gll2gl, this%IntrpMat1D_gll2gl_tr )
 
    deallocate( this%V1D_GL, this%inv_V1D_GL )
    deallocate( this%IntrpMat1D_gl2gll, this%IntrpMat1D_gl2gll_tr )
    return
  end subroutine elementoptrtensorprod3dinterp_final
 
!OCL SERIAL
  subroutine elementoptrtensorprod3dinterp_gll2gl( this, q, q_intrp )
    implicit none
    class(elementoptrtensorprod3dinterp), intent(in) :: this
    real(RP), intent(in) :: q(this%Np3D)
    real(RP), intent(out) :: q_intrp(this%Np3D_hGLvGLL)
    !----------------------------------------
 
    call intrp_core( q_intrp, &
      q, this%elem%Nnode_h1D, this%NnodeH1D_GL, this%elem%Nnode_v, &
      this%IntrpMat1D_gll2gl_tr )
    return
  end subroutine elementoptrtensorprod3dinterp_gll2gl
 
!OCL SERIAL
  subroutine elementoptrtensorprod3dinterp_gl2gll( this, q, q_intrp )
    implicit none
    class(elementoptrtensorprod3dinterp), intent(in) :: this
    real(RP), intent(in) :: q(this%Np3D_hGLvGLL)
    real(RP), intent(out) :: q_intrp(this%Np3D)
    !----------------------------------------
 
    call intrp_core( q_intrp, &
      q, this%NnodeH1D_GL, this%elem%Nnode_h1D, this%elem%Nnode_v, &
      this%IntrpMat1D_gl2gll_tr )
    return
  end subroutine elementoptrtensorprod3dinterp_gl2gll
 
!---------
!OCL SERIAL
  subroutine intrp_core( q_intrp, &
    q, Nnode_h1D, Nnode_h1D_ip, Nnode_v, IntrpMat1D_tr )
    integer, intent(in) :: Nnode_h1D
    integer, intent(in) :: Nnode_h1D_ip
    integer, intent(in) :: Nnode_v
    real(RP), intent(out) :: q_intrp(Nnode_h1D_ip,Nnode_h1D_ip,Nnode_v)
    real(RP), intent(in) :: q(Nnode_h1D,Nnode_h1D,Nnode_v)
    real(RP), intent(in) :: IntrpMat1D_tr(Nnode_h1D,Nnode_h1D_ip)
 
    integer :: p1, p2, p3
    integer :: pp
    real(RP) :: tmp
    real(RP) :: q_tmp(Nnode_h1D_ip,Nnode_h1D)
    !--------------------------------
 
    do p3=1, nnode_v
      do p2=1, nnode_h1d
      do p1=1, nnode_h1d_ip
        tmp = 0.0_rp
        do pp=1, nnode_h1d
          tmp = tmp + intrpmat1d_tr(pp,p1) * q(pp,p2,p3)
        end do
        q_tmp(p1,p2)= tmp
      end do
      end do
 
      q_intrp(:,:,p3) = 0.0_rp
      do p2=1, nnode_h1d_ip
      do pp=1, nnode_h1d
        do p1=1, nnode_h1d_ip
          q_intrp(p1,p2,p3) = q_intrp(p1,p2,p3) + intrpmat1d_tr(pp,p2) * q_tmp(p1,pp)
        end do
      end do
      end do
    end do
    return
  end subroutine intrp_core
 
end module scale_element_operation_tensorprod3d_interp