xref: /phasta/phSolver/common/asithf.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)

      subroutine asithf (y, x, strnrm, ien, fres, shgl, shp, Qwtf)

      include "common.h"

      dimension y(nshg,ndof),            fres(nshg,24)
      dimension x(numnp,nsd),            xl(npro,nenl,nsd)
      dimension ien(npro,nshl),        ycl(npro,nshl,ndof),
     &          fresl(npro,24),        WdetJ(npro),
     &          u1(npro),              u2(npro),
     &          u3(npro),              dxdxi(npro,nsd,nsd),
     &          strnrm(npro,maxsh),    dxidx(npro,nsd,nsd),
     &          shgl(nsd,nshl,maxsh),       shg(npro,nshl,nsd),
     &          shp(nshl,maxsh),
     &          fresli(npro,24),       Qwtf(ngaussf)

      dimension tmp(npro)

      call localy (y,      ycl,     ien,    5,  'gather  ')
      call localx (x,      xl,     ien,    3,  'gather  ')
c

      if(matflg(1,1).eq.0) then ! compressible
         ycl (:,:,1) = ycl(:,:,1) / (Rgas * ycl(:,:,5)) !get density
      else
         ycl(:,:,1) = one ! Even if density non unity, it would cancel out
      endif

      fresl = zero

      do intp = 1, ngaussf


c  calculate the metrics
c
c
c.... --------------------->  Element Metrics  <-----------------------
c
c.... compute the deformation gradient
c
        dxdxi = zero
c
          do n = 1, nenl
            dxdxi(:,1,1) = dxdxi(:,1,1) + xl(:,n,1) * shgl(1,n,intp)
            dxdxi(:,1,2) = dxdxi(:,1,2) + xl(:,n,1) * shgl(2,n,intp)
            dxdxi(:,1,3) = dxdxi(:,1,3) + xl(:,n,1) * shgl(3,n,intp)
            dxdxi(:,2,1) = dxdxi(:,2,1) + xl(:,n,2) * shgl(1,n,intp)
            dxdxi(:,2,2) = dxdxi(:,2,2) + xl(:,n,2) * shgl(2,n,intp)
            dxdxi(:,2,3) = dxdxi(:,2,3) + xl(:,n,2) * shgl(3,n,intp)
            dxdxi(:,3,1) = dxdxi(:,3,1) + xl(:,n,3) * shgl(1,n,intp)
            dxdxi(:,3,2) = dxdxi(:,3,2) + xl(:,n,3) * shgl(2,n,intp)
            dxdxi(:,3,3) = dxdxi(:,3,3) + xl(:,n,3) * shgl(3,n,intp)
          enddo
c
c.... compute the inverse of deformation gradient
c
        dxidx(:,1,1) =   dxdxi(:,2,2) * dxdxi(:,3,3)
     &                 - dxdxi(:,3,2) * dxdxi(:,2,3)
        dxidx(:,1,2) =   dxdxi(:,3,2) * dxdxi(:,1,3)
     &                 - dxdxi(:,1,2) * dxdxi(:,3,3)
        dxidx(:,1,3) =   dxdxi(:,1,2) * dxdxi(:,2,3)
     &                 - dxdxi(:,1,3) * dxdxi(:,2,2)
        tmp          = one / ( dxidx(:,1,1) * dxdxi(:,1,1)
     &                       + dxidx(:,1,2) * dxdxi(:,2,1)
     &                       + dxidx(:,1,3) * dxdxi(:,3,1) )
        dxidx(:,1,1) = dxidx(:,1,1) * tmp
        dxidx(:,1,2) = dxidx(:,1,2) * tmp
        dxidx(:,1,3) = dxidx(:,1,3) * tmp
        dxidx(:,2,1) = (dxdxi(:,2,3) * dxdxi(:,3,1)
     &                - dxdxi(:,2,1) * dxdxi(:,3,3)) * tmp
        dxidx(:,2,2) = (dxdxi(:,1,1) * dxdxi(:,3,3)
     &                - dxdxi(:,3,1) * dxdxi(:,1,3)) * tmp
        dxidx(:,2,3) = (dxdxi(:,2,1) * dxdxi(:,1,3)
     &                - dxdxi(:,1,1) * dxdxi(:,2,3)) * tmp
        dxidx(:,3,1) = (dxdxi(:,2,1) * dxdxi(:,3,2)
     &                - dxdxi(:,2,2) * dxdxi(:,3,1)) * tmp
        dxidx(:,3,2) = (dxdxi(:,3,1) * dxdxi(:,1,2)
     &                - dxdxi(:,1,1) * dxdxi(:,3,2)) * tmp
        dxidx(:,3,3) = (dxdxi(:,1,1) * dxdxi(:,2,2)
     &                - dxdxi(:,1,2) * dxdxi(:,2,1)) * tmp
c
c        wght=Qwt(lcsyst,intp)  ! may be different now
        wght=Qwtf(intp)
        WdetJ = wght / tmp
c
      fresli=zero
c
      if(matflg(1,1).eq.0) then ! compressible
         do i=1,nshl
            fresli(:,22) = fresli(:,22)+shp(i,intp)*ycl(:,i,1) !density at
                                !qpt
         enddo
      else   ! incompressible, set density
         fresli(:,22)= one ! reduce comp2incompr regardless of rho  datmat(1,1,1)
      endif
c
      do n = 1,nshl
        shg(:,n,1) = (shgl(1,n,intp) * dxidx(:,1,1)
     &              + shgl(2,n,intp) * dxidx(:,2,1)
     &              + shgl(3,n,intp) * dxidx(:,3,1))
        shg(:,n,2) = (shgl(1,n,intp) * dxidx(:,1,2)
     &              + shgl(2,n,intp) * dxidx(:,2,2)
     &              + shgl(3,n,intp) * dxidx(:,3,2))
        shg(:,n,3) = (shgl(1,n,intp) * dxidx(:,1,3)
     &              + shgl(2,n,intp) * dxidx(:,2,3)
     &              + shgl(3,n,intp) * dxidx(:,3,3))
      enddo

      do j=10,12  ! normal strainrate u_{i,i} no sum on i
       ig=j-9
       iv=j-8
       do i=1,nshl
        fresli(:,j) = fresli(:,j)+shg(:,i,ig)*ycl(:,i,iv)
       enddo
      enddo

c shear stresses  NOTE  there may be faster ways to do this
c                  check agains CM5 code for speed WTP

       do i=1,nshl
        fresli(:,13) = fresli(:,13)+shg(:,i,2)*ycl(:,i,2)
     &                             +shg(:,i,1)*ycl(:,i,3)
        fresli(:,14) = fresli(:,14)+shg(:,i,3)*ycl(:,i,2)
     &                             +shg(:,i,1)*ycl(:,i,4)
        fresli(:,15) = fresli(:,15)+shg(:,i,3)*ycl(:,i,3)
     &                             +shg(:,i,2)*ycl(:,i,4)
       enddo

      fresli(:,13) = pt5 * fresli(:,13)
      fresli(:,14) = pt5 * fresli(:,14)
      fresli(:,15) = pt5 * fresli(:,15)

      strnrm(:,intp) = fresli(:,22) * sqrt(
     &   two * (fresli(:,10)**2 + fresli(:,11)**2 + fresli(:,12)**2)
     &  + four * ( fresli(:,13)**2 + fresli(:,14)**2 +
     &    fresli(:,15)**2 ) )

c
c S_ij
c

      fresli(:,10) = fresli(:,10) * WdetJ ! u_{1,1}*WdetJ
      fresli(:,11) = fresli(:,11) * WdetJ ! u_{2,2}*WdetJ
      fresli(:,12) = fresli(:,12) * WdetJ ! u_{3,3}*WdetJ
      fresli(:,13) = fresli(:,13) * WdetJ ! (1/2)*(u_{1,2}+u_{2,1})*WdetJ
      fresli(:,14) = fresli(:,14) * WdetJ ! (1/2)*(u_{1,3}+u_{3,1})*WdetJ
      fresli(:,15) = fresli(:,15) * WdetJ ! (1/2)*(u_{2,3}+u_{3,2})*WdetJ

      fresli(:,22) = fresli(:,22) * WdetJ   !rho * WdetJ
c     fresli(:,24) = fresli(:,24) * WdetJ

      u1=zero
      u2=zero
      u3=zero
      do i=1,nshl
       u1 = u1 + shp(i,intp)*ycl(:,i,2)
       u2 = u2 + shp(i,intp)*ycl(:,i,3)
       u3 = u3 + shp(i,intp)*ycl(:,i,4)
      enddo

      fresli(:,1) = fresli(:,22) * u1   !rho u1 * WdetJ
      fresli(:,2) = fresli(:,22) * u2   !rho u2 * WdetJ
      fresli(:,3) = fresli(:,22) * u3   !rho u3 * WdetJ

      fresli(:,4) = fresli(:,1) * u1    !rho u1 u1 *WdetJ
      fresli(:,5) = fresli(:,2) * u2    !rho u2 u2 *WdetJ
      fresli(:,6) = fresli(:,3) * u3    !rho u3 u3 *WdetJ
      fresli(:,7) = fresli(:,1) * u2    !rho u1 u2 *WdetJ
      fresli(:,8) = fresli(:,1) * u3    !rho u1 u3 *WdetJ
      fresli(:,9) = fresli(:,2) * u3    !rho u2 u3 *WdetJ

      fresli(:,16) = strnrm(:,intp) * fresli(:,10) ! rho *|Eps| *Eps11 *WdetJ
      fresli(:,17) = strnrm(:,intp) * fresli(:,11) ! rho *|Eps| *Eps22 *WdetJ
      fresli(:,18) = strnrm(:,intp) * fresli(:,12) ! rho *|Eps| *Eps33 *WdetJ
      fresli(:,19) = strnrm(:,intp) * fresli(:,13) ! rho *|Eps| *Eps12 *WdetJ
      fresli(:,20) = strnrm(:,intp) * fresli(:,14) ! rho *|Eps| *Eps13 *WdetJ
      fresli(:,21) = strnrm(:,intp) * fresli(:,15) ! rho *|Eps| *Eps23 *WdetJ

      fresli(:,23) = WdetJ   !    Integral of 1 over the element
c
      do i = 1, 23
         fresl(:,i) = fresl(:,i) + fresli(:,i)
      enddo

      enddo !end of loop over integration points
c
      do j = 1,nshl
      do nel = 1,npro
        fres(ien(nel,j),:) = fres(ien(nel,j),:) + fresl(nel,:)
      enddo
      enddo

      return
      end