!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2007 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!     
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of 
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine umat_abaqusnl(amat,iel,iint,kode,elconloc,emec,emec0,
     &        beta,xokl,voj,xkl,vj,ithermal,t1l,dtime,time,ttime,
     &        icmd,ielas,mint_,
     &        nstate_,xstateini,xstate,stre,stiff,iorien,pgauss,orab)
!
!     calculates stiffness and stresses for a nonlinear material
!     defined by an ABAQUS umat routine
!
!     icmd=3: calcutates stress at mechanical strain
!     else: calculates stress at mechanical strain and the stiffness
!           matrix
!
!     INPUT:
!
!     amat               material name
!     iel                element number
!     iint               integration point number
!
!     kode               material type (-100-#of constants entered
!                        under *USER MATERIAL): can be used for materials
!                        with varying number of constants
!
!     elconloc(21)       user defined constants defined by the keyword
!                        card *USER MATERIAL (max. 21, actual # =
!                        -kode-100), interpolated for the
!                        actual temperature t1l
!
!     emec(6)            Lagrange mechanical strain tensor (component order:
!                        11,22,33,12,13,23) at the end of the increment
!                        (thermal strains are subtracted)
!     emec0(6)           Lagrange mechanical strain tensor at the start of the
!                        increment (thermal strains are subtracted)
!     beta(6)            residual stress tensor (the stress entered under
!                        the keyword *INITIAL CONDITIONS,TYPE=STRESS)
!
!     xokl(3,3)          deformation gradient at the start of the increment
!     voj                Jacobian at the start of the increment
!     xkl(3,3)           deformation gradient at the end of the increment
!     vj                 Jacobian at the end of the increment
!
!     ithermal           0: no thermal effects are taken into account
!                        1: thermal effects are taken into account (triggered
!                        by the keyword *INITIAL CONDITIONS,TYPE=TEMPERATURE)
!     t1l                temperature at the end of the increment
!     dtime              time length of the increment
!     time               step time at the end of the current increment
!     ttime              total time at the start of the current increment
!
!     icmd               not equal to 3: calculate stress and stiffness
!                        3: calculate only stress
!     ielas              0: no elastic iteration: irreversible effects
!                        are allowed
!                        1: elastic iteration, i.e. no irreversible
!                           deformation allowed
!
!     mint_              max. # of integration points per element in the
!                        model
!     nstate_            max. # of state variables in the model
!
!     xstateini(nstate_,mint_,# of elements)
!                        state variables at the start of the increment
!     xstate(nstate_,mint_,# of elements)
!                        state variables at the end of the increment
!
!     stre(6)            Piola-Kirchhoff stress of the second kind
!                        at the start of the increment
!
!     iorien             number of the local coordinate axis system
!                        in the integration point at stake (takes the value
!                        0 if no local system applies)
!     pgauss(3)          global coordinates of the integration point
!     orab(7,*)          description of all local coordinate systems.
!                        If a local coordinate system applies the global 
!                        tensors can be obtained by premultiplying the local
!                        tensors with skl(3,3). skl is  determined by calling
!                        the subroutine transformatrix: 
!                        call transformatrix(orab(1,iorien),pgauss,skl)
!
!
!     OUTPUT:
!
!     xstate(nstate_,mint_,# of elements)
!                        updated state variables at the end of the increment
!     stre(6)            Piola-Kirchhoff stress of the second kind at the
!                        end of the increment
!     stiff(21):         consistent tangent stiffness matrix in the material
!                        frame of reference at the end of the increment. In
!                        other words: the derivative of the PK2 stress with
!                        respect to the Lagrangian strain tensor. The matrix
!                        is supposed to be symmetric, only the upper half is
!                        to be given in the same order as for a fully
!                        anisotropic elastic material (*ELASTIC,TYPE=ANISO).
!
!     This routine allows for the use of an ABAQUS umat user subroutine
!     in CalculiX. 
!
!     Note that the following fields are not supported
!     so far: sse,spd,scd,rpl,ddsddt,drplde,drpldt,predef,
!     dpred,pnewdt,celent,layer,kspt,kstep,kinc
!
!     Furthermore, the following fields have a different meaning in
!     ABAQUS and CalculiX:
!
!     temp:  in CalculiX: temperature at the end of the increment
!              in ABAQUS: temperature at the start of the increment
!     dtemp: in CalculiX: zero
!              in ABAQUS: temperature increment
!
      implicit none
!
      character*80 amat
!
      integer ithermal,icmd,kode,ielas,iel,iint,nstate_,mint_,i,iorien,
     &  ndi,nshr,ntens,nprops,layer,kspt,kstep,kinc,kal(2,6),kel(4,21),
     &  j1,j2,j3,j4,j5,j6,j7,j8,jj,n,ier,j,matz
!
      real*8 elconloc(21),stiff(21),emec(6),emec0(6),beta(6),stre(6),
     &  vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*),
     &  time,ttime,skl(3,3),xa(3,3),ya(3,3,3,3),xstate(nstate_,mint_,*),
     &  xstateini(nstate_,mint_,*),w(3),fv1(3),fv2(3),d(6),
     &  v1,v2,v3,c(6),r(3,3),r0(3,3),eln0(6),eln(6),e(3,3),tkl(3,3),
     &  u(6),c2(6),dd,um1(3,3),z(3,3),u0(3,3)
!
      real*8 ddsdde(6,6),sse,spd,scd,rpl,ddsddt(6),drplde(6),
     &  drpldt,stran(6),dstran(6),abqtime(2),predef,temp,dtemp,
     &  dpred,drot(3,3),celent,pnewdt
!
      data kal /1,1,2,2,3,3,1,2,1,3,2,3/
!
      data kel /1,1,1,1,1,1,2,2,2,2,2,2,1,1,3,3,2,2,3,3,3,3,3,3,
     &          1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3,
     &          3,3,1,3,1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3,
     &          1,2,2,3,1,3,2,3,2,3,2,3/
!
      data d /1.d0,1.d0,1.d0,0.d0,0.d0,0.d0/
!
!     calculating the logarithmic mechanical strain at the
!     start of the increment
!
      e(1,1)=emec0(1)
      e(2,2)=emec0(4)
      e(3,3)=emec0(2)
      e(1,2)=emec0(5)
      e(1,3)=emec0(6)
      e(2,3)=emec0(3)
      e(2,1)=emec0(4)
      e(3,1)=emec0(5)
      e(3,2)=emec0(6)
!
!     calculating the eigenvalues and eigenvectors
!
      n=3
      matz=1
!
      call rs(n,n,e,w,matz,z,fv1,fv2,ier)
!
      if(ier.ne.0) then
         write(*,*) '
     &  *ERROR calculating the eigenvalues/vectors in umat_abaqusnl'
         stop
      endif
!
!     calculating the principal stretches at the start of the increment
!
      do i=1,3
         w(i)=dsqrt(2.d0*w(i)+1.d0)
      enddo
!
!     calculating the invariants at the start of the increment
!
      v1=w(1)+w(2)+w(3)
      v2=w(1)*w(2)+w(2)*w(3)+w(3)*w(1)
      v3=w(1)*w(2)*w(3)
!
!     calculating the right Cauchy-Green tensor at the start of the
!     increment
!
      do i=1,3
         c(i)=2.d0*emec0(i)+1.d0
      enddo
      do i=4,6
         c(i)=2.d0*emec0(i)
      enddo
!
!     calculating the square of the right Cauchy-Green tensor at the
!     start of the increment
!
      c2(1)=c(1)*c(1)+c(4)*c(4)+c(5)*c(5)
      c2(2)=c(4)*c(4)+c(2)*c(2)+c(6)*c(6)
      c2(3)=c(5)*c(5)+c(6)*c(6)+c(3)*c(3)
      c2(4)=c(1)*c(4)+c(4)*c(2)+c(5)*c(6)
      c2(5)=c(1)*c(5)+c(4)*c(6)+c(5)*c(3)
      c2(6)=c(4)*c(5)+c(2)*c(6)+c(6)*c(3)
!
!     calculating the right stretch tensor at the start of the increment
!     (cf. Simo and Hughes, Computational Inelasticity)
!
      dd=v1*v2-v3
      do i=1,6
         u(i)=(-c2(i)+(v1*v1-v2)*c(i)+v1*v3*d(i))/dd
      enddo
!
      u0(1,1)=u(1)
      u0(2,2)=u(2)
      u0(3,3)=u(3)
      u0(1,2)=u(4)
      u0(1,3)=u(5)
      u0(2,3)=u(6)
      u0(2,1)=u(4)
      u0(3,1)=u(5)
      u0(3,2)=u(6)
!
!     calculating the inverse of the right stretch tensor at the start
!     of the increment
!
      um1(1,1)=(c(1)-v1*u(1)+v2)/v3
      um1(2,2)=(c(2)-v1*u(2)+v2)/v3
      um1(3,3)=(c(3)-v1*u(3)+v2)/v3
      um1(1,2)=(c(4)-v1*u(4))/v3
      um1(1,3)=(c(5)-v1*u(5))/v3
      um1(2,3)=(c(6)-v1*u(6))/v3
      um1(2,1)=um1(1,2)
      um1(3,1)=um1(1,3)
      um1(3,2)=um1(2,3)
!
!     calculation of the local rotation tensor at the start of the
!     increment
!
      do i=1,3
         do j=1,3
            r0(i,j)=xokl(i,1)*um1(1,j)+xokl(i,2)*um1(2,j)+
     &              xokl(i,3)*um1(3,j)
         enddo
      enddo
!
!     calculating the logarithmic strain at the start of the increment
!
      do i=1,3
         w(i)=log(w(i))
      enddo
!
!     logarithmic strain in global coordinates at the start of the
!     increment
!
      eln0(1)=z(1,1)*z(1,1)*w(1)+z(1,2)*z(1,2)*w(2)+
     &        z(1,3)*z(1,3)*w(3)          
      eln0(2)=z(2,1)*z(2,1)*w(1)+z(2,2)*z(2,2)*w(2)+
     &        z(2,3)*z(2,3)*w(3)          
      eln0(3)=z(3,1)*z(3,1)*w(1)+z(3,2)*z(3,2)*w(2)+
     &        z(3,3)*z(3,3)*w(3)          
      eln0(4)=z(1,1)*z(2,1)*w(1)+z(1,2)*z(2,2)*w(2)+
     &        z(1,3)*z(2,3)*w(3)          
      eln0(5)=z(1,1)*z(3,1)*w(1)+z(1,2)*z(3,2)*w(2)+
     &        z(1,3)*z(3,3)*w(3)          
      eln0(6)=z(2,1)*z(3,1)*w(1)+z(2,2)*z(3,2)*w(2)+
     &        z(2,3)*z(3,3)*w(3)          
!
!     calculating the logarithmic mechanical strain at the
!     end of the increment
!
      e(1,1)=emec(1)
      e(2,2)=emec(4)
      e(3,3)=emec(2)
      e(1,2)=emec(5)
      e(1,3)=emec(6)
      e(2,3)=emec(3)
      e(2,1)=emec(4)
      e(3,1)=emec(5)
      e(3,2)=emec(6)
!
!     calculating the eigenvalues and eigenvectors
!
      call rs(n,n,e,w,matz,z,fv1,fv2,ier)
!
      if(ier.ne.0) then
         write(*,*) '
     &  *ERROR calculating the eigenvalues/vectors in umat_abaqusnl'
         stop
      endif
!
!     calculating the principal stretches at the end of the increment
!
      do i=1,3
         w(i)=dsqrt(2.d0*w(i)+1.d0)
      enddo
!
!     calculating the invariants at the end of the increment
!
      v1=w(1)+w(2)+w(3)
      v2=w(1)*w(2)+w(2)*w(3)+w(3)*w(1)
      v3=w(1)*w(2)*w(3)
!
!     calculating the right Cauchy-Green tensor at the end of the
!     increment
!
      do i=1,3
         c(i)=2.d0*emec0(i)+1.d0
      enddo
      do i=4,6
         c(i)=2.d0*emec0(i)
      enddo
!
!     calculating the square of the right Cauchy-Green tensor at the
!     end of the increment
!
      c2(1)=c(1)*c(1)+c(4)*c(4)+c(5)*c(5)
      c2(2)=c(4)*c(4)+c(2)*c(2)+c(6)*c(6)
      c2(3)=c(5)*c(5)+c(6)*c(6)+c(3)*c(3)
      c2(4)=c(1)*c(4)+c(4)*c(2)+c(5)*c(6)
      c2(5)=c(1)*c(5)+c(4)*c(6)+c(5)*c(3)
      c2(6)=c(4)*c(5)+c(2)*c(6)+c(6)*c(3)
!
!     calculating the right stretch tensor at the end of the increment
!     (cf. Simo and Hughes, Computational Inelasticity)
!
      dd=v1*v2-v3
      do i=1,6
         u(i)=(-c2(i)+(v1*v1-v2)*c(i)+v1*v3*d(i))/dd
      enddo
!
!     calculating the inverse of the right stretch tensor at the end
!     of the increment
!
      um1(1,1)=(c(1)-v1*u(1)+v2)/v3
      um1(2,2)=(c(2)-v1*u(2)+v2)/v3
      um1(3,3)=(c(3)-v1*u(3)+v2)/v3
      um1(1,2)=(c(4)-v1*u(4))/v3
      um1(1,3)=(c(5)-v1*u(5))/v3
      um1(2,3)=(c(6)-v1*u(6))/v3
      um1(2,1)=um1(1,2)
      um1(3,1)=um1(1,3)
      um1(3,2)=um1(2,3)
!
!     calculation of the local rotation tensor at the end of the
!     increment
!
      do i=1,3
         do j=1,3
            r(i,j)=xokl(i,1)*um1(1,j)+xokl(i,2)*um1(2,j)+
     &              xokl(i,3)*um1(3,j)
         enddo
      enddo
!
!     calculating the logarithmic strain at the end of the increment
!     Elog=Z.ln(w).Z^T
!
      do i=1,3
         w(i)=log(w(i))
      enddo
!
!     logarithmic strain in global coordinates at the end of the
!     increment
!
      eln(1)=z(1,1)*z(1,1)*w(1)+z(1,2)*z(1,2)*w(2)+
     &        z(1,3)*z(1,3)*w(3)          
      eln(2)=z(2,1)*z(2,1)*w(1)+z(2,2)*z(2,2)*w(2)+
     &        z(2,3)*z(2,3)*w(3)          
      eln(3)=z(3,1)*z(3,1)*w(1)+z(3,2)*z(3,2)*w(2)+
     &        z(3,3)*z(3,3)*w(3)          
      eln(4)=z(1,1)*z(2,1)*w(1)+z(1,2)*z(2,2)*w(2)+
     &        z(1,3)*z(2,3)*w(3)          
      eln(5)=z(1,1)*z(3,1)*w(1)+z(1,2)*z(3,2)*w(2)+
     &        z(1,3)*z(3,3)*w(3)          
      eln(6)=z(2,1)*z(3,1)*w(1)+z(2,2)*z(3,2)*w(2)+
     &        z(2,3)*z(3,3)*w(3)  
c      write(*,*) 'iel', iel
c      write(*,*) 'emec',(emec(i),i=1,6)
c      write(*,*) 'eln',(eln(i),i=1,6)
c      write(*,*) 'r0',((r0(i,j),j=1,3),i=1,3)
c      write(*,*) 'r',((r(i,j),j=1,3),i=1,3)
!
!     calculating the incremental rotation tensor
!     drot=r.r0^T
!
      do i=1,3
         do j=1,3
            drot(i,j)=r(i,1)*r0(j,1)+r(i,2)*r0(j,2)+r(i,3)*r0(j,3)
         enddo
      enddo
!
      ntens=6
!
      do i=1,nstate_
         xstate(i,iint,iel)=xstateini(i,iint,iel)
      enddo
!
      abqtime(1)=time-dtime
      abqtime(2)=ttime-dtime
!
      temp=t1l
      dtemp=0.d0
!
      ndi=3
      nshr=3
      ntens=ndi+nshr
!
      nprops=-kode-100
!
!     taking local material orientations into account
!
      if(iorien.ne.0) then
         call transformatrix(orab(1,iorien),pgauss,skl)
!
!        rotating the strain at the start of the increment
!        into the local system: Elog'=T^T.Elog.T
!
         xa(1,1)=eln0(1)
         xa(1,2)=eln0(4)
         xa(1,3)=eln0(5)
         xa(2,1)=eln0(4)
         xa(2,2)=eln0(2)
         xa(2,3)=eln0(6)
         xa(3,1)=eln0(5)
         xa(3,2)=eln0(6)
         xa(3,3)=eln0(3)
!
         do jj=1,6
            stran(jj)=0.d0
            j1=kal(1,jj)
            j2=kal(2,jj)
            do j3=1,3
               do j4=1,3
                  stran(jj)=stran(jj)+
     &                 xa(j3,j4)*skl(j3,j1)*skl(j4,j2)
               enddo
            enddo
         enddo
!
!        rotating the strain at the end of the increment
!        into the local system
!
         xa(1,1)=eln(1)
         xa(1,2)=eln(4)
         xa(1,3)=eln(5)
         xa(2,1)=eln(4)
         xa(2,2)=eln(2)
         xa(2,3)=eln(6)
         xa(3,1)=eln(5)
         xa(3,2)=eln(6)
         xa(3,3)=eln(3)
!
         do jj=1,6
            dstran(jj)=-stran(jj)
            j1=kal(1,jj)
            j2=kal(2,jj)
            do j3=1,3
               do j4=1,3
                  dstran(jj)=dstran(jj)+
     &                 xa(j3,j4)*skl(j3,j1)*skl(j4,j2)
               enddo
            enddo
         enddo
      else
         do jj=1,6
            stran(jj)=eln0(jj)
            dstran(jj)=eln(jj)-eln0(jj)
         enddo
      endif
!
!     rotating the stress into the local system
!     s'=J^(-1).U.S.U^T (no orientation card) or
!     s'=J^(-1).U.T^T.S.T.U^T (orientation card)
!
      if(iorien.ne.0) then
         do i=1,3
            do j=1,3
               tkl(i,j)=u0(i,1)*skl(j,1)+u0(i,2)*skl(j,2)+
     &                  u0(i,3)*skl(j,3)
            enddo
         enddo
      else
         do i=1,3
            do j=1,3
               tkl(i,j)=u0(i,j)
            enddo
         enddo
      endif
!
      xa(1,1)=stre(1)
      xa(1,2)=stre(4)
      xa(1,3)=stre(5)
      xa(2,1)=stre(4)
      xa(2,2)=stre(2)
      xa(2,3)=stre(6)
      xa(3,1)=stre(5)
      xa(3,2)=stre(6)
      xa(3,3)=stre(3)
!     
      do jj=1,6
         stre(jj)=0.d0
         j1=kal(1,jj)
         j2=kal(2,jj)
         do j3=1,3
            do j4=1,3
               stre(jj)=stre(jj)+
     &              xa(j3,j4)*tkl(j1,j3)*tkl(j2,j4)
            enddo
         enddo
         stre(jj)=stre(jj)/voj
      enddo
!
      call umat(stre,xstate(1,iint,iel),ddsdde,sse,spd,scd,rpl,ddsddt,
     &  drplde,drpldt,stran,dstran,abqtime,dtime,temp,dtemp,predef,
     &  dpred,amat,ndi,nshr,ntens,nstate_,elconloc,nprops,pgauss,drot,
     &  pnewdt,celent,xokl,xkl,iel,iint,layer,kspt,kstep,kinc)
!
!     rotating the stress into the global system
!     S=J.U^(-1).s'.U^(-T)  (no orientation card) or
!     S=J.T.U^(-1).s'.U^(-T).T^T (orientation card)
!
      if(iorien.ne.0) then
         do i=1,3
            do j=1,3
               tkl(i,j)=skl(i,1)*um1(1,j)+skl(i,2)*um1(2,j)+
     &                  skl(i,3)*um1(3,j)
            enddo
         enddo
      else
         do i=1,3
            do j=1,3
               tkl(i,j)=um1(i,j)
            enddo
         enddo
      endif
!
      xa(1,1)=stre(1)
      xa(1,2)=stre(4)
      xa(1,3)=stre(5)
      xa(2,1)=stre(4)
      xa(2,2)=stre(2)
      xa(2,3)=stre(6)
      xa(3,1)=stre(5)
      xa(3,2)=stre(6)
      xa(3,3)=stre(3)
!
      do jj=1,6
         stre(jj)=0.d0
         j1=kal(1,jj)
         j2=kal(2,jj)
         do j3=1,3
            do j4=1,3
               stre(jj)=stre(jj)+
     &              xa(j3,j4)*tkl(j1,j3)*tkl(j2,j4)
            enddo
         enddo
         stre(jj)=stre(jj)*vj
      enddo
!
!     calculate the stiffness matrix
!
      if(icmd.ne.3) then
         stiff(1)=ddsdde(1,1)
         stiff(2)=ddsdde(1,2)
         stiff(3)=ddsdde(2,2)
         stiff(4)=ddsdde(1,3)
         stiff(5)=ddsdde(2,3)
         stiff(6)=ddsdde(3,3)
         stiff(7)=ddsdde(1,4)
         stiff(8)=ddsdde(2,4)
         stiff(9)=ddsdde(3,4)
         stiff(10)=ddsdde(4,4)
         stiff(11)=ddsdde(1,5)
         stiff(12)=ddsdde(2,5)
         stiff(13)=ddsdde(3,5)
         stiff(14)=ddsdde(4,5)
         stiff(15)=ddsdde(5,5)
         stiff(16)=ddsdde(1,6)
         stiff(17)=ddsdde(2,6)
         stiff(18)=ddsdde(3,6)
         stiff(19)=ddsdde(4,6)
         stiff(20)=ddsdde(5,6)
         stiff(21)=ddsdde(6,6)
!
!        rotating the stiffness coefficients into the global system
!
         call anisotropic(stiff,ya)
!
         do jj=1,21
            j1=kel(1,jj)
            j2=kel(2,jj)
            j3=kel(3,jj)
            j4=kel(4,jj)
            stiff(jj)=0.d0
            do j5=1,3
               do j6=1,3
                  do j7=1,3
                     do j8=1,3
                        stiff(jj)=stiff(jj)+ya(j5,j6,j7,j8)*
     &                       tkl(j1,j5)*tkl(j2,j6)*tkl(j3,j7)*tkl(j4,j8)
                     enddo
                  enddo
               enddo
            enddo
         enddo
      endif
!
      return
      end