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
  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(:,:) !< Number of vertical faces of the reference element
 
    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
  end type elementbase3d
  
  public :: elementbase3d_init
  public :: elementbase3d_final
 
  !-----------------------------------------------------------------------------
  !
  !++ Private type & procedure
  ! 
 
  private :: elementbase_init
  private :: elementbase_final
 
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
    !-----------------------------------------------------------------------------
 
    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) )  
 
    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
      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 )
 
    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) )
 
    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
      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 )
 
    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) )
 
    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
      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
 
  !-- 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 )
 
    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) )
  
    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
      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
 
end module scale_element_base