c     user element for truss along the X global axis
c     the Young's modulus is made a function of temperature
      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 )

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)
         eMod1  = props(2)
         temp1    = props(3)
         eMod2  = props(4)
         temp2    = props(5)
         anu   = props(6)
         rho   = props(7)

         eDampTra    = zero         
         amassFact0  = half*area0*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)
               am0   = amassFact0*alen0
               amass(kblock,1,1) = am0
               amass(kblock,2,2) = am0
               amass(kblock,3,3) = am0
               amass(kblock,4,4) = am0
               amass(kblock,5,5) = am0
               amass(kblock,6,6) = am0

            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)
               vol0 = area0*alen0
               amElem0 = two*amassFact0*alen0

               alenX = alenX0
     *              +  (u(kblock,4)        -  u(kblock,1))
               alenY = alenY0
     *              +  (u(kblock,5)        -  u(kblock,2))
               alenZ = alenZ0
     *              +  (u(kblock,6)        -  u(kblock,3))
               alen = sqrt(alenX*alenX + alenY*alenY + alenZ*alenZ)
               area   = vol0/alen

c              linear interpolation to determine the value of E at the current value of 
c              temperature. It the temp is outside the range specified, E is constant
              
               temp_node1=predef(kblock,1,1,iPredValueNew)
               temp_node2=predef(kblock,2,1,iPredValueNew)
               temp=half*(temp_node1+temp_node2)

               if(temp.lt.temp1) then
                  eMod=eMod1
               else if(temp.gt.temp2) then
                  eMod=eMod2
               else
                  eMod=eMod1+(eMod2-eMod1)*(temp-temp1)/(temp2-temp1)                 
               end if

               ak     = area*eMod/alen

c              undamped stable time increment for translations
               dtTrialTransl = sqrt(amElem0/ak)

c              damped stable time increment; since eDampTra=0, the
c              stable time increment does not change because of damping
               critDampTransl = two*sqrt(amElem0*ak)
               csiTra = eDampTra/critDampTransl
               factDamp = sqrt(one+csiTra*csiTra) - csiTra
               dtTrialTransl = dtTrialTransl*factDamp*factorStable
               dtimeStable(kblock) = dtTrialTransl

c              force = E * logarithmic strain *current area
               strainLog = log(alen/alen0)
               fElasTra = eMod*strainLog*area
               
               forceTra = fElasTra

c              assemble internal load in RHS
               rhs(kblock,1) = -forceTra
               rhs(kblock,4) =  forceTra

c              internal energy calculation 
               alenOld  = svars(kblock,1)
               fElasTraOld = svars(kblock,2)
               
               energy(kblock, iElIe) = energy(kblock, iElIe) +
     *              half*(fElasTra+fElasTraOld)*(alen - alenOld) 

c              update state variables
               svars(kblock,1) = alen
               svars(kblock,2) = fElasTra
               
            end do
         end if
      end if
c
      return
      end