subroutine forcdt(u,v,a,dp,du,time,dtime,ndag,node,
     1ug,xord,ncrv,iacflg,jnc,ipess)
c
c CONTACT: ap@marc.de
c      
c* * * * *
c
c     input of time dependent forcing functions and boundary
c     conditions.
c
c     u            total displacements at a node
c     v            velocity at a node
c     a            acceleration at a node
c     dp           load increments at a node
c     du           displacement increments at a node
c     time         time
c     dtime        time increment
c     ndag         number of degrees of freedom per node
c     node         node number
c     ug           total displacements at node in global system
c     xord         orriginal coordinates
c     ncrv         number of coordinates
c     iacflg       acceleration flag - set to 1 if accelerations given
c     jnc          increment number
c
c* * * * *
c
c PURPOSE: Follow force type 1 for point loads.
c          This routine defines point loads at two nodes
c          in the direction along the line connecting these 
c          two nodes in the deformed configuration.
c          A positive force means the nodes attract each other.
c          A negative force means the nodes repel each other.
c
c NOTE: 1) Always make incrememnt 0 a null step.
c       2) This subroutine "forcdt" works in conjunction with
c          subroutines "ubginc", "uedinc" and "nodnum".
c       3) This subroutine does not work in conjunction with
c          the "auto increment" history definition option.
c
c USAGE: One must list all nodes which attract or repel each other
c        in the FORCDT model definition part of the MARC input.
c        Forcdt will be called for each node in this list.
c        For each node (i.e. each call) one must specify the node
c        numbers of the two nodes which interact in "nid1" and "nid2".
c        The magnitude of the interaction force at the start of each
c        increment must be defined as "ft" in subroutine ubginc.
c
c* * * * *
c
      implicit real*8 (a-h,o-z)                                               dp
c
      include '../common/space'
      include '../common/elmcom'
      include '../common/concom'
      include '../common/dimen'
      include '../common/prepro'
      include '../common/arrays'
      include '../common/heat'
      include '../common/array2'
      include '../common/strvar'
      include '../common/blnk'
      include '../common/lass'
c
      dimension u(ndag),v(ndag),a(ndag),dp(ndag),du(ndag)
      dimension ug(1),xord(1)
c
c* * * * *
c
c local varaiables:
c xcoor1(1..3) = coordinates of 1st node in deformed configuration.
c xcoor2(1..3) = coordinates of 2nd node in deformed configuration.
c rcoor1(1..3) = coordinates of 1st node in undeformed configuration.
c rcoor2(1..3) = coordinates of 2nd node in undeformed configuration.
c rdisp1(1..3) = total displacements of 1st node at start of increment.
c rdisp2(1..3) = total displacements of 2nd node at start of increment.
c ddisp1(1..3) = displacement changes of 1st node in current increment.
c ddisp2(1..3) = displacement changes of 2nd node in current increment.
c xl12(1..3)   = unit direction vector from 1st to 2nd node.
c
      dimension xcoor1(6),xcoor2(6),
     *          rcoor1(6),rcoor2(6),
     *          rdisp1(6),rdisp2(6),
     *          ddisp1(6),ddisp2(6),
     *          xl12(6)
c
c* * * * *
c
c save variables below in a new common block "force_com":
c nid1      = user node number of 1st node.
c nid2      = user node number of 2nd node.
c f01(1..3) = force vector at 1st node at start of increment.
c f02(1..3) = force vector at 2nd node at start of increment.
c f11(1..3) = force vector at 1st node in current iteration.
c f12(1..3) = force vector at 2nd node in current iteration.
c ft        = magnitude of the force at the two nodes.
c             (positive means the nodes attract each other
c              negative means the nodes repel each other).
c
      common/force_com/nid1,nid2,f01(6),f02(6),f11(6),f12(6),ft
c
c* * * * *
c
c define a flag for additional printing in output file.
      logical debug
      data debug/.false./
c 
      if (debug) write(6,*) 'ncycle:',ncycle
c
c for a 3 dimensional problem set ncrc to 3.
c for a 2 dimensional problem set ncrc to 2.
      ncrc = 3
c
c get the internal node number for 1st node in nin1.
      nid1 = 2
      call nodnum(nid1,nex1,nin1,1)
c get the internal node number for 2nd node in nin2.
      nid2 = 4
      call nodnum(nid2,nex2,nin2,1)
c
c get the undeformed coordinates of the 1st node in rcoor1.
      jrdpre=0
      call vecftc(rcoor1,vars(ixord),ncrdmx,ncrd,nin1,jrdpre,2,1)
c get the total displacements at the start of the increment of the
c 1st node in rdisp1.
      jrdpre=0
      call vecftc(rdisp1,vars(idsxts),ndegmx,ndeg,nin1,jrdpre,2,5)
c get the incremental displacements in the current increment of
c the 1st node in ddisp1.
      jrdpre=0
      call vecftc(ddisp1,vars(idsx),ndegmx,ndeg,nin1,jrdpre,2,5)
c get the undeformed coordinates of the 2nd node in rcoor2.
      jrdpre=0
      call vecftc(rcoor2,vars(ixord),ncrdmx,ncrd,nin2,jrdpre,2,1)
c get the total displacements at the start of the increment of the
c 2nd node in rdisp2.
      jrdpre=0
      call vecftc(rdisp2,vars(idsxts),ndegmx,ndeg,nin2,jrdpre,2,5)
c get the incremental displacements in the current increment of
c the 2nd node in ddisp2.
      jrdpre=0
      call vecftc(ddisp2,vars(idsx),ndegmx,ndeg,nin2,jrdpre,2,5)
c
c compute the coordinates of 1st and 2nd node in the current deformed
c configuration.
      if (ncycle.eq.0) then
       do i1=1,ncrc
        xcoor1(i1) = rcoor1(i1) + rdisp1(i1)
        xcoor2(i1) = rcoor2(i1) + rdisp2(i1)
       enddo
      else
       do i1=1,ncrc
        xcoor1(i1) = rcoor1(i1) + rdisp1(i1) + ddisp1(i1)
        xcoor2(i1) = rcoor2(i1) + rdisp2(i1) + ddisp2(i1)
       enddo
      endif
c
      if (debug) then
       write(6,20) nid1,(rcoor1(ijk),ijk=1,ncrc)
       write(6,20) nid1,(rdisp1(ijk),ijk=1,ncrc)
       write(6,20) nid1,(ddisp1(ijk),ijk=1,ncrc)
       write(6,20) nid1,(xcoor1(ijk),ijk=1,ncrc)
       write(6,20) nid2,(rcoor2(ijk),ijk=1,ncrc)
       write(6,20) nid2,(rdisp2(ijk),ijk=1,ncrc)
       write(6,20) nid2,(ddisp2(ijk),ijk=1,ncrc)
       write(6,20) nid2,(xcoor2(ijk),ijk=1,ncrc)
      endif
 20   format(i5,3e14.5)
c
c compute the unit direction vector from 1st to 2nd node.
      dlx = 0d0
      do i1=1,ncrc
       xl12(i1) = xcoor2(i1) - xcoor1(i1)
       dlx = dlx + xl12(i1)*xl12(i1)
      enddo
      dlx = dsqrt(dlx)
      if (dlx.gt.0d0) then
       do i1=1,ncrc
        xl12(i1) = xl12(i1)/dlx
       enddo
      endif
c
c initialize the incremental load vector.
      do i1=1,ndeg
       dp(i1) = 0d0
      enddo
c define the force increments at 1st node.
      if (node.eq.nid1) then
       do i1=1,ncrc
        f11(i1) = ft*xl12(i1)
        dp(i1) = f11(i1) - f01(i1)
       enddo
      endif
c define the force increments at 2nd node.
      if (node.eq.nid2) then
       do i1=1,ncrc
        f12(i1) = -ft*xl12(i1)
        dp(i1) = f12(i1) - f02(i1)
       enddo
      endif
      if (debug) write(6,20) node,(dp(ijk),ijk=1,ncrc)
c
      return
      end
      subroutine ubginc(inc,incsub)
      implicit real*8 (a-h,o-z)
c
c     user subroutine that gets called at the start of each increment.
c
      common/force_com/nid1,nid2,f01(6),f02(6),f11(6),f12(6),ft
c
c initialize vectors in increment 0.
      if (inc.eq.0) then
       do i1=1,3
        f01(i1) = 0d0
        f02(i1) = 0d0
        f11(i1) = 0d0
        f12(i1) = 0d0
       enddo
      endif
c
c define the magnitude of the force between the nodes at the start of
c a new increment.
      df = 0.1d0
      ft = dfloat(inc)*df
c
      return
      end
      subroutine uedinc(inc,incsub)
      implicit real*8 (a-h,o-z)
c
c     user subroutine that gets called at the end of each increment.
c
      common/force_com/nid1,nid2,f01(6),f02(6),f11(6),f12(6),ft
c
c update the force vectors f01 and f02.
c the converged force vectors f11 and f22 at the end of this increment
c become the initial vectors f01 and f02 of the next increment.
      do i1=1,3
       f01(i1) = f11(i1)
       f02(i1) = f12(i1)
      enddo
c
      return
      end
      subroutine nodnum(noid,noex,noin,ic)
c
      implicit real*8 (a-h,o-z)
c
c     ic=1   determine noex and noin for noid
c     ic=2   determine noin and noid for noex
c     ic=3   determine noid and noex for noin
c
c     noid = noid id (given by user)
c     noex = marc node number before optimize
c     noin = marc node number after optimize
c
      include '../common/dimen'
      include '../common/space'
      include '../common/arrays'
      include '../common/develp'
      include '../common/prepro'
c
      if(ic.eq.1) then
       noex=noid
       if(nnoids.ne.0) noex=ibsrch(noid,ints(inoids),numnp,1)
       noin=noex
       if(joptit.ne.0) noin=igetsh(ints(inpnum+noex-1),0)
      else if(ic.eq.2) then
       noid=noex
       if(nnoids.ne.0) noid=ints(inoids+noex-1)
       noin=noex
       if(joptit.ne.0) noin=igetsh(ints(inpnum+noex-1),0)
      else if(ic.eq.3) then
       noex=noin
       if(joptit.ne.0) noex=igetsh(ints(inpnum+noin-1),1)
       noid=noex
       if(nnoids.ne.0) noid=ints(inoids+noex-1)
      endif
c
      return
      end