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 (epscritical = 2.d2, scaleAcous = 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 .and.
     *    lflags(iProcedure).eq.jDynExplicit) then 

         area0  = props(1)
         rho    = props(2)
         bulkMod= props(3)
         volDrag= props(4)

         bulkModInv = one/bulkMod
         rhoInv     = one/rho

         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)

c              \int V_f H^P H^Q dV in the Theory Man
               Vol = alen0*area0
               bodyOp1 = half*Vol
               bodyOp2 = half*Vol
c              M_f^PQ in the Theory Man
               amass(kblock,1,1) = bodyOp1*bulkModInv
               amass(kblock,2,2) = bodyOp2*bulkModInv
            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
               bodyOp1 = half*Vol
               bodyOp2 = half*Vol
               dMass1  = bodyOp1*bulkModInv
               dMass2  = bodyOp2*bulkModInv

               alen0Inv = one/alen0
               
c              compute K_f^PQ = \int_V_f H^P,x H^Q,x in the theory man
               Vol   = alen0 * area0
               fact  = Vol*rhoInv
               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
               p1    = u(kblock,1)
               p2    = u(kblock,2)
               p1Dot = v(kblock,1)
               p2Dot = v(kblock,2)
               eForce1 = dKf11*p1 + dKf12*p2
               eForce2 = dKf21*p1 + dKf22*p2
               
c              assemble internal load in RHS
               rhs(kblock,1) = eForce1
               rhs(kblock,2) = eForce2

c              energies calculation - assume linear
c              using a state variable to do an average with the old value 
c              would be better
               elemDamp   = bulkModInv*volDrag*rhoInv
               energy(kblock, iElIe) = eForce1*p1 + eForce2*p2
c               energy(kblock, iElDd) =  elemDamp*
c     *                 (bodyOp1*p1*p1 + bodyOp2*p2*p2)
               energy(kblock, iElKe) = dMass1*p1*p1Dot +
     *                                      dMass2*p2*p2Dot

c              damped stable time increment; since eDampTra=0, the
c              stable time increment does not change because of damping
               eps = volDrag*sqrt(bulkModInv*rhoInv)
               epsfactor = sqrt( one + eps * eps ) - eps
               if(eps .ge. epscritical) epsfactor = half / eps
               dtimeStable(kblock) = scaleAcous*alen0*epsfactor
     *                              *sqrt(bulkModInv/rhoInv)
            end do
         end if
      end if
c
      return
      end