!======================================================================!
!                                                                      !
! Software Name : FrontISTR Ver. 3.6                                   !
!                                                                      !
!      Module Name : lib                                               !
!                                                                      !
!                    Written by Xi YUAN (AdavanceSoft)                 !
!                               K. Satoh (Advancesoft)                 !
!                                                                      !
!      Contact address :  IIS,The University of Tokyo, CISS            !
!                                                                      !
!      "Structural Analysis for Large Scale Assembly"                  !
!                                                                      !
!======================================================================!
!======================================================================!
!
!> \brief  This module provides functions for hyperelastic calculation
!
!>  \author     X.Yuan(Advancesoft), K. Sato(Advancesoft)
!>  \date       2009/1/26
!>  \version    0.00
!
!======================================================================!
module mHyperElastic

  use mMaterial
  implicit none
  integer, parameter, private :: kreal = kind(0.0d0)

  contains

  !> This subroutine calculates derivative of the invariant with respect to Cauchy-Green tensor
  ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) <
  ! subroutine cderiv( matl, sectType, ctn, itn,               &
  !             inv1b, inv2b, inv3b, dibdc, d2ibdc2, strain    )
  subroutine cderiv(matl, sectType, ctn, itn,                            &
                    inv1b, inv2b, inv3b, dibdc, d2ibdc2, strain, ctninv) 
  ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23)
    type( tMaterial ), intent(in) :: matl                  !< material rpoperties
    integer, intent(in)           :: sectType              !< not used currently
    real(kind=kreal), intent(out) :: inv1b                 !< invariants
    real(kind=kreal), intent(out) :: inv2b                 !< invariants
    real(kind=kreal), intent(out) :: inv3b                 !< invariants
    real(kind=kreal), intent(out) :: dibdc(3,3,3)          !< derivative of the invariant with respect to c(i,j)
    real(kind=kreal), intent(out) :: d2ibdc2(3,3,3,3,3)    !< derivative of the invariant with respect to c(i,j)
    real(kind=kreal), intent(in)  :: strain(6)             !< Cauchy-Lagrange strain tensor
    real(kind=kreal), intent(out) :: ctn(3,3)              !< right Cauchy-Green deformation tensor
    real(kind=kreal), intent(out) :: itn(3,3)              !< identity tensor
    ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) <
    REAL( KIND = kreal ), INTENT(OUT) :: ctninv(3, 3)              !< inverse of right Cauchy-Green deformation tensor
    ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23)
    

    integer :: i, j, k, l, m, n

    real(kind=kreal) :: inv1, inv2, inv3, inv33
    real(kind=kreal) :: delta(3,3)
    ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) <
    ! real(kind=kreal) :: didc(3,3,3), ctninv(3,3)
    real(kind=kreal) :: didc(3,3,3)
    ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23)
    real(kind=kreal) :: d2idc2(3,3,3,3,3)

!   Kronecker's delta
    delta(:,:) = 0.d0
    delta(1,1) = 1.d0
    delta(2,2) = 1.d0
    delta(3,3) = 1.d0

    ! Green-Lagrange strain to right Cauchy-Green deformation tensor
    ctn(1,1)=strain(1)*2.d0+1.d0
    ctn(2,2)=strain(2)*2.d0+1.d0
    ctn(3,3)=strain(3)*2.d0+1.d0
    ctn(1,2)=strain(4);   ctn(2,1)=ctn(1,2)
    ctn(2,3)=strain(5);   ctn(3,2)=ctn(2,3)
    ctn(3,1)=strain(6);   ctn(1,3)=ctn(3,1)

    itn(:,:) = delta(:,:)

! ----- calculate the invariant of C
    inv1 = ctn(1,1)+ctn(2,2)+ctn(3,3)
    inv2 = ctn(2,2)*ctn(3,3)+ctn(1,1)*ctn(3,3)+ctn(1,1)*ctn(2,2)                   &
           -ctn(2,3)*ctn(2,3)-ctn(1,3)*ctn(1,3)-ctn(1,2)*ctn(1,2)
    inv3 =  ctn(1,1)*ctn(2,2)*ctn(3,3)                                             &
           +ctn(2,1)*ctn(3,2)*ctn(1,3)                                             &
           +ctn(3,1)*ctn(1,2)*ctn(2,3)                                             &
           -ctn(3,1)*ctn(2,2)*ctn(1,3)                                             &
           -ctn(2,1)*ctn(1,2)*ctn(3,3)                                             &
           -ctn(1,1)*ctn(3,2)*ctn(2,3)
    inv33 = inv3**(-1.d0/3.d0)

! ----- calculate the inverse of C
    ctninv(1,1)=(ctn(2,2)*ctn(3,3)-ctn(2,3)*ctn(2,3))/inv3
    ctninv(2,2)=(ctn(1,1)*ctn(3,3)-ctn(1,3)*ctn(1,3))/inv3
    ctninv(3,3)=(ctn(1,1)*ctn(2,2)-ctn(1,2)*ctn(1,2))/inv3
    ctninv(1,2)=(ctn(1,3)*ctn(2,3)-ctn(1,2)*ctn(3,3))/inv3
    ctninv(1,3)=(ctn(1,2)*ctn(2,3)-ctn(2,2)*ctn(1,3))/inv3
    ctninv(2,3)=(ctn(1,2)*ctn(1,3)-ctn(1,1)*ctn(2,3))/inv3
    ctninv(2,1)=ctninv(1,2)
    ctninv(3,1)=ctninv(1,3)
    ctninv(3,2)=ctninv(2,3)

! ----- calculate the derivative of the C-invariant with respect to c(i,j)
    do i=1,3
      do j=1,3
        didc(i,j,1) = delta(i,j)
        didc(i,j,2) = inv1*delta(i,j)-ctn(i,j)
        didc(i,j,3) = inv3*ctninv(i,j)
      enddo
    enddo

! ----- calculate the secont derivative of the C-invariant
    do k=1,3
     do l=1,3
      do m=1,3
       do n=1,3
         d2idc2(k,l,m,n,1)=0.d0
         d2idc2(k,l,m,n,2)=delta(k,l)*delta(m,n)-             &
            (delta(k,m)*delta(l,n)+delta(k,n)*delta(l,m))/2.d0
         d2idc2(k,l,m,n,3)=inv3*(ctninv(m,n)*ctninv(k,l)-   &
            (ctninv(k,m)*ctninv(n,l)+ctninv(k,n)*ctninv(m,l))/2.d0)
       enddo
      enddo
     enddo
    enddo

! ----- derivatives for the reduced invariants
    inv1b = inv1*inv33
    inv2b = inv2*inv33*inv33
    inv3b = dsqrt(inv3)

!   --- first derivative the reduced c-invarians w.r.t. c(i,j)
    do i=1,3
      do j=1,3
        dibdc(i,j,1) = -inv33**4*inv1*didc(i,j,3)/3.d0 + inv33*didc(i,j,1)
        dibdc(i,j,2) = -2.d0*inv33**5*inv2*didc(i,j,3)/3.d0 + inv33**2*didc(i,j,2)
        dibdc(i,j,3) = didc(i,j,3)/(2.d0*dsqrt(inv3))
      enddo
    enddo

!   --- second derivative of the reduced c-invariants w.r.t. c(i,j)
    do i=1,3
     do j=1,3
      do k=1,3
       do l=1,3
         d2ibdc2(i,j,k,l,1) =  4.d0/9.d0*inv33**7*inv1*didc(i,j,3)*didc(k,l,3)                       &
                              -inv33**4/3.d0*(didc(k,l,1)*didc(i,j,3)+didc(i,j,1)*didc(k,l,3))       &
                              -inv33**4/3.d0*inv1*d2idc2(i,j,k,l,3)                                  &
                              +inv33*d2idc2(i,j,k,l,1)
         d2ibdc2(i,j,k,l,2) = 10.d0/9.d0*inv33**8*inv2*didc(i,j,3)*didc(k,l,3)                       &
                              -2.d0/3.d0*inv33**5*(didc(k,l,2)*didc(i,j,3)+didc(i,j,2)*didc(k,l,3))  &
                              -2.d0/3.d0*inv33**5*inv2*d2idc2(i,j,k,l,3)                             &
                              +inv33**2*d2idc2(i,j,k,l,2)
         d2ibdc2(i,j,k,l,3) = -didc(i,j,3)*didc(k,l,3)/(4.d0*inv3**1.5d0)                            &
                              +d2idc2(i,j,k,l,3)/(2.d0*dsqrt(inv3))
       enddo
      enddo
     enddo
    enddo

  end subroutine cderiv

!-------------------------------------------------------------------------------
!> This subroutine provides elastic tangent coefficient for Arruda-Boyce hyperelastic material
!
!-------------------------------------------------------------------------------
  subroutine calElasticArrudaBoyce( matl, sectType, cijkl, strain )
    type( tMaterial ), intent(in) :: matl             !< material rpoperties
    integer, intent(in)           :: sectType         !< not used currently
    real(kind=kreal), intent(out) :: cijkl(3,3,3,3)   !< constitutive relation
    real(kind=kreal), intent(in)  :: strain(6)        !< Cauchy-Lagrange strain tensor

    integer :: i, j, k, l
    real(kind=kreal) :: ctn(3,3), itn(3,3)
    real(kind=kreal) :: inv1b, inv2b, inv3b
    real(kind=kreal) :: dibdc(3,3,3)
    real(kind=kreal) :: d2ibdc2(3,3,3,3,3)
    real(kind=kreal) :: constant(3), coef
    ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) <
    REAL( KIND = kreal ) :: ctninv(3, 3)
    ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23)
    
    ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) <
    ! call cderiv(  matl, sectType, ctn, itn, inv1b, inv2b, inv3b,            &
    !                         dibdc, d2ibdc2, strain    )
    call cderiv(matl, sectType, ctn, itn, inv1b, inv2b, inv3b, &
                dibdc, d2ibdc2, strain, ctninv)                
    ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23)
    
    constant(1:3)=matl%variables(M_PLCONST1:M_PLCONST3)
    coef = constant(2)

    forall( i=1:3, j=1:3, k=1:3, l=1:3 )
        cijkl(i,j,k,l) =  constant(1)*(1.d0/(10.d0*coef**2)                            &
                                        +66.d0*inv1b/(1050.d0*coef**4)                 &
                                        +228.d0*inv1b**2/(7000.d0*coef**6)             &
                                        +10380.d0*inv1b**3/(673750.d0*coef**8))        &
                                      *dibdc(i,j,1)*dibdc(k,l,1)                       &
                          +constant(1)*(0.5d0+inv1b/(10.d0*coef**2)                    &
                                        +33.d0*inv1b**2/(1050.d0*coef**4)              &
                                        +76.d0*inv1b**3/(7000.d0*coef**6)              &
                                        +2595.d0*inv1b**4/(673750.d0*coef**8))         &
                                      *d2ibdc2(i,j,k,l,1)                              &
                          +(1.d0+1.d0/inv3b**2)*dibdc(i,j,3)*dibdc(k,l,3)/constant(3)  &
                          +(inv3b-1.d0/inv3b)*d2ibdc2(i,j,k,l,3)/constant(3)
    end forall
    cijkl(:,:,:,:) = 4.d0*cijkl(:,:,:,:)

  end subroutine calElasticArrudaBoyce

!-------------------------------------------------------------------------------
!> This subroutine provides to update stress and strain for Arrude-Royce material
  subroutine calUpdateElasticArrudaBoyce( matl, sectType, dstrain, dstress )
    type( tMaterial ), intent(in) :: matl         !> material properties
    integer, intent(in)           :: sectType     !> not used currently
    real(kind=kreal), intent(out) :: dstress(6)   !> Cauchy-Green strain
    real(kind=kreal), intent(in)  :: dstrain(6)   !> 2nd Piola-Kirchhoff stress

    integer :: i, j, k
    real(kind=kreal) :: ctn(3,3), itn(3,3)
    real(kind=kreal) :: inv1b, inv2b, inv3b
    real(kind=kreal) :: dibdc(3,3,3)
    real(kind=kreal) :: d2ibdc2(3,3,3,3,3)
    real(kind=kreal) :: constant(3), coef
    real(kind=kreal) :: PKstress(3,3)
    ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) <
    REAL( KIND = kreal ) :: ctninv(3, 3)
    ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23)

    constant(1:3)=matl%variables(M_PLCONST1:M_PLCONST3)
    coef = constant(2)
    ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) <
    ! call cderiv(  matl, sectType, ctn, itn,      &
    !                          inv1b, inv2b, inv3b,           &
    !                          dibdc, d2ibdc2, dstrain    )
    call cderiv(matl, sectType, ctn, itn,        &
                inv1b, inv2b, inv3b,             &
                dibdc, d2ibdc2, dstrain, ctninv) 
    ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23)

! ----- calculate the stress
    do i=1,3
     do j=1,3
       PKstress(i,j) = constant(1)*( 0.5d0+inv1b/(10.d0*coef**2)                     &
                                +33.d0*inv1b*inv1b/(1050.d0*coef**4)                 &
                                +76.d0*inv1b**3/(7000.d0*coef**6)                    &
                                +2595.d0*inv1b**4/(673750.d0*coef**8))*dibdc(i,j,1)  &
                      +(inv3b-1.d0/inv3b)*dibdc(i,j,3)/constant(3)
       enddo
    enddo

    dstress(1) = 2.d0*PKstress(1,1)
    dstress(2) = 2.d0*PKstress(2,2)
    dstress(3) = 2.d0*PKstress(3,3)
    dstress(4) = PKstress(1,2) + PKstress(2,1)
    dstress(5) = PKstress(2,3) + PKstress(3,2)
    dstress(6) = PKstress(1,3) + PKstress(3,1)

  end subroutine calUpdateElasticArrudaBoyce

!-------------------------------------------------------------------------------
!> This subroutine provides elastic tangent coefficient for Mooney-Rivlin hyperelastic material
!-------------------------------------------------------------------------------
  subroutine calElasticMooneyRivlin( matl, sectType, cijkl, strain )
    type( tMaterial ), intent(in) :: matl             !< material rpoperties
    integer, intent(in)           :: sectType         !< not used curr
    real(kind=kreal), intent(out) :: cijkl(3,3,3,3)   !< constitutive relation
    real(kind=kreal), intent(in)  :: strain(6)        !< Cauchy-Lagrange strain tensor

    integer :: i, j, k, l, m, n, jj
    real(kind=kreal) :: ctn(3,3), itn(3,3)
    real(kind=kreal) :: inv1b, inv2b, inv3b
    real(kind=kreal) :: dibdc(3,3,3)
    real(kind=kreal) :: d2ibdc2(3,3,3,3,3)
    real(kind=kreal) :: constant(3), coef
    ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) <
    REAL( KIND = kreal ) :: ctninv(3, 3)
    ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23)

    constant(1:3)=matl%variables(M_PLCONST1:M_PLCONST3)
    ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) <
    ! call cderiv( matl, sectType, ctn, itn, inv1b, inv2b, inv3b,            &
    !                          dibdc, d2ibdc2, strain    )
    call cderiv(matl, sectType, ctn, itn, inv1b, inv2b, inv3b, &
                dibdc, d2ibdc2, strain, ctninv)                
    ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23)

        forall( k=1:3, l=1:3, m=1:3, n=1:3 )
            cijkl(k,l,m,n) = d2ibdc2(k,l,m,n,1)*constant(1) +  &
                             d2ibdc2(k,l,m,n,2)*constant(2) +  &
		                     2.d0*(dibdc(k,l,3)*dibdc(m,n,3)+  &
                (inv3b-1.d0)*d2ibdc2(k,l,m,n,3))/constant(3)
        end forall
        cijkl(:,:,:,:)=4.d0*cijkl(:,:,:,:)

  end subroutine calElasticMooneyRivlin

!-------------------------------------------------------------------------------
!> This subroutine provides to update stress and strain for Mooney-Rivlin material
!-------------------------------------------------------------------------------
  subroutine calUpdateElasticMooneyRivlin( matl, sectType, strain, stress )
    type( tMaterial ), intent(in) :: matl        !< material properties
    integer, intent(in)           :: sectType    !< not used currently
    real(kind=kreal), intent(out) :: stress(6)   !< 2nd Piola-Kirchhoff stress
    real(kind=kreal), intent(in)  :: strain(6)   !< Green-Lagrangen strain

    integer :: k, l
    real(kind=kreal) :: ctn(3,3), itn(3,3)
    real(kind=kreal) :: inv1b, inv2b, inv3b
    real(kind=kreal) :: dibdc(3,3,3)
    real(kind=kreal) :: d2ibdc2(3,3,3,3,3)
    real(kind=kreal) :: constant(3)
    real(kind=kreal) :: CGstrain(3,3), dudc(3,3)
    ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) <
    REAL( KIND = kreal ) :: ctninv(3, 3)
    ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23)

    constant(1:3)=matl%variables(M_PLCONST1:M_PLCONST3)
    ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) <
    ! call cderiv( matl, sectType, ctn, itn, inv1b, inv2b, inv3b,      &
    !                          dibdc, d2ibdc2, strain    )
    call cderiv(matl, sectType, ctn, itn, inv1b, inv2b, inv3b, &
                dibdc, d2ibdc2, strain, ctninv)                
    ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23)

! ----- stress
        do l=1,3
        do k=1,3
		    dudc(k,l) = dibdc(k,l,1)*constant(1)+dibdc(k,l,2)*constant(2)  &
                            +2.d0*(inv3b-1.d0)*dibdc(k,l,3)/constant(3)
        enddo
        enddo

    stress(1)=2.d0*dudc(1,1)
    stress(2)=2.d0*dudc(2,2)
    stress(3)=2.d0*dudc(3,3)
    stress(4)=2.d0*dudc(1,2)
    stress(5)=2.d0*dudc(2,3)
    stress(6)=2.d0*dudc(1,3)

  end subroutine calUpdateElasticMooneyRivlin

      
      ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) <
!####################################################################
      SUBROUTINE calElasticMooneyRivlin_Incomp                             &
                 ( matl, sectType, cijkl, cijkl_p, strain, lambda, d2, g ) 
!####################################################################
      
      TYPE( tMaterial ), INTENT(IN) :: matl                    !< material rpoperties 
      INTEGER, INTENT(IN) :: sectType                          !< not used curr
      REAL( KIND = kreal ), INTENT(OUT) :: cijkl(3, 3, 3, 3)   !< constitutive relation
      REAL( KIND = kreal ), INTENT(OUT) :: cijkl_p(3, 3, 3, 3) !< constitutive relation
      REAL( KIND = kreal ), INTENT(IN)  :: strain(6)           !< Cauchy-Lagrange strain tensor 
      REAL( KIND = kreal ), INTENT(IN)  :: lambda
      REAL( KIND = kreal ), INTENT(OUT) :: d2(6)
      REAL( KIND = kreal ), INTENT(OUT) :: g
      
      INTEGER :: i, j, k, l, m, n, jj
      REAL( KIND = kreal ) :: ctn(3, 3), itn(3, 3)
      REAL( KIND = kreal ) :: inv1b, inv2b, inv3b
      REAL( KIND = kreal ) :: dibdc(3, 3, 3)
      REAL( KIND = kreal ) :: d2ibdc2(3, 3, 3, 3, 3)
      REAL( KIND = kreal ) :: dudc(3, 3)
      
!--------------------------------------------------------------------
      
      REAL( KIND = kreal ) :: ctninv(3, 3)
      
!--------------------------------------------------------------------
      
      REAL( KIND = kreal ) :: cijkl_tmp(3, 3, 3, 3)
      REAL( KIND = kreal ) :: c10, c01, c11, c20, c02, &
                              c21, c12, c30, c03       
      REAL( KIND = kreal ) :: det_f
      REAL( KIND = kreal ) :: i_rc, ii_rc
      REAL( KIND = kreal ) :: i_rc_minus_3, ii_rc_minus_3
      REAL( KIND = kreal ) :: dwdi_rc, dwdii_rc
      REAL( KIND = kreal ) :: ddwdidi_rc , ddwdidii_rc, &
                              ddwdiidi_rc, ddwdiidii_rc 
      
!--------------------------------------------------------------------
      
      c10 = matl%variables(M_PLCONST1)
      c01 = matl%variables(M_PLCONST2)
      c20 = matl%variables(M_PLCONST3)
      c11 = matl%variables(M_PLCONST4)
      c02 = matl%variables(M_PLCONST5)
      c30 = matl%variables(M_PLCONST6)
      c21 = matl%variables(M_PLCONST7)
      c12 = matl%variables(M_PLCONST8)
      c03 = matl%variables(M_PLCONST9)
      
!--------------------------------------------------------------------
      
      CALL cderiv                                           &
           ( matl, sectType, ctn, itn, inv1b, inv2b, inv3b, &
             dibdc, d2ibdc2, strain, ctninv )               
      
!--------------------------------------------------------------------
      
      i_rc = inv1b
      ii_rc = inv2b
      
      i_rc_minus_3  = i_rc-3.0D0
      ii_rc_minus_3 = ii_rc-3.0D0
      
!--------------------------------------------------------------------
      
      dwdi_rc                                 &
      = c10                                   &
       +c11*ii_rc_minus_3                     &
       +2.0D0*c20*i_rc_minus_3                &
       +2.0D0*c21*i_rc_minus_3*ii_rc_minus_3  &
       +c12*ii_rc_minus_3*ii_rc_minus_3       &
       +3.0D0*c30*i_rc_minus_3*i_rc_minus_3   
      dwdii_rc                                &
      = c01                                   &
       +c11*i_rc_minus_3                      &
       +2.0D0*c02*ii_rc_minus_3               &
       +c21*i_rc_minus_3*i_rc_minus_3         &
       +2.0D0*c12*i_rc_minus_3*ii_rc_minus_3  &
       +3.0D0*c03*ii_rc_minus_3*ii_rc_minus_3 
      
      ddwdidi_rc                &
      = 2.0D0*c20               &
       +2.0D0*c21*ii_rc_minus_3 &
       +6.0D0*c30*i_rc_minus_3  
      ddwdidii_rc               &
      = c11                     &
       +2.0D0*c21*i_rc_minus_3  &
       +2.0D0*c12*ii_rc_minus_3 
      ddwdiidi_rc = ddwdidii_rc 
      ddwdiidii_rc              &
      = 2.0D0*c02               &
       +2.0D0*c12*i_rc_minus_3  &
       +6.0D0*c03*ii_rc_minus_3 
      
!--------------------------------------------------------------------
      
      FORALL( k = 1:3, l = 1:3, m = 1:3, n = 1:3 )
       
       cijkl_tmp(k, l, m, n)                                    &
       = 4.0D0*( dibdc(k, l, 1)*( ddwdidi_rc *dibdc(m, n, 1)    &
                                 +ddwdidii_rc*dibdc(m, n, 2) )  &
                +dwdi_rc*d2ibdc2(k, l, m, n, 1)                 &
                +dibdc(k, l, 2)*( ddwdiidi_rc *dibdc(m, n, 1)   &
                                 +ddwdiidii_rc*dibdc(m, n, 2) ) &
                +dwdii_rc*d2ibdc2(k, l, m, n, 2) )              
       
      END FORALL
      
      FORALL( k = 1:3, l = 1:3, m = 1:3, n = 1:3 )
       
       cijkl(k, l, m, n) = 0.5D0*( cijkl_tmp(k, l, m, n)+cijkl_tmp(k, l, n, m) )
       
      END FORALL
      
!--------------------------------------------------------------------
      
      det_f = inv3b
      
      g = 2.0D0*( det_f-1.0D0 )
      
!--------------------------------------------------------------------
      
      DO l = 1, 3
       
       DO k = 1, 3
        
        dudc(k, l) = det_f*ctninv(k, l)
        
       END DO
       
      END DO
      
!--------------------------------------------------------------------
      
      FORALL( k = 1:3, l = 1:3, m = 1:3, n = 1:3 )
       
       cijkl_tmp(k, l, m, n)                         &
       = 4.0D0                                       &
         *( lambda                                   &
            *( 0.5D0*det_f*ctninv(k, l)*ctninv(m, n) &
           -det_f*ctninv(k, m)*ctninv(n, l) ) )      
       
      END FORALL
      
      FORALL( k = 1:3, l = 1:3, m = 1:3, n = 1:3 )
       
       cijkl_p(k, l, m, n) = 0.5D0*( cijkl_tmp(k, l, m, n)+cijkl_tmp(k, l, n, m) )
       
      END FORALL
      
!--------------------------------------------------------------------
      
      d2(1) = 2.0D0*dudc(1, 1)
      d2(2) = 2.0D0*dudc(2, 2)
      d2(3) = 2.0D0*dudc(3, 3)
      d2(4) = 2.0D0*dudc(1, 2)
      d2(5) = 2.0D0*dudc(2, 3)
      d2(6) = 2.0D0*dudc(3, 1)
      
!--------------------------------------------------------------------
      
      RETURN
      
!####################################################################
      END SUBROUTINE calElasticMooneyRivlin_Incomp
!####################################################################
      ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23)
      
      
      ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) <
!####################################################################
      SUBROUTINE calUpdateElasticMooneyRivlin_Incomp                     &
                 ( matl, sectType, strain, stress, stress_p, lambda, g ) 
!####################################################################
      
      TYPE( tMaterial ), INTENT(IN) :: matl            !< material properties
      INTEGER, INTENT(IN) :: sectType                  !< not used currently
      REAL( KIND = kreal ), INTENT(OUT) :: stress(6)   !< 2nd Piola-Kirchhoff stress
      REAL( KIND = kreal ), INTENT(OUT) :: stress_p(6)
      REAL( KIND = kreal ), INTENT(IN) :: strain(6)    !< Green-Lagrangen strain
      REAL( KIND = kreal ), INTENT(IN) :: lambda
      REAL( KIND = kreal ), INTENT(OUT) :: g
      
      INTEGER :: k, l
      REAL( KIND = kreal ) :: ctn(3, 3), itn(3, 3)
      REAL( KIND = kreal ) :: inv1b, inv2b, inv3b
      REAL( KIND = kreal ) :: dibdc(3, 3, 3)
      REAL( KIND = kreal ) :: d2ibdc2(3, 3, 3, 3, 3)
      REAL( KIND = kreal ) :: dudc(3, 3)
      
!--------------------------------------------------------------------
      
      REAL( KIND = kreal ) :: ctninv(3, 3)
      
!--------------------------------------------------------------------
      
      REAL( KIND = kreal ) :: c10, c01, c11, c20, c02, &
                              c21, c12, c30, c03       
      REAL( KIND = kreal ) :: i_rc, ii_rc
      REAL( KIND = kreal ) :: i_rc_minus_3, ii_rc_minus_3
      REAL( KIND = kreal ) :: det_f
      REAL( KIND = kreal ) :: lambda_2
      REAL( KIND = kreal ) :: dwdi_rc, dwdii_rc
      
!--------------------------------------------------------------------
      
      c10 = matl%variables(M_PLCONST1)
      c01 = matl%variables(M_PLCONST2)
      c20 = matl%variables(M_PLCONST3)
      c11 = matl%variables(M_PLCONST4)
      c02 = matl%variables(M_PLCONST5)
      c30 = matl%variables(M_PLCONST6)
      c21 = matl%variables(M_PLCONST7)
      c12 = matl%variables(M_PLCONST8)
      c03 = matl%variables(M_PLCONST9)
      
!--------------------------------------------------------------------
      
      CALL cderiv                                           &
           ( matl, sectType, ctn, itn, inv1b, inv2b, inv3b, &
             dibdc, d2ibdc2, strain, ctninv )               
      
!--------------------------------------------------------------------
      
      i_rc = inv1b
      ii_rc = inv2b
      
      i_rc_minus_3  = i_rc-3.0D0
      ii_rc_minus_3 = ii_rc-3.0D0
      
!--------------------------------------------------------------------
      
      dwdi_rc                                 &
      = c10                                   &
       +c11*ii_rc_minus_3                     &
       +2.0D0*c20*i_rc_minus_3                &
       +2.0D0*c21*i_rc_minus_3*ii_rc_minus_3  &
       +c12*ii_rc_minus_3*ii_rc_minus_3       &
       +3.0D0*c30*i_rc_minus_3*i_rc_minus_3   
      dwdii_rc                                &
      = c01                                   &
       +c11*i_rc_minus_3                      &
       +2.0D0*c02*ii_rc_minus_3               &
       +c21*i_rc_minus_3*i_rc_minus_3         &
       +2.0D0*c12*i_rc_minus_3*ii_rc_minus_3  &
       +3.0D0*c03*ii_rc_minus_3*ii_rc_minus_3 
      
!--------------------------------------------------------------------
      
      DO l = 1, 3
       
       DO k = 1, 3
        
        dudc(k, l) = dwdi_rc*dibdc(k, l, 1)+dwdii_rc*dibdc(k, l, 2)
        
       END DO
       
      END DO
      
      stress(1) = 2.0D0*dudc(1, 1)
      stress(2) = 2.0D0*dudc(2, 2)
      stress(3) = 2.0D0*dudc(3, 3)
      stress(4) = 2.0D0*dudc(1, 2)
      stress(5) = 2.0D0*dudc(2, 3)
      stress(6) = 2.0D0*dudc(3, 1)
      
!--------------------------------------------------------------------
      
      det_f = inv3b
      
      g = 2.0D0*( det_f-1.0D0 )
      
!--------------------------------------------------------------------
      
      DO l = 1, 3
       
       DO k = 1, 3
        
        dudc(k, l) = det_f*ctninv(k, l)
        
       END DO
       
      END DO
      
      lambda_2 = 2.0D0*lambda
      
      ! { S } vector for pressure
      stress_p(1) = lambda_2*dudc(1, 1)
      stress_p(2) = lambda_2*dudc(2, 2)
      stress_p(3) = lambda_2*dudc(3, 3)
      stress_p(4) = lambda_2*dudc(1, 2)
      stress_p(5) = lambda_2*dudc(2, 3)
      stress_p(6) = lambda_2*dudc(3, 1)
      
!--------------------------------------------------------------------
      
      RETURN
      
!####################################################################
      END SUBROUTINE calUpdateElasticMooneyRivlin_Incomp
!####################################################################
      ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23)
      
end module
