static_LIB_Fbar.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 m_static_lib_fbar
 
  use hecmw, only : kint, kreal
  use elementinfo
 
  implicit none
 
  real(kind=kreal), private, parameter :: i33(3,3) = &
    &  reshape( (/1.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 1.d0/), (/3,3/) )
 
contains
 
 
  !----------------------------------------------------------------------*
  subroutine stf_c3d8fbar &
      (etype, nn, ecoord, gausses, stiff, cdsys_id, coords, &
      time, tincr, u, temperature)
    !----------------------------------------------------------------------*
 
    use mmechgauss
    use m_matmatrix
    use m_common_struct
    use m_static_lib_3d, only: geomat_c3
    use m_utilities
 
    !---------------------------------------------------------------------
 
    integer(kind=kint), intent(in)  :: etype
    integer(kind=kint), intent(in)  :: nn
    real(kind=kreal), intent(in)    :: ecoord(3, nn)          
    type(tgaussstatus), intent(in)  :: gausses(:)
    real(kind=kreal), intent(out)   :: stiff(:,:)             
    integer(kind=kint), intent(in)  :: cdsys_id
    real(kind=kreal), intent(inout) :: coords(3, 3)           
    real(kind=kreal), intent(in)    :: time                   
    real(kind=kreal), intent(in)    :: tincr                  
    real(kind=kreal), intent(in), optional :: u(:, :)         
    real(kind=kreal), intent(in)    :: temperature(nn)        
 
    !---------------------------------------------------------------------
 
    integer(kind=kint) :: flag
    integer(kind=kint), parameter :: ndof = 3
    real(kind=kreal) :: d(6, 6),b(6, ndof*nn),db(6, ndof*nn)
    real(kind=kreal) :: gderiv(nn, 3),stress(6),mat(6, 6)
    real(kind=kreal) :: det, wg, temp, spfunc(nn)
    integer(kind=kint) :: i, j, p, q, lx, serr
    real(kind=kreal) :: naturalcoord(3)
    real(kind=kreal) :: gdispderiv(3, 3)
    real(kind=kreal) :: b1(6, ndof*nn)
    real(kind=kreal) :: smat(9, 9), elem(3, nn)
    real(kind=kreal) :: bn(9, ndof*nn), sbn(9, ndof*nn)
    real(kind=kreal) :: coordsys(3, 3)
 
    real(kind=kreal) :: elem0(3,nn), elem1(3, nn), gderiv1(nn,ndof), b2(6, ndof*nn), z1(3,2)
    real(kind=kreal) :: v0, jacob, jacob_ave, gderiv1_ave(nn,ndof)
    real(kind=kreal) :: gderiv2_ave(ndof*nn,ndof*nn)
    real(kind=kreal) :: fbar(3,3), jratio(8), coeff, sff, dstrain(6), ddfs, fs(3,3), gfs(3,2)
 
    !---------------------------------------------------------------------
 
    stiff(:, :) = 0.0d0
    ! we suppose the same material type in the element
    flag = gausses(1)%pMaterial%nlgeom_flag
    if( .not. present(u) ) flag = infinitesimal    ! enforce to infinitesimal deformation analysis
    elem(:, :) = ecoord(:, :)
    elem0(:, :) = ecoord(:, :)
    if( flag == updatelag ) elem(:, :) = ecoord(:, :)+u(:, :)
    elem1(:, :) = ecoord(:, :)+u(:, :)
 
    !cal volumetric average of J=detF and dN/dx
    v0 = 0.d0
    jacob_ave = 0.d0
    gderiv1_ave(1:nn,1:ndof) = 0.d0
    gderiv2_ave(1:ndof*nn,1:ndof*nn) = 0.d0
    do lx = 1, numofquadpoints(etype)
      call getquadpoint( etype, lx, naturalcoord(:) )
      call getglobalderiv( etype, nn, naturalcoord, elem0, det, gderiv)
      wg = getweight(etype, lx)*det
      if( flag == infinitesimal ) then
        jacob = 1.d0
        gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + jacob*wg*gderiv(1:nn, 1:ndof)
      else
        gdispderiv(1:3, 1:3) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
        jacob = determinant33( i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof) )
        jratio(lx) = dsign(dabs(jacob)**(-1.d0/3.d0),jacob)
        call getglobalderiv( etype, nn, naturalcoord, elem1, det, gderiv1)
        gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + jacob*wg*gderiv1(1:nn, 1:ndof)
        do p=1,nn
          do q=1,nn
            do i=1,ndof
              do j=1,ndof
                gderiv2_ave(3*p-3+i,3*q-3+j) = gderiv2_ave(3*p-3+i,3*q-3+j)  &
                  & + jacob*wg*(gderiv1(p,i)*gderiv1(q,j)-gderiv1(q,i)*gderiv1(p,j))
              end do
            end do
          end do
        end do
      endif
      v0 = v0 + wg
      jacob_ave = jacob_ave + jacob*wg
    enddo
    if( dabs(v0) > 1.d-12 ) then
      if( dabs(jacob_ave) < 1.d-12 ) stop 'Error in Update_C3D8Fbar: too small average jacob'
      jacob_ave = jacob_ave/v0
      jratio(1:8) = (dsign(dabs(jacob_ave)**(1.d0/3.d0),jacob_ave))*jratio(1:8) !Jratio(1:8) = (jacob_ave**(1.d0/3.d0))*Jratio(1:8)
      gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof)/(v0*jacob_ave)
      gderiv2_ave(1:ndof*nn,1:ndof*nn) = gderiv2_ave(1:ndof*nn,1:ndof*nn)/(v0*jacob_ave)
    else
      stop 'Error in Update_C3D8Fbar: too small element volume'
    end if
 
    do lx = 1, numofquadpoints(etype)
 
      call getquadpoint( etype, lx, naturalcoord(:) )
      call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
 
      if( cdsys_id > 0 ) then
        call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
        if( serr == -1 ) stop "Fail to setup local coordinate"
        if( serr == -2 ) then
          write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
        end if
      end if
 
      call getshapefunc( etype, naturalcoord, spfunc )
      temp = dot_product( temperature, spfunc )
      call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
 
      if( flag == updatelag ) then
        call geomat_c3( gausses(lx)%stress, mat )
        d(:, :) = d(:, :)-mat
      endif
 
      wg = getweight(etype, lx)*det
      b(1:6, 1:nn*ndof) = 0.0d0
      do j = 1, nn
        b(1, 3*j-2) = gderiv(j, 1)
        b(2, 3*j-1) = gderiv(j, 2)
        b(3, 3*j  ) = gderiv(j, 3)
        b(4, 3*j-2) = gderiv(j, 2)
        b(4, 3*j-1) = gderiv(j, 1)
        b(5, 3*j-1) = gderiv(j, 3)
        b(5, 3*j  ) = gderiv(j, 2)
        b(6, 3*j-2) = gderiv(j, 3)
        b(6, 3*j  ) = gderiv(j, 1)
      end do
 
      ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
      if( flag == infinitesimal ) then
        b2(1:6, 1:nn*ndof) = 0.0d0
        do j = 1, nn
          z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
          b2(1,3*j-2:3*j) = z1(1:3,1)
          b2(2,3*j-2:3*j) = z1(1:3,1)
          b2(3,3*j-2:3*j) = z1(1:3,1)
        end do
 
        ! BL = Jratio*(BL0 + BL1)+BL2
        do j = 1, nn*ndof
          b(1:3, j) = b(1:3,j)+b2(1:3,j)
        end do
 
      else if( flag == totallag ) then
        gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
        fbar(1:ndof, 1:ndof) = jratio(lx)*(i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof))
        call getglobalderiv( etype, nn, naturalcoord, elem1, det, gderiv1)
 
        ! ---dudx(i,j) ==> gdispderiv(i,j)
        b1(1:6, 1:nn*ndof) = 0.0d0
        do j = 1, nn
          b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
          b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
          b1(1, 3*j  ) = gdispderiv(3, 1)*gderiv(j, 1)
          b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
          b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
          b1(2, 3*j  ) = gdispderiv(3, 2)*gderiv(j, 2)
          b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
          b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
          b1(3, 3*j  ) = gdispderiv(3, 3)*gderiv(j, 3)
          b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
          b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
          b1(4, 3*j  ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
          b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
          b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
          b1(5, 3*j  ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
          b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
          b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
          b1(6, 3*j  ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
        end do
 
        dstrain(1) = 0.5d0*(dot_product(fbar(1:3,1),fbar(1:3,1))-1.d0)
        dstrain(2) = 0.5d0*(dot_product(fbar(1:3,2),fbar(1:3,2))-1.d0)
        dstrain(3) = 0.5d0*(dot_product(fbar(1:3,3),fbar(1:3,3))-1.d0)
        dstrain(4) = dot_product(fbar(1:3,1),fbar(1:3,2))
        dstrain(5) = dot_product(fbar(1:3,2),fbar(1:3,3))
        dstrain(6) = dot_product(fbar(1:3,3),fbar(1:3,1))
 
        b2(1:6, 1:nn*ndof) = 0.0d0
        do j = 1, nn
          z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
          b2(1,3*j-2:3*j) = (2.d0*dstrain(1)+1.d0)*z1(1:3,1)
          b2(2,3*j-2:3*j) = (2.d0*dstrain(2)+1.d0)*z1(1:3,1)
          b2(3,3*j-2:3*j) = (2.d0*dstrain(3)+1.d0)*z1(1:3,1)
          b2(4,3*j-2:3*j) = 2.d0*dstrain(4)*z1(1:3,1)
          b2(5,3*j-2:3*j) = 2.d0*dstrain(5)*z1(1:3,1)
          b2(6,3*j-2:3*j) = 2.d0*dstrain(6)*z1(1:3,1)
        end do
 
        ! BL = Jratio*(BL0 + BL1)+BL2
        do j = 1, nn*ndof
          b(:, j) = jratio(lx)*jratio(lx)*(b(:,j)+b1(:,j))+b2(:,j)
        end do
 
      else if( flag == updatelag ) then
        wg = (jratio(lx)**3.d0)*getweight(etype, lx)*det
 
        b2(1:3, 1:nn*ndof) = 0.0d0
        do j = 1, nn
          z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
          b2(1, 3*j-2:3*j) = z1(1:3,1)
          b2(2, 3*j-2:3*j) = z1(1:3,1)
          b2(3, 3*j-2:3*j) = z1(1:3,1)
        end do
 
        do j = 1, nn*ndof
          b(1:3, j) = b(1:3,j)+b2(1:3,j)
        end do
 
      end if
 
      db(1:6, 1:nn*ndof) = matmul( d, b(1:6, 1:nn*ndof) )
      do j=1,nn*ndof 
        do i=1,nn*ndof
          stiff(i, j) = stiff(i, j)+dot_product( b(:, i), db(:, j) )*wg
        end do
      end do
 
      ! calculate the initial stress matrix(1): dFbar*dFbar*Stress
      if( flag == totallag .or. flag == updatelag ) then
        stress(1:6) = gausses(lx)%stress
        if( flag == totallag ) then
          coeff = jratio(lx)
          sff = dot_product(stress(1:6),dstrain(1:6))
        else if( flag == updatelag ) then
          coeff = 1.d0
          sff = stress(1)+stress(2)+stress(3)
          gderiv1 = gderiv
          fbar(1:3,1:3) = i33(1:3,1:3)
        end if
 
        bn(1:9, 1:nn*ndof) = 0.0d0
        do j = 1, nn
          bn(1, 3*j-2) = coeff*gderiv(j, 1)
          bn(2, 3*j-1) = coeff*gderiv(j, 1)
          bn(3, 3*j  ) = coeff*gderiv(j, 1)
          bn(4, 3*j-2) = coeff*gderiv(j, 2)
          bn(5, 3*j-1) = coeff*gderiv(j, 2)
          bn(6, 3*j  ) = coeff*gderiv(j, 2)
          bn(7, 3*j-2) = coeff*gderiv(j, 3)
          bn(8, 3*j-1) = coeff*gderiv(j, 3)
          bn(9, 3*j  ) = coeff*gderiv(j, 3)
          z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
          bn(1, 3*j-2:3*j) = bn(1, 3*j-2:3*j) + fbar(1,1)*z1(1:3,1)
          bn(2, 3*j-2:3*j) = bn(2, 3*j-2:3*j) + fbar(2,1)*z1(1:3,1)
          bn(3, 3*j-2:3*j) = bn(3, 3*j-2:3*j) + fbar(3,1)*z1(1:3,1)
          bn(4, 3*j-2:3*j) = bn(4, 3*j-2:3*j) + fbar(1,2)*z1(1:3,1)
          bn(5, 3*j-2:3*j) = bn(5, 3*j-2:3*j) + fbar(2,2)*z1(1:3,1)
          bn(6, 3*j-2:3*j) = bn(6, 3*j-2:3*j) + fbar(3,2)*z1(1:3,1)
          bn(7, 3*j-2:3*j) = bn(7, 3*j-2:3*j) + fbar(1,3)*z1(1:3,1)
          bn(8, 3*j-2:3*j) = bn(8, 3*j-2:3*j) + fbar(2,3)*z1(1:3,1)
          bn(9, 3*j-2:3*j) = bn(9, 3*j-2:3*j) + fbar(3,3)*z1(1:3,1)
        end do
        smat(:, :) = 0.0d0
        do j= 1, 3
          smat(j  , j  ) = stress(1)
          smat(j  , j+3) = stress(4)
          smat(j  , j+6) = stress(6)
          smat(j+3, j  ) = stress(4)
          smat(j+3, j+3) = stress(2)
          smat(j+3, j+6) = stress(5)
          smat(j+6, j  ) = stress(6)
          smat(j+6, j+3) = stress(5)
          smat(j+6, j+6) = stress(3)
        end do
        sbn(1:9, 1:nn*ndof) = matmul( smat(1:9, 1:9), bn(1:9, 1:nn*ndof) )
        do j=1,nn*ndof 
          do i=1,nn*ndof
            stiff(i, j) = stiff(i, j)+dot_product( bn(:, i), sbn(:, j) )*wg
          end do
        end do
 
        ! calculate the initial stress matrix(2): d(dFbar)*Stress
        fs(1,1) = fbar(1,1)*stress(1)+fbar(1,2)*stress(4)+fbar(1,3)*stress(6)
        fs(1,2) = fbar(1,1)*stress(4)+fbar(1,2)*stress(2)+fbar(1,3)*stress(5)
        fs(1,3) = fbar(1,1)*stress(6)+fbar(1,2)*stress(5)+fbar(1,3)*stress(3)
        fs(2,1) = fbar(2,1)*stress(1)+fbar(2,2)*stress(4)+fbar(2,3)*stress(6)
        fs(2,2) = fbar(2,1)*stress(4)+fbar(2,2)*stress(2)+fbar(2,3)*stress(5)
        fs(2,3) = fbar(2,1)*stress(6)+fbar(2,2)*stress(5)+fbar(2,3)*stress(3)
        fs(3,1) = fbar(3,1)*stress(1)+fbar(3,2)*stress(4)+fbar(3,3)*stress(6)
        fs(3,2) = fbar(3,1)*stress(4)+fbar(3,2)*stress(2)+fbar(3,3)*stress(5)
        fs(3,3) = fbar(3,1)*stress(6)+fbar(3,2)*stress(5)+fbar(3,3)*stress(3)
        do i=1,nn
          z1(1:3,1) = (gderiv1_ave(i,1:3)-gderiv1(i,1:3))/3.d0
          gfs(1:3,1) = coeff*matmul(fs,gderiv(i,1:3))
          do j=1,nn
            z1(1:3,2) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
            gfs(1:3,2) = coeff*matmul(fs,gderiv(j,1:3))
            do p=1,ndof
              do q=1,ndof
                ddfs = z1(p,1)*z1(q,2)
                ddfs = ddfs + (gderiv2_ave(3*i-3+p,3*j-3+q)-gderiv1_ave(i,p)*gderiv1_ave(j,q))/3.d0
                ddfs = ddfs + gderiv1(i,q)*gderiv1(j,p)/3.d0
                ddfs = sff*ddfs + z1(p,1)*gfs(q,2)+z1(q,2)*gfs(p,1)
                stiff(3*i-3+p, 3*j-3+q) = stiff(3*i-3+p, 3*j-3+q) + ddfs*wg
              end do
            end do
          end do
        end do
      end if
 
 
    end do ! gauss roop
 
  end subroutine stf_c3d8fbar
 
 
  !----------------------------------------------------------------------*
  subroutine update_c3d8fbar                              &
      (etype, nn, ecoord, u, du, cdsys_id, coords, &
      qf, gausses, iter, time, tincr, tt, t0, tn  )
    !----------------------------------------------------------------------*
 
    use m_fstr
    use mmaterial
    use mmechgauss
    use m_matmatrix
    use m_elastoplastic
    use mhyperelastic
    use m_utilities
    use m_static_lib_3d
 
    !---------------------------------------------------------------------
 
    integer(kind=kint), intent(in)    :: etype
    integer(kind=kint), intent(in)    :: nn
    real(kind=kreal), intent(in)      :: ecoord(3, nn) 
    real(kind=kreal), intent(in)      :: u(3, nn)      
    real(kind=kreal), intent(in)      :: du(3, nn)     
    integer(kind=kint), intent(in)    :: cdsys_id
    real(kind=kreal), intent(inout)   :: coords(3, 3)  
    real(kind=kreal), intent(out)     :: qf(nn*3)      
    type(tgaussstatus), intent(inout) :: gausses(:)
    integer, intent(in) :: iter
    real(kind=kreal), intent(in)      :: time          
    real(kind=kreal), intent(in)      :: tincr         
    real(kind=kreal), intent(in)      :: tt(nn)        
    real(kind=kreal), intent(in)      :: t0(nn)        
    real(kind=kreal), intent(in)      :: tn(nn)        
 
    !---------------------------------------------------------------------
 
    integer(kind=kint) :: flag
    integer(kind=kint), parameter :: ndof = 3
    real(kind=kreal) :: b(6, ndof*nn), b1(6, ndof*nn)
    real(kind=kreal) :: gderiv(nn, 3), gdispderiv(3, 3), det, wg
    integer(kind=kint) :: j, lx, serr
    real(kind=kreal) :: naturalcoord(3), rot(3, 3), spfunc(nn)
    real(kind=kreal) :: totaldisp(3, nn), elem(3, nn), elem1(3, nn), coordsys(3, 3)
    real(kind=kreal) :: dstrain(6)
    real(kind=kreal) :: dvol
    real(kind=kreal) :: ttc, tt0, ttn, alpo(3), ina(1), epsth(6)
    logical :: ierr, matlaniso
 
    real(kind=kreal) :: elem0(3,nn), gderiv1(nn,ndof), b2(6, ndof*nn), z1(3)
    real(kind=kreal) :: v0, jacob, jacob_ave, gderiv1_ave(nn,ndof)
    real(kind=kreal) :: fbar(3,3), jratio(8)
    real(kind=kreal) :: jacob_ave05, gderiv05_ave(nn,ndof)
 
    !---------------------------------------------------------------------
 
    qf(:) = 0.0d0
    ! we suppose the same material type in the element
    flag = gausses(1)%pMaterial%nlgeom_flag
    elem0(1:ndof,1:nn) = ecoord(1:ndof,1:nn)
    totaldisp(:, :) = u(:, :)+du(:, :)
    if( flag == infinitesimal ) then
      elem(:, :) = ecoord(:, :)
      elem1(:, :) = ecoord(:, :)
    else if( flag == totallag ) then
      elem(:, :) = ecoord(:, :)
      elem1(:, :) = ( du(:, :)+u(:, :) )+ecoord(:, :)
    else if( flag == updatelag ) then
      elem(:, :) = ( 0.5d0*du(:, :)+u(:, :) )+ecoord(:, :)
      elem1(:, :) = ( du(:, :)+u(:, :) )+ecoord(:, :)
      totaldisp(:, :) = du(:, :)
    end if
 
    matlaniso = .false.
    ina = tt(1)
    call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
    if( .not. ierr ) matlaniso = .true.
 
    !cal volumetric average of J=detF and dN/dx
    v0 = 0.d0
    jacob_ave = 0.d0
    gderiv1_ave(1:nn,1:ndof) = 0.d0
    if( flag == updatelag ) then
      jacob_ave05 = 0.d0
      gderiv05_ave(1:nn,1:ndof) = 0.d0
    endif
    do lx = 1, numofquadpoints(etype)
      call getquadpoint( etype, lx, naturalcoord(:) )
      call getglobalderiv( etype, nn, naturalcoord, elem0, det, gderiv)
      wg = getweight(etype, lx)*det
      if( flag == infinitesimal ) then
        jacob = 1.d0
        gderiv1(1:nn, 1:ndof) = gderiv(1:nn, 1:ndof)
      else
        gdispderiv(1:3, 1:3) = matmul( du(1:ndof, 1:nn)+u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
        jacob = determinant33( i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof) )
        jratio(lx) = dsign(dabs(jacob)**(-1.d0/3.d0),jacob) ! Jratio(LX) = jacob**(-1.d0/3.d0)
 
        call getglobalderiv( etype, nn, naturalcoord, elem1, det, gderiv1)
      endif
      v0 = v0 + wg
      jacob_ave = jacob_ave + jacob*wg
      gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + jacob*wg*gderiv1(1:nn, 1:ndof)
      if( flag == updatelag ) then
        call getglobalderiv( etype, nn, naturalcoord, elem, det, gderiv1)
        wg = getweight(etype, lx)*det
        jacob_ave05 = jacob_ave05 + wg
        gderiv05_ave(1:nn,1:ndof) = gderiv05_ave(1:nn,1:ndof) + wg*gderiv1(1:nn, 1:ndof)
      endif
    enddo
    if( dabs(v0) > 1.d-12 ) then
      if( dabs(jacob_ave) < 1.d-12 ) stop 'Error in Update_C3D8Fbar: too small average jacob'
      jacob_ave = jacob_ave/v0
      jratio(1:8) = (dsign(dabs(jacob_ave)**(1.d0/3.d0),jacob_ave))*jratio(1:8) !Jratio(1:8) = (jacob_ave**(1.d0/3.d0))*Jratio(1:8)
      gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof)/(v0*jacob_ave)
      if( flag == updatelag ) gderiv05_ave(1:nn,1:ndof) = gderiv05_ave(1:nn,1:ndof)/jacob_ave05
    else
      stop 'Error in Update_C3D8Fbar: too small element volume'
    end if
 
    do lx = 1, numofquadpoints(etype)
 
      call getquadpoint( etype, lx, naturalcoord(:) )
      call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
 
      if( cdsys_id > 0 ) then
        call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
        if( serr == -1 ) stop "Fail to setup local coordinate"
        if( serr == -2 ) then
          write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
        end if
      end if
 
      gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
 
      ! ========================================================
      !     UPDATE STRAIN and STRESS
      ! ========================================================
 
      ! Thermal Strain
      call getshapefunc(etype, naturalcoord, spfunc)
      ttc = dot_product(tt, spfunc)
      tt0 = dot_product(t0, spfunc)
      ttn = dot_product(tn, spfunc)
      call cal_thermal_expansion_c3( tt0, ttc, gausses(lx)%pMaterial, coordsys, matlaniso, epsth )
 
      ! Update strain
      if( flag == infinitesimal ) then
        dvol = dot_product( totaldisp(1,1:nn), gderiv1_ave(1:nn,1) ) !du1/dx1
        dvol = dvol + dot_product( totaldisp(2,1:nn), gderiv1_ave(1:nn,2) ) !du2/dx2
        dvol = dvol + dot_product( totaldisp(3,1:nn), gderiv1_ave(1:nn,3) ) !du3/dx3
        dvol = (dvol-(gdispderiv(1,1)+gdispderiv(2,2)+gdispderiv(3,3)))/3.d0
        dstrain(1) = gdispderiv(1, 1) + dvol
        dstrain(2) = gdispderiv(2, 2) + dvol
        dstrain(3) = gdispderiv(3, 3) + dvol
        dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
        dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
        dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
        dstrain(:) = dstrain(:)-epsth(:)
 
        gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
 
      else if( flag == totallag ) then
        fbar(1:ndof, 1:ndof) = jratio(lx)*(i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof))
 
        ! Green-Lagrange strain
        dstrain(1) = 0.5d0*(dot_product(fbar(1:3,1),fbar(1:3,1))-1.d0)
        dstrain(2) = 0.5d0*(dot_product(fbar(1:3,2),fbar(1:3,2))-1.d0)
        dstrain(3) = 0.5d0*(dot_product(fbar(1:3,3),fbar(1:3,3))-1.d0)
        dstrain(4) = dot_product(fbar(1:3,1),fbar(1:3,2))
        dstrain(5) = dot_product(fbar(1:3,2),fbar(1:3,3))
        dstrain(6) = dot_product(fbar(1:3,3),fbar(1:3,1))
        dstrain(:) = dstrain(:)-epsth(:)
 
        gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
 
      else if( flag == updatelag ) then
        dvol = dot_product( totaldisp(1,1:nn), gderiv05_ave(1:nn,1) ) !du1/dx1
        dvol = dvol + dot_product( totaldisp(2,1:nn), gderiv05_ave(1:nn,2) ) !du2/dx2
        dvol = dvol + dot_product( totaldisp(3,1:nn), gderiv05_ave(1:nn,3) ) !du3/dx3
        dvol = (dvol-(gdispderiv(1,1)+gdispderiv(2,2)+gdispderiv(3,3)))/3.d0
        dstrain(1) = gdispderiv(1, 1) + dvol
        dstrain(2) = gdispderiv(2, 2) + dvol
        dstrain(3) = gdispderiv(3, 3) + dvol
        dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
        dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
        dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
        dstrain(:) = dstrain(:)-epsth(:)
 
        rot = 0.0d0
        rot(1, 2)= 0.5d0*( gdispderiv(1, 2)-gdispderiv(2, 1) );  rot(2, 1) = -rot(1, 2)
        rot(2, 3)= 0.5d0*( gdispderiv(2, 3)-gdispderiv(3, 2) );  rot(3, 2) = -rot(2, 3)
        rot(1, 3)= 0.5d0*( gdispderiv(1, 3)-gdispderiv(3, 1) );  rot(3, 1) = -rot(1, 3)
 
        gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+dstrain(1:6)+epsth(:)
 
        call getglobalderiv( etype, nn, naturalcoord, elem0, det, gderiv1)
        gdispderiv(1:ndof, 1:ndof) = matmul( du(1:ndof, 1:nn)+u(1:ndof, 1:nn), gderiv1(1:nn, 1:ndof) )
        fbar(1:ndof, 1:ndof) = jratio(lx)*(i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof))
 
      end if
 
      ! Update stress
      call update_stress3d( flag, gausses(lx), rot, dstrain, fbar, coordsys, time, tincr, ttc, tt0, ttn )
 
      ! ========================================================
      ! calculate the internal force ( equivalent nodal force )
      ! ========================================================
      ! Small strain
      b(1:6, 1:nn*ndof) = 0.0d0
      do j = 1, nn
        b(1,3*j-2) = gderiv(j, 1)
        b(2,3*j-1) = gderiv(j, 2)
        b(3,3*j  ) = gderiv(j, 3)
        b(4,3*j-2) = gderiv(j, 2)
        b(4,3*j-1) = gderiv(j, 1)
        b(5,3*j-1) = gderiv(j, 3)
        b(5,3*j  ) = gderiv(j, 2)
        b(6,3*j-2) = gderiv(j, 3)
        b(6,3*j  ) = gderiv(j, 1)
      end do
 
      wg=getweight( etype, lx )*det
      if( flag == infinitesimal ) then
        gderiv1(1:nn, 1:ndof) = gderiv(1:nn, 1:ndof)
 
        b2(1:6, 1:nn*ndof) = 0.0d0
        do j = 1, nn
          z1(1:3) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
          b2(1,3*j-2:3*j) = z1(1:3)
          b2(2,3*j-2:3*j) = z1(1:3)
          b2(3,3*j-2:3*j) = z1(1:3)
        end do
 
        ! BL = BL0 + BL2
        do j = 1, nn*ndof
          b(:, j) = b(:,j)+b2(:,j)
        end do
 
      else if( flag == totallag ) then
 
        b1(1:6, 1:nn*ndof) = 0.0d0
        do j = 1, nn
          b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
          b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
          b1(1, 3*j  ) = gdispderiv(3, 1)*gderiv(j, 1)
          b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
          b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
          b1(2, 3*j  ) = gdispderiv(3, 2)*gderiv(j, 2)
          b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
          b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
          b1(3, 3*j  ) = gdispderiv(3, 3)*gderiv(j, 3)
          b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
          b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
          b1(4, 3*j  ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
          b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
          b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
          b1(5, 3*j  ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
          b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
          b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
          b1(6, 3*j  ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
        end do
 
        call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv1)
 
        b2(1:6, 1:nn*ndof) = 0.0d0
        do j = 1, nn
          z1(1:3) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
          b2(1,3*j-2:3*j) = (2.d0*dstrain(1)+1.d0)*z1(1:3)
          b2(2,3*j-2:3*j) = (2.d0*dstrain(2)+1.d0)*z1(1:3)
          b2(3,3*j-2:3*j) = (2.d0*dstrain(3)+1.d0)*z1(1:3)
          b2(4,3*j-2:3*j) = 2.d0*dstrain(4)*z1(1:3)
          b2(5,3*j-2:3*j) = 2.d0*dstrain(5)*z1(1:3)
          b2(6,3*j-2:3*j) = 2.d0*dstrain(6)*z1(1:3)
        end do
 
        ! BL = R^2*(BL0 + BL1)+BL2
        do j = 1, nn*ndof
          b(:, j) = jratio(lx)*jratio(lx)*(b(:,j)+b1(:,j))+b2(:,j)
        end do
 
      else if( flag == updatelag ) then
 
        call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv)
        wg = (jratio(lx)**3.d0)*getweight(etype, lx)*det
        b(1:6, 1:nn*ndof) = 0.0d0
        do j = 1, nn
          b(1, 3*j-2) = gderiv(j, 1)
          b(2, 3*j-1) = gderiv(j, 2)
          b(3, 3*j  ) = gderiv(j, 3)
          b(4, 3*j-2) = gderiv(j, 2)
          b(4, 3*j-1) = gderiv(j, 1)
          b(5, 3*j-1) = gderiv(j, 3)
          b(5, 3*j  ) = gderiv(j, 2)
          b(6, 3*j-2) = gderiv(j, 3)
          b(6, 3*j  ) = gderiv(j, 1)
        end do
 
        b2(1:3, 1:nn*ndof) = 0.0d0
        do j = 1, nn
          z1(1:3) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
          b2(1, 3*j-2:3*j) = z1(1:3)
          b2(2, 3*j-2:3*j) = z1(1:3)
          b2(3, 3*j-2:3*j) = z1(1:3)
        end do
 
        do j = 1, nn*ndof
          b(1:3, j) = b(1:3,j)+b2(1:3,j)
        end do
 
      end if
 
      qf(1:nn*ndof)                                                          &
        = qf(1:nn*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,1:nn*ndof) )*wg
 
      ! integrate strain energy
      gausses(lx)%strain_energy = gausses(lx)%strain_energy*wg
    end do
 
  end subroutine update_c3d8fbar
 
  !----------------------------------------------------------------------*
  subroutine tload_c3d8fbar                    &
      (etype, nn, xx, yy, zz, tt, t0,   &
      gausses, vect, cdsys_id, coords)
    !----------------------------------------------------------------------*
 
    use m_fstr
    use mmechgauss
    use m_matmatrix
    use m_utilities
 
    !---------------------------------------------------------------------
 
    integer(kind=kint), parameter   :: ndof = 3
    integer(kind=kint), intent(in)  :: etype, nn
    type(tgaussstatus), intent(in)  :: gausses(:)
    real(kind=kreal), intent(in)    :: xx(nn), yy(nn), zz(nn)
    real(kind=kreal), intent(in)    :: tt(nn), t0(nn)
    real(kind=kreal), intent(out)   :: vect(nn*ndof)
    integer(kind=kint), intent(in)  :: cdsys_ID
    real(kind=kreal), intent(inout) :: coords(3, 3)           
 
    !---------------------------------------------------------------------
 
    real(kind=kreal) :: alp, alp0, d(6, 6), b(6, ndof*nn)
    real(kind=kreal) :: z1(3), det, ecoord(3, nn)
    integer(kind=kint) :: j, LX, serr
    real(kind=kreal) :: estrain(6), sgm(6), h(nn)
    real(kind=kreal) :: naturalcoord(3), gderiv(nn, 3)
    real(kind=kreal) :: wg, outa(1), ina(1), gderiv1_ave(nn, 3), alpo(3), alpo0(3), coordsys(3, 3)
    real(kind=kreal) :: tempc, temp0, v0, jacob_ave, thermal_eps, tm(6,6)
    logical :: ierr, matlaniso
 
    !---------------------------------------------------------------------
 
    matlaniso = .false.
 
    if( cdsys_id > 0 ) then   ! cannot define aniso expansion when no local coord defined
      ina = tt(1)
      call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
      if( .not. ierr ) matlaniso = .true.
    end if
 
    vect(:) = 0.0d0
 
    ecoord(1, :) = xx(:)
    ecoord(2, :) = yy(:)
    ecoord(3, :) = zz(:)
 
    !cal volumetric average of J=detF and dN/dx
    v0 = 0.d0
    jacob_ave = 0.d0
    gderiv1_ave(1:nn,1:ndof) = 0.d0
    do lx = 1, numofquadpoints(etype)
      call getquadpoint( etype, lx, naturalcoord(:) )
      call getglobalderiv( etype, nn, naturalcoord, ecoord, det, gderiv)
      wg = getweight(etype, lx)*det
      v0 = v0 + wg
      jacob_ave = jacob_ave + wg
      gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + wg*gderiv(1:nn, 1:ndof)
    enddo
    if( dabs(v0) > 1.d-12 ) then
      if( dabs(jacob_ave) < 1.d-12 ) stop 'Error in TLOAD_C3D8Fbar: too small average jacob'
      jacob_ave = jacob_ave/v0
      gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof)/(v0*jacob_ave)
    else
      stop 'Error in TLOAD_C3D8Fbar: too small element volume'
    end if
 
    ! LOOP FOR INTEGRATION POINTS
    do lx = 1, numofquadpoints(etype)
 
      call getquadpoint( etype, lx, naturalcoord(:) )
      call getshapefunc( etype, naturalcoord, h(1:nn) )
      call getglobalderiv(etype, nn, naturalcoord, ecoord, det, gderiv)
 
      if( cdsys_id > 0 ) then
        call set_localcoordsys(coords, g_localcoordsys(cdsys_id), coordsys, serr)
        if( serr == -1 ) stop "Fail to setup local coordinate"
        if( serr == -2 ) then
          write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
        end if
      end if
 
      !  WEIGHT VALUE AT GAUSSIAN POINT
      wg = getweight(etype, lx)*det
      b(1:6, 1:nn*ndof) = 0.0d0
      do j = 1, nn
        b(1,3*j-2) = gderiv(j, 1)
        b(2,3*j-1) = gderiv(j, 2)
        b(3,3*j  ) = gderiv(j, 3)
        b(4,3*j-2) = gderiv(j, 2)
        b(4,3*j-1) = gderiv(j, 1)
        b(5,3*j-1) = gderiv(j, 3)
        b(5,3*j  ) = gderiv(j, 2)
        b(6,3*j-2) = gderiv(j, 3)
        b(6,3*j  ) = gderiv(j, 1)
      end do
 
      do j = 1, nn
        z1(1:3) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
        b(1,3*j-2:3*j) = b(1,3*j-2:3*j)+z1(1:3)
        b(2,3*j-2:3*j) = b(2,3*j-2:3*j)+z1(1:3)
        b(3,3*j-2:3*j) = b(3,3*j-2:3*j)+z1(1:3)
      end do
 
      tempc = dot_product( h(1:nn),tt(1:nn) )
      temp0 = dot_product( h(1:nn),t0(1:nn) )
 
      call matlmatrix( gausses(lx), d3, d, 1.d0, 0.0d0, coordsys, tempc )
 
      ina(1) = tempc
      if( matlaniso ) then
        call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo(:), ierr, ina )
        if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
      else
        call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
        if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
        alp = outa(1)
      end if
      ina(1) = temp0
      if( matlaniso ) then
        call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo0(:), ierr, ina )
        if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
      else
        call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
        if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
        alp0 = outa(1)
      end if
 
      !**
      !** THERMAL strain
      !**
      if( matlaniso ) then
        do j = 1, 3
          estrain(j) = alpo(j)*(tempc-ref_temp)-alpo0(j)*(temp0-ref_temp)
        end do
        estrain(4:6) = 0.0d0
        call transformation(coordsys, tm)
        estrain(:) = matmul( estrain(:), tm  )      ! to global coord
        estrain(4:6) = estrain(4:6)*2.0d0
      else
        thermal_eps = alp*(tempc-ref_temp)-alp0*(temp0-ref_temp)
        estrain(1:3) = thermal_eps
        estrain(4:6) = 0.0d0
      end if
 
      !**
      !** SET SGM  {s}=[D]{e}
      !**
      sgm(:) = matmul( d(:, :), estrain(:) )
 
      !**
      !** CALCULATE LOAD {F}=[B]T{e}
      !**
      vect(1:nn*ndof) = vect(1:nn*ndof)+matmul( sgm(:), b(:, :) )*wg
 
    end do
 
  end subroutine tload_c3d8fbar
 
 
end module m_static_lib_fbar