C-------------------------------------------------------------
      subroutine vuanisohyper_inv (nb, nFiber, nInv, nElement, nIntPt,
     $     nLayer, nSecPt, cmname, nstatev, nfieldv, numprops,
     $     props, tempOld, tempNew, fieldOld, fieldNew, stateOld,
     $     sInvariant, zeta, uDev, duDi,d2uDiDi, stateNew)
C
*unUsedArgument zeta pvw
c
      include 'vaba_param.inc'
C
      character*80 cmname
      dimension  props(numprops),
     $     tempOld(nb), fieldOld(nb,nfieldv),stateOld(nb,nstatev),
     $     tempNew(nb), fieldNew(nb,nfieldv),stateNew(nb,nstatev),
     $     sInvariant(nb,nInv), zeta(nb,nFiber*(nFiber-1)/2),
     $     uDev(nb), duDi(nb,nInv), d2uDiDi(nb,nInv*(nInv+1)/2)
C
C
C
      if (cmname(1:11) .eq. 'VUANISO_HGO') then
         call VUANISOHYPER_INVHGO(sInvariant, uDev, nFiber, ninv,
     $     duDi, d2uDiDi, nb, numprops, props)
      else if(cmname(1:14) .eq. 'VUANISO_INVISO') then
         call VUANISOHYPER_INVISO(sInvariant, uDev, nFiber, ninv,
     $     duDi, d2uDiDi, nb, numprops, props)
      else if(cmname(1:17) .eq. 'VUANISO_FUNGINV44') then
         call VUANISOHYPER_FUNGINV44(sInvariant, uDev, nFiber,
     $     ninv, duDi, d2uDiDi, nb, numprops, props)
      else if(cmname(1:17) .eq. 'VUANISO_FUNGINV45') then
         call VUANISOHYPER_FUNGINV45(sInvariant, uDev, nFiber, 
     $     ninv, duDi, d2uDiDi, nb, numprops, props)
      else
         call xplb_abqerr(-2,'User subroutine VUANISOHYPER_INV missing!'
     *        ,intv,zero,' ')
         call xplb_exit
      end if
C
C
C
      return
      end
c------------------------------------------------------------------
c
c     HGO model
c
      subroutine vuanisohyper_invhgo (ainv, ua, nfibers, ninv,
     $     ui1, ui2, nb, numprops, props)
C
      include 'vaba_param.inc'
C
      dimension ua(nb), ainv(nb,ninv), ui1(nb,ninv),
     $     ui2(nb,ninv*(ninv+1)/2), props(numprops)
C
c     ainv: invariants
c     ua  : udev
c     ui1 : dUdI
c     ui2 : d2U/dIdJ
C
      parameter ( half = 0.5d0,
     *            zero = 0.d0, 
     *            one  = 1.d0, 
     *            two  = 2.d0, 
     *            three= 3.d0, 
     *            four = 4.d0, 
     *            five = 5.d0, 
     *            six  = 6.d0,
c
     *            index_I1 = 1,
     *            index_J  = 3,
     *            asmall   = 2.d-16  )
C
C     HGO model
C
      C10 = props(1)
      rk1 = props(3)
      rk2 = props(4)
      rkp = props(5)
c
      do kb = 1,nb
        ua(kb) = zero
        om3kp = one - three * rkp
        do k1 = 1, nfibers
          index_i4 = indxInv4(k1,k1)
          E_alpha1 = rkp  * (ainv(kb,index_i1) - three) 
     *            + om3kp * (ainv(kb,index_i4) - one  )
          E_alpha = max(E_alpha1, zero)
          ht4a    = half + sign(half,E_alpha1 + asmall)
          aux     = exp(rk2*E_alpha*E_alpha)
c energy
          ua(kb) = ua(kb) +  aux - one
c ui1
          ui1(kb,index_i1) = ui1(kb,index_i1) + aux * E_alpha
          ui1(kb,index_i4) = rk1 * om3kp * aux * E_alpha
c ui2
          aux2 = ht4a + two * rk2 * E_alpha * E_alpha
          ui2(kb,indx(index_I1,index_I1))=
     *           ui2(kb,indx(index_I1,index_I1))
     *         + aux * aux2
          ui2(kb,indx(index_I1,index_i4)) = rk1*rkp*om3kp * aux * aux2
          ui2(kb,indx(index_i4,index_i4)) = rk1*om3kp*om3kp*aux * aux2
        end do
c
c     deviatoric energy
c
        ua(kb) = ua(kb) * rk1 / (two * rk2)
        ua(kb) = ua(kb) + C10 * (ainv(kb,index_i1) - three)
c
c     compute derivatives
c
        ui1(kb,index_i1) = rk1 * rkp * ui1(kb,index_i1) + C10
        ui2(kb,indx(index_I1,index_I1))= ui2(kb,indx(index_I1,index_I1)) 
     *                            * rk1 * rkp * rkp
      end do
c     
c     compressible case
      if(props(2).gt.zero) then
        do kb = 1,nb
          Dinv = one / props(2)
          det = ainv(kb,index_J)
          ui1(kb,index_J) = Dinv * (det - one/det)
          ui2(kb,indx(index_J,index_J))= Dinv * (one + one / det / det)
        end do
      end if
c
      return
      end
C-------------------------------------------------------------
C     Function to map index from Square to Triangular storage 
C 		 of symmetric matrix
C
      integer function indx( i, j )
      include 'vaba_param.inc'
      ii = min(i,j)
      jj = max(i,j)
      indx = ii + jj*(jj-1)/2
      return
      end
C-------------------------------------------------------------
C
C     Function to generate enumeration of scalar
C     Pseudo-Invariants of type 4

      integer function indxInv4( i, j )
      include 'vaba_param.inc'
      ii = min(i,j)
      jj = max(i,j)
      indxInv4 = 4 + jj*(jj-1) + 2*(ii-1)
      return
      end
C-------------------------------------------------------------
C
C     Function to generate enumeration of scalar
C     Pseudo-Invariants of type 5
C
      integer function indxInv5( i, j )
      include 'vaba_param.inc'
      ii = min(i,j)
      jj = max(i,j)
      indxInv5 = 5 + jj*(jj-1) + 2*(ii-1)
      return
      end
C-------------------------------------------------------------
c
c    generalized Fung Anisotropic Model 
c    reformulated in terms of J and I4ab
c
c
      subroutine vuanisohyper_funginv44(ainv, ua, nfibers, ninv,
     $     ui1, ui2, nb, numprops, props)
C
      include 'vaba_param.inc'
C
      dimension ua(nb), ainv(nb,ninv), ui1(nb,ninv),
     $     ui2(nb,ninv*(ninv+1)/2), props(numprops)
C
c     ainv: invariants
c     ua  : udev
c     ui1 : dUdI
c     ui2 : d2U/dIdJ
C
      parameter ( half = 0.5d0,
     *            zero = 0.d0, 
     *            one  = 1.d0, 
     *            two  = 2.d0, 
     *            three= 3.d0, 
     *            four = 4.d0, 
     *            five = 5.d0, 
     *            six  = 6.d0,
c
     *            indx_J = 3 )
C
      dimension bmx(6,6),Cg(6),unit(6),index4(6)
      dimension dQdC(6),d2QdC2(6,6),ind4C2I(6),ind4I2C(6)
      dimension eg(6)
C
C     Fung model in invariant form --- in terms of J and I4ab
C
c     active invariants:
c
c     1    
c     2
c     3    yes   J
c     4    yes   I4(11)
c     5
c     6    yes   I4(12)
c     7
c     8    yes   I4(22)
c     9
c     10   yes   I4(13)
c     11
c     12   yes   I4(23)
c     13
c     14   yes   I4(33)
c     15
c
c     bmx 11, 22, 33, 12, 13, 23
c

c      write(*,*) 'inside 44 44 44'

      do k1=1,6
         do k2=1,k1
            ktmp = k2 + (k1-1)*k1/2
            if(k1.le.3) then
               bmx(k2,k1) = props(ktmp)
            else 
               if(k2.le.3) then
                  bmx(k2,k1) = two * props(ktmp)
               else 
                  bmx(k2,k1) = four * props(ktmp)
               end if
            end if
         end do
      end do
      do k1=1,6
         do k2=k1+1,6
            bmx(k2,k1) = bmx(k1,k2)
         end do
      end do
c
      Cval = props(22)
c
c      indx_411 = 4
c      indx_412 = 6
c      indx_422 = 8
c      indx_413 = 10
c      indx_423 = 12
c      indx_433 = 14
c
c      indxC_411 = 1
c      indxC_422 = 2
c      indxC_433 = 3
c      indxC_412 = 4  to be consistent with Fung input
c      indxC_413 = 5  to be consistent with Fung input
c      indxC_423 = 6  to be consistent with Fung input

c
c     index of I4ab (I4_11,12,22,13,23,33) in ainv
c
      index4(1) = 4
      index4(2) = 6
      index4(3) = 8
      index4(4) = 10
      index4(5) = 12
      index4(6) = 14
c
c     index of components in C to invariants
c
      ind4C2I(1) = 4
      ind4C2I(2) = 8
      ind4C2I(3) = 14
      ind4C2I(4) = 6
      ind4C2I(5) = 10
      ind4C2I(6) = 12
c
c     index of invariant number to components in C
c
      ind4I2C(1) = 1
      ind4I2C(2) = 4
      ind4I2C(3) = 2
      ind4I2C(4) = 5
      ind4I2C(5) = 6
      ind4I2C(6) = 3
c
      unit(1) = one
      unit(2) = one
      unit(3) = one
      unit(4) = zero
      unit(5) = zero
      unit(6) = zero
c
      do kb = 1,nb
        do k1=1,6
          Cg(k1) = ainv(kb,ind4C2I(k1))
          eg(k1) = half * (Cg(k1) - unit(k1))
        end do
c
        Qval = zero
        do k1=1,6
          do k2=1,6
            Qval = Qval + eg(k1)*bmx(k1,k2)*eg(k2)
          end do
        end do
        ua(kb) = half * Cval * (exp(Qval) - one)
c
        do k1=1,6
          dQdC(k1) = zero
          do k2=1,6
            d2QdC2(k1,k2) = half * bmx(k1,k2)
            dQdC(k1) = dQdC(k1) + bmx(k1,k2)*eg(k2)
          end do
        end do
        dUdQ = half * Cval * exp(Qval)
        d2UdQ2 = dUdQ
c
        do k1=1,6
          ui1(kb,index4(k1)) = dUdQ * dQdC(ind4I2C(k1))
          do k2=1,k1
            ui2(kb,indx(index4(k1),index4(k2))) = 
     *           dUdQ   * d2QdC2(ind4I2C(k1),ind4I2C(k2))
     *         + d2UdQ2 * dQdC(ind4I2C(k1)) * dQdC(ind4I2C(k2))
          end do
        end do
      end do
c
      if(props(23).gt.zero) then
        do kb = 1,nb
          det = ainv(kb,indx_J)
          Dinv = one / props(23)
          ui1(kb,indx_J) = Dinv * (det - one/det)
          ui2(kb,indx(indx_J,indx_J)) = Dinv * (one + one / det / det)
        end do
      end if
c
      return
      end
C----------------------------------------------------------------------
c
c     isotropic hyperelasticity through VUANISOHYPER_INV
c
      subroutine vuanisohyper_inviso (ainv, ua, nfibers, ninv,
     $     ui1, ui2, nb, numprops, props)
C
      include 'vaba_param.inc'
C
      dimension ua(nb), ainv(nb,ninv), ui1(nb,ninv),
     $     ui2(nb,ninv*(ninv+1)/2), props(numprops)
C
c     ainv: invariants
c     ua  : udev
c     ui1 : dUdI
c     ui2 : d2U/DiDi
C
      parameter ( half = 0.5d0,
     *            zero = 0.d0, 
     *            one  = 1.d0, 
     *            two  = 2.d0, 
     *            three= 3.d0, 
     *            four = 4.d0, 
     *            five = 5.d0, 
     *            six  = 6.d0,
     *            twt4 = 24.d0,
c
     *            index_I1 = 1,
     *            index_I2 = 2,
     *            index_J  = 3  )
C
C     -- Polynomial with N=2 --
C
      C10 = props(1)
      C01 = props(2)
      C20 = props(3)
      C11 = props(4)
      C02 = props(5)
      D1  = props(6)
      D2  = props(7)
c
      do kb = 1,nb
        rI1 = ainv(kb,index_I1)
        rI2 = ainv(kb,index_I2)
        rJ  = ainv(kb,index_J )
c
        rI1m3 = rI1 - three
        rI2m3 = rI2 - three
        rJm1  = rJ  - one
c
        ua(kb) =   C10 * rI1m3 + C01 * rI2m3
     *           + C20 * rI1m3 * rI1m3 + C11 * rI1m3 * rI2m3 
     *           + C02 * rI2m3 * rI2m3
        ui1(kb,index_I1) = C10 + two * C20 * rI1m3 + C11 * rI2m3
        ui1(kb,index_I2) = C01 + two * C02 * rI2m3 + C11 * rI1m3
c
        ui2(kb,indx(index_I1,index_I1)) = two * C20
        ui2(kb,indx(index_I1,index_I2)) = C11
        ui2(kb,indx(index_I2,index_I2)) = two * C02
      end do
*
c     compressible case
*
      if(D1.gt.zero) then
        do kb = 1,nb
          Dinv1 = one / D1
          rJ  = ainv(kb,index_J )
          rJm1  = rJ  - one
          ui1(kb,index_J) = two * Dinv1 * rJm1
          ui2(kb,indx(index_J,index_J)) = two * Dinv1
        end do
      end if
c     
      if(D2.gt.zero) then
        do kb = 1,nb
          Dinv2 = one / D2
          rJ  = ainv(kb,index_J )
          rJm1  = rJ  - one
          ui1(kb,index_J) = ui1(kb,index_J) + four * Dinv2 * rJm1**3
          ui2(kb,indx(index_J,index_J)) = ui2(kb,indx(index_J,index_J))
     *                              + three * four * Dinv2 * rJm1**2
        end do
      end if
c
      return
      end
c-------------------------------------------------------------------------
c
c    generalized Fung Anisotropic Model 
c    reformulated in terms of J, I4ab, and I5aa
c
      subroutine vuanisohyper_funginv45(ainv, ua, nfibers, ninv,
     $     ui1, ui2, nb, numprops, props)
C
      include 'vaba_param.inc'
C
      dimension ua(nb), ainv(nb,ninv), ui1(nb,ninv),
     $     ui2(nb,ninv*(ninv+1)/2), props(numprops)
C
c     ainv: invariants
c     ua  : udev
c     ui1 : dUdI
c     ui2 : d2U/dIdJ
C
      parameter ( half = 0.5d0,
     *            zero = 0.d0, 
     *            one  = 1.d0, 
     *            two  = 2.d0, 
     *            three= 3.d0, 
     *            four = 4.d0, 
     *            five = 5.d0, 
     *            six  = 6.d0,
     *            quat = 0.25d0,
c
     *            indx_J = 3 )
C
      dimension bmx(6,6),Cg(6),unit(6),indI(6),ind4C2I(6),eg(6)
      dimension dQdC(6),d2QdC2(6,6),aI(6),d2QdCdI(6,6)
      dimension dQdI(6),d2QdI2(6,6),dCdI(6,6),d2CdI2(6,6,6)
C
C     Fung model in invariant form --- in terms of J and I4ab
C
c     active invariants:
c
c     1    
c     2
c     3    yes   J
c     4    
c     5    yes   I5(11)
c     6    yes   I4(12)
c     7
c     8    
c     9    yes   I5(22)
c     10   yes   I4(13)
c     11
c     12   yes   I4(23)
c     13
c     14   
c     15   yes   I5(33)
c
c     bmx 11, 22, 33, 12, 13, 23
c
c      write(*,*) 'inside 45 45 45'

      do k1=1,6
         do k2=1,k1
            ktmp = k2 + (k1-1)*k1/2
            if(k1.le.3) then
               bmx(k2,k1) = props(ktmp)
            else 
               if(k2.le.3) then
                  bmx(k2,k1) = two * props(ktmp)
               else 
                  bmx(k2,k1) = four * props(ktmp)
               end if
            end if
         end do
      end do
      do k1=1,6
         do k2=k1+1,6
            bmx(k2,k1) = bmx(k1,k2)
         end do
      end do
c
      Cval = props(22)
c
c      indx_411 = 4
c      indx_412 = 6
c      indx_422 = 8
c      indx_413 = 10
c      indx_423 = 12
c      indx_433 = 14
c
c      indxC_411 = 1
c      indxC_422 = 2
c      indxC_433 = 3
c      indxC_412 = 4  to be consistent with Fung input
c      indxC_413 = 5  to be consistent with Fung input
c      indxC_423 = 6  to be consistent with Fung input

c
c     index of invariants 
c
      indI(1) = 5
      indI(2) = 6
      indI(3) = 9
      indI(4) = 10
      indI(5) = 12
      indI(6) = 15
c
      ind4C2I(1) = 4
      ind4C2I(2) = 8
      ind4C2I(3) = 14
      ind4C2I(4) = 6
      ind4C2I(5) = 10
      ind4C2I(6) = 12
c
      unit(1) = one
      unit(2) = one
      unit(3) = one
      unit(4) = zero
      unit(5) = zero
      unit(6) = zero
c
      do kb = 1,nb
        do k1=1,6
          aI(k1) = ainv(kb,indI(k1))
          Cg(k1) = ainv(kb,ind4C2I(k1))
          eg(k1) = half * (Cg(k1) - unit(k1))
        end do
c
        Qval = zero
        do k1=1,6
          do k2=1,6
            Qval = Qval + eg(k1)*bmx(k1,k2)*eg(k2)
          end do
        end do
        dUdQ = half * Cval * exp(Qval)
        d2UdQ2 = dUdQ
        ua(kb) = half * Cval * (exp(Qval) - one)
c
c     dQdC and d2QdC2
        do k1=1,6
          dQdC(k1) = zero
          do k2=1,6
            d2QdC2(k1,k2) = half * bmx(k1,k2)
            dQdC(k1) = dQdC(k1) + bmx(k1,k2)*eg(k2)
          end do
        end do
c
c     compute dCdI(6,6) and d2CdI2(6,6,6) 
c     not the best way to do.......
        do k3=1,6
          do k2=1,6
            dCdI(k2,k3) = zero
            do k1=1,6
              d2CdI2(k1,k2,k3) = zero
            end do
          end do
        end do
c
c     dC1dI(6)
        dCdI(1,1) = one / (two * Cg(1))
        dCdI(1,2) = - aI(2) / Cg(1)
        dCdI(1,4) = - aI(4) / Cg(1)
c     d2C1dI2(1,6,6)
        aux = - one / (two * Cg(1) * Cg(1))
        d2CdI2(1,1,1) = aux * dCdI(1,1)
        d2CdI2(1,1,2) = aux * dCdI(1,2)
        d2CdI2(1,1,4) = aux * dCdI(1,4)
        d2CdI2(1,2,1) = d2CdI2(1,1,2)
        d2CdI2(1,4,1) = d2CdI2(1,1,4)

        aux1 = - one / Cg(1)
        aux  = one / (Cg(1) * Cg(1))
        d2CdI2(1,2,2) = aux1 + aux * aI(2) * dCdI(1,2)
        d2CdI2(1,2,4) =        aux * aI(2) * dCdI(1,4)
        d2CdI2(1,4,2) = d2CdI2(1,2,4)
c
        d2CdI2(1,4,4) = aux1 + aux * aI(4) * dCdI(1,4)
c
c
c     dC2dI(6)
        dCdI(2,2) = - aI(2) / Cg(2)
        dCdI(2,3) = one / (two * Cg(2))
        dCdI(2,5) = - aI(5) / Cg(2)
c     d2C2dI(2,6,6)
        aux1 = - one / Cg(2)
        aux  = one / (Cg(2) * Cg(2))
        d2CdI2(2,2,2) = aux1 + aux * aI(2) * dCdI(2,2)
        d2CdI2(2,2,3) =        aux * aI(2) * dCdI(2,3)
        d2CdI2(2,2,5) =        aux * aI(2) * dCdI(2,5)
        d2CdI2(2,3,2) = d2CdI2(2,2,3)
        d2CdI2(2,5,2) = d2CdI2(2,2,5)

        d2CdI2(2,3,3) =  - half * aux * dCdI(2,3)
        d2CdI2(2,3,5) =  - half * aux * dCdI(2,5)
        d2CdI2(2,5,3) = d2CdI2(2,3,5)

        d2CdI2(2,5,5) = aux1 + aux * aI(5) * dCdI(2,5)
c
c
c     dC3dI(6)
        dCdI(3,4) = - aI(4) / Cg(3)
        dCdI(3,5) = - aI(5) / Cg(3)
        dCdI(3,6) = one / (two * Cg(3))
c     d2C3dI(3,6,6)
        aux1 = - one / Cg(3)
        aux  = one / (Cg(3) * Cg(3))
        d2CdI2(3,4,4) = aux1 + aux * aI(4) * dCdI(3,4)
        d2CdI2(3,4,5) =        aux * aI(4) * dCdI(3,5)
        d2CdI2(3,4,6) =        aux * aI(4) * dCdI(3,6)
        d2CdI2(3,5,4) = d2CdI2(3,4,5)
        d2CdI2(3,6,4) = d2CdI2(3,4,6)

        d2CdI2(3,5,5) = aux1 + aux * aI(5) * dCdI(3,5)
        d2CdI2(3,5,6) =        aux * aI(5) * dCdI(3,6)
        d2CdI2(3,6,5) = d2CdI2(3,5,6)

        d2CdI2(3,6,6) =   - half * aux * dCdI(3,6)
c
        dCdI(4,2) = one
        dCdI(5,4) = one
        dCdI(6,5) = one
c
        do k1=1,6
          dQdI(k1) = zero
          do k2=1,6
            dQdI(k1) = dQdI(k1) + dQdC(k2) * dCdI(k2,k1)
          end do
        end do
        call vaprd(d2QdC2,dCdI,d2QdCdI,6,6,6)
        do k1=1,6
          do k2=1,6
            d2QdI2(k1,k2) = zero
            do k3=1,6
              d2QdI2(k1,k2) = d2QdI2(k1,k2) 
     *             + dCdI(k3,k1)*d2QdCdI(k3,k2)
     *             + dQdC(k3)*d2CdI2(k3,k1,k2)
            end do
          end do
        end do
c
        do k1=1,6
          ui1(kb,indI(k1)) = dUdQ * dQdI(k1)
        end do
        do k1=1,6
          do k2=1,6
            ui2(kb,indx(indI(k1),indI(k2))) = d2UdQ2*dQdI(k1)*dQdI(k2)
     *           + dUdQ*d2QdI2(k1,k2)
          end do
        end do
      end do
*
      if(props(23).gt.zero) then
        do kb = 1,nb
          det = ainv(kb,indx_J)
          Dinv = one / props(23)
          ui1(kb,indx_J) = Dinv * (det - one/det)
          ui2(kb,indx(indx_J,indx_J)) = Dinv * (one + one / det / det)
        end do
      end if
c
      return
      end
C----------------------------------------------------------------------
      subroutine vaprd(A,B,C,n,m,k)
c
      include 'vaba_param.inc'
c
      parameter (zero = 0.d0)
      dimension A(n,m),B(m,k),C(n,k)
c
      do k1=1,n
        do k2=1,k
          C(k1,k2) = zero
          do k3=1,m
            C(k1,k2)=C(k1,k2)+A(k1,k3)*B(k3,k2)
          end do
        end do
      end do
c
      return
      end
C----------------------------------------------------------------------
      subroutine aTprd(A,B,C,n,m,k)
c
      include 'vaba_param.inc'
c
      parameter (zero = 0.d0)
      dimension A(n,m),B(m,k),C(n,k)
c
      do k1=1,n
         do k2=1,k
            C(k1,k2) = zero
            do k3=1,m
               C(k1,k2)=C(k1,k2)+A(k1,k3)*B(k2,k3)
            end do
         end do
      end do
c
      return
      end