hecmw_Jacob_341.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
 
module hecmw_jacob341
contains
 
  subroutine hecmw_jacob_341 ( hecMESH, iElem, DET, W, N, NX, NY, NZ)
    use hecmw_util
    implicit none
 
    type(hecmwst_local_mesh):: hecmesh
    integer(kind=kint)::       ielem
    real(kind=kreal)::         det
    real(kind=kreal)::         w(4),n(4,4),nx(4,4),ny(4,4),nz(4,4)
 
    integer(kind=kint):: i, ii, j, ilocal
    real(kind=kreal):: dum
    real(kind=kreal):: xx(4), yy(4), zz(4)
    real(kind=kreal):: nr(4), ns(4), nt(4)
    real(kind=kreal):: jacob(3,3), jinv(3,3)
    real(kind=kreal),parameter:: alpha = 0.58541020, beta = 0.13819660
 
    !!
    !! ******** Set values of coordinates********
    !!
    do i = 1, 4
      ii = 4 - i
      ilocal = hecmesh%elem_node_item( 4*ielem -ii )
      xx(i) = hecmesh%node( ilocal*3 -2 )
      yy(i) = hecmesh%node( ilocal*3 -1 )
      zz(i) = hecmesh%node( ilocal*3    )
      w(i)  = 0.25d0
    end do
 
    !!
    !! ******** Set values to shape functions ********
    !!
    n(1,1) = alpha
    n(2,1) = beta
    n(3,1) = beta
    n(4,1) = beta
 
    n(1,2) = beta
    n(2,2) = alpha
    n(3,2) = beta
    n(4,2) = beta
 
    n(1,3) = beta
    n(2,3) = beta
    n(3,3) = alpha
    n(4,3) = beta
 
    n(1,4) = beta
    n(2,4) = beta
    n(3,4) = beta
    n(4,4) = alpha
    !!
    !! ******** Derivative of shape functions ********
    !!
    !! ----------- For R-Coordinate -------------
    nr(1) = 1.0
    nr(2) = 0.0
    nr(3) = 0.0
    nr(4) =-1.0
    !! ----------- For S-Coordinate -------------
    ns(1) = 0.0
    ns(2) = 1.0
    ns(3) = 0.0
    ns(4) =-1.0
    !! ----------- For T-Coordinate -------------
    nt(1) = 0.0
    nt(2) = 0.0
    nt(3) = 1.0
    nt(4) =-1.0
    !!
    !! ******** Jacobi matrix calculation********
    !!
    jacob(1,1) = nr(1)*xx(1)+nr(2)*xx(2)+nr(3)*xx(3)+nr(4)*xx(4)
    jacob(2,1) = ns(1)*xx(1)+ns(2)*xx(2)+ns(3)*xx(3)+ns(4)*xx(4)
    jacob(3,1) = nt(1)*xx(1)+nt(2)*xx(2)+nt(3)*xx(3)+nt(4)*xx(4)
    jacob(1,2) = nr(1)*yy(1)+nr(2)*yy(2)+nr(3)*yy(3)+nr(4)*yy(4)
    jacob(2,2) = ns(1)*yy(1)+ns(2)*yy(2)+ns(3)*yy(3)+ns(4)*yy(4)
    jacob(3,2) = nt(1)*yy(1)+nt(2)*yy(2)+nt(3)*yy(3)+nt(4)*yy(4)
    jacob(1,3) = nr(1)*zz(1)+nr(2)*zz(2)+nr(3)*zz(3)+nr(4)*zz(4)
    jacob(2,3) = ns(1)*zz(1)+ns(2)*zz(2)+ns(3)*zz(3)+ns(4)*zz(4)
    jacob(3,3) = nt(1)*zz(1)+nt(2)*zz(2)+nt(3)*zz(3)+nt(4)*zz(4)
 
    det = jacob(1,1) * jacob(2,2) * jacob(3,3)   &
      &         + jacob(1,2) * jacob(2,3) * jacob(3,1)   &
      &         + jacob(1,3) * jacob(2,1) * jacob(3,2)   &
      &         - jacob(1,3) * jacob(2,2) * jacob(3,1)   &
      &         - jacob(1,2) * jacob(2,1) * jacob(3,3)   &
      &         - jacob(1,1) * jacob(2,3) * jacob(3,2)
    !!
    !! ******** Inverse Jacobi matrix calculation ********
    !!
    dum = 1.0 / det
    jinv(1,1) = dum*(  jacob(2,2)*jacob(3,3) -jacob(2,3)*jacob(3,2) )
    jinv(2,1) = dum*( -jacob(2,1)*jacob(3,3) +jacob(2,3)*jacob(3,1) )
    jinv(3,1) = dum*(  jacob(2,1)*jacob(3,2) -jacob(2,2)*jacob(3,1) )
    jinv(1,2) = dum*( -jacob(1,2)*jacob(3,3) +jacob(1,3)*jacob(3,2) )
    jinv(2,2) = dum*(  jacob(1,1)*jacob(3,3) -jacob(1,3)*jacob(3,1) )
    jinv(3,2) = dum*( -jacob(1,1)*jacob(3,2) +jacob(1,2)*jacob(3,1) )
    jinv(1,3) = dum*(  jacob(1,2)*jacob(2,3) -jacob(1,3)*jacob(2,2) )
    jinv(2,3) = dum*( -jacob(1,1)*jacob(2,3) +jacob(1,3)*jacob(2,1) )
    jinv(3,3) = dum*(  jacob(1,1)*jacob(2,2) -jacob(1,2)*jacob(2,1) )
 
    do i = 1, 4
      do j = 1, 4
 
        nx(i,j) = jinv(1,1)*nr(i) + jinv(1,2)*ns(i) +jinv(1,3)*nt(i)
        ny(i,j) = jinv(2,1)*nr(i) + jinv(2,2)*ns(i) +jinv(2,3)*nt(i)
        nz(i,j) = jinv(3,1)*nr(i) + jinv(3,2)*ns(i) +jinv(3,3)*nt(i)
 
      end do
    end do
 
  end subroutine hecmw_jacob_341
 
end module hecmw_jacob341