scale_element_quadrilateral.F90

Go to the documentation of this file.

!> module FElib / Element / Quadrilateral
!!
!! @par Description
!!           A module for a quadrilateral finite element
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
#include "scaleFElib.h"
module scale_element_quadrilateral
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
 
  use scale_element_base, only: &
    elementbase2d, &
    elementbase2d_init, elementbase2d_final
  
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !  
 
  !> Derived type representing a quadrilateral element
  type, public, extends(elementbase2d) :: quadrilateralelement
  contains
    procedure :: init => quadrilateralelement_init
    procedure :: final => quadrilateralelement_final
    procedure :: genintgausslegendreintrpmat => quadrilateralelement_gen_intgausslegendreintrpmat
  end type quadrilateralelement
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !   
  private :: construct_element
 
contains
 
!> Initialize an object to manage a hexahedral element
!!
!! @param elem Object of finite element
!! @param elemOrder Polynomial order with 1D direction
!! @param LumpedMassMatFlag Flag whether mass lumping is considered
!OCL SERIAL
  subroutine quadrilateralelement_init( &
      elem, elemOrder,  &
      LumpedMassMatFlag )
    
    implicit none
 
    class(quadrilateralelement), intent(inout) :: elem
    integer, intent(in) :: elemOrder
    logical, intent(in) :: LumpedMassMatFlag
 
    !-----------------------------------------------------------------------------
    
    elem%PolyOrder = elemorder
    elem%Nv = 4
    elem%Np = (elemorder + 1)**2
    elem%Nfp = elemorder + 1
    elem%Nfaces = 4
    elem%NfpTot = elem%Nfp*elem%Nfaces
       
    call elementbase2d_init(elem, lumpedmassmatflag)
    call construct_element(elem)
 
    return
  end subroutine quadrilateralelement_init
 
!> Finalize an object to manage a quadrilateral element
!!
!! @param elem Object of finite element
!OCL SERIAL
  subroutine quadrilateralelement_final(elem)
    implicit none
 
    class(quadrilateralelement), intent(inout) :: elem
    !-----------------------------------------------------------------------------
 
    call elementbase2d_final(elem)
 
    return
  end subroutine quadrilateralelement_final
 
!OCL SERIAL
  subroutine construct_element(elem)
 
    use scale_linalgebra, only: linalgebra_inv
    use scale_polynomial, only: &
      polynomial_gengausslobattopt, polynomial_gengausslobattoptintweight,   &
      polynomial_genlegendrepoly, polynomial_gendlegendrepoly,               &
      polynomial_genlagrangepoly, polynomial_gendlagrangepoly_lglpt
    use scale_element_base, only: & 
      elementbase_construct_massmat, elementbase_construct_stiffmat, &
      elementbase_construct_liftmat
 
    implicit none
 
    type(quadrilateralelement), intent(inout) :: elem
 
    integer :: nodes_ij(elem%Nfp, elem%Nfp)
 
    real(RP) :: lglPts1D(elem%Nfp)
    real(DP) :: intWeight_lgl1DPts(elem%Nfp)
 
    real(RP) :: P1D_ori(elem%Nfp, elem%Nfp)
    real(RP) :: DP1D_ori(elem%Nfp, elem%Nfp)
    real(RP) :: DLagr1D(elem%Nfp, elem%Nfp)
    real(RP) :: V1D(elem%Nfp, elem%Nfp)
    real(RP) :: Emat(elem%Np, elem%NfpTot)
    real(RP) :: MassEdge(elem%Nfp, elem%Nfp)
 
    real(RP) :: eta, etac
    real(RP) :: filter1D(elem%Nfp), filter2D(elem%Np)
 
    integer :: i, j
    integer :: p1, p2
    integer :: n, l, f
    integer :: Nord
    !-----------------------------------------------------------------------------
 
    lglpts1d(:)      = polynomial_gengausslobattopt( elem%PolyOrder )
    p1d_ori(:,:)     = polynomial_genlegendrepoly( elem%PolyOrder, lglpts1d )
    dp1d_ori(:,:) = polynomial_gendlegendrepoly( elem%PolyOrder, lglpts1d, p1d_ori )
    dlagr1d(:,:) = polynomial_gendlagrangepoly_lglpt(elem%PolyOrder, lglpts1d)
    
    !* Preparation 
    
    do j=1, elem%Nfp
    do i=1, elem%Nfp
        nodes_ij(i,j) = i + (j-1)*elem%Nfp
    end do
    end do
 
    ! Set the mask to extract the values at faces
    
    elem%Fmask(:,1) = nodes_ij(:,1)
    elem%Fmask(:,2) = nodes_ij(elem%Nfp,:)
    elem%Fmask(:,3) = nodes_ij(:,elem%Nfp)
    elem%Fmask(:,4) = nodes_ij(1,:)
    !$acc update device(elem%Fmask)
 
    !* Set the coordinates of LGL points, and the Vandermonde and differential matricies
 
    elem%Dx1(:,:) = 0.0_rp
    elem%Dx2(:,:) = 0.0_rp
 
    do j=1, elem%Nfp
    do i=1, elem%Nfp
      n = i + (j-1)*elem%Nfp
      
      !* Set the coordinates of LGL points
      elem%x1(n) = lglpts1d(i)
      elem%x2(n) = lglpts1d(j)
 
      !* Set the Vandermonde and differential matricies
      do p2=1, elem%Nfp
      do p1=1, elem%Nfp
        l = p1 + (p2-1)*elem%Nfp
        elem%V(n,l) = (p1d_ori(i,p1)*p1d_ori(j,p2))                       &
                      * sqrt((dble(p1-1) + 0.5_dp)*(dble(p2-1) + 0.5_dp))
  
        if(p2==j) elem%Dx1(n,l) = dlagr1d(p1,i)
        if(p1==i) elem%Dx2(n,l) = dlagr1d(p2,j)
      end do
      end do
    end do
    end do
    elem%invV(:,:) = linalgebra_inv(elem%V)
    !$acc update device(elem%x1, elem%x2, elem%V, elem%Dx1, elem%Dx2, elem%invV)
    
    !* Set the weights at LGL points to integrate over element
 
    intweight_lgl1dpts(:) = polynomial_gengausslobattoptintweight(elem%PolyOrder)
        
    do j=1, elem%Nfp
    do i=1, elem%Nfp
        l = i + (j - 1)*elem%Nfp
        elem%IntWeight_lgl(l) = &
          intweight_lgl1dpts(i) * intweight_lgl1dpts(j)
    end do
    end do
    !$acc update device(elem%IntWeight_lgl)
 
    !* Set the mass matrix
 
    if (elem%IsLumpedMatrix()) then
      elem%invM(:,:) = 0.0_rp
      elem%M(:,:)    = 0.0_rp
      do j=1, elem%Nfp
      do i=1, elem%Nfp
        l = i + (j - 1)*elem%Nfp
        elem%M(l,l) = elem%IntWeight_lgl(l)
        elem%invM(l,l) = 1.0_rp/elem%IntWeight_lgl(l)
      end do
      end do      
    else
      call elementbase_construct_massmat( elem%V, elem%Np, & ! (in)
        elem%M, elem%invM )                                  ! (out)
    end if
    !$acc update device(elem%M, elem%invM)
 
    !* Set the stiffness matrix
 
    call elementbase_construct_stiffmat( elem%M, elem%invM, elem%Dx1, elem%Np, & ! (in)
      elem%Sx1 )                                                                 ! (out)
    call elementbase_construct_stiffmat( elem%M, elem%invM, elem%Dx2, elem%Np, & ! (in)
      elem%Sx2 )                                                                 ! (out)
    !$acc update device(elem%Sx1, elem%Sx2)
 
    !* Set the lift matrix
 
    do p1=1, elem%Nfp
      v1d(:,p1) = p1d_ori(:,p1)*sqrt(dble(p1-1) + 0.5_dp)
    end do
 
    emat(:,:) = 0.0_rp
    do f=1, elem%Nfaces
 
      if (elem%IsLumpedMatrix()) then
        massedge = 0.0_rp
        do l=1, elem%Nfp
          massedge(l,l) = intweight_lgl1dpts(l)
        end do  
      else
        call elementbase_construct_massmat( v1d, elem%Nfp, & ! (in)
          massedge )                                         ! (out)
      end if
 
      emat(elem%Fmask(:,f), (f-1)*elem%Nfp+1:f*elem%Nfp) = massedge
    end do
    call elementbase_construct_liftmat( elem%invM, emat, elem%Np, elem%NfpTot, & ! (in)
      elem%Lift )                                                                ! (out)
    !$acc update device(elem%Lift)
  
    return
  end subroutine construct_element
 
!OCL SERIAL
  function quadrilateralelement_gen_intgausslegendreintrpmat( this, IntrpPolyOrder, &
    intw_intrp, x_intrp, y_intrp ) result(IntrpMat)
 
    use scale_polynomial, only: &
      polynomial_genlegendrepoly,    &
      polynomial_gengausslegendrept, &
      polynomial_gengausslegendreptintweight
    
    implicit none
 
    class(quadrilateralelement), 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)
 
    real(RP) :: r_int1D_i(IntrpPolyOrder)
    real(RP) :: r_int1Dw_i(IntrpPolyOrder)
    real(RP) :: P_int1D_ori(IntrpPolyOrder,this%PolyOrder+1)
    real(RP) :: Vint(IntrpPolyOrder**2,(this%PolyOrder+1)**2)
 
    integer :: p1, p2, p1_, p2_
    integer :: n_, l_
    !-----------------------------------------------------
 
    r_int1d_i(:) = polynomial_gengausslegendrept( intrppolyorder )
    r_int1dw_i(:) = polynomial_gengausslegendreptintweight( intrppolyorder )
    p_int1d_ori(:,:) = polynomial_genlegendrepoly( this%PolyOrder, r_int1d_i)
 
    do p2_=1, intrppolyorder
    do p1_=1, intrppolyorder
      n_= p1_ + (p2_-1)*intrppolyorder
      if (present(intw_intrp)) intw_intrp(n_) = r_int1dw_i(p1_) * r_int1dw_i(p2_)
      if (present(x_intrp)) x_intrp(n_) = r_int1d_i(p1_)
      if (present(y_intrp)) y_intrp(n_) = r_int1d_i(p2_)
      
      do p2=1, this%Nfp
      do p1=1, this%Nfp
        l_ = p1 + (p2-1)*this%Nfp
        vint(n_,l_) =  p_int1d_ori(p1_,p1) * sqrt(dble(p1-1) + 0.5_dp) &
                     * p_int1d_ori(p2_,p2) * sqrt(dble(p2-1) + 0.5_dp)
      end do
      end do
    end do
    end do
    intrpmat(:,:) = matmul(vint, this%invV)
 
    return
  end function quadrilateralelement_gen_intgausslegendreintrpmat
 
  !-------------------
 
end module scale_element_quadrilateral