C Infinitesimal damage model 
C by Ravi Kiran Yellavajjala
C Date: 12/23/2015
      subroutine vumat(
C Read only (unmodifiable)variables -
     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
     2  stepTime, totalTime, dt, cmname, coordMp, charLength,
     3  props, density, strainInc, relSpinInc,
     4  tempOld, stretchOld, defgradOld, fieldOld,
     5  stressOld, stateOld, enerInternOld, enerInelasOld,
     6  tempNew, stretchNew, defgradNew, fieldNew,
C Write only (modifiable) variables -
     7  stressNew, stateNew, enerInternNew, enerInelasNew )
C      include 'vaba_param.inc'
C
      dimension props(nprops), density(nblock), coordMp(nblock,*),
     1  charLength(nblock), strainInc(nblock,ndir+nshr),
     2  relSpinInc(nblock,nshr), tempOld(nblock),
     3  stretchOld(nblock,ndir+nshr),
     4  defgradOld(nblock,ndir+nshr+nshr),
     5  fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
     6  stateOld(nblock,nstatev), enerInternOld(nblock),
     7  enerInelasOld(nblock), tempNew(nblock),
     8  stretchNew(nblock,ndir+nshr),
     8  defgradNew(nblock,ndir+nshr+nshr),
     9  fieldNew(nblock,nfieldv),
     1  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     2  enerInternNew(nblock), enerInelasNew(nblock)
C
      character*80 cmname
C
      REAL(8) :: E_mod, v, sig_0, eps_e_tr(1,6), zz, H, Q, b
      REAL(8) :: Ra, p_strain, fac, dev_str(1,6), s_tri(1,6)
      REAL(8) :: sig_y, phi, pre, tot, kb, mu, vol_strain
      REAL(8) :: cons, K_tp1, K_tp2, sig_v, K_tp, Dt1, dg
      REAL(8) :: p_strain1, q_tri, q_tri1, dir(1,6)
      REAL(8) :: s_new(1,6), cr_strain, K_sig, plas_strain(1,6)
      REAL(8) :: trans_strain(1,6), C2, red_flag, hh
      REAL(8) :: H_a, H_m, Q_a, Q_m, A_con, B_con, sig_max, p_max 
      REAL(8) :: vol_exp, K_vol, tot_strain(1,6), sig_a, sig_m
      REAL(8) :: dil_strain(1,3), K_tv, tr_dil, aa_con, stress_thresold
      REAL(8) :: flag_stress, mmm, cons1, phidot, cons3, m, n1, dam
      REAL(8) :: d_dam, A1, B1, C1, D1, F1, F2, der(2,2), res(2,1), dam1
      REAL(8) :: max_iter, T_sig, det_der, der_inv(2,2), Xa(2,1), conv
      REAL(8) :: P1, P2, cri_dam, T1, lambda, lambda_cr, num, qtri_new 
      REAL(8) :: T_avg, E1, G1
      COMPLEX(8) :: dgc, damc, F1_dg, F1_dam, F2_dg, F2_dam, cons_c


      Parameter (Zero = 0.d0, One = 1.d0, Two = 2.d0, Three = 3.d0,
     1 Tole = 1.0D-8)
      hh=10.0d0**-12.0d0

      max_iter=30.0
C------material properties---------
C E_mod - Youngs modulus
C v - Poissons Ratio
      fac=sqrt(Three/Two)
      E_mod=200000.0d0
      v=0.30d0
	  

C--sig_a,sig_m - yield stress of austenite, martensite  
      sig_0=props(1)

C-- H, Q, b - isotropic hardening parameters
      H=props(2)
      Q=props(3)
      b=props(4)
C--- power hard mod
      n1=props(5)
C----Strain rate parameters
      phidot=props(6)
      m=props(7)
C-- cr_strain-strain at which damage initiates	  	  
      cr_strain=props(8)
      cri_dam=props(9)
      lambda_cr=props(10)
C----kb-bulk mod. and mu-shear mod.	  
      kb=E_mod/(3.0d0*(1.0d0-(2.0d0*v))) 
      mu=E_mod/(2.0d0*(1.0d0+v))

	  
      do i=1,nblock   
C------initializing state variables---------
C----stateOld(i,1:6)-elastic strain
C----stateOld(i,7)-p strain
C----stateOld(i,8)-martensitic vol fraction
C----stateOld(i,9)-old Von Mises stress
         tot=totalTime
	     if (tot.eq.0.0d0) then
              stateOld(i,1)=0.0d0
              stateOld(i,2)=0.0d0
              stateOld(i,3)=0.0d0
              stateOld(i,4)=0.0d0
              stateOld(i,5)=0.0d0
              stateOld(i,6)=0.0d0
			  
              stateOld(i,7)=0.0d0
              stateOld(i,8)=0.0d0
              stateOld(i,9)=0.0d0
              stateOld(i,10)=0.0d0
              stateOld(i,11)=0.0d0
              stateOld(i,12)=0.0d0
			  
              stateOld(i,13)=0.0d0
              stateOld(i,14)=0.0d0
			  
              stateOld(i,15)=1.0d0
              stateOld(i,16)=1
              stateOld(i,17)=0.0d0



         end if 
C----finding total strain-----

         tot_strain(1,1)=stateOld(i,1)+strainInc(i,1)
         tot_strain(1,2)=stateOld(i,2)+strainInc(i,2)
         tot_strain(1,3)=stateOld(i,3)+strainInc(i,3)
         tot_strain(1,4)=stateOld(i,4)+strainInc(i,4)
         tot_strain(1,5)=stateOld(i,5)+strainInc(i,5)
         tot_strain(1,6)=stateOld(i,6)+strainInc(i,6)
		 
C----finding elastic strain-----
         eps_e_tr(1,1)=tot_strain(1,1)-stateOld(i,7)
         eps_e_tr(1,2)=tot_strain(1,2)-stateOld(i,8)
         eps_e_tr(1,3)=tot_strain(1,3)-stateOld(i,9)
         eps_e_tr(1,4)=tot_strain(1,4)-stateOld(i,10)
         eps_e_tr(1,5)=tot_strain(1,5)-stateOld(i,11)
         eps_e_tr(1,6)=tot_strain(1,6)-stateOld(i,12)
		 
         p_strain=stateOld(i,13)
		 

		 
         dam=stateOld(i,14)
         lambda=stateOld(i,15)
         T_sig=stateOld(i,17)
         num=stateOld(i,18)

C-----volumetric component of eps_e_tr
         vol_strain=(eps_e_tr(1,1)+eps_e_tr(1,2)
     1 +eps_e_tr(1,3))
C-----deviatoric component of eps_e_tr--------------------
         dev_str(1,1)=eps_e_tr(1,1)-(One/Three)*vol_strain
         dev_str(1,2)=eps_e_tr(1,2)-(One/Three)*vol_strain
         dev_str(1,3)=eps_e_tr(1,3)-(One/Three)*vol_strain
         dev_str(1,4)=eps_e_tr(1,4)
         dev_str(1,5)=eps_e_tr(1,5)
         dev_str(1,6)=eps_e_tr(1,6)
		 
C-----Determination of pressure------------------------------
         pre=kb*(tot_strain(1,1)+tot_strain(1,2)+
     1 tot_strain(1,3))
C-----Determination of trial deviatoric stress (s_tri)--------
         s_tri(1,1)=Two*mu*dev_str(1,1)
         s_tri(1,2)=Two*mu*dev_str(1,2)
         s_tri(1,3)=Two*mu*dev_str(1,3)
         s_tri(1,4)=Two*mu*dev_str(1,4)
         s_tri(1,5)=Two*mu*dev_str(1,5)
         s_tri(1,6)=Two*mu*dev_str(1,6)
		 
C-----Determination of trial J2 stress-----------------
         q_tri=fac*sqrt(s_tri(1,1)**2.0d0+s_tri(1,2)**2.0d0
     1 +s_tri(1,3)**2.0d0+Two*(s_tri(1,4)**2.0d0+s_tri(1,5)**2.0d0
     2 +s_tri(1,6)**2.0d0))
         q_tri1=q_tri/(fac)
C-----direction vector
         dir=s_tri/((q_tri1))

C-----isotropic hardening law
         Ra=(H*p_strain**n1)+Q*(One-exp(-b*p_strain))
C         call mprint(Ra,1,1,1,'Ra',6)
C-----Determination of current yield strength
         sig_y=sig_0+Ra
         phi=((q_tri)-sig_y)*(1.0d0-dam)
C         call mprint(q_tri,1,1,1,'q_tri',6)
C         call mprint(phi,1,1,1,'phi',6)
      if (phi.LT.1.0D-6) then
C-----Update the state variables
             stateNew(i,1)=tot_strain(1,1)
             stateNew(i,2)=tot_strain(1,2)
             stateNew(i,3)=tot_strain(1,3)
             stateNew(i,4)=tot_strain(1,4)
             stateNew(i,5)=tot_strain(1,5)
             stateNew(i,6)=tot_strain(1,6)	
			 
             stateNew(i,7)=stateOld(i,7)
             stateNew(i,8)=stateOld(i,8)
             stateNew(i,9)=stateOld(i,9)
             stateNew(i,10)=stateOld(i,10)
             stateNew(i,11)=stateOld(i,11)
             stateNew(i,12)=stateOld(i,12)
			 
             stateNew(i,13)=stateOld(i,13)
             stateNew(i,14)=stateOld(i,14)
			 
             stateNew(i,15)=stateOld(i,15)
             stateNew(i,16)=stateOld(i,16)
			 stateNew(i,17)=(dabs(pre)/q_tri)

			 
             



C-----Update the stress
             stressNew(i,1)=(s_tri(1,1)+pre)
             stressNew(i,2)=(s_tri(1,2)+pre)
             stressNew(i,3)=(s_tri(1,3)+pre)
             stressNew(i,4)=s_tri(1,4)
             stressNew(i,5)=s_tri(1,5)
             stressNew(i,6)=s_tri(1,6)
         end if
		 
C------transformation & Plasticity step--------
         if (phi.GE.1.0D-6) then

		     red_flag=0.0d0
		     dg=0.0d0
             d_dam=0.0d0
             dam1=dam
             Dt1=dt

             conv=1.0d0
             
            
			
             A1=0.313d0
             B1=2.013d0
             G1=0.010d0
             C1=-1.93d0
             D1=3.8d0
			 
             P1=1.50d0
             P2=1.50d0
			 
             
             T1=0.0d0
             if (((dam1).GT.(cri_dam)).or.((lambda).GT.(lambda_cr))) 
     1 then
			     T1=1.0d0
             end if
			 

             do while ((conv.EQ.1.0d0).AND.(red_flag<max_iter))
			     red_flag=red_flag+1.0d0
				 
                 cons1=sig_0+(H*(p_strain+dg)**n1)+Q*(One-
     1 exp(-b*(p_strain+dg)))
                 F1=dg-Dt1*phidot*(((q_tri-Three*mu*dg/(1.0d0-dam1)
     1 )/(cons1)))**(1.0d0/m)
				 
                 E1=(A1*EXP(B1*T_sig))
                 F2=dam1-dam-dg*G1*E1*EXP(E1*T_sig)
     1 -dg*(P1*exp(P2*(p_strain+dg))*T1)
                 dgc=CMPLX(dg,hh)
                 damc=CMPLX(dam1,hh)
				 
                 cons_c=sig_0+(H*(p_strain+dgc)**n1)+Q*(One-
     1 exp(-b*(p_strain+dgc)))
                 F1_dg=dgc-Dt1*phidot*(((q_tri-Three*mu*dgc/(1.0d0-dam1)
     1 )/(cons_c)))**(1.0d0/m)
                 F1_dam=dg-Dt1*phidot*(((q_tri-Three*mu*dg/(1.0d0-damc)
     1 )/(cons1)))**(1.0d0/m)


                 F2_dg=dam1-dam-dgc*G1*E1*EXP(E1*T_sig)
     1 -dgc*(P1*exp(P2*(p_strain+dgc))*T1)
                 F2_dam=damc-dam-dg*G1*E1*EXP(E1*T_sig)
     1 -dg*(P1*exp(P2*(p_strain+dg))*T1)
	 

                 der(1,1)=AIMAG(F1_dg)/hh
                 der(1,2)=AIMAG(F1_dam)/hh
                 der(2,1)=AIMAG(F2_dg)/hh
                 der(2,2)=AIMAG(F2_dam)/hh

				 


                 res(1,1)=F1
                 res(2,1)=F2
				 
                 call two_inverse(der,der_inv)
                 det_der=der(1,1)*der(2,2)-der(1,2)*der(2,1)
                 if (isnan(det_der)) then
			         call mprint(iter,1,1,1,'isnan encountered',7)
				     call mprint(F1,1,1,1,'F1',7)
				     call mprint(F2,1,1,1,'F2',7)
				     call mprint(der,2,2,2,'der',7)
				     call mprint(dam1,1,1,1,'dam1',7)
				     call mprint(dg,1,1,1,'dg',7)
				     STOP
                 end if
                 Xa=matmul(der_inv, res)
                 dg=dg-Xa(1,1)
                 dam1=dam1-Xa(2,1)

				 
                 if ((dabs(Xa(1,1)).lt.tole*dabs(dg)).and.
     1 (dabs(F1).lt.tole*sig_0)) then
                     if ((dabs(Xa(2,1)).lt.tole*dabs(dam1)).and.
     1 (dabs(F2).lt.tole*sig_0)) then
	                     conv=0.0d0
                     end if
                 end if

             end do
             lambda=lambda+(C1*T_sig+D1)*dg



C------update eff plastic strain--------
             p_strain=p_strain+dg
			 
             if ((p_strain).LT.(cr_strain)) then
		 
                 dam1=0.0d0
             end if
			 


C-----updating plastic strain
C             plas_strain=(fac*dg*dir)/(1.0d0-dam1)
             plas_strain=(fac*dg*dir)
	 
C----updating the deviatoric stress
             s_new=(s_tri-(Two*mu*fac*dg*dir*(1.0d0-dam1)))
             q_tri_new=fac*sqrt(s_new(1,1)**2.0d0+s_new(1,2)**2.0d0
     1 +s_new(1,3)**2.0d0+Two*(s_new(1,4)**2.0d0+s_new(1,5)**2.0d0
     2 +s_new(1,6)**2.0d0))
C             T_avg=(T_sig*num+(dabs(pre)/q_tri_new))/(num+1.0d0)
             T_sig=(dabs(pre)/q_tri_new)
C-----Update the elastic strain

             stateNew(i,1)=tot_strain(1,1)
             stateNew(i,2)=tot_strain(1,2)
             stateNew(i,3)=tot_strain(1,3)
             stateNew(i,4)=tot_strain(1,4)
             stateNew(i,5)=tot_strain(1,5)
             stateNew(i,6)=tot_strain(1,6)	
			 
             stateNew(i,7)=stateOld(i,7)+plas_strain(1,1)
             stateNew(i,8)=stateOld(i,8)+plas_strain(1,2)
             stateNew(i,9)=stateOld(i,9)+plas_strain(1,3)
             stateNew(i,10)=stateOld(i,10)+plas_strain(1,4)
             stateNew(i,11)=stateOld(i,11)+plas_strain(1,5)
             stateNew(i,12)=stateOld(i,12)+plas_strain(1,6)
			 
             stateNew(i,13)=p_strain
             stateNew(i,14)=dam1
			 
             stateNew(i,15)=lambda
             stateNew(i,16)=stateOld(i,16)
			 
             if (((dam1).GT.(cri_dam)).or.((lambda).GT.(lambda_cr))) 
     1 then
		 
                 stateNew(i,16)=0
             end if

			 stateNew(i,17)=T_sig



C-----updating the stress-------
             stressNew(i,1)=s_new(1,1)+pre*(1.0d0-dam1)
             stressNew(i,2)=s_new(1,2)+pre*(1.0d0-dam1)
             stressNew(i,3)=s_new(1,3)+pre*(1.0d0-dam1)
             stressNew(i,4)=s_new(1,4)
             stressNew(i,5)=s_new(1,5)
             stressNew(i,6)=s_new(1,6)

			 
         end if 
      end do 
      return
      end subroutine VUMAT 
      
      
      SUBROUTINE two_inverse(jac,jac_inv)
      REAL(8) :: jac(2,2), jac1(2,2), jac_inv(2,2)	  
	  
      jac1(1,1)=jac(2,2)
      jac1(1,2)=-jac(1,2)
      jac1(2,1)=-jac(2,1)
      jac1(2,2)=jac(1,1)			
      jac_inv=(1.0d0/(jac(1,1)*jac(2,2)-jac(1,2)*jac(2,1)))*jac1
      end SUBROUTINE two_inverse
	  
      subroutine mprint(a,ii,jj,mm,name,iow)
c-----[--.----+----.----+----.-----------------------------------------]
c      Purpose: Output array of integer values

c      Inputs:
c         a(mm,*)  - Array to output
c         ii       - Number of rows to output
c         jj       - Number of columns to output
c         mm       - Dimension of array
c         name     - Name of array to appear with outputs

c      Outputs:
c         none
c-----[--.----+----.----+----.-----------------------------------------]
C      include 'vaba_param.inc'

      character name*(*), aname*30
      integer   ii,jj,mm,nn, ja,jb,j,i,n,irow,iow
      real*8    a(mm,*)

      
c     Print an ii x jj array whose dimension is mm for first subscript

      aname = name
      nn    = (jj+5)/6
      jb    = 0

      do n = 1,nn

        ja = jb + 1
        jb = min(jj,ja + 5)

c       Output to O file

        if(n.eq.1 .or. ii.gt.1) then
          write(iow,2000) aname,(j,j=ja,jb)
        endif
        do i = 1,ii
          if(ii.eq.1) then
            irow = n
          else
            irow = i
          endif
          write(iow,2001) irow,(a(i,j),j=ja,jb)
        end do ! i

      end do ! n

c     Flush buffer to force outputs (fix later)

c     Formats

 2000 format(/4x,'Matrix: ',a30/4x,'row/col',i6,5i12)
 2001 format(i8,1p,6e12.4)

      end
