c     user element for truss along the X global axis
      subroutine vuel(
     *     nblock,
c          to be defined
     *     rhs,amass,dtimeStable,
     *     svars,nsvars,
     *     energy,
c          
     *     nnode,ndofel,
     *     props,nprops,
     *     jprops,njprops,
     *     coords,ncrd,
     *     u,du,v,a,
     *     jtype,jelem,
     *     time,period,dtimeCur,dtimePrev,kstep,kinc,lflags,
     *     dMassScaleFactor,
     *     predef,npredef,
     *     jdltyp,adlmag)

C     
      include 'vaba_param.inc'

      parameter ( zero = 0.d0, half = 0.5d0, one = 1.d0, two=2.d0 )
      parameter (scaleTemp = 0.9d0)

c     operation code
      parameter ( jMassCalc            = 1,
     *            jIntForceAndDtStable = 2,
     *            jExternForce         = 3)

c     flags
      parameter (iProcedure = 1,
     *           iNlgeom    = 2,
     *           iOpCode    = 3,
     *           nFlags     = 3)

c     time
      parameter (iStepTime  = 1,
     *           iTotalTime = 2,
     *           nTime      = 2)

c     procedure flags
      parameter ( jDynExplicit = 17 )

c     energies 
      parameter ( iElPd = 1,
     *            iElCd = 2,
     *            iElIe = 3,
     *            iElTs = 4,
     *            iElDd = 5,
     *            iElBv = 6,
     *            iElDe = 7,
     *            iElHe = 8,
     *            iElKe = 9,
     *            iElTh = 10,
     *            iElDmd = 11,
     *            iElDc = 12,
     *            nElEnergy = 12)


c     predefined variables
      parameter ( iPredValueNew = 1,
     *            iPredValueOld = 2,
     *            nPred         = 2)    

c     indexing in a 3-long vector

      parameter (factorStable = 0.99d0)
**------------------------------------
C
      dimension rhs(nblock,ndofel), amass(nblock,ndofel,ndofel),
     *     dtimeStable(nblock),
     *     svars(nblock,nsvars), energy(nblock,nElEnergy),
     *     props(nprops), jprops(njprops),
     *     jelem(nblock), time(nTime), lflags(nFlags),
     *     coords(nblock,nnode,ncrd), u(nblock,ndofel),
     *     du(nblock,ndofel), v(nblock,ndofel), a(nblock, ndofel),
     *     dMassScaleFactor(nblock),
     *     predef(nblock, nnode, npredef, nPred), adlmag(nblock)

c     VUEL subroutine for a truss element along the global X-axis only
c     superimposed with rotations at both ends
c     Limitations: 
c       It will not work if motions along another direction other than X
c       Nodes must be noncoincident
c     Notes: 
c        Define only nonzero entries; the arrays to be defined have
c        been zeroed out just before this call
c        for jInternForceOnly, one can use state vars to 
c        compute initial forces at time zero (prestress)

      if (jtype .eq. 2) then 

         area0    = props(1)
         rho      = props(2)
         conduct  = props(3)
         specHeat = props(4)
         diffusivity = conduct/(rho*specHeat)

         if ( lflags(iOpCode).eq.jMassCalc ) then
            do kblock = 1, nblock
c              use original distance to compute mass
               alenX0 = (coords(kblock,2,1) - coords(kblock,1,1))
               alenY0 = (coords(kblock,2,2) - coords(kblock,1,2))
               alenZ0 = (coords(kblock,2,3) - coords(kblock,1,3))
               alen0 = sqrt(alenX0*alenX0 + alenY0*alenY0 +
     *                      alenZ0*alenZ0)
               Vol = alen0*area0   
               bodyOp1 = half*Vol
               bodyOp2 = half*Vol
               amass(kblock,4,4) = bodyOp1*specHeat*rho
               amass(kblock,8,8) = bodyOp2*specHeat*rho

c              dummy structural mass
               amass(kblock,1,1) = one
               amass(kblock,2,2) = one
               amass(kblock,3,3) = one
               amass(kblock,5,5) = one
               amass(kblock,6,6) = one
               amass(kblock,7,7) = one
            end do
         else if ( lflags(iOpCode) .eq.
     *             jIntForceAndDtStable) then
            do kblock = 1, nblock
               alenX0 = (coords(kblock,2,1) - coords(kblock,1,1))
               alenY0 = (coords(kblock,2,2) - coords(kblock,1,2))
               alenZ0 = (coords(kblock,2,3) - coords(kblock,1,3))
               alen0 = sqrt(alenX0*alenX0 + alenY0*alenY0 +
     *                      alenZ0*alenZ0)
               Vol     = alen0*area0   

               alen0Inv = one/alen0
               
c              compute "stiffness" matrix
               Vol   = alen0*area0
               fact  = Vol*conduct
               dKf11 =  (alen0Inv*one)* (alen0Inv*one)*fact
               dKf22 = (-alen0Inv*one)*(-alen0Inv*one)*fact
               dKf12 = (-alen0Inv*one)* (alen0Inv*one)*fact
               dKf21 =  dKf12
               
c              compute element force K_f^PQ*P^Q in the Theory Man
               t1    = u(kblock,4)
               t2    = u(kblock,8)
               eForce1 = dKf11*t1 + dKf12*t2
               eForce2 = dKf21*t1 + dKf22*t2
               
c              assemble internal load in RHS
               rhs(kblock,1) = -eForce1
               rhs(kblock,2) = -eForce2

c              stable time increment
               dtimeStable(kblock) = scaleTemp*
     *              alen0*alen0/(two*diffusivity)

c              dummy energy assigments
               energy(kblock, iElKe) = zero
               energy(kblock, iElIe) = zero
            end do
         end if
      end if
c
      return
      end