c*****************************************************************
      subroutine uelmat(rhs,amatrx,svars,energy,ndofel,nrhs,nsvars,
     1  props, nprops, coords, mcrd, nnode, u, du, v, a, jtype, time,
     2  dtime, kstep, kinc, jelem, params, ndload, jdltyp, adlmag,
     3  predef, npredf, lflags, mlvarx, ddlmag, mdload, pnewdt, jprops, 
     4  njpro, period,materiallib)
c
      include 'aba_param.inc'
C
      dimension rhs(mlvarx,*), amatrx(ndofel, ndofel), props(*),
     1  svars(*), energy(*), coords(mcrd, nnode), u(ndofel),
     2  du(mlvarx,*), v(ndofel), a(ndofel), time(2), params(*),
     3  jdltyp(mdload,*), adlmag(mdload,*), ddlmag(mdload,*),
     4  predef(2, npredf, nnode), lflags(*), jprops(*)
      parameter (zero=0.d0, one=1.d0)

      parameter (ndim=2, ndof=2, ndi=3, nshr=1, nnodemax=4,
     1     ntens=4, ninpt=4, nsvint=8)
c
c      ndim  ... number of spatial dimensions
c      ndof  ... number of degrees of freedom per node
c      ndi   ... number of direct stress components
c      nshr  ... number of shear stress component
c      ntens ... total number of stress tensor components
c                (=ndi+nshr)
c      ninpt ... number of integration points
c      nsvint... number of state variables per integration pt 
c                (including stress and strain)
c
      dimension  stiff(ndof*nnodemax,ndof*nnodemax),
     1 force(ndof*nnodemax), shape(nnodemax), dshape(ndim,nnodemax),
     2 xjaci(ndim,ndim), bmat(nnodemax*ndim), statevLocal(nsvint),
     3 stress(ntens), ddsdde(ntens, ntens),
     4 stran(ntens), dstran(ntens), wght(ninpt)
  
      dimension coords_ip(3),dfgrd0(3,3),dfgrd1(3,3),
     1 drot(3,3), ddsddt(ntens), drplde(ntens)
      dimension coords_new(mcrd,nnodemax)
c
      dimension predef_loc(npredf),dpredef_loc(npredf),xx1(3,3),
     *     xx1Old(3,3)
      dimension xjaci_new(ndim,ndim),bmat_new(nnodemax*ndim)
c
      data wght /one, one, one, one/
c
c********************************************************************
c
c     U1 = first-order, plane strain, full integration
c
c     State variables: each integration point has nsvint SDVs
c
c       isvinc=(npt-1)*nsvint    ... integration point counter
c       statev(1+isvinc        ) ... stress         
c       statev(1+isvinc+  ntens) ... strain         
c
c********************************************************************
      if (lflags(3).eq.4) then
        do i=1, ntens
          do j=1, ntens
            ddsdde(i,j) = zero
          end do
          ddsdde(i,j) = one
        enddo
        do i=1, ndofel
          do j=1, ndofel
            amatrx(i,j) = zero
          end do
          amatrx(i,i) = one
        end do
        goto 999
      end if
c
c     Preliminaries
c
      pnewdtLocal = pnewdt
      if(jtype .ne. 1) then
        write(7,*)'Incorrect element type'
        call xit
      endif 
      if(nsvars .lt. ninpt*nsvint) then
        write(7,*)'Increase the number of SDVs to', ninpt*nsvint
        call xit
      endif 
c
c     initialize rhs and lhs
c
      do k1=1, ndof*nnode
        rhs(k1, 1)= zero
        do k2=1, ndof*nnode 
          amatrx(k1, k2)= zero
        end do
      end do
c
      do k1=1,nnode
         do k2=1,mcrd
            kk = (k1-1)*mcrd + k2
            coords_new(k2,k1) = coords(k2,k1) + u(kk)
         end do
      end do

c
c     loop over integration points
c

      do kintk = 1, ninpt
c
c       evaluate shape functions and derivatives
c
        call shapefcn(kintk,ninpt,nnode,ndim,shape,dshape)
c
c       compute coordinates at the integration point
c
        do k1=1, 3
          coords_ip(k1) = zero
        end do
        do k1=1,nnode
          do k2=1,mcrd
            if (lflags(2).eq.0) then
              coords_ip(k2)=coords_ip(k2)+shape(k1)*coords(k2,k1)
            else
              coords_ip(k2)=coords_ip(k2)+shape(k1)*coords_new(k2,k1)
            end if
          end do
        end do       
c
        if(npredf.gt.0) then
           call tempfv(kintk,ninpt,nnode,ndim,shape,predef,
     *         npredf,predef_loc,dpredef_loc)
        end if
c
c       form B-matrix
c
        djac = one
        djac_new = one
        call jacobian(jelem,mcrd,ndim,nnode,coords,dshape,
     1   djac,xjaci,pnewdt,coords_new,xjaci_new,djac_new)
        if (pnewdt .lt. pnewdtLocal) pnewdtLocal = pnewdt
c
        call bmatrix(xjaci,dshape,nnode,ndim,bmat,xjaci_new,
     *       bmat_new)
c
c       calculate incremental strains
c 
        call straininc(ntens,ndof,ndim,nnode,mlvarx,bmat,du,dstran,
     *     u,xx1,xx1Old)
c        
c       call constitutive routine
c
        call statevar(kintk,nsvint,svars,statevLocal,1)
c
c       material routine for computing state variables
c
        sse = zero
        do k1=1,ntens
           stress(k1) = statevLocal(k1)
           stran (k1) = statevLocal(k1+ntens)
        end do

        dvmat = zero
        celent = one
        call material_lib_mech(materiallib,stress,ddsdde,stran,dstran,
     *       kintk,djac,dvmat,xx1,predef_loc,dpredef_loc,npredf,celent,
     *       coords_ip)
c
        do k1=1,ntens
           statevLocal(k1) = stress(k1)
           statevLocal(ntens+k1) = stran(k1) + dstran(k1)
        end do

c
        call statevar(kintk,nsvint,svars,statevLocal,0)
c        
c       form stiffness matrix and internal force vector
c
        if(lflags(2).eq.0) then
           call stiffmatrix(ntens,nnode,ndim,ndof,
     1          wght(kintk),djac,ddsdde,stress,bmat,stiff,force)
        else
           call stiffmatrix(ntens,nnode,ndim,ndof,
     1          wght(kintk),djac_new,ddsdde,stress,bmat_new,stiff,force)
        end if
c
c       assemble rhs and lhs
c
        do k1=1, ndof*nnode
            rhs(k1, 1) = rhs(k1, 1) - force(k1) 
          do k2=1, ndof*nnode
            amatrx(k1, k2) = amatrx(k1, k2) + stiff(k1,k2)
          end do
        end do
      end do       ! end loop on material integration points
      pnewdt = pnewdtLocal
c
 999  continue
c
      return
      end
c*****************************************************************
c
c shape functions
c
      subroutine shapefcn(kintk,ninpt,nnode,ndim,dN,dNdz)
c
      include 'aba_param.inc'
c
      parameter (dmone=-1.0d0,one=1.0d0,four=4.0d0,eight=8.0d0,
     *     gaussCoord=0.577350269d0)
      parameter (maxElemNode=8,maxDof=3,i2d4node=24,i3d8node=38)
      dimension dN(*),dNdz(ndim,*),coord24(2,4),coord38(3,8)
c     
      data  coord24 /dmone, dmone,
     2                 one, dmone,
     3                 one,  one,
     4               dmone,  one/
c     
      data  coord38 /dmone, dmone, dmone,
     2                 one, dmone, dmone,
     3                 one,   one, dmone,
     4               dmone,   one, dmone,
     5               dmone, dmone,   one,
     6                 one, dmone,   one,
     7                 one,   one,   one,
     8               dmone,   one,   one/
C
      iCode = 0
      if (ninpt.eq.4.and.nnode.eq.4.and.ndim.eq.2) then
        iCode = 24
      else if (ninpt.eq.8.and.nnode.eq.8.and.ndim.eq.3) then
        iCode = 38
      else
        write (6,*) '***ERROR: The shape fuctions cannot be found'
      end if
C
C  3D 8-nodes
C
      if (iCode.eq.i3d8node) then
c
c       determine (g,h,r)
c
        g = coord38(1,kintk)*gaussCoord
        h = coord38(2,kintk)*gaussCoord
        r = coord38(3,kintk)*gaussCoord
c
c       shape functions
        dN(1) = (one - g)*(one - h)*(one - r)/eight
        dN(2) = (one + g)*(one - h)*(one - r)/eight
        dN(3) = (one + g)*(one + h)*(one - r)/eight
        dN(4) = (one - g)*(one + h)*(one - r)/eight
        dN(5) = (one - g)*(one - h)*(one + r)/eight
        dN(6) = (one + g)*(one - h)*(one + r)/eight
        dN(7) = (one + g)*(one + h)*(one + r)/eight
        dN(8) = (one - g)*(one + h)*(one + r)/eight
c     
c       derivative d(Ni)/d(g)
        dNdz(1,1) = -(one - h)*(one - r)/eight
        dNdz(1,2) =  (one - h)*(one - r)/eight
        dNdz(1,3) =  (one + h)*(one - r)/eight
        dNdz(1,4) = -(one + h)*(one - r)/eight
        dNdz(1,5) = -(one - h)*(one + r)/eight
        dNdz(1,6) =  (one - h)*(one + r)/eight
        dNdz(1,7) =  (one + h)*(one + r)/eight
        dNdz(1,8) = -(one + h)*(one + r)/eight
c
c       derivative d(Ni)/d(h)
        dNdz(2,1) = -(one - g)*(one - r)/eight
        dNdz(2,2) = -(one + g)*(one - r)/eight
        dNdz(2,3) =  (one + g)*(one - r)/eight
        dNdz(2,4) =  (one - g)*(one - r)/eight
        dNdz(2,5) = -(one - g)*(one + r)/eight
        dNdz(2,6) = -(one + g)*(one + r)/eight
        dNdz(2,7) =  (one + g)*(one + r)/eight
        dNdz(2,8) =  (one - g)*(one + r)/eight
c
c       derivative d(Ni)/d(r)
        dNdz(3,1) = -(one - g)*(one - h)/eight
        dNdz(3,2) = -(one + g)*(one - h)/eight
        dNdz(3,3) = -(one + g)*(one + h)/eight
        dNdz(3,4) = -(one - g)*(one + h)/eight
        dNdz(3,5) = (one - g)*(one - h)/eight
        dNdz(3,6) = (one + g)*(one - h)/eight
        dNdz(3,7) = (one + g)*(one + h)/eight
        dNdz(3,8) = (one - g)*(one + h)/eight
C     
C  2D 4-nodes
C
      else if (iCode.eq.i2d4node) then
c
c         determine (g,h)
c
        g = coord24(1,kintk)*gaussCoord
        h = coord24(2,kintk)*gaussCoord
c
c       shape functions
        dN(1) = (one - g)*(one - h)/four;
        dN(2) = (one + g)*(one - h)/four;
        dN(3) = (one + g)*(one + h)/four;
        dN(4) = (one - g)*(one + h)/four;
c
c       derivative d(Ni)/d(g)
        dNdz(1,1) = -(one - h)/four;
        dNdz(1,2) =  (one - h)/four;
        dNdz(1,3) =  (one + h)/four;
        dNdz(1,4) = -(one + h)/four;
c
c       derivative d(Ni)/d(h)
        dNdz(2,1) = -(one - g)/four;
        dNdz(2,2) = -(one + g)/four;
        dNdz(2,3) =  (one + g)/four;
        dNdz(2,4) =  (one - g)/four;
      end if
c      
      return
      end
c********************************************************
      subroutine tempfv(kintk,ninpt,nnode,ndim,shape,predef,
     *         npredf,predef_loc,dpredef_loc)
c
      include 'aba_param.inc'
c
      dimension shape(nnode),predef(2,npredf,nnode)
      dimension predef_loc(npredf),dpredef_loc(npredf)
      parameter (zero=0.d0)
c
      do k1=1,npredf
         predef_loc(k1) = zero
         dpredef_loc(k1) = zero
         do k2=1,nnode
             predef_loc(k1) =  predef_loc(k1)+
     &          (predef(1,k1,k2)-predef(2,k1,k2))*shape(k2)
            dpredef_loc(k1) = dpredef_loc(k1)+predef(2,k1,k2)*shape(k2)
         end do
      end do
c
      return
      end

c*****************************************************************
      subroutine jacobian(jelem,mcrd,ndim,nnode,
     1 coords,dshape,djac,xjaci,pnewdt,coords_new,xjaci_new,
     2 djac_new)
c
c     Notation:  ndim ....... element dimension
c                nnode  ..... number of nodes
c                coords ..... coordinates of nodes
c                dshape ..... derivs of shape fcn
c                djac ....... determinant of Jacobian
c                xjaci ...... inverse of Jacobian matrix
c
      include 'aba_param.inc'
      parameter(zero=0.d0, fourth=0.25d0, maxDof=3)

      dimension  xjac(maxDof,maxDof), xjaci(ndim,*), coords(mcrd,*)
      dimension  dshape(ndim,*),coords_new(mcrd,*)
      dimension  xjac_new(maxDof,maxDof), xjaci_new(ndim,*)

c
      do i = 1, ndim
        do j = 1, ndim
           xjac(i,j)  = zero
           xjaci(i,j) = zero
           xjac_new(i,j)  = zero
           xjaci_new(i,j) = zero
        end do
      end do
c
      do inod= 1, nnode
          do idim = 1, ndim
            do jdim = 1, ndim
              xjac(jdim,idim) = xjac(jdim,idim) + 
     1        dshape(jdim,inod)*coords(idim,inod)    
    
              xjac_new(jdim,idim) = xjac_new(jdim,idim) + 
     1        dshape(jdim,inod)*coords_new(idim,inod)        
            end do
          end do 
      end do
C
C  ndim == 3
C
      if (ndim.eq.3) then
        djac = xjac(1,1)*xjac(2,2)*xjac(3,3) +
     &       xjac(2,1)*xjac(3,2)*xjac(1,3) +
     &       xjac(3,1)*xjac(2,3)*xjac(1,2) -
     &       xjac(3,1)*xjac(2,2)*xjac(1,3) -
     &       xjac(2,1)*xjac(1,2)*xjac(3,3) -
     &       xjac(1,1)*xjac(2,3)*xjac(3,2)
        djac_new = xjac_new(1,1)*xjac_new(2,2)*xjac_new(3,3) +
     &       xjac_new(2,1)*xjac_new(3,2)*xjac_new(1,3) +
     &       xjac_new(3,1)*xjac_new(2,3)*xjac_new(1,2) -
     &       xjac_new(3,1)*xjac_new(2,2)*xjac_new(1,3) -
     &       xjac_new(2,1)*xjac_new(1,2)*xjac_new(3,3) -
     &       xjac_new(1,1)*xjac_new(2,3)*xjac_new(3,2)
        if (djac .gt. zero) then
          ! jacobian is positive - o.k.
          xjaci(1,1) = (xjac(2,2)*xjac(3,3)-xjac(2,3)*xjac(3,2))/djac
          xjaci(1,2) = (xjac(1,3)*xjac(3,2)-xjac(1,2)*xjac(3,3))/djac
          xjaci(1,3) = (xjac(1,2)*xjac(2,3)-xjac(1,3)*xjac(2,2))/djac
          ! 
          xjaci(2,1) = (xjac(2,3)*xjac(3,1)-xjac(2,1)*xjac(3,3))/djac
          xjaci(2,2) = (xjac(1,1)*xjac(3,3)-xjac(1,3)*xjac(3,1))/djac
          xjaci(2,3) = (xjac(1,3)*xjac(2,1)-xjac(1,1)*xjac(2,3))/djac
          !
          xjaci(3,1) = (xjac(2,1)*xjac(3,2)-xjac(2,2)*xjac(3,1))/djac
          xjaci(3,2) = (xjac(1,2)*xjac(3,1)-xjac(1,1)*xjac(3,2))/djac
          xjaci(3,3) = (xjac(1,1)*xjac(2,2)-xjac(1,2)*xjac(2,1))/djac
        else
          ! negative or zero jacobian
          write(7,*)'WARNING: element',jelem,'has neg. 
     1               Jacobian'
          pnewdt = fourth
        endif
        
        if (djac_new .gt. zero) then
          ! jacobian is positive - o.k.
          xjaci_new(1,1) = (xjac_new(2,2)*xjac_new(3,3)-
     &         xjac_new(2,3)*xjac_new(3,2))/djac_new
          xjaci_new(1,2) = (xjac_new(1,3)*xjac_new(3,2)-
     &         xjac_new(1,2)*xjac_new(3,3))/djac_new
          xjaci_new(1,3) = (xjac_new(1,2)*xjac_new(2,3)-
     &         xjac_new(1,3)*xjac_new(2,2))/djac_new
          ! 
          xjaci_new(2,1) = (xjac_new(2,3)*xjac_new(3,1)-
     &         xjac_new(2,1)*xjac_new(3,3))/djac_new
          xjaci_new(2,2) = (xjac_new(1,1)*xjac_new(3,3)-
     &         xjac_new(1,3)*xjac_new(3,1))/djac_new
          xjaci_new(2,3) = (xjac_new(1,3)*xjac_new(2,1)-
     &         xjac_new(1,1)*xjac_new(2,3))/djac_new
          !
          xjaci_new(3,1) = (xjac_new(2,1)*xjac_new(3,2)-
     &         xjac_new(2,2)*xjac_new(3,1))/djac_new
          xjaci_new(3,2) = (xjac_new(1,2)*xjac_new(3,1)-
     &         xjac_new(1,1)*xjac_new(3,2))/djac_new
          xjaci_new(3,3) = (xjac_new(1,1)*xjac_new(2,2)-
     &         xjac_new(1,2)*xjac_new(2,1))/djac_new
        else
          ! negative or zero jacobian
          write(7,*)'WARNING: element',jelem,'has neg. 
     1         Jacobian'
          pnewdt = fourth
        endif
C
C  ndim == 2
C
      else if (ndim.eq.2) then
        djac = xjac(1,1)*xjac(2,2) - xjac(1,2)*xjac(2,1)
        djac_new =   xjac_new(1,1)*xjac_new(2,2)
     *       - xjac_new(1,2)*xjac_new(2,1)
        if (djac .gt. zero) then
          ! jacobian is positive - o.k.
          xjaci(1,1) =  xjac(2,2)/djac
          xjaci(2,2) =  xjac(1,1)/djac
          xjaci(1,2) = -xjac(1,2)/djac
          xjaci(2,1) = -xjac(2,1)/djac
        else
          ! negative or zero jacobian
          write(7,*)'WARNING: element',jelem,'has neg. Jacobian'
          pnewdt = fourth
        endif
        
        if (djac_new .gt. zero) then
          ! jacobian is positive - o.k.
          xjaci_new(1,1) =  xjac_new(2,2)/djac_new
          xjaci_new(2,2) =  xjac_new(1,1)/djac_new
          xjaci_new(1,2) = -xjac_new(1,2)/djac_new
          xjaci_new(2,1) = -xjac_new(2,1)/djac_new
        else
          ! negative or zero jacobian
          write(7,*)'WARNING: element',jelem,'has neg. Jacobian'
          pnewdt = fourth
        endif
      end if

      return
      end
c*****************************************************************
      subroutine bmatrix(xjaci,dshape,nnode,ndim,bmat,
     *       xjaci_new,bmat_new)
c
c     Notation:  
c                bmat(i) .....dN1/dx, dN1/dy, dN2/dx, dN2/dy..
c                xjaci ...... inverse Jabobian matrix
c                dshape ......derivative of shape functions
c
      include 'aba_param.inc'
      parameter (zero=0.d0)
      dimension bmat(*),  dshape(ndim,*)
      dimension xjaci(ndim,*)
      dimension xjaci_new(ndim,*),bmat_new(*)


      do i = 1, nnode*ndim
          bmat(i) = zero
          bmat_new(i) = zero
      end do

      do inod = 1, nnode
        do ider = 1, ndim
          do idim = 1, ndim
             irow = idim + (inod - 1)*ndim
             bmat(irow) = bmat(irow) + 
     1       xjaci(idim,ider)*dshape(ider,inod)      

             bmat_new(irow) = bmat_new(irow) + 
     1       xjaci_new(idim,ider)*dshape(ider,inod)      

          end do
        end do
      end do 

      return
      end

c*****************************************************************
      subroutine straininc(ntens,ndof,ndim,nnode,
     1 mlvarx,bmat,du,dstran,u,xx1,xx1Old)
c
c     Notation:  
c       dstran(i)  incremental strain component 
c       note:      i = 1   xx direction 
c                    = 2   yy direction 
c                    = 3   zz direction
c                    = 4   xy direction
c    u() - displacement
c   du() - increment of displacement in the last inc.
c   
c
      include 'aba_param.inc'
      parameter(zero=0.d0, one=1.d0)

      dimension dstran(*), bmat(ndim,*), 
     1     du(mlvarx, *), xdu(3), xx1(3,*), 
     2      u(mlvarx, *), utmp(3),
     3      utmpOld(3),xx1Old(3,*),eps(3,3),dInvFold(3,3)
      do i = 1, ntens
        dstran(i) = zero
      end do
      !
      ! set xx1 to Identity matrix
      do k1=1,3
         do k2=1,3
            xx1(k1,k2) = zero
            xx1Old(k1,k2) = zero
         end do
         xx1(k1,k1) = one
         xx1Old(k1,k1) = one
      end do
c
c************************************
c    Compute incremental strains
c************************************
c
      do nodi = 1, nnode
           
         incr_row = (nodi - 1)*ndof

         do i = 1, ndof
           xdu(i)= du(i + incr_row,1)
           utmp(i) = u(i + incr_row,1)
           utmpOld(i) = utmp(i)-xdu(i)
         end do

         dNidx = bmat(1,nodi)
         dNidy = bmat(2,nodi)
         if (ndim.eq.3) then
           dNidz = bmat(3,nodi)

           dstran(1) = dstran(1) + dNidx*xdu(1)
           dstran(2) = dstran(2) + dNidy*xdu(2)
           dstran(3) = dstran(3) + dNidz*xdu(3)
           ! (1-2)
           dstran(4) = dstran(4) + 
     1          dNidy*xdu(1) + 
     2          dNidx*xdu(2)  
           ! (1-3)
           dstran(5) = dstran(5) + 
     1          dNidz*xdu(1) + 
     2          dNidx*xdu(3)  
           ! (2-3)
           dstran(6) = dstran(6) + 
     1          dNidy*xdu(3) + 
     2          dNidz*xdu(2)  

c          deformation gradient

           xx1(1,1) = xx1(1,1) + dNidx*utmp(1)
           xx1(1,2) = xx1(1,2) + dNidy*utmp(1)
           xx1(1,3) = xx1(1,3) + dNidz*utmp(1)
           !
           xx1(2,1) = xx1(2,1) + dNidx*utmp(2)
           xx1(2,2) = xx1(2,2) + dNidy*utmp(2)
           xx1(2,3) = xx1(2,3) + dNidz*utmp(2)
           !
           xx1(3,1) = xx1(3,1) + dNidx*utmp(3)
           xx1(3,2) = xx1(3,2) + dNidy*utmp(3)
           xx1(3,3) = xx1(3,3) + dNidz*utmp(3)


         else if (ndim.eq.2) then
           dstran(1) = dstran(1) + dNidx*xdu(1)
           dstran(2) = dstran(2) + dNidy*xdu(2)
           dstran(4) = dstran(4) + 
     1          dNidy*xdu(1) + 
     2          dNidx*xdu(2)  

c        deformation gradient

           xx1(1,1) = xx1(1,1) + dNidx*utmp(1)
           xx1(1,2) = xx1(1,2) + dNidy*utmp(1)
           xx1(2,1) = xx1(2,1) + dNidx*utmp(2)
           xx1(2,2) = xx1(2,2) + dNidy*utmp(2)
c
           xx1Old(1,1) = xx1Old(1,1) + dNidx*utmpOld(1)
           xx1Old(1,2) = xx1Old(1,2) + dNidy*utmpOld(1)
           xx1Old(2,1) = xx1Old(2,1) + dNidx*utmpOld(2)
           xx1Old(2,2) = xx1Old(2,2) + dNidy*utmpOld(2)
         
       end if

      end do
c
      return
      end

c*****************************************************************
      subroutine statevar(npt,nsvint,statev,statev_ip,icopy)
c
c     Transfer data to/from element-level state variable array from/to
c     material-point level state variable array.
c
      include 'aba_param.inc'

      dimension statev(*),statev_ip(*)

      isvinc= (npt-1)*nsvint     ! integration point increment

      if (icopy .eq. 1) then
c
c     Prepare arrays for entry into umat
c
        do i = 1, nsvint
          statev_ip(i)=statev(i+isvinc)
        end do
c
      else
c
c     Update element state variables upon return from umat
c
        do i = 1, nsvint
          statev(i+isvinc)=statev_ip(i)
        end do
      end if

      return
      end

c*****************************************************************
      subroutine stiffmatrix(ntens,nnode,ndim,ndof,
     1 weight,djac,ddsdde,stress,bmat,stiff,force)
c
c     Stiffness matrix and internal force contributions at 
c     material integration point
c
      include 'aba_param.inc'

      parameter(zero=0.d0,maxDof=3)

      dimension stiff(ndof*nnode,*), stiff_p(maxDof,maxDof)
      dimension force(*), force_p(maxDof)
      dimension stress(*),bmat(ndim,*),ddsdde(ntens,*)

      dNjdx = zero
      dNjdy = zero
      dNjdz = zero
      do i = 1, ndof*nnode
          force(i) = zero
        do j = 1, ndof*nnode
          stiff(j,i) = zero
        end do
      end do

      dvol= weight*djac
      do nodj = 1, nnode

          incr_col = (nodj - 1)*ndof

          dNjdx = bmat(1,nodj)
          dNjdy = bmat(2,nodj)
          if (ndim.eq.2) then
            force_p(1) = dNjdx*stress(1) + dNjdy*stress(4)
            force_p(2) = dNjdy*stress(2) + dNjdx*stress(4)
          else if (ndim.eq.3) then
            dNjdz = bmat(3,nodj)
            force_p(1) = dNjdx*stress(1) +
     &           dNjdy*stress(4) +
     &           dNjdz*stress(5)
            
            force_p(2) = dNjdy*stress(2) +
     &           dNjdx*stress(4) +
     &           dNjdz*stress(6)
            
            force_p(3) = dNjdy*stress(6) +
     &           dNjdx*stress(5) +
     &           dNjdz*stress(3)
          end if            
          do jdof = 1, ndof

            jcol = jdof + incr_col

            force(jcol) = force(jcol) +
     &           force_p(jdof)*dvol
            
          end do

          do nodi = 1, nnode

            incr_row = (nodi -1)*ndof

            dNidx = bmat(1,nodi)
            dNidy = bmat(2,nodi)
            if (ndim.eq.2) then
              stiff_p(1,1) = dNidx*ddsdde(1,1)*dNjdx
     &             + dNidy*ddsdde(4,4)*dNjdy
     &             + dNidx*ddsdde(1,4)*dNjdy
     &             + dNidy*ddsdde(4,1)*dNjdx

              stiff_p(1,2) = dNidx*ddsdde(1,2)*dNjdy
     &             + dNidy*ddsdde(4,4)*dNjdx
     &             + dNidx*ddsdde(1,4)*dNjdx
     &             + dNidy*ddsdde(4,2)*dNjdy
              
              stiff_p(2,1) = dNidy*ddsdde(2,1)*dNjdx
     &             + dNidx*ddsdde(4,4)*dNjdy
     &             + dNidy*ddsdde(2,4)*dNjdy
     &             + dNidx*ddsdde(4,1)*dNjdx

              stiff_p(2,2) = dNidy*ddsdde(2,2)*dNjdy
     &             + dNidx*ddsdde(4,4)*dNjdx
     &             + dNidy*ddsdde(2,4)*dNjdx
     &             + dNidx*ddsdde(4,2)*dNjdy
              
            else if (ndim.eq.3) then
              dNidz = bmat(3,nodi)
              stiff_p(1,1) = dNjdx*ddsdde(1,1)*dNidx
     &             + dNjdx*ddsdde(1,4)*dNidy
     &             + dNjdx*ddsdde(1,5)*dNidz
              
     &             + dNjdy*ddsdde(4,1)*dNidx 
     &             + dNjdy*ddsdde(4,4)*dNidy
     &             + dNjdy*ddsdde(4,5)*dNidz
              
     &             + dNjdz*ddsdde(5,1)*dNidx
     &             + dNjdz*ddsdde(5,4)*dNidy
     &             + dNjdz*ddsdde(5,5)*dNidz
c     
              stiff_p(1,2) = dNjdx*ddsdde(1,2)*dNidy
     &             + dNjdx*ddsdde(1,4)*dNidx
     &             + dNjdx*ddsdde(1,6)*dNidz
              
     &             + dNjdy*ddsdde(4,2)*dNidy 
     &             + dNjdy*ddsdde(4,4)*dNidx
     &             + dNjdy*ddsdde(4,6)*dNidz
              
     &             + dNjdz*ddsdde(5,2)*dNidy
     &             + dNjdz*ddsdde(5,4)*dNidx
     &             + dNjdz*ddsdde(5,6)*dNidz
c     
              stiff_p(1,3) = dNjdx*ddsdde(1,3)*dNidz
     &             + dNjdx*ddsdde(1,5)*dNidx
     &             + dNjdx*ddsdde(1,6)*dNidy

     &             + dNjdy*ddsdde(4,3)*dNidz 
     &             + dNjdy*ddsdde(4,5)*dNidx
     &             + dNjdy*ddsdde(4,6)*dNidy
              
     &             + dNjdz*ddsdde(5,3)*dNidz
     &             + dNjdz*ddsdde(5,5)*dNidx
     &             + dNjdz*ddsdde(5,6)*dNidy
c
              stiff_p(2,1) = dNjdy*ddsdde(2,1)*dNidx
     &             + dNjdy*ddsdde(2,4)*dNidy
     &             + dNjdy*ddsdde(2,5)*dNidz

     &             + dNjdx*ddsdde(4,1)*dNidx 
     &             + dNjdx*ddsdde(4,4)*dNidy
     &             + dNjdx*ddsdde(4,5)*dNidz
              
     &             + dNjdz*ddsdde(6,1)*dNidx
     &             + dNjdz*ddsdde(6,4)*dNidy
     &             + dNjdz*ddsdde(6,5)*dNidz
c
              stiff_p(2,2) = dNjdy*ddsdde(2,2)*dNidy
     &             + dNjdy*ddsdde(2,4)*dNidx
     &             + dNjdy*ddsdde(2,6)*dNidz

     &             + dNjdx*ddsdde(4,2)*dNidy 
     &             + dNjdx*ddsdde(4,4)*dNidx
     &             + dNjdx*ddsdde(4,6)*dNidz

     &             + dNjdz*ddsdde(6,2)*dNidy
     &             + dNjdz*ddsdde(6,4)*dNidx
     &             + dNjdz*ddsdde(6,6)*dNidz
c
              stiff_p(2,3) = dNjdy*ddsdde(2,3)*dNidz
     &             + dNjdy*ddsdde(2,5)*dNidx
     &             + dNjdy*ddsdde(2,6)*dNidy
              
     &             + dNjdx*ddsdde(4,3)*dNidz 
     &             + dNjdx*ddsdde(4,5)*dNidx
     &             + dNjdx*ddsdde(4,6)*dNidy
              
     &             + dNjdz*ddsdde(6,3)*dNidz
     &             + dNjdz*ddsdde(6,5)*dNidx
     &             + dNjdz*ddsdde(6,6)*dNidy
c
              stiff_p(3,1) = dNjdz*ddsdde(3,1)*dNidx
     &             + dNjdz*ddsdde(3,4)*dNidy
     &             + dNjdz*ddsdde(3,5)*dNidz
              
     &             + dNjdx*ddsdde(5,1)*dNidx 
     &             + dNjdx*ddsdde(5,4)*dNidy
     &             + dNjdx*ddsdde(5,5)*dNidz
              
     &             + dNjdy*ddsdde(6,1)*dNidx
     &             + dNjdy*ddsdde(6,4)*dNidy
     &             + dNjdy*ddsdde(6,5)*dNidz
c     
              stiff_p(3,2) = dNjdz*ddsdde(3,2)*dNidy
     &             + dNjdz*ddsdde(3,4)*dNidx
     &             + dNjdz*ddsdde(3,6)*dNidz
              
     &             + dNjdx*ddsdde(5,2)*dNidy 
     &             + dNjdx*ddsdde(5,4)*dNidx
     &             + dNjdx*ddsdde(5,6)*dNidz

     &             + dNjdy*ddsdde(6,2)*dNidy
     &             + dNjdy*ddsdde(6,4)*dNidx
     &             + dNjdy*ddsdde(6,6)*dNidz
c
              stiff_p(3,3) = dNjdz*ddsdde(3,3)*dNidz
     &             + dNjdz*ddsdde(3,5)*dNidx
     &             + dNjdz*ddsdde(3,6)*dNidy

     &             + dNjdx*ddsdde(5,3)*dNidz 
     &             + dNjdx*ddsdde(5,5)*dNidx
     &             + dNjdx*ddsdde(5,6)*dNidy

     &             + dNjdy*ddsdde(6,3)*dNidz
     &             + dNjdy*ddsdde(6,5)*dNidx
     &             + dNjdy*ddsdde(6,6)*dNidy
            end if

            do jdof = 1, ndof
              icol = jdof + incr_col
              do idof = 1, ndof
                irow = idof + incr_row
                stiff(irow,icol) = stiff(irow,icol) +
     &          stiff_p(idof,jdof)*dvol
              end do
            end do
          end do
      end do
c
      return
      end

