c     user element for 3D AXIAL connector
      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
      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)

      if (jtype .eq. 1001 .and.
     *    lflags(iProcedure).eq.jDynExplicit) then 

         conn_stiffness    = props(1)
         conn_dampingcoeff = props(2)
         conn_lref         = props(3)
         conn_mass         = props(4)

         if (dTimeCur.eq.zero) then 
            dtInv   = one
         else
            dtInv = one/dTimeCur
         end if

c        mass at each node
         am0=half*conn_mass

         if ( lflags(iOpCode).eq.jMassCalc ) then
            
            do kblock = 1, nblock
               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
               
c              initial length of the connector
               aLengthX0=coords(kblock,2,1)-coords(kblock,1,1)
               aLengthY0=coords(kblock,2,2)-coords(kblock,1,2)
               aLengthZ0=coords(kblock,2,3)-coords(kblock,1,3)

               aLength0=sqrt(aLengthX0*aLengthX0+aLengthY0*aLengthY0+
     *                       aLengthZ0*aLengthZ0)

c              deformed length

               aLengthX=aLengthX0+(u(kblock,4)-u(kblock,1))
               aLengthY=aLengthY0+(u(kblock,5)-u(kblock,2))
               aLengthZ=aLengthZ0+(u(kblock,6)-u(kblock,3))

               aLength=sqrt(aLengthX*aLengthX+aLengthY*aLengthY+
     *                      aLengthZ*aLengthZ)

c              internal force calcualtions

               if(conn_lref.eq.zero) then
                    aLengthRef=aLength0
               else
                    aLengthRef=conn_Lref
               endif

               u1_mat = aLength-aLengthRef
               
               aLengthOld = svars(kblock,6)
               rel_vel1 = (aLength - aLengthOld)*dtInv

c              force in the spring
               force1_el=conn_stiffness*u1_mat
c              force in the damper
               force1_vis=conn_dampingcoeff*rel_vel1
c              total force in the connector
               force1=force1_el+force1_vis

c              component of force1
               force1X=force1*(aLengthX/aLength)
               force1Y=force1*(aLengthY/aLength)
               force1Z=force1*(aLengthZ/aLength)

c              assemble internal load in RHS
               rhs(kblock,1) = -force1X
               rhs(kblock,2) = -force1Y
               rhs(kblock,3) = -force1Z
               rhs(kblock,4) =  force1X
               rhs(kblock,5) =  force1Y
               rhs(kblock,6) =  force1Z

c              energy calculation 

               u1_mat_old     =   svars(kblock,1)      
               force1_el_old  =   svars(kblock,2)
               force1_vis_old =   svars(kblock,3)
      

c              internal energy
               energy(kblock, iElIe) = energy(kblock, iElIe) +
     *          half*(force1_el+force1_el_old)*(u1_mat-u1_mat_old)
c              viscous dissipation
c               energyViscOld = svars(kblock,7)
c               energyViscNew = energyViscOld +
c     $              half*(force1_vis+force1_vis_old)*rel_vel1*dTimeCur
c               energy(kblock, iElDd) = energy(kblock, iElDd) +
c     $              energyViscNew
               energy(kblock, iElDd) = energy(kblock, iElDd) +
     *          half*(force1_vis+force1_vis_old)*(u1_mat-u1_mat_old)     
c               energy(kblock, iElDd) = energy(kblock, iElDd) +
c     *          half*(force1_vis+force1_vis_old)*(aLength-aLengthOld) 
               
c              update state variables
               svars(kblock,1) = u1_mat     
               svars(kblock,2) = force1_el
               svars(kblock,3) = force1_vis


c              for output of connector constitutive displacement (CCU)
c                   in the 1-direction (along conn axis)
               svars(kblock,4) = force1
               svars(kblock,5) = u1_mat

c              lenght is a state variable (CP); so is the dissip energy
               svars(kblock,6) = aLength
c               svars(kblock,7) = energyViscNew
               

c              undamped freq
               f_undamp=sqrt(conn_stiffness/conn_mass)
c              stable time increment in the element (2*1/frequency)
c              the actual stable time increment would also depend
c              on damping; however it is ignored here
               dtimeStable(kblock) = two/f_undamp

            end do
         end if

      end if
c
      return
      end