c$Id:$
      subroutine trussnd(d,ul,xl,ix,tl,s,p,ndf,ndm,nst,isw)

c      * * F E A P * * A Finite Element Analysis Program

c....  Copyright (c) 1984-2005: Robert L. Taylor
c                               All rights reserved

c-----[--.----+----.----+----.-----------------------------------------]
c     1, 2 and 3 dimensional truss element routine

c     Outputs: (isw = 4)

c         xx  - Coordinates at mid-length of truss bar
c         sig - Force on truss bar
c         eps - Strain on truss bar

c-----[--.----+----.----+----.-----------------------------------------]
      implicit  none

      include  'augdat.h'
      include  'bdata.h'
      include  'cdata.h'
      include  'eldata.h'
      include  'elplot.h'
      include  'eltran.h'
      include  'hdata.h'
      include  'iofile.h'
      include  'pmod2d.h'
      include  'prld1.h'
      include  'prstrs.h'
      include  'ptdat6.h'
      include  'comblk.h'

      logical   nonlin
      integer   i, j, ii, jj, i1, j1, nhv,nhi, tdof
      integer   ndf, ndm, nst, isw, istrt
      real*8    xlen, xlen0, xlen2, sig(2), eps(2), ta, body, dd(2)
      real*8    lm,cmd,cmo,hmd,hmo,geo,forc, alphar, one3

      real*8    db(3,2),bl(3,2),br(3,2),xx(3),th(2)

      integer   ix(*)
      real*8    d(*),ul(ndf,nen,*),xl(ndm,nen),tl(*),s(nst,nst),p(ndf,*)

      save

      data      one3 / 0.3333333333333333d0 /

c     Set Analysis Type

      if(nint(d(39)).eq.1) then
        dtype  =  1
        nonlin = .true.
      else
        dtype  =  nint(d(18))
        nonlin =  dtype.lt.0
      endif
      hflag  = nint(d(30)).eq.1
      nhv    = d(15)

c     Set nodal temperatures: Can be specified or computed

      if(isw.gt.1) then
        tdof = d(19)
        if(tdof.le.0) then
          do i = 1,nel ! {
            th(i) = tl(i)
          end do ! i     }
        else
          do i = 1,nel ! {
            th(i) = ul(tdof,i,1)
          end do ! i     }
        endif
      endif

c     INPUT MATERIAL PROPERTIES

      if(isw.eq.1) then

        write(iow,2000)
        if(ior.lt.0) write(*,2000)

c       Input material parameters

        call inmate(d,tdof,ndm*2,2)
        if(nint(d(160)).eq.2) then
          nh1    = nh1 + 1 ! Augmented storage
          d(166) = 1
        else
          d(166) = 0
        endif

c       Set tdof to zero if 1, 2, or larger than ndf

        if(tdof.gt.ndf) then
          write(iow,3003)
          if(ior.lt.0) write(*,3003)
          tdof = 0
        elseif(tdof.gt.0 .and. tdof.le.ndm) then
          write(iow,3004)
          if(ior.lt.0) write(*,3004)
          tdof = 0
        endif

c       Deactivate dof in element for dof > ndm

        do i = ndm+1,ndf
          ix(i) = 0
        end do

c       If temperature dof is specified activate dof

        if(tdof.gt.0) then
          ix(tdof) = 1
        endif

c       Set plot sequence

        call pltln2(iel)

c     CHECK FOR ZERO LENGTH ELEMENTS

      elseif(isw.eq.2) then

        if(ix(1).eq.0 .or. ix(2).eq.0) then
          write(iow,4000) n,ix(1),ix(2)
          if(ior.lt.0) then
            write(*,4000) n,ix(1),ix(2)
          endif
        else
          xlen = 0.0d0
          eps(1)  = 0.0d0
          do i = 1,ndm
            eps(1)  = max(eps(1) , abs(xl(i,1)), abs(xl(i,2)))
            xlen = max(xlen, abs(xl(i,2) - xl(i,1)))
          end do
          if(xlen.eq.1.0d-10*eps(1)) then
            write(iow,4001) n
            if(ior.lt.0) then
              write(*,4001) n
            endif
          endif
        endif

c     COMPUTE ELEMENT STIFFNESS AND RESIDUAL ARRAYS

      elseif(isw.eq.3 .or. isw.eq.6) then

c       Compute element stress, strain, and tangent modulus

        call strn1d(d,xl,ul,th,ndm,ndf,nen,xlen,xlen0,xx,eps,ta,
     &              nonlin,isw)
        nhi    = nint(d(166))
        if(nint(d(160)).eq.2) then ! Set augmented value
          sig(1) = hr(nh2)
        else
          sig(1) = 0.0d0
        endif
        istrt = nint(d(84))
        call modl1d(d,ta,eps,hr(nhi+nh1),hr(nhi+nh2),nhv,
     &              1,istrt,dd,sig,isw)

c       Save stress and strain for tplot

        tt(1) = sig(1)
        tt(2) = eps(1)

c       Multiply tangent modulus by length and area

        if(dtype.lt.0) then
          dd(1) = dd(1) - sig(1)
        endif

        dd(1) = (dd(1)*ctan(1) + dd(2)*d(78)*ctan(2))*d(32)*xlen

c       Compute strain-displacement matrix

        xlen2  = 1.d0/(xlen*xlen)
        alphar = xlen2/ctan(1) - xlen2
        do i = 1,ndm

          bl(i,1) = (xl(i,1) - xl(i,2))*xlen2

c         Set non-linear terms in strain-displacement matrix

          if(nonlin) then
            bl(i,1) = bl(i,1) + (ul(i,1,1) - ul(i,2,1))*xlen2
          endif
          br(i,1) = bl(i,1)

c         Set remaining terms and multiply by elastic modulus

          bl(i,2) = -bl(i,1)
          br(i,2) = -br(i,1)

          db(i,1) =  bl(i,1)*dd(1)
          db(i,2) = -db(i,1)

        end do

c       Compute mass terms

        cmd = d(4)*d(32)*xlen0*one3
        cmo = cmd*0.5d0
        lm  = cmd*1.5d0
        hmd = d(7)*cmd + (1.d0 - d(7))*lm
        hmo = d(7)*cmo

c       Form internal, body  and inertia force vector
c       Include stiffness proportional Rayleigh damping effect in force
c       Include mass proportional Rayleigh damping in residual

        forc = (sig(1) + d(78)*sig(2))*xlen*d(32)

        do i = 1,ndm

c         Set body loading factors

          if(int(d(73+i)).gt.0) then
            body = 0.5*xlen0*(d(10+i) + prldv(int(d(73+i)))*d(70+i))
          else
            body = 0.5*xlen0*d(10+i)*dm
          endif

          p(i,1) = body   - bl(i,1)*forc
     &           - hmd*(ul(i,1,5) + d(77)*ul(i,1,4))
     &           - hmo*(ul(i,2,5) + d(77)*ul(i,2,4))

          p(i,2) = body   - bl(i,2)*forc
     &           - hmo*(ul(i,1,5) + d(77)*ul(i,1,4))
     &           - hmd*(ul(i,2,5) + d(77)*ul(i,2,4))

        end do

c       Compute stiffness terms

        if(isw.eq.3) then

          i1 = 0
          do ii = 1,2
            j1 = 0
            do jj = 1,2
              do i = 1,ndm
                do j = 1,ndm
                  s(i+i1,j+j1) = db(i,ii)*br(j,jj)
                end do
              end do
              j1 = j1 + ndf
            end do
            i1 = i1 + ndf
          end do

c         Correct for geometric and inertial tangent effects

          if(nonlin .and. gflag) then
            geo = sig(1)*d(32)/xlen*ctan(1)
          else
            geo = 0.0d0
          endif

c         Set diagonal and off-diagonal terms
c         Include mass proportional rayleigh damping

          hmd = hmd*(ctan(3) + d(77)*ctan(2)) + geo
          hmo = hmo*(ctan(3) + d(77)*ctan(2)) - geo

          do i = 1,ndm
            j      = i + ndf
            s(i,i) = s(i,i) + hmd
            s(i,j) = s(i,j) + hmo
            s(j,i) = s(j,i) + hmo
            s(j,j) = s(j,j) + hmd
          end do

        endif

c     OUTPUT STRESS AND STRAIN IN ELEMENT

      elseif(isw.eq.4 .or. isw.eq.8) then

c       Form strain and stress

        call strn1d(d,xl,ul,th,ndm,ndf,nen,xlen,xlen0,xx,eps,ta,
     &              nonlin,isw)
        nhi  = nint(d(166))
        if(nint(d(160)).eq.2) then ! Set augmented value
          sig(1) = hr(nh2)
        else
          sig(1) = 0.0d0
        endif
        istrt = nint(d(84))
        call modl1d(d,ta,eps,hr(nhi+nh1),hr(nhi+nh2),nhv,
     &              1,istrt,dd,sig,isw)

c       Form truss force: multiply stress by area

        sig(1) = d(32)*sig(1)

c       Output element force/strains

        if(isw.eq.4) then
          mct = mct - 1
          if(mct.le.0) then
            write(iow,2001) o,head
            if(ior.lt.0) write(*,2001) o,head
            mct = 50
          endif
          write(iow,2002) n,ma,xx,sig(1),eps(1)
          if(ior.lt.0) then
            write(*,2002) n,ma,xx,sig(1),eps(1)
          endif
        else
          call trcnnd(ix,sig(1),hr(nph),hr(nph+numnp))
        endif

c     COMPUTE ELEMENT MASS MATRICES

      elseif(isw.eq.5) then

        xlen0 = 0.0d0
        do i = 1,ndm
          xlen0 = xlen0 + (xl(i,2)-xl(i,1))**2
        end do

c       Compute mass matrix terms

        cmd = d(4)*d(32)*sqrt(xlen0)*one3
        lm  = 1.5d0*cmd
        cmo = 0.5d0*cmd
        hmd = d(7)*cmd + (1.d0 - d(7))*lm
        hmo = d(7)*cmo

        do i = 1,ndm

c         Higher order mass

          j      = i + ndf
          s(i,i) = hmd
          s(j,j) = hmd
          s(j,i) = hmo
          s(i,j) = hmo

c         Lumped mass

          p(i,1) = lm
          p(i,2) = p(i,1)

        end do

c     Augmented update

      elseif(isw.eq.10) then

c       Form strain and stress

        if(nint(d(160)).eq.2) then ! Set augmented value
          call strn1d(d,xl,ul,th,ndm,ndf,nen,xlen,xlen0,xx,eps,ta,
     &                nonlin,isw)
          nhi    = nint(d(166))
          sig(1) = hr(nh2)
          istrt = nint(d(84))
          call modl1d(d,ta,eps,hr(nhi+nh1),hr(nhi+nh2),nhv,
     &                1,istrt,dd,sig,isw)

          hr(nh2) = hr(nh2) + augf*(sig(1) - hr(nh2))
        endif

c     COMPUTE ELEMENT ENERGY

      elseif(isw.eq.13) then

c       Compute element stress, strain, and tangent modulus

        call strn1d(d,xl,ul,th,ndm,ndf,nen,xlen,xlen0,xx,eps,ta,
     &              nonlin,isw)
        istrt = nint(d(84))
        call modl1d(d,ta,eps,hr(nh1),hr(nh2),nhv,
     &              1,istrt,dd,sig,isw)

c       Stored energy

        epl(8) = epl(8) + 0.5d0*eps(1)*sig(1)*d(32)*xlen

c       Compute mass terms

        cmd = d(4)*d(32)*xlen0*one3
        cmo = cmd*0.5d0
        lm  = cmd*1.5d0
        hmd =(d(7)*cmd + (1.d0 - d(7))*lm)*0.5d0
        hmo = d(7)*cmo

c       Kinetic energy

        do i = 1,ndm
          epl(7) = epl(7) + hmd*(ul(i,1,4)**2 + ul(i,2,4)**2)
     &                    + hmo* ul(i,1,4)*ul(i,2,4)
        end do

      endif

c     FORMATS

2000  format(9x,'T r u s s    E l e m e n t'/1x)
2001  format(a1,20a4//9x,'Truss Element'//' Elmt Matl    ',
     & '1-coord    2-coord    3-coord',5x,'Force',9x,'Strain')
2002  format(2i5,1p3e11.3,1p2e14.5)

3003  format(' *WARNING* Thermal d.o.f. > active d.o.f.s : Set to 0')
3004  format(' *WARNING* Thermal d.o.f. can not be <= ndm: Set to 0')

4000  format(' *ERROR* Element',i7,' has nodes',2i8)
4001  format(' *ERROR* Element',i7,' has zero length')

      end

      subroutine strn1d(d,xl,ul,tl,ndm,ndf,nen,xlen,xlen0,xx,eps,ta,
     &                  nonlin,isw)

c     Compute constitutive equation

      implicit  none

      include  'iofile.h'
      include  'eltran.h'
      include  'pmod2d.h'

      logical   nonlin
      integer   i,ndm,ndf,nen, isw
      real*8    xlen,xlen0,xlenn,xlen1,xlena,dx,du,dd, ta
      real*8    alpha, alphar
      real*8    xx(3),d(*),xl(ndm,*),ul(ndf,nen,*),tl(*),eps(*)

      save

c     Set integration parameter: For energy = 1.0

      if(isw.eq.13) then
        alpha = 1.d0
      else
        alpha = ctan(1)
      endif
      alphar = 1.d0/alpha - 1.d0

c     Compute length and strain terms

      xlen0  = 0.0d0
      xlenn  = 0.0d0
      xlen1  = 0.0d0
      xlena  = 0.0d0
      eps(1) = 0.0d0
      eps(2) = 0.0d0
      do i = 1,ndm
        dx    = xl(i,2)   - xl(i,1)
        du    = ul(i,2,1) - ul(i,1,1)
        dd    = ul(i,2,2) - ul(i,1,2)
        xlen0 = xlen0 +  dx**2
        xlen1 = xlen1 + (dx + du + dd*alphar)**2
        xlenn = xlenn + (dx + du - dd       )**2
        xlena = xlena + (dx + du            )**2
        eps(1)= eps(1) +  dx*du
        eps(2)= eps(2) +  dx*(ul(i,2,4) - ul(i,1,4))
        xx(i) = (xl(i,2) + xl(i,1))*0.5d0
      end do

c     Compute temperature change

      ta   = 0.5d0*(tl(1) + tl(2)) - d(9)

c     Compute strain forms

      if(dtype.lt.0) then

c       Logarithmic stretch strain and deformed length

        eps(1) = 0.5d0*log(xlena/xlen0)
        xlen0 = sqrt(xlen0)
        xlen  = sqrt(xlena)

      else

c       Green or linear strain

        if(nonlin) then
          eps(1) = 0.5d0*xlena/xlen0 - 0.5d0
        else
          eps(1) = eps(1)/xlen0
          eps(2) = eps(2)/xlen0
        endif
        xlen0 = sqrt(xlen0)
        xlen  = xlen0
      endif

      end

      subroutine trcnnd(ix,sig,dt,st)

      implicit  none

      integer   i,ll

      integer   ix(*)
      real*8    dt(*),st(*),sig(*)

      save

      do i = 1,2

        ll = ix(i)
        if(ll.gt.0) then

          dt(ll) = dt(ll) + 1.d0

c         Stress projections

          st(ll) = st(ll) + sig(1)

        endif
      end do

      end