scale_element_hexahedral.F90

Go to the documentation of this file.

!> module FElib / Element / hexahedron
!!
!! @par Description
!!           A module for a hexahedral finite element
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
#include "scaleFElib.h"
module scale_element_hexahedral
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
 
  use scale_element_base, only: &
    elementbase3d, &
    elementbase3d_init, elementbase3d_final
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !  
  
  !> Derived type representing a hexahedral element
  type, public, extends(elementbase3d) :: hexahedralelement
  contains
    procedure :: init => hexhedralelement_init
    procedure :: final => hexhedralelement_final
    procedure :: genintgausslegendreintrpmat => hexhedralelement_gen_intgausslegendreintrpmat
  end type hexahedralelement
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !   
  private :: construct_element
 
contains
!> Initialize an object to manage a hexahedral element
!!
!! @param elem Object of finite element
!! @param elemOrder_h Polynomial order with 1D horizontal direction
!! @param elemOrder_v Polynomial order with vertical direction
!! @param LumpedMassMatFlag Flag whether mass lumping is considered
!OCL SERIAL
  subroutine hexhedralelement_init( &
      elem, elemOrder_h, elemOrder_v,      &
      LumpedMassMatFlag )
    implicit none
    
    class(hexahedralelement), intent(inout) :: elem
    integer, intent(in) :: elemOrder_h
    integer, intent(in) :: elemOrder_v
    logical, intent(in) :: LumpedMassMatFlag
 
    !-----------------------------------------------------------------------------
    
    elem%PolyOrder_h = elemorder_h
    elem%PolyOrder_v = elemorder_v
    elem%Nnode_h1D = elemorder_h + 1
    elem%Nnode_v = elemorder_v + 1
 
    elem%Nv = 8
    elem%Nfaces_h = 4
    elem%Nfaces_v = 2  
    elem%Nfaces = elem%Nfaces_h + elem%Nfaces_v
 
    elem%Nfp_h = elem%Nnode_h1D*elem%Nnode_v
    elem%Nfp_v = elem%Nnode_h1D**2
    elem%NfpTot = elem%Nfp_h*elem%Nfaces_h + elem%Nfp_v*elem%Nfaces_v
 
    elem%Np = elem%Nfp_v * elem%Nnode_v
    
    call elementbase3d_init(elem, lumpedmassmatflag)
    call construct_element(elem)
 
    return
  end subroutine hexhedralelement_init
 
!> Finalize an object to manage a hexahedral element
!!
!! @param elem Object of finite element
!OCL SERIAL  
  subroutine hexhedralelement_final(elem)
    implicit none
 
    class(hexahedralelement), intent(inout) :: elem
    !-----------------------------------------------------------------------------
 
    call elementbase3d_final(elem)
 
    return
  end subroutine hexhedralelement_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
    use scale_element_quadrilateral, only: quadrilateralelement
 
    implicit none
    
    type(hexahedralelement), intent(inout) :: elem
 
    integer :: nodes_ijk(elem%Nnode_h1D, elem%Nnode_h1D, elem%Nnode_v)
 
    real(RP) :: lglPts1D_h(elem%Nnode_h1D)
    real(RP) :: lglPts1D_v(elem%Nnode_v)
 
    real(DP) :: intWeight_lgl1DPts_h(elem%Nnode_h1D)
    real(DP) :: intWeight_lgl1DPts_v(elem%Nnode_v)
 
    real(RP) :: P1D_ori_h(elem%Nnode_h1D, elem%Nnode_h1D)
    real(RP) :: P1D_ori_v(elem%Nnode_v, elem%Nnode_v)
    real(RP) :: DP1D_ori_h(elem%Nnode_h1D, elem%Nnode_h1D)
    real(RP) :: DP1D_ori_v(elem%Nnode_v, elem%Nnode_v)
    real(RP) :: DLagr1D_h(elem%Nnode_h1D, elem%Nnode_h1D)
    real(RP) :: DLagr1D_v(elem%Nnode_v, elem%Nnode_v)
    real(RP) :: V2D_h(elem%Nfp_h, elem%Nfp_h)
    real(RP) :: V2D_v(elem%Nfp_v, elem%Nfp_v)
    real(RP) :: Emat(elem%Np, elem%NfpTot)
    real(RP) :: MassEdge_h(elem%Nfp_h, elem%Nfp_h)
    real(RP) :: MassEdge_v(elem%Nfp_v, elem%Nfp_v)
 
    integer :: i, j, k
    integer :: p1, p2, p3
    integer :: n, l, f
    integer :: Nord
    integer :: is, ie
 
    integer :: f_h, f_v
    integer :: fp, fp_h1, fp_h2, fp_v
 
    type(quadrilateralelement) :: elem2D
    !-----------------------------------------------------------------------------
 
    lglpts1d_h(:)      = polynomial_gengausslobattopt( elem%PolyOrder_h )
    p1d_ori_h(:,:)     = polynomial_genlegendrepoly( elem%PolyOrder_h, lglpts1d_h )
    dp1d_ori_h(:,:) = polynomial_gendlegendrepoly( elem%PolyOrder_h, lglpts1d_h, p1d_ori_h )
    dlagr1d_h(:,:) = polynomial_gendlagrangepoly_lglpt( elem%PolyOrder_h, lglpts1d_h )
 
    lglpts1d_v(:)      = polynomial_gengausslobattopt( elem%PolyOrder_v )
    p1d_ori_v(:,:)     = polynomial_genlegendrepoly( elem%PolyOrder_v, lglpts1d_v )
    dp1d_ori_v(:,:) = polynomial_gendlegendrepoly( elem%PolyOrder_v, lglpts1d_v, p1d_ori_v )
    dlagr1d_v(:,:) = polynomial_gendlagrangepoly_lglpt( elem%PolyOrder_v, lglpts1d_v )
    
    !* Preparation 
 
    do k=1, elem%Nnode_v   
    do j=1, elem%Nnode_h1D
    do i=1, elem%Nnode_h1D
        nodes_ijk(i,j,k) = i + (j-1)*elem%Nnode_h1D + (k-1)*elem%Nnode_h1D**2
    end do
    end do
    end do
    
    ! Set the mask to extract the values at faces
    
    elem%Fmask_h(:,1) = reshape(nodes_ijk(:,1,:), (/ elem%Nfp_h /))
    elem%Fmask_h(:,2) = reshape(nodes_ijk(elem%Nnode_h1D,:,:), (/ elem%Nfp_h /))
    elem%Fmask_h(:,3) = reshape(nodes_ijk(:,elem%Nnode_h1D,:), (/ elem%Nfp_h /))
    elem%Fmask_h(:,4) = reshape(nodes_ijk(1,:,:), (/ elem%Nfp_h /))
 
    elem%Fmask_v(:,1) = reshape(nodes_ijk(:,:,1), (/ elem%Nfp_v /))
    elem%Fmask_v(:,2) = reshape(nodes_ijk(:,:,elem%Nnode_v), (/ elem%Nfp_v /))
 
    !$acc update device(elem%Fmask_h, elem%Fmask_v)
    
    !- ColMask
 
    do j=1, elem%Nnode_h1D
    do i=1, elem%Nnode_h1D
      n = i + (j-1)*elem%Nnode_h1D
      elem%Colmask(:,n) = nodes_ijk(i,j,:)
    end do
    end do
    !$acc update device(elem%Colmask)
    
    != Hslice
 
    do k=1, elem%Nnode_v
      elem%Hslice(:,k) = reshape(nodes_ijk(:,:,k), (/ elem%Nfp_v /))
    end do
    !$acc update device(elem%Hslice)
 
    !- IndexH2Dto3D
 
    do k=1, elem%Nnode_v    
    do j=1, elem%Nnode_h1D
    do i=1, elem%Nnode_h1D
      n = i + (j-1)*elem%Nnode_h1D + (k-1)*elem%Nnode_h1D**2
      elem%IndexH2Dto3D(n) = nodes_ijk(i,j,1)
    end do
    end do    
    end do
    !$acc update device(elem%IndexH2Dto3D)
 
    !- IndexH2Dto3D_bnd
 
    call elem2d%Init( elem%PolyOrder_h, .false. )
 
    do f_h=1, 4
      do fp_v=1, elem%Nnode_v
      do fp_h1=1, elem%Nnode_h1D
        fp = fp_h1 + (fp_v-1)*elem%Nnode_h1D + (f_h-1)*elem%Nfp_h
        elem%IndexH2Dto3D_bnd(fp) = elem2d%Fmask(fp_h1,f_h)
      end do  
      end do
    end do
    do f_v=1, 2
      do fp_h2=1, elem%Nnode_h1D
      do fp_h1=1, elem%Nnode_h1D
        fp = fp_h1 + (fp_h2-1)*elem%Nnode_h1D    &
           + (f_v-1) * elem%Nfp_v                &
           + 4 * elem%Nnode_h1D * elem%Nnode_v
          elem%IndexH2Dto3D_bnd(fp) = fp_h1 + (fp_h2-1)*elem%Nnode_h1D
      end do  
      end do
    end do
    !$acc update device(elem%IndexH2Dto3D_bnd)
 
    call elem2d%Final()
 
    !- IndexZ1Dto3D
 
    do k=1, elem%Nnode_v    
    do j=1, elem%Nnode_h1D
    do i=1, elem%Nnode_h1D
      n = i + (j-1)*elem%Nnode_h1D + (k-1)*elem%Nnode_h1D**2
      elem%IndexZ1Dto3D(n) = k
    end do
    end do    
    end do
    !$acc update device(elem%IndexZ1Dto3D)
    
    !* Set the coordinates of LGL points, and the Vandermonde and differential matricies
 
    elem%Dx1(:,:) = 0.0_rp
    elem%Dx2(:,:) = 0.0_rp
    elem%Dx3(:,:) = 0.0_rp
 
    do k=1, elem%Nnode_v
    do j=1, elem%Nnode_h1D
    do i=1, elem%Nnode_h1D
      n = i + (j-1)*elem%Nnode_h1D + (k-1)*elem%Nnode_h1D**2
      
      !* Set the coordinates of LGL points
      elem%x1(n) = lglpts1d_h(i)
      elem%x2(n) = lglpts1d_h(j)
      elem%x3(n) = lglpts1d_v(k)
 
      !* Set the Vandermonde and differential matricies
      do p3=1, elem%Nnode_v
      do p2=1, elem%Nnode_h1D
      do p1=1, elem%Nnode_h1D
        l = p1 + (p2-1)*elem%Nnode_h1D + (p3-1)*elem%Nnode_h1D**2
        elem%V(n,l) = (p1d_ori_h(i,p1)*p1d_ori_h(j,p2)*p1d_ori_v(k,p3))                         &
                      * sqrt((dble(p1-1) + 0.5_dp)*(dble(p2-1) + 0.5_dp)*(dble(p3-1) + 0.5_dp))
  
        if(p2==j .and. p3==k) elem%Dx1(n,l) = dlagr1d_h(p1,i)
        if(p1==i .and. p3==k) elem%Dx2(n,l) = dlagr1d_h(p2,j)
        if(p1==i .and. p2==j) elem%Dx3(n,l) = dlagr1d_v(p3,k)
      end do
      end do
      end do
    end do
    end do
    end do
    elem%invV(:,:) = linalgebra_inv(elem%V)
    !$acc update device(elem%x1, elem%x2, elem%x3, elem%V, elem%Dx1, elem%Dx2, elem%Dx3, elem%invV)
    
    !* Set the weights at LGL points to integrate over element
 
    intweight_lgl1dpts_h(:) = polynomial_gengausslobattoptintweight(elem%PolyOrder_h)
    intweight_lgl1dpts_v(:) = polynomial_gengausslobattoptintweight(elem%PolyOrder_v)
 
    do k=1, elem%Nnode_v
    do j=1, elem%Nnode_h1D
    do i=1, elem%Nnode_h1D
        l = i + (j-1)*elem%Nnode_h1D + (k-1)*elem%Nnode_h1D**2
        elem%IntWeight_lgl(l) = &
          intweight_lgl1dpts_h(i) * intweight_lgl1dpts_h(j) * intweight_lgl1dpts_v(k)
    end do
    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 k=1, elem%Nnode_v
      do j=1, elem%Nnode_h1D
      do i=1, elem%Nnode_h1D
        l = i + (j-1)*elem%Nnode_h1D + (k-1)*elem%Nnode_h1D**2
        elem%M(l,l) = elem%IntWeight_lgl(l)
        elem%invM(l,l) = 1.0_dp/elem%IntWeight_lgl(l)
      end do
      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)
    call elementbase_construct_stiffmat( elem%M, elem%invM, elem%Dx3, elem%Np, & ! (in)
      elem%Sx3 )                                                                 ! (out)
    !$acc update device(elem%Sx1, elem%Sx2, elem%Sx3)
    
    !* Set the lift matrix
 
    do k=1, elem%Nnode_v
    do i=1, elem%Nnode_h1D
      n = i + (k-1)*elem%Nnode_h1D
      do p3=1, elem%Nnode_v
      do p1=1, elem%Nnode_h1D
        l = p1 + (p3-1)*elem%Nnode_h1D
        v2d_h(n,l) =   p1d_ori_h(i,p1)*p1d_ori_v(k,p3) &
                     * sqrt( (dble(p1-1) + 0.5_dp)*(dble(p3-1) + 0.5_dp) )
      end do
      end do
    end do
    end do
    do j=1, elem%Nnode_h1D
    do i=1, elem%Nnode_h1D
      n = i + (j-1)*elem%Nnode_h1D
      do p2=1, elem%Nnode_h1D
      do p1=1, elem%Nnode_h1D
        l = p1 + (p2-1)*elem%Nnode_h1D
        v2d_v(n,l) =   p1d_ori_h(i,p1)*p1d_ori_h(j,p2) &
                     * sqrt( (dble(p1-1) + 0.5_dp)*(dble(p2-1) + 0.5_dp) )
      end do
      end do
    end do
    end do
  
    !--
 
    emat(:,:) = 0.0_rp
    do f=1, elem%Nfaces_h
      if (elem%IsLumpedMatrix()) then
        massedge_h(:,:) = 0.0_rp
        do k=1, elem%Nnode_v
        do i=1, elem%Nnode_h1D
          l = i + (k-1)*elem%Nnode_h1D
          massedge_h(l,l) = intweight_lgl1dpts_h(i) * intweight_lgl1dpts_v(k)
        end do
        end do
      else
        call elementbase_construct_massmat( v2d_h, elem%Nfp_h, & ! (in)
          massedge_h )                                           ! (out)
      end if
 
      is = (f-1)*elem%Nfp_h + 1
      ie = is + elem%Nfp_h - 1
      emat(elem%Fmask_h(:,f), is:ie) = massedge_h
    end do
 
    do f=1, elem%Nfaces_v
      if (elem%IsLumpedMatrix()) then
        massedge_v(:,:) = 0.0_rp
        do j=1, elem%Nnode_h1D
        do i=1, elem%Nnode_h1D
          l = i + (j-1)*elem%Nnode_h1D
          massedge_v(l,l) = intweight_lgl1dpts_h(i) * intweight_lgl1dpts_h(j)
        end do
        end do
      else
        call elementbase_construct_massmat( v2d_v, elem%Nfp_v, & ! (in)
          massedge_v )                                           ! (out)
      end if
 
      is = elem%Nfaces_h*elem%Nfp_h + (f-1)*elem%Nfp_v + 1
      ie = is + elem%Nfp_v - 1
      emat(elem%Fmask_v(:,f), is:ie) = massedge_v
    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 hexhedralelement_gen_intgausslegendreintrpmat( this, IntrpPolyOrder, &
    intw_intrp, x_intrp, y_intrp, z_intrp ) result(IntrpMat)
 
    use scale_polynomial, only: &
      polynomial_genlegendrepoly,    &
      polynomial_gengausslegendrept, &
      polynomial_gengausslegendreptintweight
    
    implicit none
 
    class(hexahedralelement), intent(in) :: this
    integer, intent(in) :: IntrpPolyOrder
    real(RP), intent(out), optional :: intw_intrp(IntrpPolyOrder**3)
    real(RP), intent(out), optional :: x_intrp(IntrpPolyOrder**3)
    real(RP), intent(out), optional :: y_intrp(IntrpPolyOrder**3)
    real(RP), intent(out), optional :: z_intrp(IntrpPolyOrder**3)
    real(RP) :: IntrpMat(IntrpPolyOrder**3,this%Np)
 
    real(RP) :: r_int1D_i(IntrpPolyOrder)
    real(RP) :: r_int1Dw_i(IntrpPolyOrder)
    real(RP) :: P_int1D_ori_h(IntrpPolyOrder,this%Nnode_h1D)
    real(RP) :: P_int1D_ori_v(IntrpPolyOrder,this%Nnode_v)
    real(RP) :: Vint(IntrpPolyOrder**3,this%Np)
 
    integer :: p1, p2, p3, p1_, p2_, p3_
    integer :: n_, l_, m_
    !-----------------------------------------------------
 
    r_int1d_i(:) = polynomial_gengausslegendrept( intrppolyorder )
    r_int1dw_i(:) = polynomial_gengausslegendreptintweight( intrppolyorder )
    p_int1d_ori_h(:,:) = polynomial_genlegendrepoly( this%PolyOrder_h, r_int1d_i)
    p_int1d_ori_v(:,:) = polynomial_genlegendrepoly( this%PolyOrder_v, r_int1d_i)
 
    do p3_=1, intrppolyorder
    do p2_=1, intrppolyorder
    do p1_=1, intrppolyorder
      n_= p1_ + (p2_-1)*intrppolyorder + (p3_-1)*intrppolyorder**2
      if (present(intw_intrp)) intw_intrp(n_) = r_int1dw_i(p1_) * r_int1dw_i(p2_) * r_int1dw_i(p3_)
      if (present(x_intrp)) x_intrp(n_) = r_int1d_i(p1_)
      if (present(y_intrp)) y_intrp(n_) = r_int1d_i(p2_)
      if (present(z_intrp)) z_intrp(n_) = r_int1d_i(p3_)
 
      do p3=1, this%Nnode_v
      do p2=1, this%Nnode_h1D
      do p1=1, this%Nnode_h1D
        l_ = p1 + (p2-1)*this%Nnode_h1D + (p3-1)*this%Nnode_h1D**2
        vint(n_,l_) =  p_int1d_ori_h(p1_,p1) * sqrt(dble(p1-1) + 0.5_dp) &
                     * p_int1d_ori_h(p2_,p2) * sqrt(dble(p2-1) + 0.5_dp) &
                     * p_int1d_ori_v(p3_,p3) * sqrt(dble(p3-1) + 0.5_dp)
      end do
      end do
      end do
    end do
    end do
    end do
    intrpmat(:,:) = matmul(vint, this%invV)
 
    return
  end function hexhedralelement_gen_intgausslegendreintrpmat
 
end module scale_element_hexahedral