subroutine vuanisohyper_inv (
C Read only -
     *     nblock, nFiber, nInv, 
     *     jElem, kIntPt, kLayer, kSecPt, 
     *     cmname,
     *     nstatev, nfieldv, nprops,
     *     props, tempOld, tempNew, fieldOld, fieldNew,
     *     stateOld, sInvariant, zeta,  
C     Write only -
     *     uDev, duDi, d2uDiDi, 
     *     stateNew )
         implicit real (a-h,o-z)
      parameter (j_sys_Dimension = 1)
      parameter( n_vec_Length = 136  )
      parameter( maxblk = n_vec_Length  )
C
      dimension props(nprops), 
     *  tempOld(nblock),
     *  fieldOld(nblock,nfieldv), 
     *  stateOld(nblock,nstatev), 
     *  tempNew(nblock),
     *  fieldNew(nblock,nfieldv),
     *  sInvariant(nblock,nInv),
     *  uDev(nblock), duDi(nblock,*), 
     *  d2uDiDi(nblock,*),
     *  stateNew(nblock,nstatev)
C
      character*80 cmname
C
      parameter ( zero = 0.d0, one = 1.d0, two = 2.d0, 
     *     three = 3.d0, four = 4.d0, five = 5.d0, six = 6.d0 )
C
C     Tissue Strain energy function (3D)
C
C     Read material properties
      nu=props(1)
      nu1=props(1)/two
	nu2=props(1)/three
      kappa=props(2)
	kappa1=props(2)/four
	kapp=props(2)/two
C
      do k = 1, nblock
         uDev(k) = kappa1 - (three/two)*nu
C     Compute Udev and 1st and 2nd derivatives w.r.t invariants
C     - I1
         bi1 = sInvariant(k,1)
      	bi3 = sInvariant(k,3)
             uDev(k) = uDev(k) 
     *        + nu1*bi1*bi3**(2/3)
         duDi(k,1) = nu1*bi3**(2/3)
         d2uDiDi(k,indx(1,1)) = zero
C     - I2
         bi2 = sInvariant(k,2)
         uDev(k) = uDev(k) 
     *        + zero
         d2uDiDi(k,indx(2,2)) = zero      

C     - I3 (=J)
              uDev(k) = uDev(k) 
     *        + kappa1*((bi3**2)-2*log(bi3))
         duDi(k,3) = nu2*bi1*bi3**(-1/3)-kappa1*(two*bi3-(two/bi3))
          d2uDiDi(k,indx(3,3)) = -1*nu2*(1/three)*bi1*bi3**(-4/3)
     *     -kapp-kapp/((bi3)**2)

C
      end do
C
      return
      end
C
C     Function to map index from Square to Triangular storage 
C     of symmetric matrix
C
      integer function indx( i, j )
               implicit real (a-h,o-z)
      parameter (j_sys_Dimension = 1)
      parameter( n_vec_Length = 136  )
      parameter( maxblk = n_vec_Length  )
      ii = min(i,j)
      jj = max(i,j)
      indx = ii + jj*(jj-1)/2
      return
      end