C
C  User element accessing Abaqus materials 
C  Heat Transfer -- conduction
C
c*****************************************************************
      subroutine uelmat(rhs, amatrx, svars, energy, ndofel, nrhs, 
     1  nsvars, props, nprops, coords, mcrd, nnode, u, du, v, a, jtype, 
     2  time, 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(*)
c
c     local arrays
c
      parameter (zero=0.d0, one=1.d0)

      parameter (ndim=2, ndof=1, ninpt=4, nnodemax=4)
c
c      ndim  ... number of spatial dimensions
c      ndof  ... number of degrees of freedom per node
c      ninpt ... number of integration points
c
      dimension  stiff(ndof*nnodemax,ndof*nnodemax),
     1 force(ndof*nnodemax), shape(nnodemax), dshape(ndim,nnodemax),
     2 xjaci(ndim,ndim), bmat(nnodemax*ndim), wght(ninpt)
  
      dimension coords_ip(3),dfgrd0(3,3),dfgrd1(3,3),
     1 drot(3,3)
      dimension coords_new(mcrd,nnodemax)
c
      dimension predef_loc(npredf),dpredef_loc(npredf),xx1(3,3),
     1     xx1Old(3,3)
      dimension xjaci_new(ndim,ndim),bmat_new(nnodemax*ndim)
      dimension dtemdx(ndim),rhoUdotdg(3),flux(ndim),dfdt(ndim),
     1     dfdg(ndim,ndim)
c
      data wght /one, one, one, one/

c
c********************************************************************
c
c     U1 = first-order, plane strain, full integration
c
c********************************************************************
      if (lflags(3).eq.4) goto 999
c
c     Preliminaries
c
      pnewdtLocal = pnewdt
      if(jtype .ne. 1) then
        write(7,*)'Incorrect element type'
        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      initialization
c
        rho = zero
        rhoUdot   = zero
        rhoUdotdt = zero
        rhoUdotdg = zero
        do i=1, 3
          rhoUdotdg(i) = zero
        end do
        do i=1, ndim
          flux(i)   = zero
          dfdt(i)   = zero
        end do
        do i=1, ndim
          do j=1, ndim
            dfdg(i,j) = zero
          end do
        end do
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
            coords_ip(k2)=coords_ip(k2)+shape(k1)*coords(k2,k1)
          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)
c
        if (pnewdt .lt. pnewdtLocal) pnewdtLocal = pnewdt
c
        call bmatrix(xjaci,dshape,nnode,ndim,bmat,xjaci_new,
     1       bmat_new)
c
c       compute temp. and temp. gradient
c 
        temp = zero
        dtemp = zero
        call settemp(ndofel,ndof,ndim,nnode,mlvarx,bmat,du,
     *       dstran,u,xx1,xx1Old,temp,dtemp,dtemdx,shape)
c
c       get Abaqus material
c
        rpl = zero
        drpldt = zero
        celent = one
        call material_lib_ht(materiallib,rhoUdot,rhoUdotdt,rhoUdotdg,
     *       flux,dfdt,dfdg,rpl,drpldt,kintk,djac,predef_loc,
     *       dpredef_loc,npredf,temp,dtemp,dtemdx,celent,coords_ip)
c
c        
c       form stiffness matrix and internal force vector
c
        call rhsjacobian(nnode,ndim,ndof,
     1       wght(kintk),djac,rhoUdot,rhoUdotdt,rhoUdotdg,flux,
     2       dfdt,dfdg,shape,bmat,stiff,force,dtime,lflags)

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   Compute shape fuctions
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,
     1     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*****************************************************************
c   Get local predefined fileds
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*****************************************************************
c   Compute jacobian matrix
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
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(idim,jdim) = xjac(idim,jdim) + 
     1             dshape(jdim,inod)*coords(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)
        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
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. 
     1               Jacobian'
          pnewdt = fourth
        endif        
      end if

      return
      end
c*****************************************************************
c
c   Compute the B matrix
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'
c
      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)      

          end do
        end do
      end do 

      return
      end
c*****************************************************************
c
c   Set temperatures
c
      subroutine settemp(ndofel,ndof,ndim,nnode,
     1 mlvarx,bmat,du,dstran,u,xx1,xx1Old,temp,dtemp,dtemdx,dN)
c
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(ndofel), utmp(3),
     3      utmpOld(3),xx1Old(3,*),eps(3,3),dInvFold(3,3)
      dimension dtemdx(*),dN(*)
C      
c
c****************************************************************
c    Compute temp, dtemp, and temp gradient at the material point
c****************************************************************
c

      temp  = zero
      dtemp = zero
      do iNode=1, nnode
        temp  = temp + dN(iNode)*u(iNode)
        dtemp = dtemp + dN(iNode)*du(iNode,1)
      end do
      do iDof = 1, ndim
        dtemdx(iDof) = zero
        do iNode=1, nnode
          dtemdx(iDof) = dtemdx(iDof) + bmat(idof,iNode)*u(iNode)
        end do
      end do
c
      return
      end
c*****************************************************************
c
c   Compute element jacobian and nodal forces
c
      subroutine rhsjacobian(nnode,ndim,ndof,
     1     weight,djac,rhoUdot,rhoUdotdt,rhoUdotdg,flux,dfdt,
     2     dfdg,dN,bmat,stiff,force,dtime,lflags)
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,*)
      dimension force(*)
      dimension flux(*),dfdt(*),dfdg(ndim,*),rhoUdotdg(*)
      dimension dN(*),bmat(ndim,*),lflags(*)

      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
        if (lflags(1).eq.32.or.lflags(1).eq.33) 
     &       force(nodj) = dN(nodj)*rhoUdot*dvol
ccc        force(nodj) = dN(nodj)*(rhoUdot-rpl)*dvol
        do jDof=1, ndim
          force(nodj) = force(nodj)+bmat(jDof,nodj)*flux(jDof)*dvol
          do nodi=1, nnode
ccc            stiff(nodj,nodi) = stiff(nodj,nodi) + 
ccc     *           bmat(jDof,nodj)*dN(nodi)*dfdt(jDof)*dvol
            do iDof=1, ndim
              stiff(nodj,nodi) = stiff(nodj,nodi) + 
     *             bmat(jDof,nodj)*bmat(iDof,nodi)*dfdg(jDof,iDof)*dvol
            end do
          end do
        end do
      end do
      do nodj=1, nnode
        do nodi=1, nnode
          do iDof=1, ndim
ccc            stiff(nodj,nodi) = stiff(nodj,nodi) + 
ccc     *           dN(nodj)*bmat(iDof,nodi)*rhoUdotdg(iDof)*dvol
          end do
        end do
      end do
c
      if (lflags(1).eq.32.or.lflags(1).eq.33) then
        do nodj=1, nnode
          do nodi=1, nnode
            stiff(nodj,nodi) = stiff(nodj,nodi) + 
     *           rhoUdotdt*dN(nodj)*dN(nodi)*dvol/dtime
          end do
        end do
      end if
c
      return
      end