      subroutine elmt04(d,u,x,ix,t,s,p,ndf,ndm,nst,isw)

c     Purpose: Any dimensional truss element routine

c     Inputs:
c        d(*)      - Material set parameters
c        ul(ndf,*) - Solution variables for element
c        xl(ndf,*) - Nodal coordinates for element
c        ix(*)     - Nodal connections for element
c        tl(*)     - Nodal temperatures for element
c        ndf       - Number dof/node
c        ndm       - Spatial dimension of FEM mesh
c        nst       - Size of element matrix
c        isw       - Switch parameter to control action taken by routine

c     Outputs:
c        s(nst,*)  - Element matrix
c        p(*)      - Element vector

      implicit real*8 (a-h,o-z)
      real*8 x(ndm,*),t(*),dm
      real*8 d(*),u(ndf,*),s(nst,*),p(*),db(3),dx(3),xx(3)
      character*4 yyy*80
      common /cdata/ numnp,numel,nummat,nen,neq
      common /eldata/ dm,n,ma,mct,iel,nel
      common /hdata/  nh1,nh2
      common /iofile/ ior,iow
      go to (1,2,3,3,3,3,2,2),isw
c.... input material properties
1     if(ior.lt.0) write(*,3000) 
      call pintio(yyy,10)
      read(yyy,1000,err=110) (d(i),i=1,6)
      d( 7) = d(1)*d(2)
      d( 8) = d(4)*d(2)
      d( 9) = d(5)*d(2)
      d(10) = d(6)*d(2)
      d(11) = d(3)*d(2)
      call pconsd(xx,3,0.0d0)
      if(d(4).gt.0.0d0) then
        write(iow,2000) (d(i),i=1,6)
        if(ior.lt.0) write(*,2000) (d(i),i=1,6)
      else
        write(iow,2001) (d(i),i=1,3)
        if(ior.lt.0) write(*,2001) (d(i),i=1,3)
      endif
      nh1 = 3
      return
110   call perror('PCELM4',yyy)
      go to 1
2     return
c.... compute element arrays
3     xl  = 0.0
      eps = 0.0
      do 30 i = 1,ndm
        dx(i) = x(i,2) - x(i,1)
        xl = xl + dx(i)**2
        eps = eps + dx(i)*(u(i,2)-u(i,1))
        xx(i) = (x(i,2) + x(i,1))/2.
30    continue
      eps = eps/xl
      if(mod(isw,3).eq.0) then
        call modl04(d,eps, sig,ad)
c.... form a residual
        sig = sig/sqrt(xl)
        do 31 i = 1,ndf
          p(i) = dx(i)*sig
          p(i+ndf) = -p(i)
31      continue
c.... compute tangent stiffness
        if(isw.eq.3) then
          xl = xl*sqrt(xl)
          do 32 i = 1,ndm
            db(i) = ad*dx(i)
            dx(i) = dx(i)/xl
32        continue
          i1 = 0
          do 36 ii = 1,2
            j1 = i1
            do 35 jj = ii,2
              do 33 i = 1,ndm
              do 33 j = 1,ndm
                s(i+i1,j+j1) = db(i)*dx(j)
33            continue
              j1 = j1 + ndf
              do 34 j = 1,ndm
                dx(j) = -dx(j)
34            continue
35          continue
            i1 = i1 + ndf
36        continue
          do 37 i = 1,ndm
          do 37 j = 1,ndm
            s(i+ndf,j) = s(j,i+ndf)
37        continue
        endif
c.... output stress and strain in element
      elseif(isw.eq.4) then
        call modl04(d,eps, sig,ad)
        mct = mct - 1
        if(mct.le.0) then
          call prthed(iow)
          write(iow,2002) 
          if(ior.lt.0) write(*,2002)
          mct = 50
        endif
        write(iow,2003) n,ma,xx,sig,eps
        if(ior.lt.0) write(*,2003) n,ma,xx,sig,eps
c.... compute element lumped mass matrix
      elseif(isw.eq.5) then
        sm = d(11)*sqrt(xl)/2.0
        do 52 i = 1,ndm
          p(i    ) = sm
          p(i+ndf) = sm
52      continue
      endif
c.... formats
1000  format(8f10.0)
2000  format(5x,'T r u s s    E l e m e n t'//
     1   10x,'Modulus  =',e12.5/10x,'Area     =',e12.5/
     2   10x,'Density  =',e12.5/10x,'Yield    =',e12.5/
     3   10x,'Iso. Hard=',e12.5/10x,'Kin. Hard=',e12.5/)
2001  format(5x,'T r u s s    E l e m e n t'//
     1   10x,'Modulus  =',e12.5/10x,'Area     =',e12.5/
     2   10x,'Density  =',e12.5/)
2002  format(5x,'T r u s s    E l e m e n t'//' elem mate',
     1  4x,'1-coord',4x,'2-coord',4x,'3-coord',5x,'force',7x,'strain')
2003  format(2i5,3f11.4,2e13.5)
3000  format(' Input  El/Plas: E, A, rho, Y, H-iso, H-Kin'/
     1       '        Elastic: E, A, rho'/3x,'>',$)
      end
c
      subroutine modl04(d,eps, sig,ad)
      real*8 d(10), eps, sig, ad, yld, sum, gam, h
      common /hdata/  nh1,nh2
      common h(1)
c.... trial stress
      ad  = d(7)
      sig = d(7)*(eps - h(nh1))
      if(d(8).gt.0.0d0) then
c.... compute plastic corrections
        yld = d(8) + d(9)*h(nh1+2)
        phi = abs(sig - h(nh1+1))
        if(phi.gt.yld) then
          sum      = d(7) + d(9) + d(10)
          gam      = (phi-yld-1.d-08*d(8))/sum
          sig      = sig - d(7)*gam
          ad       = ad -  d(7)**2/sum
c.... update the history terms
          sum      = gam*(sig - h(nh1+1))/phi
          h(nh1  ) = h(nh1  ) + gam*sum
          h(nh1+1) = h(nh1+1) + gam*sum*d(10)
          h(nh1+2) = h(nh1+2) + gam
        endif
      endif
      end
