scale_element_line.F90

Go to the documentation of this file.

!> module FElib / Element / line
!!
!! @par Description
!!           A module for a line finite element
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
#include "scaleFElib.h"
module scale_element_line
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
 
  use scale_element_base, only: &
    elementbase1d, &
    elementbase1d_init, elementbase1d_final
  
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !  
 
  !> Derived type representing a line element
  type, public, extends(elementbase1d) :: lineelement
  contains
    procedure :: init => lineelement_init
    procedure :: final => lineelement_final
    procedure :: genintgausslegendreintrpmat => lineelement_gen_intgausslegendreintrpmat    
  end type lineelement
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !   
  private :: construct_element
  
contains
 
!> Initialize an object to manage a line element
!!
!! @param elem Object of finite element
!! @param elemOrder Polynomial order
!! @param LumpedMassMatFlag Flag whether mass lumping is considered
!OCL SERIAL
  subroutine lineelement_init( &
    elem, elemOrder,           &
    LumpedMassMatFlag )
    
    implicit none
 
    class(lineelement), intent(inout) :: elem
    integer, intent(in) :: elemOrder
    logical, intent(in) :: LumpedMassMatFlag
 
    !-----------------------------------------------------------------------------
    
    elem%PolyOrder = elemorder
    elem%Nv = 2
    elem%Np = elemorder + 1
    elem%Nfp = 1
    elem%Nfaces = 2
    elem%NfpTot = elem%Nfp*elem%Nfaces
       
    call elementbase1d_init(elem, lumpedmassmatflag)
    call construct_element(elem)
 
    return
  end subroutine lineelement_init
 
!> Finalize an object to manage a line element
!!
!! @param elem Object of finite element
!OCL SERIAL  
  subroutine lineelement_final(elem)
    implicit none
 
    class(lineelement), intent(inout) :: elem
    !-----------------------------------------------------------------------------
 
    call elementbase1d_final(elem)
 
    return
  end subroutine lineelement_final
 
  !> Construct the element matrices and coordinates of LGL points
!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(lineelement), intent(inout) :: elem
 
    integer :: nodes(elem%Np)
 
    real(RP) :: lglPts1D(elem%Np)
    real(DP) :: intWeight_lgl1DPts(elem%Np)
 
    real(RP) :: P1D_ori(elem%Np, elem%Np)
    real(RP) :: DP1D_ori(elem%Np, elem%Np)
    real(RP) :: DLagr1D(elem%Np, elem%Np)
    real(RP) :: Emat(elem%Np, elem%Nfp*elem%Nfaces)
    real(RP) :: MassEdge(elem%Nfp, elem%Nfp)
 
    integer :: i
    integer :: p1
    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 i=1, elem%Np
        nodes(i) = i
    end do
 
    ! Set the mask to extract the values at faces
    
    elem%Fmask(:,1) = 1
    elem%Fmask(:,2) = elem%Np
    !$acc update device(elem%Fmask)
 
    !* Set the coordinates of LGL points, and the Vandermonde and differential matricies
 
    elem%Dx1(:,:) = 0.0_rp
  
    do n=1, elem%Np      
      !* Set the coordinates of LGL points
      elem%x1(n) = lglpts1d(n)
 
      !* Set the Vandermonde and differential matricies
      do l=1, elem%Np
        elem%V(n,l) = p1d_ori(n,l) * sqrt(dble(l-1) + 0.5_rp)
        elem%Dx1(n,l) = dlagr1d(l,n)
      end do
    end do
    elem%invV(:,:) = linalgebra_inv(elem%V)
    !$acc update device(elem%x1, elem%V, elem%Dx1, elem%invV)
    
    !* Set the weights at LGL points to integrate over element
 
    elem%IntWeight_lgl(:) = polynomial_gengausslobattoptintweight(elem%PolyOrder)
    !$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 i=1, elem%Np
        elem%M(i,i) = elem%IntWeight_lgl(i)
        elem%invM(i,i) = 1.0_rp/elem%IntWeight_lgl(i)
      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)
    !$acc update device(elem%Sx1)
 
    !* Set the lift matrix
 
    emat(:,:) = 0.0_rp  
    do f=1, elem%Nfaces
      massedge(:,:) = 0.0_rp
      do l=1, elem%Nfp
        massedge(l,l) = 1.0_rp
      end do  
      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 lineelement_gen_intgausslegendreintrpmat( this, IntrpPolyOrder, &
    intw_intrp, x_intrp ) result(IntrpMat)
 
    use scale_polynomial, only: &
      polynomial_genlegendrepoly,    &
      polynomial_gengausslegendrept, &
      polynomial_gengausslegendreptintweight
    
    implicit none
 
    class(lineelement), intent(in) :: this
    integer, intent(in) :: IntrpPolyOrder
    real(RP), intent(out), optional :: intw_intrp(IntrpPolyOrder)
    real(RP), intent(out), optional :: x_intrp(IntrpPolyOrder)
    real(RP) :: IntrpMat(IntrpPolyOrder,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,this%PolyOrder+1)
 
    integer :: p1, p1_
    !-----------------------------------------------------
 
    r_int1d_i(:) = polynomial_gengausslegendrept( intrppolyorder )
    r_int1dw_i(:) = polynomial_gengausslegendreptintweight( intrppolyorder )
    p_int1d_ori(:,:) = polynomial_genlegendrepoly( this%PolyOrder, r_int1d_i)
 
    do p1_=1, intrppolyorder
      if (present(intw_intrp)) intw_intrp(p1_) = r_int1dw_i(p1_)
      if (present(x_intrp)) x_intrp(p1_) = r_int1d_i(p1_)
      do p1=1, this%Np
        vint(p1_,p1) =  p_int1d_ori(p1_,p1) * sqrt(real(p1-1,kind=rp) + 0.5_rp)
      end do
    end do
    intrpmat(:,:) = matmul(vint, this%invV)
 
    return
  end function lineelement_gen_intgausslegendreintrpmat
 
end module scale_element_line