scale_element_base.F90

Go to the documentation of this file.

 
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_element_base
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision  
  
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
 
  !- Base
 
  type, public :: elementbase
    integer :: np
    integer :: nfaces
    integer :: nfptot
    integer :: nv
    logical, private :: lumpedmatflag
 
    real(rp), allocatable :: v(:,:)
    real(rp), allocatable :: invv(:,:)
    real(rp), allocatable :: m(:,:)
    real(rp), allocatable :: invm(:,:)
    real(rp), allocatable :: lift(:,:)
    real(rp), allocatable :: intweight_lgl(:)
  contains
    procedure :: islumpedmatrix => elementbase_islumpedmatrix
  end type elementbase
  public :: elementbase_construct_massmat
  public :: elementbase_construct_stiffmat
  public :: elementbase_construct_liftmat
 
  !- 1D
 
  type, public, extends(elementbase) :: elementbase1d
    integer :: polyorder
    integer :: nfp
    integer, allocatable :: fmask(:,:)
 
    real(rp), allocatable :: x1(:) 
 
    real(rp), allocatable :: dx1(:,:)
 
    real(rp), allocatable :: sx1(:,:) 
  end type elementbase1d
  
  public :: elementbase1d_init
  public :: elementbase1d_final
 
  !- 2D
 
  type, public, extends(elementbase) :: elementbase2d
    integer :: polyorder
    integer :: nfp
    integer, allocatable :: fmask(:,:)
 
    real(rp), allocatable :: x1(:)
    real(rp), allocatable :: x2(:)    
 
    real(rp), allocatable :: dx1(:,:)
    real(rp), allocatable :: dx2(:,:)
 
    real(rp), allocatable :: sx1(:,:)
    real(rp), allocatable :: sx2(:,:) 
  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
 
  type, public, extends(elementbase) :: elementbase3d    
    integer :: polyorder_h
    integer :: nnode_h1d
    integer :: nfaces_h
    integer :: nfp_h
    integer, allocatable :: fmask_h(:,:)
 
    integer :: polyorder_v
    integer :: nnode_v
    integer :: nfaces_v
    integer :: nfp_v
    integer, allocatable :: fmask_v(:,:)
 
    integer, allocatable :: colmask(:,:)
    integer, allocatable :: hslice(:,:)
    integer, allocatable :: indexh2dto3d(:)
    integer, allocatable :: indexh2dto3d_bnd(:)    
    integer, allocatable :: indexz1dto3d(:)
 
    real(rp), allocatable :: x1(:)
    real(rp), allocatable :: x2(:)
    real(rp), allocatable :: x3(:)
    
    real(rp), allocatable :: dx1(:,:)
    real(rp), allocatable :: dx2(:,:)    
    real(rp), allocatable :: dx3(:,:)
  
    real(rp), allocatable :: sx1(:,:)
    real(rp), allocatable :: sx2(:,:)    
    real(rp), allocatable :: sx3(:,:)
  end type elementbase3d
  
  public :: elementbase3d_init
  public :: elementbase3d_final
 
  !-----------------------------------------------------------------------------
  !
  !++ Private type & procedure
  ! 
 
  private :: elementbase_init
  private :: elementbase_final
 
contains
  !-- Base 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
 
!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
 
!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
 
 
!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
 
!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
 
!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 ------------------------------------------------------------------------------
 
!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
 
!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 ------------------------------------------------------------------------------
 
!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
 
!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 ------------------------------------------------------------------------------
  
!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
 
!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