hecmw_Jacob_241.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_jacob241
contains
 
  subroutine hecmw_jacob_241 ( hecMESH, iElem, DET, W, N, NX, NY )
    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)
 
    integer(kind=kint):: i, j, jj, ilocal
    integer(kind=kint):: lx, ly
    real(kind=kreal)::   dum
    real(kind=kreal)::   xx(4), yy(4)
    real(kind=kreal)::   xg(2), hr(2,2,4), hs(2,2,4), ht(2,2,4)
    real(kind=kreal)::   h(2,2,4),bx(2,2,4),by(2,2,4),bz(2,2,4)
    real(kind=kreal)::   ri, si, rp, sp, rm, sm
    real(kind=kreal)::   xj11,xj12,xj21,xj22
    real(kind=kreal)::   xji11,xji12,xji21,xji22
 
    data xg/-0.5773502691896d0, 0.5773502691896d0/
 
    do j = 1, 4
      jj = 4 - j
      ilocal = hecmesh%elem_node_item  ( 4*ielem -jj )
      xx(j)  = hecmesh%node( ilocal*2 -1 )
      yy(j)  = hecmesh%node( ilocal*2    )
      w(j)  = 1.0d0
    end do
    !C
    !C*LOOP OVER ALL INTEGRATION POINTS
    do lx=1,2
      ri=xg(lx)
      do ly=1,2
        si=xg(ly)
        rp=1.0+ri
        sp=1.0+si
        rm=1.0-ri
        sm=1.0-si
        !!
        !! ******** Shape functions ********
        !!
        h(lx,ly,1)=0.25*rm*sm
        h(lx,ly,2)=0.25*rp*sm
        h(lx,ly,3)=0.25*rp*sp
        h(lx,ly,4)=0.25*rm*sp
        !!
        !! ******** Derivative of shape functions ********
        !!
        !! ----------- For R-Coordinate -------------
        hr(lx,ly,1)=-.25*sm
        hr(lx,ly,2)= .25*sm
        hr(lx,ly,3)= .25*sp
        hr(lx,ly,4)=-.25*sp
        !! ----------- For S-Coordinate -------------
        hs(lx,ly,1)=-.25*rm
        hs(lx,ly,2)=-.25*rp
        hs(lx,ly,3)= .25*rp
        hs(lx,ly,4)= .25*rm
        !!
        !! ******** Jacobi matrix calculation********
        !!
        xj11=0.0
        xj21=0.0
        xj12=0.0
        xj22=0.0
        do i=1,4
          xj11=xj11+hr(lx,ly,i)*xx(i)
          xj21=xj21+hs(lx,ly,i)*xx(i)
          xj12=xj12+hr(lx,ly,i)*yy(i)
          xj22=xj22+hs(lx,ly,i)*yy(i)
        end do
        det=xj11*xj22-xj21*xj12
        !!
        !! ******** Inverse Jacobi matrix calculation ********
        !!
        dum=1.0/det
        xji11= xj22*dum
        xji12=-xj12*dum
        xji21=-xj21*dum
        xji22= xj11*dum
        do j=1, 4
          bx(lx,ly,j)=xji11*hr(lx,ly,j)+xji12*hs(lx,ly,j)
          by(lx,ly,j)=xji21*hr(lx,ly,j)+xji22*hs(lx,ly,j)
        end do
        !C
      end do
    end do
    !C
 
    j = 1
    do lx = 1, 2
      do ly = 1, 2
 
        n(j,1)  = h(lx,ly,1)
        n(j,2)  = h(lx,ly,2)
        n(j,3)  = h(lx,ly,3)
        n(j,4)  = h(lx,ly,4)
 
        nx(j,1) = bx(lx,ly,1)
        nx(j,2) = bx(lx,ly,2)
        nx(j,3) = bx(lx,ly,3)
        nx(j,4) = bx(lx,ly,4)
 
        ny(j,1) = by(lx,ly,1)
        ny(j,2) = by(lx,ly,2)
        ny(j,3) = by(lx,ly,3)
        ny(j,4) = by(lx,ly,4)
 
        j = j + 1
      end do
    end do
 
    !C
    !C
  end subroutine hecmw_jacob_241
 
end module hecmw_jacob241