c***********************************************************************
c    A STATIC, VISCO AND DYNAMIC PROCEDURE VERSION OF THE
c    COHESIVE ELEMENT LIBRARY
c
c    Version 4.0 
c    ------------
c   
c    Date:
c    -----
c     3 March, 2000
c
c    Developed by:
c    ------------
c    Pakal Rahul Kumar
c    Oak Ridge National Laboratory
c
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c    NOTES:
c    -----
c
c    The integration for the stiffness matirx and rhs load vector
c    for all the cohesive elements are done over the original area
c
c    All calculations for these elements are done over the mid-surface
c    of the cohesive element in the deformed state.
c
c    The mode I form of the formulation is implemented for all
c    the cohesive elements. It is assumed in this implementaion that
c    the bottom side/face of the cohesive element lies along the plane
c    of symmetry.
c
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c    THE SAVED VARIBLES REQUIREMENT FOR VARAIOUS ELEMENT TYPES
c    AS FOLLOWS:
c
c    STATIC and VISCO PROCEDURES
c    ---------------------------
c    4-node 2D Element Plane Strain/Stress  TYPE=U1    VARIABLES = 12
c    4-node 2D Element Axisymmetric         TYPE=U2    VARIABLES = 18
c    6-node 2D Element plane Strain/Stress  TYPE=U3    VARIABLES = 24
c    6-node 2D Element Axisymmetric         TYPE=U4    VARIABLES = 30
c    8-node 3D Element                      TYPE=U5    VARIABLES = 81
c    12-node 3D Element                     TYPE=U6    VARIABLES = 138
c
c    DYNAMIC PROCEDURE
c    -----------------
c    4-node 2D Element Plane Strain/Stress  TYPE=U1    VARIABLES = 28
c    4-node 2D Element Axisymmetric         TYPE=U2    VARIABLES = 34
c    6-node 2D Element Plane Strain/Stress  TYPE=U3    VARIABLES = 48
c    6-node 2D Element Axisymmetric         TYPE=U4    VARIABLES = 54
c    8-node 3D Element                      TYPE=U5    VARIABLES = 138
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c MATERIAL PROPERTIES FOR THE COHESIVE ELEMENTS ARE INPUT
c USING THE REAL PROPERTIES OF THE USER ELEMENT.
c THE FOLLOWING FORMS OF POTENTIALS ARE IMPLEMENTED HERE.
c
c --------------------------------------------------------------
c Xu & Needleman Cohesive Zone Model.
c Ref: Numerical Simulations of Fast Crack Growth in
c      Brittle Solids. JMPS VOl. 42, No. 9, pp. 1397-1434, 1994
c
c   2D Elements: (5 real properties + 4 integer properties)
c   --------
c
c    props(1) = phi_n
c    props(2) = phi_t
c    props(3) = del_cr_n
c    props(4) = del_cr_t
c    props(5) = r
c
c   3D Elements: (4 real properties + 4 integer properties)
c   --------
c    props(1) = phi_n
c    props(2) = del_cr_n
c    props(3) = del_cr_t1
c    props(4) = del_cr_t2
c
c --------------------------------------------------------------
c Needleman's Exponential form of Cohesive Zone Model.
c Ref: An Analysis of Decohesionalong an imperfect interface
c Int. J. of Fracture, VOl. 42, pp. 21-40, 1990.
c Equation: 3.7
c
c NOTE: NOT IMPLEMENTED FOR 3D.
c
c
c   2D Elements:  (4 real properties + 4 integer properties)
c   -------------
c      props(1) = phi_n
c      props(2) = phi_t
c      props(3) = del_n
c      props(4) = del_n
c      Note variable ALPHA in the equation is set to 1.0
c      alpha = 1.d0
c--------------------------------------------------------------
c
c Leonard-Jone's Potential. (-A/u^2 + B/u^8)*exp(-v^2)
c
c Where: u = delta_n/del_cr, v = delta_t/del_cr
c        A, and B are constants. 
c
c In this implementation user can define the following 
c properties:
c
c   fracture energy, phi
c   critical opening, del_cr
c   factor, A.
c
c   factor B is evaluated in terms of phi, and A.
c
c   2D and 3D elements: (3 real properties + 4 integer properties)
c   ------------------
c      props(1) = phi_n (fracture energy)
c      props(2) = del_cr (critical opening displacement)
c      props(3) = factor A.
c
c------------------------------------------------------------------
c
c  THE 4 INTEGER PROPERTIES FOR EACH MODEL TO BE SPECIFIED
c  ARE DESCRIBED BELOW:
c  -------------------------------------------------
c
c    jprops(1) 
c        1 = Invokes penalty augmentation in compression
c        0 = No penatly augmentation in compression.
c
c    jprops(2) = Model flag (1, 2, or 3)
c
c        1 = Xu & Needleman
c        2 = Exponential potential (Implemented for 2d elements only)
c        3 = L-J Potential
c
c    jprops(3) = 1
c        Invokes the mode I crack formulation.
c        The integration is performed over original area
c        and the normal does not change with the deforming
c        cohesive element. The normal to the coheive element
c        is the normal to the symmetry plane.
c        NOTE : THE Bottom face of cohesive element is taken
c               as the plane of symmetry.
c
c    jprops(4) = 1
c        Prints the element stiffness and internal force 
c        vector for every element for every iteration to the
c        jobname.dat file.
c
c
c***********************************************************************
c
       subroutine uel(rhs,amatrx,svars,energy,ndofel,nrhs,nsvars,
     x props,nprops,coords,mcrd,nnode,u,du,v,a,jtype,time,dtime,
     x kstep,kinc,jelem,params,ndload,jdltyp,adlmag,predef,npredf,
     x lflags,mlvarx,ddlmag,mdload,pnewdt,jprops,njprop,period)
c
c Functions:
c ---------
c   Main driver routine for the Cohesive elements for 2D and
c   3D modeling. Call the respective element formulation routines
c
c
      include 'ABA_PARAM.INC'
c
      dimension rhs(mlvarx,*),amatrx(ndofel,ndofel),props(*),
     x svars(*),energy(8),coords(mcrd,nnode),u(ndofel),
     x du(mlvarx,*),v(ndofel),a(ndofel),time(2),params(*),
     x jdltyp(mdload,*),adlmag(mdload,*),ddlmag(mdload,*),
     x predef(2,npredf,nnode),lflags(*),jprops(*)
c
      logical procflags(6),mattype(2)
c
      parameter(ne=9600)
c
      common /kwrite/kfirst1,kincold,kcount2,
     x gpdelt(ne*2),xloc(ne*2),gptract(ne*2)
c
c --- Check for the procedure type being run currrently
      iproc=0
c -- STATIC
      if(lflags(1).eq.1.or.lflags(1).eq.2)then
c       --  NOT A STATIC, RIKS procedure
        if(nrhs.ne.2)iproc=1
      endif
c -- DYNAMIC
      if(lflags(1).eq.11.or.lflags(1).eq.12)iproc=1
c -- VISCO, Implicit
      if(lflags(1).eq.22)iproc=1
c -- FREQENCY, analysis
      if(lflags(1).eq.41)iproc=1
c
      if(iproc.eq.0)then
        write(7,*)'**** ERROR The user element is not written'
        write(7,*)'     for this procedure type'
        stop
      endif
c
c
c --- Based on the solution procedure set some variables

      do i = 1, 6
         procflags(i) = .false.
      enddo
      mattype(1) = .false.
      mattype(2) = .false.
c
c     -- STATIC and VISCO
      if (lflags(3) .eq. 1 .and.
     x        (lflags(1) .eq. 1 .or. lflags(1) .eq. 2 .or.
     x            lflags(1) .eq. 22) ) then
         procflags(1) = .true.
c     -- DYANAMIC Normal time incrementaion.
      elseif (lflags(3) .eq. 1 .and.
     x        (lflags(1) .eq. 11 .or. lflags(1) .eq. 12) ) then
         procflags(2) = .true.
c     -- DYNAMIC Half step residual calculation.
      elseif ( lflags(3) .eq. 5) then
         procflags(3) = .true.
c     -- DYNAMIC Initial accelration calculation
      elseif ( lflags(3) .eq. 6) then
         procflags(4) = .true.
c     -- DYNAMIC Velocity jump calculation
      elseif ( lflags(3) .eq. 4 ) then
         procflags(5) = .true.
         call kvzero (amatrx(1,1),ndofel*ndofel)
         goto 999
      elseif (lflags(3) .eq. 2) then
c     --- just stiffness matrix
c         write(7,101)lflags(1), jelem
c 101     format('stiff: lflags 1 and jelem = ',2I6)
         mattype(1) = .true.
      elseif (lflags(3) .eq. 4) then
c     --- just mass matrix (zero for all cohesive elements)
c         write(7,102)lflags(1), jelem
c 102     format('mass: lflags 1 and jelem = ',2I6)
         call kvzero (amatrx(1,1),ndofel*ndofel)
         goto 999
      endif
c
c --- Find the element type and branch out the relevent element
      if (jtype .eq. 1 .or. jtype .eq. 2) then
c         --- 4 Noded 2D Linear Element...
          kax = 0
          if (jtype .eq. 2) kax = 1
          call kelem2D4Node (kax,procflags,mattype,ndofel,nnode,
     x njprop,nrhs,nsvars,mlvarx,mcrd,jprops(1),props(1),coords(1,1),
     x params(1),svars(1),u(1),amatrx(1,1),rhs(1,1))
c
      elseif (jtype .eq. 3 .or. jtype .eq. 4) then
c         --- 6 Noded 2D Quadratic Element...
          kax = 0
          if (jtype .eq. 4) kax = 1
          call kelem2D6Node (kax,procflags,mattype,ndofel,nnode,njprop,
     x nrhs,nsvars,mlvarx,mcrd,jprops(1),props(1),coords(1,1),params(1),
     x svars(1),u(1),amatrx(1,1),rhs(1,1))
      elseif (jtype .eq. 5) then
c         --- 8 Noded 3D Linear Element...BRick-to-Brick
          call kelem3D8Node (procflags,mattype,ndofel,nnode,njprop,nrhs,
     x nsvars,mlvarx,mcrd,jprops(1),props(1),coords(1,1),svars(1),u(1),
     x params(1),amatrx(1,1),rhs(1,1))
      elseif (jtype .eq. 6) then
c         --- 12 Noded 3D Quadratic Element..Tet-to-Tet 
          call kelem3D12Node (procflags,mattype,ndofel,nnode,njprop,nrhs,
     x nsvars,mlvarx,mcrd,jprops(1),props(1),coords(1,1),svars(1),u(1),
     x params(1),amatrx(1,1),rhs(1,1))
      endif
c
c --- Perform stiffness matrix and internal force vector output
c     if jprops(4)=1
c
      if (jprops(4) .eq. 1) then

         write (6,900)  jelem
 900     format (/6x,'STIFFNESS MATRIX FOR ELEMENT ',i8)
         do i = 1,ndofel
            n = ndofel
            if (ndofel.gt.6) n = 6
            write (6,1000) i,(amatrx(i,j),j=1,n)
 1000       format (1x,i3,1x,6e15.7)
            if (ndofel.gt.6) write (6,1100) (amatrx(i,j),j=7,ndofel)
 1100       format (5x,6e15.7)
         enddo
c
         write (6,1200) jelem,(rhs(j,1),j=1,ndofel)
 1200    format (/6x,'INTERNAL FORCE VECTOR FOR ELEM',i8/
     x           (5x,6e15.7))
      endif
c
 999  return
      end
c
c ==================================================================
      subroutine kelem2D4Node (kax,procflags,mattype,ndofel,
     x numnode,njprop,nrhs,nsvars,mlvarx,mcrd,jprops,props,coords,
     x params,svars,u,amatrx,rhs)
c
c --- 2D 4 Node Cohesive Element---------------------------------
c
      include 'ABA_PARAM.INC'
c
c Function:
c --------
c   Compute element matrices, residual load vectors,
c    --- Arguments
      dimension
     x rhs(*),amatrx(ndofel,ndofel),props(*),
     x svars(*),coords(mcrd,numnode),u(ndofel),
     x jprops(*),params(*)
c
      logical procflags(*),mattype(*)
c
c     --- Local Variables
      dimension
     x TransMat(8,8), Shapes(4),
     x xcur(2,4), dnx(2),dny(2),dnz(2),
     x work(64),xcurRot(2,4),stiff1(8,8),
     x residual1(8),GPdelta(2),delta_n(2),delta_t(2),
     x AMat(2,8),Trac(2),Trac_Jacob(2,2),
     x stiff2(8,8),residual2(8),GP_Trac(2,2),
     x ShapeDerv(2),xcurMid(2,2),
     x x2cor(2),x3cor(3),wt2(2),wt3(3),
     x GPts(1),rhssave(8),x1cor(1),wt1(1),r8temp1(8),r8temp2(8)
c
      logical staticproc,dynaproc,hsresd,veljump,initialacel,
     x stiffmat
c
      data x1cor /0.d0/
      data wt1 /2.d0/
c
      data x2cor /-0.5773502691896d0, 0.5773502691896d0/
      data wt2  /1.0d0, 1.0d0/
c
      data x3cor  / -0.774596669241483d0, 0.0d0, 0.774596669241483d0 /
c
      data wt3  / 0.555555555555556d0, 0.8888888888888889d0,
     x            0.555555555555556d0 /
c
      data dPI /3.1415927d0/
c
c
c            4   side 2   3
c            +------------+
c            |            |
c            |            |
c            +------------+
c            1   side 1   2
c
      keshp = 0
      kemid = 0
      kmodeI = 0
c
      if (jprops(3) .eq. 1) kmodeI = 1
c
c --- Update element nodal coordinates, x = X + dx operation
      do inode = 1, 4
         do i = 1, 2
            xcur(i,inode) = coords(i,inode) + u(2*(inode-1)+i)
         enddo
      enddo
c
c --- Compute the coordinates of the mid surface...
      if (kmodeI .eq. 0) then
         do i = 1, 2
            xcurMid(i,1) = 0.5d0*(xcur(i,1) + xcur(i,4))
            xcurMid(i,2) = 0.5d0*(xcur(i,2) + xcur(i,3))
         enddo
      elseif (kmodeI .eq. 1) then
         do i = 1, 2
            xcurMid(i,1) = coords(i,1)
            xcurMid(i,2) = coords(i,2)
         enddo
      endif
c
c -- Compute the deformed length
c
c
      dlength = (xcurMid(1,1) - xcurMid(1,2))**2 +
     x            (xcurMid(2,1) - xcurMid(2,2))**2
      dlength = dlength**0.5
c
c --- Compute the original length
      dlength0 = (coords(1,1)-coords(1,2))**2 +
     x              (coords(2,1)-coords(2,2))**2
      dlength0 = dlength0**0.5
c
c -- Compute the direction cosines for the mid-plane

      dnx1 = (xcurMid(1,2) - xcurMid(1,1))/dlength
      dny1 = (xcurMid(2,2) - xcurMid(2,1))/dlength
c
      dnx(1) = dnx1
      dny(1) = dny1
      dnz(1) = 0.d0
      dnx(2) = -dny(1)
      dny(2) = dnx(1)
      dnz(1) = 0.d0
c
c --- Form the current nodal coordinates in the element mid-plane
c     normal-tangential coordinate system.
c
      do inode = 1, 4
         xcurRot(1,inode) = dnx(1)*xcur(1,inode) + dny(1)*xcur(2,inode)
         xcurRot(2,inode) = dnx(2)*xcur(1,inode) + dny(2)*xcur(2,inode)
      enddo
c
c --- Compute the normal and tamgential openings between nodes
c     1 & 4 and 2 & 3
      delta_n(1) = xcurRot(2,4) - xcurRot(2,1)
      delta_n(2) = xcurRot(2,3) - xcurRot(2,2)
      delta_t(1) = xcurRot(1,4) - xcurRot(1,1)
      delta_t(2) = xcurRot(1,3) - xcurRot(1,2)
c
c --- Get the Transformation matrix from global to the
c     element normal-tangential coordinate system.
      call kFormTransMatrix (keshp,kemid,dnx(1),dny(1),dnz(1),8,
     x TransMat(1,1))
c
c --- Iniitlaizations....
      nini = 2
      ipt1 = 5
      ipt2 = 9
      ipt3 = 13
      ipt4 = 21
      if (kax .eq. 1) then
         nini = 3
         ipt1 = 7
         ipt2 = 13
         ipt3 = 19
         ipt4 = 27
      endif
c
      alpha = params(1)
      staticproc = procflags(1)
      dynaproc = procflags(2)
      hsresd = procflags(3)
      initialacel = procflags(4)
c
      stiffmat = .false.
      if (mattype(1)) stiffmat = .true.
c
      if (staticproc .or. stiffmat) then
         aconst = 1.d0
      elseif (dynaproc .or. hsresd) then
         aconst = 1.d0 + alpha
      endif
c
      call kvzero (amatrx(1,1),ndofel*ndofel)
      call kvzero (rhs(1),ndofel)
c
c --- Begin numerical integration....
c
      do ini = 1, nini
         if (nini .eq. 2) then
            GPts(1) = x2cor(ini)
            GaussWt = wt2(ini)
         elseif (nini .eq. 3) then
            GPts(1) = x3cor(ini)
            GaussWt = wt3(ini)
         endif
         call kFormShapes (keshp,kemid,2,GPts(1),Shapes(1),
     x      ShapeDerv(1))
         if (kax .eq. 1) then
            if (kmodeI .eq. 0) then
               xx = Shapes(1)*xcurMid(1,1) + Shapes(2)*xcurMid(1,2)
            elseif (kmodeI .eq. 1) then
               xx = Shapes(1)*coords(1,1) + Shapes(2)*coords(1,2)
            endif
         endif
c
         GPdelta(1) = delta_n(1)*Shapes(1) + delta_n(2)*Shapes(2)
         GPdelta(2) = delta_t(1)*Shapes(1) + delta_t(2)*Shapes(2)
c
c        --- Compute the tractions based on the supplied traction law.
         call kComputeTractions (2,props(1),GPdelta(1),jprops(1),
     x      jprops(2),dlength0,Trac(1),Trac_Jacob(1,1),
     x      GPdelcr)
c
         call kFormAMatrix (keshp,kemid,2,Shapes(1),AMat(1,1))
c
         DefGrad = dlength/dlength0
         DetJ = 0.5d0*dlength
c
         if (kax .ne. 1) then
            deltaV = GaussWt*DetJ/DefGrad
         elseif (kax .eq. 1) then
            deltaV = GaussWt*DetJ*(2.d0*dPI*xx)/DefGrad
         endif
c
         if (staticproc .or. dynaproc .or. stiffmat) then
c
c        --- Compute the Gauss Point element stiffness.
c        --- Perform [A]T*[del_T]*[A] operation in element local
c             coordinate system
            call kmaxb (Trac_Jacob(1,1),AMat(1,1),work(1),2,2,2,
     x            2,8,2)
            call kmatxb (AMat(1,1),work(1),stiff1(1,1), 2,2,8, 8,8,2)
            call kvmult1 (stiff1(1,1),64,deltaV)
c
c           --- Transform element stiffness into the global coordinate
c            system
            call kmaxb (stiff1(1,1),TransMat(1,1),work(1),8,8,8,
     x            8,8,8)
            call kmatxb (TransMat(1,1),work(1),stiff2(1,1), 8,8,8,
     x            8,8,8)
c           --- Add the gauss point contribution to the total stiffness
c             matrix.
            do i = 1, 8
               do j = 1, 8
                  amatrx(i,j) = amatrx(i,j) + aconst*stiff2(i,j)
               enddo
            enddo
c
         endif
c
c      --- Compute the Gauss Point residual load vector
         if (.not.stiffmat) then
            Trac(1) = -Trac(1)
            Trac(2) = -Trac(2)
            call kmatxv (AMat(1,1),Trac(1),residual1(1),2,8)
            call kvmult1 (residual1(1),8,deltaV)
c        --- Transform the force vector into global system
            call kmatxv (TransMat(1,1),residual1(1),residual2(1),8,8)
            call kvapb1 (rhs(1),residual2(1),8)
         endif
c
c        -- Save the gp tractions and openings etc...
         if (staticproc .or. dynaproc .or. initialacel) then
            i = 2*(ini-1)
            svars(i+1) = -Trac(1)
            svars(i+2) = -Trac(2)
            svars(i+ipt1) = GPdelta(1)
            svars(i+ipt1+1) = GPdelta(2)
            goto 123
            if (kmodeI .eq. 0) then
               svars(i+ipt2) = xcurMid(1,1)*Shapes(1) +
     x                           xcurMid(1,2)*Shapes(2)
               svars(i+ipt2+1) = xcurMid(2,1)*Shapes(1) +
     x                           xcurMid(2,2)*Shapes(2)
            elseif (kmodeI .eq. 1) then
               svars(i+ipt2) = xcur(1,1)*Shapes(1) +
     x                           xcur(1,2)*Shapes(2)
               svars(i+ipt2+1) = xcur(2,1)*Shapes(1) +
     x                           xcur(2,2)*Shapes(2)
            endif
123         continue
            svars(i+ipt2) = coords(1,1)*Shapes(1) +
     x                        coords(1,2)*Shapes(2)
            svars(i+ipt2+1) = coords(2,1)*Shapes(1) +
     x                           coords(2,2)*Shapes(2)
c            svars(nini*2*3 + ini) = Trac_Jacob(2,2)
         endif
c
      enddo
c
      if (dynaproc .or. initialacel) then
c      if (dynaproc) then
         call kvmove (rhs(1),rhssave(1),8)
      endif
c
      if (dynaproc) then
         call kvmult1(rhs(1),8,aconst)
         call kvmove (svars(ipt3),r8temp1(1),8)
         call kvmult1(r8temp1(1),8,-alpha)
         call kvapb1(rhs(1),r8temp1(1),8)
      endif
c
c     --- Half step residual LV
      if (hsresd) then
         call kvmult1(rhs(1),8,aconst)
         xx = -0.5d0*alpha
         call kvmove(svars(ipt3),r8temp1(1),8)
         call kvmult1(r8temp1(1),8,xx)
         call kvmove(svars(ipt4),r8temp2(1),8)
         call kvmult1(r8temp2(1),8,xx)
         call kvapb1(rhs(1),r8temp1(1),8)
         call kvapb1(rhs(1),r8temp2(1),8)
      endif
c
      if (dynaproc .or. initialacel) then
c      if (dynaproc) then
         call kvmove(svars(ipt3),svars(ipt4),8)
         call kvmove(rhssave(1),svars(ipt3),8)
      endif
c
 999  return
      end
c
c=======================================================================
      subroutine kelem2D6Node (kax,procflags,mattype,ndofel,numnode,
     x njprop,nrhs,nsvars,mlvarx,mcrd,jprops,props,coords,params,
     x svars,u,amatrx,rhs)
c
c --- 2D 6-Noded Cohesive Element----------------------------------
c
      include 'ABA_PARAM.INC'
c
c Function:
c --------
c   Compute element matrices, residual load vectors,
c    --- Arguments
      dimension
     x rhs(*),amatrx(ndofel,ndofel),props(*),
     x svars(*),coords(mcrd,numnode),u(ndofel),
     x jprops(*),params(*)
c
      logical procflags(*),mattype(*)
c
      dimension
     x TransMat(12,12), Shapes(3),
     x xcur(2,6),xcurMid(2,3),
     x work(144), dnxx(2),dnyy(2),dnzz(2),
     x xcurRot(2,6),stiff1(12,12),residual1(12),
     x GPdelta(2),delta_n(3),delta_t(3),AMat(2,12),
     x Trac(2),Trac_Jacob(2,2),stiff2(12,12),residual2(12),
     x ShapeDerv(3),x5cor(5),wt5(5),x4cor(4), wt4(4),
     x GPts(1), rhssave(12)
c
      logical staticproc,dynaproc,hsresd,initialacel,
     x          stiffmat
c
      data x4cor /-0.861136311594953d0,-0.339981043584856d0,
     x             0.339981043584856d0, 0.861136311594953d0/
      data wt4   / 0.347854845137454d0, 0.652145154862546d0,
     x             0.652145154862546d0, 0.347854845137454d0/
c
      data x5cor /-0.906179845938664d0,-0.538469310105683d0,
     x             0.000000000000000d0, 0.538469310105683d0,
     x             0.906179845938664d0/

      data wt5  /  0.236926885056189d0, 0.478628670499366d0,
     x             0.568888888888889d0, 0.478628670499366d0,
     x             0.236926885056189d0/
c
      data dPI /3.1415927d0/
c
c
c                side 2
c
c            4      6     3
c            +------+-----+
c            |            |
c            |            |
c            +------+-----+
c            1      5     2
c
c               side 1
c
c
      keshp = 0
      kemid = 1
      kmodeI = 0
c
      if (jprops(3) .eq. 1) kmodeI = 1
c
c --- Update element nodal coordinates, x = X + dx operation
      do inode = 1, 6
         do i = 1, 2
            xcur(i,inode) = coords(i,inode) + u(2*(inode-1)+i)
         enddo
      enddo
c
c --- Compute the coordinates of the mid surface...
      if (kmodeI .eq. 0) then
         do i = 1, 2
            xcurMid(i,1) = 0.5d0*(xcur(i,1) + xcur(i,4))
            xcurMid(i,2) = 0.5d0*(xcur(i,2) + xcur(i,3))
            xcurMid(i,3) = 0.5d0*(xcur(i,5) + xcur(i,6))
         enddo
      elseif (kmodeI .eq. 1) then
         do i = 1, 2
            xcurMid(i,1) = coords(i,1)
            xcurMid(i,2) = coords(i,2)
            xcurMid(i,3) = coords(i,5)
         enddo
      endif
c
c --- Compute the tangent at the left end of the mid-surface

      dn1s = -1.5d0
      dn2s = -0.5d0
      dn3s = 2.d0
      dnx = xcurMid(1,1)*dn1s + xcurMid(1,2)*dn2s + xcurMid(1,3)*dn3s
      dny = xcurMid(2,1)*dn1s + xcurMid(2,2)*dn2s + xcurMid(2,3)*dn3s
      r8temp = sqrt(dnx*dnx + dny*dny)
      dnx = dnx/r8temp
      dny = dny/r8temp
c
c --- Rotate the current node 1 & node 4 coordinates into the
c     the above computed coordinate system.
c
      inode = 4
      xcurRot(1,inode) = dnx*xcur(1,inode) + dny*xcur(2,inode)
      xcurRot(2,inode) = -dny*xcur(1,inode) + dnx*xcur(2,inode)
      inode = 1
      xcurRot(1,inode) = dnx*xcur(1,inode) + dny*xcur(2,inode)
      xcurRot(2,inode) = -dny*xcur(1,inode) + dnx*xcur(2,inode)
c     -- Get the normal and tangential openings at left end....
      delta_n(1) = xcurRot(2,4) - xcurRot(2,1)
      delta_t(1) = xcurRot(1,4) - xcurRot(1,1)
c
c --- Compute the tangent at the right end of the mid-surface

      dn1s = 0.5d0
      dn2s = 1.5d0
      dn3s = -2.d0
      dnx = xcurMid(1,1)*dn1s + xcurMid(1,2)*dn2s + xcurMid(1,3)*dn3s
      dny = xcurMid(2,1)*dn1s + xcurMid(2,2)*dn2s + xcurMid(2,3)*dn3s
      r8temp = sqrt(dnx*dnx + dny*dny)
      dnx = dnx/r8temp
      dny = dny/r8temp
c
c --- Rotate the current node 1 & node 4 coordinates into the
c     the above computed coordinate system.
c
      inode = 2
      xcurRot(1,inode) = dnx*xcur(1,inode) + dny*xcur(2,inode)
      xcurRot(2,inode) = -dny*xcur(1,inode) + dnx*xcur(2,inode)
      inode = 3
      xcurRot(1,inode) = dnx*xcur(1,inode) + dny*xcur(2,inode)
      xcurRot(2,inode) = -dny*xcur(1,inode) + dnx*xcur(2,inode)
c     -- Get the normal and tangential openings at right end...
      delta_n(2) = xcurRot(2,3) - xcurRot(2,2)
      delta_t(2) = xcurRot(1,3) - xcurRot(1,2)

c
c --- Compute the tangent at the middle of the mid-surface

      dn1s = -0.5d0
      dn2s = 0.5d0
      dn3s = 0.d0
      dnx = xcurMid(1,1)*dn1s + xcurMid(1,2)*dn2s + xcurMid(1,3)*dn3s
      dny = xcurMid(2,1)*dn1s + xcurMid(2,2)*dn2s + xcurMid(2,3)*dn3s
      r8temp = sqrt(dnx*dnx + dny*dny)
      dnx = dnx/r8temp
      dny = dny/r8temp
c
c --- Rotate the current node 5 & node 6 coordinates into the
c     the above computed coordinate system.
c
      inode = 5
      xcurRot(1,inode) = dnx*xcur(1,inode) + dny*xcur(2,inode)
      xcurRot(2,inode) = -dny*xcur(1,inode) + dnx*xcur(2,inode)
      inode = 6
      xcurRot(1,inode) = dnx*xcur(1,inode) + dny*xcur(2,inode)
      xcurRot(2,inode) = -dny*xcur(1,inode) + dnx*xcur(2,inode)
c     -- Get the normal and tangential openings at center
      delta_n(3) = xcurRot(2,6) - xcurRot(2,5)
      delta_t(3) = xcurRot(1,6) - xcurRot(1,5)
c
c --- Perform some initializations
      call kvzero (amatrx(1,1),ndofel*ndofel)
      call kvzero (rhs(1),ndofel)
c
      nini = 4
      ipt1 = 9
      ipt2 = 17
      ipt3 = 25
      ipt4 = 37
      if (kax .eq. 1) then
         nini = 5
         ipt1 = 11
         ipt2 = 21
         ipt3 = 31
         ipt4 = 43
      endif
c
c --- Define procedure based flags...
      alpha = params(1)
      staticproc = procflags(1)
      dynaproc = procflags(2)
      hsresd = procflags(3)
      initialacel = procflags(4)
c
      stiffmat = .false.
      if (mattype(1)) stiffmat = .true.
c
      if (staticproc .or. stiffmat) then
         aconst = 1.d0
      elseif (dynaproc .or. hsresd) then
         aconst = 1.d0 + alpha
      endif
c
c --- Beging numerical integration....
c
      do ini = 1, nini

         if (nini .eq. 4) then
            GPts(1) = x4cor(ini)
            GaussWt = wt4(ini)
         elseif (nini .eq. 5) then
            GPts(1) = x5cor(ini)
            GaussWt = wt5(ini)
         endif

         call kFormShapes (keshp,kemid,3,GPts(1),Shapes(1),
     x         ShapeDerv(1))
c
c        --- Get the normal to side 1 and 2 at the current gauss point
c
          dn1s = ShapeDerv(1)
          dn2s = ShapeDerv(2)
          dn3s = ShapeDerv(3)
          dnxi = xcurMid(1,1)*dn1s + xcurMid(1,2)*dn2s + xcurMid(1,3)*dn3s
          dnyi = xcurMid(2,1)*dn1s + xcurMid(2,2)*dn2s + xcurMid(2,3)*dn3s
          detj = sqrt(dnxi*dnxi + dnyi*dnyi)

          dnxi = dnxi/detj
          dnyi = dnyi/detj
          dnxx(1) = dnxi
          dnyy(1) = dnyi
          dnxx(2) = -dnyi
          dnyy(2) = dnxi
          dnzz(1) = 0.d0
          dnzz(2) = 0.d0

          dnxi = coords(1,1)*dn1s + coords(1,2)*dn2s + coords(1,5)*dn3s
          dnyi = coords(2,1)*dn1s + coords(2,2)*dn2s + coords(2,5)*dn3s
          detj0 = sqrt(dnxi*dnxi + dnyi*dnyi)
c
c    --- Get the Transformation matrix from global to the
c        element normal-tangential coordinate system.
         call kFormTransMatrix (keshp,kemid,dnxx(1),dnyy(1),dnzz(1),
     x         12,TransMat(1,1))
c
         if (kax .eq. 1) then
            if (kmodeI .eq. 0) then
               xx = Shapes(1)*xcurMid(1,1) + Shapes(2)*xcurMid(1,2) +
     x             Shapes(3)*xcurMid(1,3)
            elseif (kmodeI .eq. 1) then
               xx = Shapes(1)*coords(1,1) + Shapes(2)*coords(1,2) +
     x             Shapes(3)*coords(1,3)
            endif
         endif

         GPdelta(1) = delta_n(1)*Shapes(1) + delta_n(2)*Shapes(2) +
     x         delta_n(3)*Shapes(3)
         GPdelta(2) = delta_t(1)*Shapes(1) + delta_t(2)*Shapes(2) +
     x         delta_t(3)*Shapes(3)
c
c           --- Compute the tractions based on the supplied traction law.
c

         call kComputeTractions (2,props(1),GPdelta(1),jprops(1),
     x      jprops(2),detj0,Trac(1),Trac_Jacob(1,1),GPdelcr)
c
         call kFormAMatrix (keshp,kemid,2,Shapes(1),AMat(1,1))
c
         DefGrad = detj/detj0
         if (kax .ne. 1) then
            deltaV = GaussWt*detj/DefGrad
         elseif (kax .eq. 1) then
            deltaV = GaussWt*detj*(2.d0*dPI*xx)/DefGrad
         endif

         if (staticproc .or. dynaproc .or. stiffmat) then
c
c           --- Compute the Gauss Point element stiffness.
c           --- Perform [A]T*[del_T]*[A] operation in element local
c                coordinate system

             call kmaxb (Trac_Jacob(1,1),AMat(1,1),work(1),2,2,2,
     x                2,12,2)
             call kmatxb (AMat(1,1),work(1),stiff1(1,1), 2,2,12,
     x                12,12,2)
             call kvmult1 (stiff1(1,1),144,deltaV)
c
c            --- Transform the element stiffness into the global
c                coordinate system
             call kmaxb (stiff1(1,1),TransMat(1,1),work(1),
     x              12,12,12,12,12,12)
             call kmatxb (TransMat(1,1),work(1),stiff2(1,1),
     x            12,12,12, 12,12,12)
c            --- Add the gauss point contribution to the total stiffness
c                matrix.
             do i = 1, 12
                do j = 1, 12
                   amatrx(i,j) = amatrx(i,j) + aconst*stiff2(i,j)
                enddo
             enddo

         endif
c        --- Compute the Gauss Point element load vector
         if (.not.stiffmat) then
            Trac(1) = -Trac(1)
            Trac(2) = -Trac(2)
            call kmatxv (AMat(1,1),Trac(1),residual1(1),2,12)
            call kvmult1 (residual1(1),12,deltaV)
c          --- Transform the force vector into global system
            call kmatxv (TransMat(1,1),residual1(1),residual2(1),12,12)
            call kvapb1 (rhs(1),residual2(1),12)
         endif
c
         i = 2*(ini-1)
         svars(i+1) = -Trac(1)
         svars(i+2) = -Trac(2)
         svars(i+ipt1) = GPdelta(1)
         svars(i+ipt1+1) = GPdelta(2)
         if (kmodeI .eq. 0) then
            svars(i+ipt2) = xcurMid(1,1)*Shapes(1) +
     x            xcurMid(1,2)*Shapes(2)+xcurMid(1,3)*Shapes(3)
            svars(i+ipt2+1) = xcurMid(2,1)*Shapes(1) +
     x            xcurMid(2,2)*Shapes(2)+xcurMid(2,3)*Shapes(3)
         elseif (kmodeI .eq. 1) then
            svars(i+ipt2) = xcur(1,1)*Shapes(1) +
     x                  xcur(1,2)*Shapes(2)+xcur(1,3)*Shapes(3)
            svars(i+ipt2+1) = xcur(2,1)*Shapes(1) +
     x                  xcur(2,2)*Shapes(2)+xcur(2,3)*Shapes(3)
         endif
c
      enddo
c
      if (dynaproc .or. initialacel) then
         call kvmove (rhs(1),rhssave(1),12)
      endif
c
      if (dynaproc) then
         call kvmult1(rhs(1),12,aconst)
         call kvmult1(svars(ipt3),12,-alpha)
         call kvapb1(rhs(1),svars(ipt3),12)
      endif
c
      if (hsresd) then
         call kvmult1(rhs(1),12,aconst)
         xx = -0.5d0*alpha
         call kvmult1(svars(ipt3),12,xx)
         call kvmult1(svars(ipt4),12,xx)
         call kvapb1(rhs(1),svars(ipt3),12)
         call kvapb1(rhs(1),svars(ipt4),12)
      endif
c
      if (dynaproc .or. initialacel) then
         call kvmove(svars(ipt3),svars(ipt4),12)
         call kvmove(rhssave(1),svars(ipt3),12)
      endif
c
999   return
      end
c====================================================================
      subroutine kelem3D8Node(procflags,mattype,ndofel,numnode,njprop,
     x nrhs,nsvars,mlvarx,mcrd,jprops,props,coords,svars,u,params,
     x amatrx,rhs)
c
c --- 3D-8 Noded Brick Linear Cohesive Element------------------
c
      include 'ABA_PARAM.INC'
c
c Function:
c --------
c   Compute element matrices, residual load vectors, for 3d
c   8-node brick to 8-node brick cohesive element
c
c    --- Arguments
      dimension
     x rhs(*),amatrx(ndofel,ndofel),props(*),
     x svars(*),coords(mcrd,numnode),u(ndofel),
     x jprops(*),params(*)
c
      logical procflags(*),mattype(*)
c
      parameter (intg=3)
      parameter (nrow=24)
      parameter (nnode=8)
      parameter (nglodim=3)
      parameter (ngausspts=intg*intg)
c
      dimension
     x TransMat(nrow,nrow),Shapes(4),xcur(nglodim,nnode),
     x xcurMid(nglodim,4),work(nrow*nrow),dnxt(3),dnyt(3),dnzt(3),
     x dnxb(3),dnyb(3),dnzb(3),
     x xcurRot(nglodim,nnode),stiff1(nrow,nrow),residual1(nrow),
     x GPdelta(nglodim),delta_n(4),delta_t1(4),
     x AMat(nglodim,nrow),Trac(nglodim),Trac_Jacob(nglodim,nglodim),
     x dnx(3),dny(3),dnz(3),delta_t2(4),
     x stiff2(nrow,nrow),residual2(nrow),GP_Trac(nglodim,ngausspts),
     x ShapeDerv(2,4), GaussPts(2,81),GaussWts(81), dNodGPts(2,4),
     x GPts(2),GaussPtst(2,9),GaussWtst(9),coordst(3,4),
     x rhssave(nrow)
c
      logical staticproc,dynaproc,hsresd,initialacel,stiffmat
c
c       data GaussPts /-0.577350269189626d0, -0.577350269189626d0,
c     x  -0.577350269189626d0, 0.577350269189626d0,
c     x  0.577350269189626d0, 0.577350269189626d0,
c     x  0.577350269189626d0, -0.577350269189626d0/
c
c      data GaussWts /1.d0, 1.d0, 1.d0, 1.d0/
c
      data dNodGPts /-1.d0, -1.d0, 1.d0, -1.d0,
     x              1.d0, 1.d0, -1.d0, 1.d0/
c
      keshp = 1
      kemid = 0
      kmodeI = jprops(3)
c
      do i = 1, 4
         call kvmove (coords(1,i),coordst(1,i),3)
      enddo
c
c --- Update element nodal coordinates, x = X + dx operation
      do inode = 1, nnode
         do i = 1, nglodim
            xcur(i,inode) = coords(i,inode) + u(nglodim*(inode-1)+i)
         enddo
      enddo
c
c --- Get the element mid-plane in the defomred configuration
      do i = 1, 4
         do j = 1, 3
            if (kmodeI .eq. 0) then
               xcurMid(j,i) = 0.5d0*(xcur(j,i) + xcur(j,i+4))
            elseif (kmodeI .eq. 1) then
               xcurMid(j,i) = coords(j,i)
            endif
         enddo
      enddo
c
c --- Form the current nodal coordinates in the
c     normal-tangential coordinate system at the element corner
c     nodes.
c
      do inode = 1, 4
         call kFormNodTangents (keshp,kemid,nglodim,dNodGPts(1,inode),
     x      xcurMid(1,1),dnx(1),dny(1),dnz(1),dArea)
c
         ii = inode
         xcurRot(1,ii) = dnx(1)*xcur(1,ii) + dny(1)*xcur(2,ii)+
     x                         dnz(1)*xcur(3,ii)
         xcurRot(2,ii) = dnx(2)*xcur(1,ii) + dny(2)*xcur(2,ii)+
     x                         dnz(2)*xcur(3,ii)
         xcurRot(3,ii) = dnx(3)*xcur(1,ii) + dny(3)*xcur(2,ii)+
     x                         dnz(3)*xcur(3,ii)
         ii = 4 + inode
         xcurRot(1,ii) = dnx(1)*xcur(1,ii) + dny(1)*xcur(2,ii)+
     x                         dnz(1)*xcur(3,ii)
         xcurRot(2,ii) = dnx(2)*xcur(1,ii) + dny(2)*xcur(2,ii)+
     x                         dnz(2)*xcur(3,ii)
         xcurRot(3,ii) = dnx(3)*xcur(1,ii) + dny(3)*xcur(2,ii)+
     x                         dnz(3)*xcur(3,ii)
      enddo
c
c -- Compute the normal openings and tangential slidings
      do inode = 1, 4
         delta_n(inode) = xcurRot(2,inode+4) - xcurRot(2,inode)
         delta_t1(inode) = xcurRot(1,inode+4) - xcurRot(1,inode)
         delta_t2(inode) = xcurRot(3,inode+4) - xcurRot(3,inode)
      enddo
c
c --- Perform some initializations
      call kvzero (amatrx(1,1),ndofel*ndofel)
      call kvzero (rhs(1),ndofel)
c
c --- Iniitlaizations....
      ipt1 = 1 + nglodim*ngausspts
      ipt2 = ipt1 + nglodim*ngausspts
      ipt3 = ipt2 + nglodim*ngausspts
      ipt4 = ipt3 + nrow
      ipt5 = ipt4 + nrow
c
c --- Define procedure based flags...
      alpha = params(1)
      staticproc = procflags(1)
      dynaproc = procflags(2)
      hsresd = procflags(3)
      initialacel = procflags(4)
c
      stiffmat = .false.
      if (mattype(1)) stiffmat = .true.
c
      if (staticproc .or. stiffmat) then
         aconst = 1.d0
      elseif (dynaproc .or. hsresd) then
         aconst = 1.d0 + alpha
      endif
c
c --- Get the gauss points and wts
      call kQuadIntPts (intg,i4temp,GaussPtst(1,1),GaussWtst(1))
      ii = 0
      do i = 1, intg
         do j = 1, intg
            ii = ii + 1
            GaussPts(1,ii) = GaussPtst(1,i)
            GaussPts(2,ii) = GaussPtst(1,j)
            GaussWts(ii) = GaussWtst(i)*GaussWtst(j)
         enddo
      enddo
       
c       --- Begin numerical integration....
c
      do ini = 1, ngausspts

         GPts(1) = GaussPts(1,ini)
         GPts(2) = GaussPts(2,ini)

         call kFormShapes (keshp,kemid,2,GPts(1),Shapes(1),
     x         ShapeDerv(1,1))

         call kFormNodTangents (keshp,kemid,nglodim,GPts(1),
     x         xcurMid(1,1),dnx(1),dny(1),dnz(1),dArea)
c
c        --- Get the Transformation matrix from global to the
c            element normal-tangential coordinate system.

         call kFormTransMatrix (keshp,kemid,dnx(1),dny(1),dnz(1),
     x         nrow,TransMat(1,1))
c
c        --- form the GP undef. area
         call kFormNodTangents (keshp,kemid,nglodim,GPts(1),
     x         coordst(1,1),dnx(1),dny(1),dnz(1),dArea0)
c
         call kvzero (GPdelta(1),3)
         do i = 1, 4
            GPdelta(1) = GPdelta(1) + delta_n(i)*Shapes(i)
            GPdelta(2) = GPdelta(2) + delta_t1(i)*Shapes(i)
            GPdelta(3) = GPdelta(3) + delta_t2(i)*Shapes(i)
         enddo
c
c           --- Compute the tractions based on the supplied traction law.
c

         call kComputeTractions (nglodim,props(1),GPdelta(1),
     x      jprops(1),jprops(2),dArea0,Trac(1),Trac_Jacob(1,1),GPdelcr)

c
         svars(ipt5+ini-1) = GPdelcr
c
         call kFormAMatrix (keshp,kemid,nglodim,Shapes(1),AMat(1,1))
c
         DefGrad = dArea/dArea0
         DetJ = dArea
         deltaV = GaussWts(ini)*DetJ/DefGrad

         if (staticproc .or. dynaproc .or. stiffmat) then
c
c          --- Compute the Gauss Point element stiffness.
c          --- Perform [A]T*[del_T]*[A] operation in element local
c                coordinate system

            call kmaxb (Trac_Jacob(1,1),AMat(1,1),work(1),
     x                  nglodim,nglodim,nglodim,nglodim,nrow,nglodim)
            call kmatxb (AMat(1,1),work(1),stiff1(1,1),
     x                   nglodim,nglodim,nrow,nrow,nrow,nglodim)
            call kvmult1 (stiff1(1,1),nrow*nrow,deltaV)
c
c           --- Transform the element stiffness into the global
c            coordinate system, [T]t*[Ke]*[T]

            call kmaxb (stiff1(1,1),TransMat(1,1),work(1),nrow,
     x                  nrow,nrow,nrow,nrow,nrow)
            call kmatxb (TransMat(1,1),work(1),stiff2(1,1),nrow,
     x                   nrow,nrow,nrow,nrow,nrow)

c           --- Add the gauss point contribution to the total stiffness
c               matrix.
            do i = 1, nrow
               do j = 1, nrow
                  amatrx(i,j) = amatrx(i,j) + aconst*stiff2(i,j)
               enddo
            enddo
c
         endif
c
c        --- Compute the rhs vector...
         if (.not.stiffmat) then
c           --- Compute the Gauss Point residual load vector
            call kvmult1 (Trac(1),nglodim,-1.d0)
            call kmatxv (AMat(1,1),Trac(1),residual1(1),nglodim,nrow)
            call kvmult1 (residual1(1),nrow,deltaV)
c           --- Transform the force vector into global system
            call kmatxv (TransMat(1,1),residual1(1),residual2(1),
     x         nrow,nrow)
            call kvapb1 (rhs(1),residual2(1),nrow)
            call kvmult1 (Trac(1),nglodim,-1.d0)
         endif
c
         i = nglodim*(ini-1)
         call kvmove(Trac(1),svars(i+1),nglodim)
         call kvmove(GPdelta(1),svars(ipt1+i),nglodim)
         if (kmodeI .eq. 0) then
            do j = 1, 3
               svars(i+ipt2+j-1) = 0.d0
               do k = 1, 4
                  svars(i+ipt2+j-1) = svars(i+ipt2+j-1) + 
     x                 xcurMid(j,k)*Shapes(k)
               enddo
            enddo
         elseif (kmodeI .eq. 1) then
            do j = 1, 3
               svars(i+ipt2+j-1) = 0.d0
               do k = 1, 4
                  svars(i+ipt2+j-1) = svars(i+ipt2+j-1) + 
     x                 xcur(j,k)*Shapes(k)
               enddo
            enddo
         endif
c
      enddo
c
      if (dynaproc .or. initialacel) then
         call kvmove(rhs(1),rhssave(1),nrow)
      endif
c
      if (dynaproc) then
         call kvmult1(rhs(1),nrow,aconst)
         call kvmult1(svars(ipt3),nrow,-alpha)
         call kvapb1(rhs(1),svars(ipt3),nrow)
      endif
c
      if (hsresd) then
         call kvmult1(rhs(1),nrow,aconst)
         xx = -0.5d0*alpha
         call kvmult1(svars(ipt3),nrow,xx)
         call kvmult1(svars(ipt4),nrow,xx)
         call kvapb1(rhs(1),svars(ipt3),nrow)
         call kvapb1(rhs(1),svars(ipt4),nrow)
      endif
c
      if (dynaproc .or. initialacel) then
         call kvmove(svars(ipt3),svars(ipt4),nrow)
         call kvmove(rhssave(1),svars(ipt3),nrow)
      endif
c
999   return
      end
c=======================================================================
      subroutine kelem3D12Node (procflags,mattype,ndofel,numnode,
     x njprop,nrhs,nsvars,mlvarx,mcrd,jprops,props,coords,svars,u,
     x params,amatrx,rhs)
c
c --- 3D-12 Noded Tet-to-Tet Quadratic Cohesive Element--------
c
      include 'ABA_PARAM.INC'
c
c Function:
c --------
c   Compute element matrices, residual load vectors,

c    --- Arguments

      dimension
     x rhs(*),amatrx(ndofel,ndofel),props(*),
     x svars(*),coords(mcrd,numnode),u(ndofel),
     x jprops(*),params(*)
c
      logical procflags(*),mattype(*)
c
      parameter (nrow=36)
      parameter (nnode=12)
      parameter (nglodim=3)
      parameter (ngausspts=3)
c
      dimension 
     x TransMat(nrow,nrow), xtop(nglodim,6),
     x xbot(nglodim,6),Shapes(6),xcur(nglodim,nnode), 
     x work(nrow*nrow),dnxt(3),dnyt(3),dnzt(3),dnxb(3),dnyb(3),dnzb(3),
     x xcurRot(nglodim,nnode),stiff1(nrow,nrow),residual1(nrow),
     x GPdelta(nglodim),delta_n(6),delta_t1(6),
     x AMat(nglodim,nrow),Trac(nglodim),Trac_Jacob(nglodim,nglodim),
     x dnx(3),dny(3),dnz(3),delta_t2(6),
     x stiff2(nrow,nrow),residual2(nrow),GP_Trac(nglodim,ngausspts),
     x ShapeDerv(2,6), dNodGPts(2,6),xcurMid(3,6),
     x GPts(2),dl1n3(3),dl2n3(3),wt3(3),coordst(3,6),
     x rhssave(nrow)
c
      logical staticproc,dynaproc,hsresd,initialacel,stiffmat
c
      data dl1n3
     x  /0.666666666666667d0,0.166666666666667d0,0.166666666666667d0/
      data dl2n3
     x  /0.166666666666667d0,0.666666666666667d0,0.166666666666667d0/
      data wt3
     x  /0.333333333333333d0,0.333333333333333d0,0.333333333333333d0/
c
      data dNodGPts /0.d0,0.d0,1.d0,0.d0,0.d0,1.d0,
     x              0.5d0,0.d0,0.5d0,0.5d0,0.d0,0.5d0/
c
      keshp = 2
      kemid = 1
      kmodeI = jprops(3)
c
      do i = 1, 6
         do j = 1, 3
            coordst(j,i) = coords(j,i)
         enddo
      enddo
c
c --- Update element nodal coordinates, x = X + dx operation

      do inode = 1, nnode
         do i = 1, nglodim
            xcur(i,inode) = coords(i,inode) + u(nglodim*(inode-1)+i)
         enddo
      enddo
c
c --- Get the element Mid surface nodal coordinates
      if (kmodeI .eq. 0) then
         do i = 1, 6
            do j = 1, 3
               xcurMid(j,i) = 0.5d0*(xcur(j,i) + xcur(j,i+6))
            enddo
         enddo
      elseif (kmodeI .eq. 1) then
         do i = 1, 6
            do j = 1, 3
               xcurMid(j,i) = coords(j,i) 
            enddo
         enddo
      endif
c
c --- Form the current nodal coordinates in the 
c     normal-tangential coordinate system at the element corner
c     nodes.
c
      do inode = 1, 6

         call kFormNodTangents (keshp,kemid,nglodim,dNodGPts(1,inode),
     x      xcurMid(1,1),dnx(1),dny(1),dnz(1),dArea)

         ii = inode
         xcurRot(1,ii) = dnx(1)*xcur(1,ii) + dny(1)*xcur(2,ii)+
     x                         dnz(1)*xcur(3,ii)
         xcurRot(2,ii) = dnx(2)*xcur(1,ii) + dny(2)*xcur(2,ii)+
     x                         dnz(2)*xcur(3,ii)
         xcurRot(3,ii) = dnx(3)*xcur(1,ii) + dny(3)*xcur(2,ii)+
     x                         dnz(3)*xcur(3,ii)
         ii = 6 + inode
         xcurRot(1,ii) = dnx(1)*xcur(1,ii) + dny(1)*xcur(2,ii)+
     x                         dnz(1)*xcur(3,ii)
         xcurRot(2,ii) = dnx(2)*xcur(1,ii) + dny(2)*xcur(2,ii)+
     x                         dnz(2)*xcur(3,ii)
         xcurRot(3,ii) = dnx(3)*xcur(1,ii) + dny(3)*xcur(2,ii)+
     x                         dnz(3)*xcur(3,ii)
      enddo
c
c -- Compute the normal openings and tangential slidings 

      do inode = 1, 6
         delta_n(inode) = xcurRot(2,inode+6) - xcurRot(2,inode)
         delta_t1(inode) = xcurRot(1,inode+6) - xcurRot(1,inode)
         delta_t2(inode) = xcurRot(3,inode+6) - xcurRot(3,inode)
      enddo

c
c --- materix zero initializations
      call kvzero (amatrx(1,1),ndofel*ndofel)
      call kvzero (rhs(1),ndofel)

c --- pointer Iniitlaizations....
      ipt1 = 1 + nglodim*ngausspts
      ipt2 = ipt1 + nglodim*ngausspts
      ipt3 = ipt2 + nglodim*ngausspts
      ipt4 = ipt3 + nrow
      ipt5 = ipt4 + nrow
c
c --- Define procedure based flags...
      alpha = params(1)
      staticproc = procflags(1)
      dynaproc = procflags(2)
      hsresd = procflags(3)
      initialacel = procflags(4)
c
      stiffmat = .false.
      if (mattype(1)) stiffmat = .true.
c
      if (staticproc .or. stiffmat) then
         aconst = 1.d0
      elseif (dynaproc .or. hsresd) then
         aconst = 1.d0 + alpha
      endif
c
c --- Begin numerical integration....
c 
      do ini = 1, ngausspts

         GPts(1) = dl1n3(ini)
         GPts(2) = dl2n3(ini)

         call kFormShapes (keshp,kemid,2,GPts(1),Shapes(1),
     x         ShapeDerv(1,1))

         call kFormNodTangents (keshp,kemid,nglodim,GPts(1),
     x         xcurMid(1,1),dnx(1),dny(1),dnz(1),dArea)

c        --- Get the Transformation matrix from global to the 
c            element normal-tangential coordinate system.

         call kFormTransMatrix (keshp,kemid,dnx(1),dny(1),dnz(1),
     x         nrow,TransMat(1,1))
c
c        --- form the GP undef. area
         call kFormNodTangents (keshp,kemid,nglodim,GPts(1),
     x         coordst(1,1),dnx(1),dny(1),dnz(1),dArea0)

         call kvzero (GPdelta(1),3)
         do i = 1, 6
            GPdelta(1) = GPdelta(1) + delta_n(i)*Shapes(i) 
            GPdelta(2) = GPdelta(2) + delta_t1(i)*Shapes(i) 
            GPdelta(3) = GPdelta(3) + delta_t2(i)*Shapes(i) 
         enddo
c
c        --- Compute the tractions based on the supplied traction law.
c

         call kComputeTractions (nglodim,props(1),GPdelta(1),
     x      jprops(1),jprops(2),dArea0,Trac(1),Trac_Jacob(1,1),GPdelcr)
c
         DefGrad = dArea/dArea0
         DetJ = dArea
         deltaV = wt3(ini)*DetJ/DefGrad
c
         call kFormAMatrix (keshp,kemid,nglodim,Shapes(1),AMat(1,1))
c
c        --- Compute the Gauss Point element stiffness.

         if (staticproc .or. dynaproc .or. stiffmat) then

c         --- Perform [A]T*[del_T]*[A] operation in element local 
c             coordinate system
             call kmaxb (Trac_Jacob(1,1),AMat(1,1),work(1),
     x                  nglodim,nglodim,nglodim,nglodim,nrow,nglodim)
             call kmatxb (AMat(1,1),work(1),stiff1(1,1), nglodim,
     x                   nglodim,nrow,nrow,nrow,nglodim)
             call kvmult1 (stiff1(1,1),nrow*nrow,deltaV)

c            --- Transform the element stiffness into the global 
c                coordinate system, [T]t*[Ke]*[T]
             call kmaxb (stiff1(1,1),TransMat(1,1),work(1),nrow,
     x               nrow,nrow,nrow,nrow,nrow)
             call kmatxb (TransMat(1,1),work(1),stiff2(1,1),nrow, 
     x                nrow,nrow,nrow,nrow,nrow)

c            --- Add the gauss point contribution to the total stiffness
c                matrix.
             do i = 1, nrow
                do j = 1, nrow
                   amatrx(i,j) = amatrx(i,j) + aconst*stiff2(i,j)
                enddo
             enddo
c
         endif
c
c        --- compute the rhs vector...
         if (.not.stiffmat) then
            call kvmult1 (Trac(1),nglodim,-1.d0)
            call kmatxv (AMat(1,1),Trac(1),residual1(1),nglodim,nrow)
            call kvmult1 (residual1(1),nrow,deltaV)
c           --- Transform the force vector into global system
            call kmatxv (TransMat(1,1),residual1(1),residual2(1),
     x         nrow,nrow)
            call kvapb1 (rhs(1),residual2(1),nrow)
            call kvmult1 (Trac(1),nglodim,-1.d0)
         endif
c
c        --- store history variables...
         i = nglodim*(ini-1)
         call kvmove(Trac(1),svars(i+1),nglodim)
         call kvmove(GPdelta(1),svars(ipt1+i),nglodim)
         if (kmodeI .eq. 0) then
            do j = 1, 3
               svars(i+ipt2+j-1) = 0.d0
               do k = 1, 6
                     svars(i+ipt2+j-1) = svars(i+ipt2+j-1) +
     x              xcurMid(j,k)*Shapes(k)
               enddo
            enddo
         elseif (kmodeI .eq. 1) then
            do j = 1, 3
               svars(i+ipt2+j-1) = 0.d0
               do k = 1, 6
                  svars(i+ipt2+j-1) = svars(i+ipt2+j-1) +
     x              xcur(j,k)*Shapes(k)
               enddo
            enddo
         endif
c
      enddo
c
      if (dynaproc .or. initialacel) then
         call kvmove(rhs(1),rhssave(1),nrow)
      endif
c
      if (dynaproc) then
         call kvmult1(rhs(1),nrow,aconst)
         call kvmult1(svars(ipt3),nrow,-alpha)
         call kvapb1(rhs(1),svars(ipt3),nrow)
      endif
c
      if (hsresd) then
         call kvmult1(rhs(1),nrow,aconst)
         xx = -0.5d0*alpha
         call kvmult1(svars(ipt3),nrow,xx)
         call kvmult1(svars(ipt4),nrow,xx)
         call kvapb1(rhs(1),svars(ipt3),nrow)
         call kvapb1(rhs(1),svars(ipt4),nrow)
      endif
c
      if (dynaproc .or. initialacel) then
         call kvmove(svars(ipt3),svars(ipt4),nrow)
         call kvmove(rhssave(1),svars(ipt3),nrow)
      endif
c
999   return
      end
c======================================================================
      subroutine kFormAMatrix (keshp,kemid,nglodim,Shapes,AMat)
c
      include 'ABA_PARAM.INC'
c
c Inputs:
c ------
c   keshp           - 0-line, 1-quad, 2-tri
c   kemid           - 0-no mid node, 1-mide node
c   nglodim         - Global dimensionality
c   Shapes          - Shape functions.
c
c Outputs:
c -------
c  AMat             - Matrix relation nodal dispalcements to
c                   - the normal/tangential openings.
c
c
      dimension Shapes(*),AMat(nglodim,*)
c
      if (nglodim .eq. 2) then
         if (keshp .eq. 0 .and. kemid .eq. 0) then
            call kvzero (AMat(1,1),16)
         elseif (keshp .eq. 0 .and. kemid .eq. 1) then
            call kvzero (AMat(1,1),24)
         endif
         AMat(1,2) = -Shapes(1)
         AMat(1,4) = -Shapes(2)
         AMat(1,6) = Shapes(2)
         AMat(1,8) = Shapes(1)
         AMat(2,1) = -Shapes(1)
         AMat(2,3) = -Shapes(2)
         AMat(2,5) = Shapes(2)
         AMat(2,7) = Shapes(1)
        if (keshp .eq. 0 .and. kemid .eq. 1) then
            AMat(1,10) = -Shapes(3)
            AMat(1,12) = Shapes(3)
            AMat(2,9) = -Shapes(3)
            AMat(2,11) = Shapes(3)
        endif
      elseif (nglodim .eq. 3) then
         if (keshp .eq. 1 .and. kemid .eq. 0) then
            call kvzero (AMat(1,1),3*24)
            do i = 1, 4
               AMat(2,3*i-2) = -Shapes(i)
               AMat(1,3*i-1) = -Shapes(i)
               AMat(3,3*i) = -Shapes(i)
               AMat(2,12+3*i-2) = Shapes(i)
               AMat(1,12+3*i-1) = Shapes(i)
               AMat(3,12+3*i) = Shapes(i)
            enddo
         elseif (keshp .eq. 2 .and. kemid .eq. 1) then
            call kvzero (AMat(1,1),3*36)
            do i = 1, 6
               AMat(2,3*i-2) = -Shapes(i)
               AMat(1,3*i-1) = -Shapes(i)
               AMat(3,3*i) = -Shapes(i)
               AMat(2,18+3*i-2) = Shapes(i)
               AMat(1,18+3*i-1) = Shapes(i)
               AMat(3,18+3*i) = Shapes(i)
            enddo
         endif
      endif
c
      return
      end
c======================================================================
      subroutine kFormNodTangents (keshp,kemid,nglodim,GPts,xyz,dnx,
     x   dny,dnz,dArea)
c
      include 'ABA_PARAM.INC'
c
c  Inputs:
c  ------
c   keshp          - Element shape key (0-Line, 1-quad, 2-tri)
c   kemid          - Element mid node key (0-no mid nodes,
c                                                     1-yes mid nodes)
c   nglodim         - Dimensionality
c   GPts           - Input gauss points
c   xyz            - Array of current nodal cooridinates
c
c Outputs:
c --------
c   dnx        - Direction x cosines of the x,y,z axes
c   dny        - Direction y cosines of the x,y,z axes
c   dnz        - Direction z cosines of the x,y,z axes
c                                 AT EACH OF THE NODES
c   dArea      - Cross prod of vectors in the plane of
c                of the cohesive element area of intg.
c
      dimension xyz(nglodim,*), dnx(*), dny(*), dnz(*),
     x          GPts(*),Shapes(6),ShapeDerv(2,6)
c
      call kFormShapes (keshp,kemid,2,GPts(1),Shapes(1),ShapeDerv(1,1))
c
      if (nglodim .eq. 3) then
c        --- Linear Quad. face
         if (keshp .eq. 1 .and. kemid .eq. 0) then
               dnxs = xyz(1,1)*ShapeDerv(1,1) +
     x               xyz(1,2)*ShapeDerv(1,2) +
     x               xyz(1,3)*ShapeDerv(1,3) +
     x               xyz(1,4)*ShapeDerv(1,4)
               dnys = xyz(2,1)*ShapeDerv(1,1) +
     x               xyz(2,2)*ShapeDerv(1,2) +
     x               xyz(2,3)*ShapeDerv(1,3) +
     x               xyz(2,4)*ShapeDerv(1,4)
               dnzs = xyz(3,1)*ShapeDerv(1,1) +
     x               xyz(3,2)*ShapeDerv(1,2) +
     x               xyz(3,3)*ShapeDerv(1,3) +
     x               xyz(3,4)*ShapeDerv(1,4)
               dnxt = xyz(1,1)*ShapeDerv(2,1) +
     x               xyz(1,2)*ShapeDerv(2,2) +
     x               xyz(1,3)*ShapeDerv(2,3) +
     x               xyz(1,4)*ShapeDerv(2,4)
               dnyt = xyz(2,1)*ShapeDerv(2,1) +
     x               xyz(2,2)*ShapeDerv(2,2) +
     x               xyz(2,3)*ShapeDerv(2,3) +
     x               xyz(2,4)*ShapeDerv(2,4)
               dnzt = xyz(3,1)*ShapeDerv(2,1) +
     x               xyz(3,2)*ShapeDerv(2,2) +
     x               xyz(3,3)*ShapeDerv(2,3) +
     x               xyz(3,4)*ShapeDerv(2,4)
c        --- Quadratic Tri. face
         elseif (keshp .eq. 2 .and. kemid .eq. 1) then
               dnxs = xyz(1,1)*ShapeDerv(1,1) +
     x               xyz(1,2)*ShapeDerv(1,2) +
     x               xyz(1,3)*ShapeDerv(1,3) +
     x               xyz(1,4)*ShapeDerv(1,4) +
     x               xyz(1,5)*ShapeDerv(1,5) +
     x               xyz(1,6)*ShapeDerv(1,6)
               dnys = xyz(2,1)*ShapeDerv(1,1) +
     x               xyz(2,2)*ShapeDerv(1,2) +
     x               xyz(2,3)*ShapeDerv(1,3) +
     x               xyz(2,4)*ShapeDerv(1,4) +
     x               xyz(2,5)*ShapeDerv(1,5) +
     x               xyz(2,6)*ShapeDerv(1,6)
               dnzs = xyz(3,1)*ShapeDerv(1,1) +
     x               xyz(3,2)*ShapeDerv(1,2) +
     x               xyz(3,3)*ShapeDerv(1,3) +
     x               xyz(3,4)*ShapeDerv(1,4) +
     x               xyz(3,5)*ShapeDerv(1,5) +
     x               xyz(3,6)*ShapeDerv(1,6)
               dnxt = xyz(1,1)*ShapeDerv(2,1) +
     x               xyz(1,2)*ShapeDerv(2,2) +
     x               xyz(1,3)*ShapeDerv(2,3) +
     x               xyz(1,4)*ShapeDerv(2,4) +
     x               xyz(1,5)*ShapeDerv(2,5) +
     x               xyz(1,6)*ShapeDerv(2,6)
               dnyt = xyz(2,1)*ShapeDerv(2,1) +
     x               xyz(2,2)*ShapeDerv(2,2) +
     x               xyz(2,3)*ShapeDerv(2,3) +
     x               xyz(2,4)*ShapeDerv(2,4) +
     x               xyz(2,5)*ShapeDerv(2,5) +
     x               xyz(2,6)*ShapeDerv(2,6)
               dnzt = xyz(3,1)*ShapeDerv(2,1) +
     x               xyz(3,2)*ShapeDerv(2,2) +
     x               xyz(3,3)*ShapeDerv(2,3) +
     x               xyz(3,4)*ShapeDerv(2,4) +
     x               xyz(3,5)*ShapeDerv(2,5) +
     x               xyz(3,6)*ShapeDerv(2,6)
         endif
c
         r8temp = sqrt(dnxs*dnxs + dnys*dnys + dnzs*dnzs)
         dnx(1) = dnxs/r8temp
         dny(1) = dnys/r8temp
         dnz(1) = dnzs/r8temp
         dnxr = dnys*dnzt - dnyt*dnzs
         dnyr = dnxt*dnzs - dnxs*dnzt
         dnzr = dnxs*dnyt - dnxt*dnys
         r8temp = sqrt(dnxr*dnxr + dnyr*dnyr + dnzr*dnzr)
         dArea = r8temp
         dnx(2) = dnxr/r8temp
         dny(2) = dnyr/r8temp
         dnz(2) = dnzr/r8temp
c
         dnxt = dnys*dnzr - dnyr*dnzs
         dnyt = dnxr*dnzs - dnxs*dnzr
         dnzt = dnxs*dnyr - dnxr*dnys
         r8temp = sqrt(dnxt*dnxt + dnyt*dnyt + dnzt*dnzt)
         dnx(3) = dnxt/r8temp
         dny(3) = dnyt/r8temp
         dnz(3) = dnzt/r8temp
      endif
c
      return
      end
c======================================================================
      subroutine kFormShapes (keshp,kemid,nn,GaussPts,Shapes,ShapeDerv)
c
      include 'ABA_PARAM.INC'
c
c Function:
c ---------
c   Shapes Function Calculations.
c
c Inputs:
c -------
c  keshp        - 0-Line, 1-Quad, 2-Tri
c  kemid        - 0-no mid node, 1-mid node
c  nn           - Row dimension of array ShapeDerv
c
c Outputs:
c --------
c   Shapes            - Shape functions at current (ini) gp.
c   ShapeDerv         - Shape func. derivatives...
c
c
      dimension GaussPts(*),Shapes(*),ShapeDerv(nn,*)
c
c --- 1D linear element
      if (keshp .eq. 0. and. kemid .eq. 0) then
         chi = GaussPts(1)
         Shapes(1) = 0.5d0*(1.d0 - chi)
         Shapes(2) = 0.5d0*(1.d0 + chi)
         ShapeDerv(1,1) = -0.5d0
         ShapeDerv(2,1) = 0.5d0
c
c --- 1D Quadratic element
      elseif (keshp .eq. 0 .and. kemid .eq. 1) then
         chi = GaussPts(1)
         Shapes(1) = -0.5d0*chi*(1.d0 - chi)
         Shapes(2) = 0.5d0*chi*(1.d0 + chi)
         Shapes(3) = 1.d0 - chi*chi
c
         ShapeDerv(1,1) = -0.5d0*(1.d0 - 2.d0*chi)
         ShapeDerv(2,1) = 0.5d0*(1.d0 + 2.d0*chi)
         ShapeDerv(3,1) = -2.d0*chi
c
c --- 2D Linear Quad Element
      elseif (keshp .eq. 1 .and. kemid .eq. 0) then
         chi = GaussPts(1)
         eta = GaussPts(2)
         Shapes(1) = 0.25d0*(1-chi)*(1-eta)
         Shapes(2) = 0.25d0*(1+chi)*(1-eta)
         Shapes(3) = 0.25d0*(1+chi)*(1+eta)
         Shapes(4) = 0.25d0*(1-chi)*(1+eta)
         ShapeDerv(1,1) = -0.25d0*(1-eta)
         ShapeDerv(2,1) = -0.25d0*(1-chi)
         ShapeDerv(1,2) = 0.25d0*(1-eta)
         ShapeDerv(2,2) = -0.25d0*(1+chi)
         ShapeDerv(1,3) = 0.25d0*(1+eta)
         ShapeDerv(2,3) = 0.25d0*(1+chi)
         ShapeDerv(1,4) = -0.25d0*(1+eta)
         ShapeDerv(2,4) = 0.25d0*(1-chi)
c
c --- 3D Quadatric Tet. Element
      elseif (keshp .eq. 2 .and. kemid .eq. 1) then
         chi=GaussPts(1)
         eta=GaussPts(2)
         Shapes(1)=2.d0*(1-chi-eta)*(0.5d0-chi-eta)
         Shapes(2)=2.d0*chi*(chi-0.5d0)
         Shapes(3)=2.d0*eta*(eta-0.5d0)
         Shapes(4)=4.d0*chi*(1.d0-chi-eta)
         Shapes(5)=4.d0*chi*eta
         Shapes(6)=4.d0*eta*(1.d0-chi-eta)
         ShapeDerv(1,1)=-2.d0*(0.5d0-chi-eta)-2.d0*(1.d0-chi-eta)
         ShapeDerv(2,1)=-2.d0*(0.5d0-chi-eta)-2.d0*(1.d0-chi-eta)
         ShapeDerv(1,2)=2.d0*(2.d0*chi-0.5d0)
         ShapeDerv(2,2)=0.d0
         ShapeDerv(1,3)=0.d0
         ShapeDerv(2,3)=2.d0*(2.d0*eta-0.5d0)
         ShapeDerv(1,4)=4.d0*(1.d0-chi-eta)-4.d0*chi
         ShapeDerv(2,4)=-4.d0*chi
         ShapeDerv(1,5)=4.d0*eta
         ShapeDerv(2,5)=4.d0*chi
         ShapeDerv(1,6)=-4.d0*eta
         ShapeDerv(2,6)=4.d0*(1.d0-chi-eta)-4.d0*eta
c
      endif
c
      return
      end
c======================================================================
      subroutine kFormTransMatrix (keshp,kemid,dnx,dny,dnz,nrow,
     x TransMat)
c
      include 'ABA_PARAM.INC'
c
c Inputs:
c -------
c     keshp           - 0-Line,1 1-Quad, 2-Tri
c     kemid           - 0-No mide node, 1-Mid node
c     dnx, dny, dnz   - direction cosines of the element
c                       tangential-normal coordinate system.
c
c Outputs:
c --------
c     TransMat        - Transformation matrix.
c
      dimension dnx(*),dny(*),dnz(*),TransMat(nrow,*)
c
      call kvzero (TransMat(1,1),nrow*nrow)
c
      if (keshp .eq. 0 .and. kemid .eq. 0) then
         ndim = 2
         nnode = 4
      elseif (keshp .eq. 0 .and. kemid .eq. 1) then
         ndim = 2
         nnode = 6
      elseif (keshp .eq. 1 .and. kemid .eq. 0) then
        ndim = 3
        nnode = 8
      elseif (keshp .eq. 2 .and. kemid .eq. 1) then
        ndim = 3
        nnode = 12
      endif
c
      do inode = 1, nnode
         i2temp = ndim*(inode-1)
         do i = 1, ndim
           TransMat(i2temp+i,i2temp+1) = dnx(i)
           TransMat(i2temp+i,i2temp+2) = dny(i)
           if (ndim .eq. 3) then
              TransMat(i2temp+i,i2temp+3) = dnz(i)
           endif
         enddo
      enddo
c
      return
      end
c================================================================
      subroutine kComputeTractions (nglodim,props,GPdelta,
     x kcomp,ModelFlag,dlength0,Trac,Trac_Jacob,GPdelcr)
c
      include 'ABA_PARAM.INC'
c
c
c Function:
c ---------
c  This Rotuine gives the Gauss point tractions and the
c  Traction Material Jacobian, ie., the Tensor that relates the
c  changes in tractions to the changes in the normal and
c  tangential relative displacements.
c
c Inputs:
c -------
c   nglodim              - Dimensionality of the the prob.
c   GPdelta              - The Gauss point normal and
c                          tangential displacements
c
c   kcomp                - 1 - Add penalty stiffness in compression
c                          0 - no penalty augmentation
c
c   ModelFlag            - 1 = Xu & Needleman 
c                          2 = Exponential potential
c                          3 = L-J Potential
c
c   dlength0             - Original length/area of cohesive 
c                          element
c
c Outputs:
c -------
c   Trac                 - The output normal and tangential
c                          traction Values
c
c   Trac_Jacob
c                        - The material tensor relating the
c                          normal and tangential change in
c                          tractions to the changes in deltas.
c  GPDelcr               - Gauss point delta critical 
c
c
      dimension
     x   GPdelta(nglodim), Trac(nglodim),
     x   Trac_Jacob(nglodim,nglodim),
     x   props(*)
c
      parameter(del_max=15.d0)
      parameter(del_maxt=50.d0)
c
      youngsmod = 1.0e8
c
      if (nglodim .eq. 3) goto 500

      if (ModelFlag .eq. 1) goto 100
      if (ModelFlag .eq. 2) goto 200
      if (ModelFlag .eq. 3) goto 300
      if (ModelFlag .eq. 4) goto 4000
c
c     --- 2D Xu & Needleman potential....
 100  continue

      phi_n = props(1)
      phi_t = props(2)
      del_n = props(3)
      del_t = props(4)
      r = props(5)
      q = phi_t/phi_n
c
      GPdelcr = del_n
c
      delta_n = GPdelta(1)
      delta_t = GPdelta(2)
c
      if ((delta_n/del_n) .gt. del_max) delta_n=del_max*del_n
      if ((delta_n/del_n) .lt. -del_max) delta_n=-del_max*del_n
c
      if ((delta_t/del_t) .gt. sqrt(del_maxt)) 
     x    delta_t=sqrt(del_maxt)*del_t
      if ((delta_t/del_t) .lt. -sqrt(del_maxt)) 
     x    delta_t=-sqrt(del_maxt)*del_t
c
c      if (abs(delta_n) .le. 1.0e-6*del_n) delta_n = 1.0e-6*del_n
c      if (abs(delta_t) .le. 1.0e-6*del_t) delta_t = 1.0e-6*del_t
c
c     -- Normal traction
      Trac(1) =
     x -phi_n/del_n*exp(-delta_n/del_n)*((1-r+delta_n/del_n)*(1-q)/(
     #r-1)-(q+(r-q)/(r-1)*delta_n/del_n)*exp(-delta_t**2/del_t**2))+phi_
     #n*exp(-delta_n/del_n)*(1/del_n*(1-q)/(r-1)-(r-q)/(r-1)/del_n*exp(-
     #delta_t**2/del_t**2))
c
c     --- Tangential traction
      Trac(2) =
     x 2*phi_n*exp(-delta_n/del_n)*(q+(r-q)/(r-1)*delta_n/del_n)*del
     #ta_t/del_t**2*exp(-delta_t**2/del_t**2)
c
c     --- Normal stiffness
      Trac_Jacob(1,1) =
     x phi_n/del_n**2*exp(-delta_n/del_n)*((1-r+delta_n/del_n)*(1-q)
     #/(r-1)-(q+(r-q)/(r-1)*delta_n/del_n)*exp(-delta_t**2/del_t**2))-2*
     #phi_n/del_n*exp(-delta_n/del_n)*(1/del_n*(1-q)/(r-1)-(r-q)/(r-1)/d
     #el_n*exp(-delta_t**2/del_t**2))
c
      Trac_Jacob(1,2) =
     x -2*phi_n/del_n*exp(-delta_n/del_n)*(q+(r-q)/(r-1)*delta_n/del
     #_n)*delta_t/del_t**2*exp(-delta_t**2/del_t**2)+2*phi_n*exp(-delta_
     #n/del_n)*(r-q)/(r-1)/del_n*delta_t/del_t**2*exp(-delta_t**2/del_t*
     #*2)

      Trac_Jacob(2,1) = Trac_Jacob(1,2)

c     --- Tangential stiffness
      Trac_Jacob(2,2) =
     x 2*phi_n*exp(-delta_n/del_n)*(q+(r-q)/(r-1)*delta_n/del_n)/del
     #_t**2*exp(-delta_t**2/del_t**2)-4*phi_n*exp(-delta_n/del_n)*(q+(r-
     #q)/(r-1)*delta_n/del_n)*delta_t**2/del_t**4*exp(-delta_t**2/del_t*
     #*2)
c
      goto 400
c
c     --- 2D Exponential potential 
200   continue
      phi_n = props(1)
      phi_t = props(2)
      del_n = props(3)
      del_t = props(4)
c
      alpha = 1.d0
c
      sigmax = 16.d0*phi_n/(9.d0*del_cr)
      del_cr = del_n
c
      GPdelcr = del_cr
c
      delta_n = GPdelta(1)
      delta_t = GPdelta(2)
c
      if ((delta_n/del_n) .gt. del_max) delta_n=del_max*del_n
      if ((delta_n/del_n) .lt. -del_max) delta_n=-del_max*del_n
      if ((delta_t/del_t) .gt. del_max) delta_t=del_max*del_t
      if ((delta_t/del_t) .lt. -del_max) delta_t=-del_max*del_t
c
      z = 16.d0*exp(1.d0)/9.d0
c
      Trac(1) =
     x 9.E0/16.E0*sigmax*del_cr*(-z/del_cr*exp(-z*delta_n/del_cr)+(1
     #+z*delta_n/del_cr-alpha*z**2*delta_t**2/del_cr**2/2)*z/del_cr*exp(
     #-z*delta_n/del_cr))

      Trac(2) =
     x 9.E0/16.E0*sigmax/del_cr*alpha*z**2*delta_t*exp(-z*delta_n/de
     #l_cr)

      Trac_Jacob(1,1) =
     x 9.E0/16.E0*sigmax*del_cr*(2*z**2/del_cr**2*exp(-z*delta_n/del
     #_cr)-(1+z*delta_n/del_cr-alpha*z**2*delta_t**2/del_cr**2/2)*z**2/d
     #el_cr**2*exp(-z*delta_n/del_cr))

      Trac_Jacob(1,2) =
     x -9.E0/16.E0*sigmax/del_cr**2*alpha*z**3*delta_t*exp(-z*delta_
     #n/del_cr)

      Trac_Jacob(2,1) = Trac_Jacob(1,2)

      Trac_Jacob(2,2) =
     x 9.E0/16.E0*sigmax/del_cr*alpha*z**2*exp(-z*delta_n/del_cr)
c
      goto 400
c
c     --- 2D L-J Potential...
 300  continue

      phi_n = props(1)
      del_n = props(2)
      del_t = del_n
      a = props(3)

      b = (27.d0/256.d0)*(a**4/phi_n**3)
c
      del_cr = (12.d0**(1.d0/6.d0) - 4.d0**(1.d0/6.d0))*
     x               del_n*((b/a)**(1.0/6.d0))
c
      GPdelcr = del_cr
      sigmax = (4.d0/(3.d0*sqrt(12.d0)))
      sigmax = sigmax*(a**1.5)/sqrt(b)
      sigmax = sigmax/del_n
c
      delta_star = del_n*((4.d0*b/a)**(1.0/6.d0))

      delta_n = GPdelta(1)
      delta_t = GPdelta(2)
      delta_n = delta_star + delta_n
c
      if ((delta_t/del_t) .gt. del_max) delta_t=del_max*del_t
      if ((delta_t/del_t) .lt. -del_max) delta_t=-del_max*del_t
c
      Trac(1) =
     x (2*a/delta_n**3*del_n**2-27.E0/32.E0*a**4/phi_n**3/delta_n**9
     x*del_n**8)*exp(-delta_t**2/del_t**2)

      Trac(2) =
     x -2*(-a/delta_n**2*del_n**2+27.E0/256.E0*a**4/phi_n**3/delta_n
     x**8*del_n**8)*delta_t/del_t**2*exp(-delta_t**2/del_t**2)

      Trac_Jacob(1,1) =
     x (-6*a/delta_n**4*del_n**2+243.E0/32.E0*a**4/phi_n**3/delta_n*
     x*10*del_n**8)*exp(-delta_t**2/del_t**2)

      Trac_Jacob(1,2) =
     x -2*(2*a/delta_n**3*del_n**2-27.E0/32.E0*a**4/phi_n**3/delta_n
     x**9*del_n**8)*delta_t/del_t**2*exp(-delta_t**2/del_t**2)

      Trac_Jacob(2,1) = Trac_Jacob(1,2)

      Trac_Jacob(2,2) =
     x -2*(-a/delta_n**2*del_n**2+27.E0/256.E0*a**4/phi_n**3/delta_n
     x**8*del_n**8)/del_t**2*exp(-delta_t**2/del_t**2)+4*(-a/delta_n**2*
     xdel_n**2+27.E0/256.E0*a**4/phi_n**3/delta_n**8*del_n**8)*delta_t**
     x2/del_t**4*exp(-delta_t**2/del_t**2)

      delta_n = delta_n - delta_star
c
      goto 400
c
c ---- Dugdale like model
4000  continue
c
      phi_n = props(1)
      del_cr = props(2)
      f = props(3)
      f1 = props(4)
c
      del1 = f*del_cr
      del2 = f1*del_cr
c
      sigmax = 2.d0*phi_n/(del_cr + del2 - del1) 
c
      delta_n = GPdelta(1)
      delta_t = GPdelta(2)
c
      if (delta_n .lt. del1 .and. delta_n .gt. -del_cr) then
         if (delta_n .ge. 0.d0) Trac(1) = sigmax*(delta_n/del1)
         if (delta_n .lt. 0.d0) then 
            ss = sigmax/del1
            Trac(1) = -ss*delta_n
         endif
         Trac_Jacob(1,1) = sigmax/del1
      elseif (delta_n .ge. del1 .and. delta_n .le. del2) then
         Trac(1) = sigmax
         Trac_Jacob(1,1) = 0.d0
      elseif (delta_n .gt. del2 .and. delta_n .le. del_cr) then
         Trac(1) = sigmax*( (del_cr - delta_n)/(del_cr - del2) )
         Trac_Jacob(1,1) = -sigmax/(del_cr - del2)
      elseif (delta_n .gt. del_cr) then
         Trac(1) = 0.d0
         Trac_Jacob(1,1) = 0.d0
      endif
c
      if (delta_t .lt. 0.d0) goto 4001
      if (delta_t .lt. del1 .and. delta_t .ge. 0.d0) then
         Trac(2) = sigmax*(delta_t/del1)
         Trac_Jacob(2,2) = sigmax/del1
      elseif (delta_t .ge. del1 .and. delta_t .le. del2) then
         Trac(2) = sigmax
         Trac_Jacob(2,2) = 0.d0
      elseif (delta_t .gt. del2 .and. delta_t .le. del_cr) then
         Trac(2) = sigmax*( (del_cr - delta_t)/(del_cr - del2) )
         Trac_Jacob(2,2) = -sigmax/(del_cr - del2)
      elseif (delta_t .gt. del_cr) then
         Trac(2) = 0.d0
         Trac_Jacob(2,2) = 0.d0
      endif
c
4001  continue
      if (delta_t .gt. -del1 .and. delta_t .lt. 0.d0) then
         dt = -delta_t
         Trac(2) = -sigmax*(dt/del1)
         Trac_Jacob(2,2) = sigmax/del1
      elseif (delta_t .ge. -del2 .and. delta_t .le. -del1) then
         Trac(2) = -sigmax
         Trac_Jacob(2,2) = 0.d0
      elseif (delta_t .ge. -del_cr .and. delta_t .lt. -del2) then
         dt = -delta_t
         Trac(2) = -sigmax*( (del_cr - dt)/(del_cr - del2) )
         Trac_Jacob(2,2) = -sigmax/(del_cr - del2)
      elseif (delta_t .lt. -del_cr) then
         Trac(2) = 0.d0
         Trac_Jacob(2,2) = 0.d0
      endif
c
      Trac_Jacob(1,2) = 0.d0
      Trac_Jacob(2,1) = 0.d0
c
c     --- Penalty stiffness augmentation in compression
 400  continue
c
c     --- Add the penalty stiffness in compression if required
c         and requested
      if (delta_n .lt. 0.d0 .and. kcomp .eq. 1) then
         const = youngsmod/dlength0**2
         Trac(1) = Trac(1) - const*3.d0*delta_n**2
         Trac_Jacob(1,1) = Trac_Jacob(1,1) - const*6.d0*delta_n
      endif
c
c     --- Jump to end of subroutine to exit 
      goto 999

c
c     ---- 3D Model computations....
 500  continue
c
      if (ModelFlag .eq. 1) goto 600
      if (ModelFlag .eq. 2) goto 999
      if (ModelFlag .eq. 3) goto 700
c
c     --- Xu and Needleman 3D potential...
c
600   continue

      phi_n = props(1)

      del_n = props(2)
      del_t1 = props(3)
      del_t2 = props(4)

      delta_n = GPdelta(1)
      delta_t1 = GPdelta(2)
      delta_t2 = GPdelta(3)

      if ((delta_n/del_n) .gt. del_max) delta_n=del_max*del_n
      if ((delta_n/del_n) .lt. -del_max) delta_n=-del_max*del_n
      if ((delta_t1/del_t1) .gt. del_max) delta_t=del_max*del_t1
      if ((delta_t1/del_t1) .lt. -del_max) delta_t=-del_max*del_t1
      if ((delta_t2/del_t2) .gt. del_max) delta_t=del_max*del_t2
      if ((delta_t2/del_t2) .lt. -del_max) delta_t=-del_max*del_t2
c
      Trac(1) =
     x phi_n/del_n*exp(-delta_n/del_n)*(1+delta_n/del_n)*exp(-delta_
     #t1**2/del_t1**2-delta_t2**2/del_t2**2)-phi_n*exp(-delta_n/del_n)/d
     #el_n*exp(-delta_t1**2/del_t1**2-delta_t2**2/del_t2**2)

      Trac(2) =
     x 2*phi_n*exp(-delta_n/del_n)*(1+delta_n/del_n)*delta_t1/del_t1
     #**2*exp(-delta_t1**2/del_t1**2-delta_t2**2/del_t2**2)

      Trac(3) =
     x 2*phi_n*exp(-delta_n/del_n)*(1+delta_n/del_n)*delta_t2/del_t2
     #**2*exp(-delta_t1**2/del_t1**2-delta_t2**2/del_t2**2)

      Trac_Jacob(1,1) =
     x -phi_n/del_n**2*exp(-delta_n/del_n)*(1+delta_n/del_n)*exp(-de
     #lta_t1**2/del_t1**2-delta_t2**2/del_t2**2)+2*phi_n/del_n**2*exp(-d
     #elta_n/del_n)*exp(-delta_t1**2/del_t1**2-delta_t2**2/del_t2**2)

      Trac_Jacob(1,2) =
     x -2*phi_n/del_n*exp(-delta_n/del_n)*(1+delta_n/del_n)*delta_t1
     #/del_t1**2*exp(-delta_t1**2/del_t1**2-delta_t2**2/del_t2**2)+2*phi
     #_n*exp(-delta_n/del_n)/del_n*delta_t1/del_t1**2*exp(-delta_t1**2/d
     #el_t1**2-delta_t2**2/del_t2**2)

      Trac_Jacob(1,3) =
     x -2*phi_n/del_n*exp(-delta_n/del_n)*(1+delta_n/del_n)*delta_t2
     #/del_t2**2*exp(-delta_t1**2/del_t1**2-delta_t2**2/del_t2**2)+2*phi
     #_n*exp(-delta_n/del_n)/del_n*delta_t2/del_t2**2*exp(-delta_t1**2/d
     #el_t1**2-delta_t2**2/del_t2**2)

      Trac_Jacob(2,2) =
     x 2*phi_n*exp(-delta_n/del_n)*(1+delta_n/del_n)/del_t1**2*exp(-
     #delta_t1**2/del_t1**2-delta_t2**2/del_t2**2)-4*phi_n*exp(-delta_n/
     #del_n)*(1+delta_n/del_n)*delta_t1**2/del_t1**4*exp(-delta_t1**2/de
     #l_t1**2-delta_t2**2/del_t2**2)

      Trac_Jacob(2,3) =
     x -4*phi_n*exp(-delta_n/del_n)*(1+delta_n/del_n)*delta_t1/del_t
     #1**2*delta_t2/del_t2**2*exp(-delta_t1**2/del_t1**2-delta_t2**2/del
     #_t2**2)

      Trac_Jacob(3,3) =
     x 2*phi_n*exp(-delta_n/del_n)*(1+delta_n/del_n)/del_t2**2*exp(-
     #delta_t1**2/del_t1**2-delta_t2**2/del_t2**2)-4*phi_n*exp(-delta_n/
     #del_n)*(1+delta_n/del_n)*delta_t2**2/del_t2**4*exp(-delta_t1**2/de
     #l_t1**2-delta_t2**2/del_t2**2)
c
      Trac_Jacob(2,1) = Trac_Jacob(1,2)
      Trac_Jacob(3,1) = Trac_Jacob(1,3)
      Trac_Jacob(3,2) = Trac_Jacob(2,3)
c
      goto 900
c
c     --- 3D L-J Potential...
 700  continue

      phi_n = props(1)
      del_n = props(2)
      del_t1 = del_n
      del_t2 = del_n
      a = props(3)

      b = (27.d0/256.d0)*(a**4/phi_n**3)
c
      del_cr = (12.d0**(1.d0/6.d0) - 4.d0**(1.d0/6.d0))*
     x               del_n*((b/a)**(1.0/6.d0))
c
      GPdelcr = del_cr
      sigmax = (4.d0/(3.d0*sqrt(12.d0)))
      sigmax = sigmax*(a**1.5)/sqrt(b)
      sigmax = sigmax/del_n
c
      delta_star = del_n*((4.d0*b/a)**(1.0/6.d0))

      delta_n = GPdelta(1)
      delta_t1 = GPdelta(2)
      delta_t2 = GPdelta(3)

      delta_n = delta_star + delta_n

      if ((delta_t1/del_t1) .gt. del_max) delta_t=del_max*del_t1
      if ((delta_t1/del_t1) .lt. -del_max) delta_t=-del_max*del_t1
      if ((delta_t2/del_t2) .gt. del_max) delta_t=del_max*del_t2
      if ((delta_t2/del_t2) .lt. -del_max) delta_t=-del_max*del_t2
c

      Trac(1) = 
     x (2*a/delta_n**3*del_n**2-8*b/delta_n**9*del_n**8)*exp(-delta_
     xt1**2/del_t1**2-delta_t2**2/del_t2**2)

      Trac(2) = 
     x -2*(-a/delta_n**2*del_n**2+b/delta_n**8*del_n**8)*delta_t1/de
     xl_t1**2*exp(-delta_t1**2/del_t1**2-delta_t2**2/del_t2**2)

      Trac(3) = 
     x -2*(-a/delta_n**2*del_n**2+b/delta_n**8*del_n**8)*delta_t2/de
     xl_t2**2*exp(-delta_t1**2/del_t1**2-delta_t2**2/del_t2**2)

      Trac_Jacob(1,1) = 
     x (-6*a/delta_n**4*del_n**2+72*b/delta_n**10*del_n**8)*exp(-del
     xta_t1**2/del_t1**2-delta_t2**2/del_t2**2)

      Trac_Jacob(1,2) = 
     x -2*(2*a/delta_n**3*del_n**2-8*b/delta_n**9*del_n**8)*delta_t1
     x/del_t1**2*exp(-delta_t1**2/del_t1**2-delta_t2**2/del_t2**2)

      Trac_Jacob(1,3) = 
     x -2*(2*a/delta_n**3*del_n**2-8*b/delta_n**9*del_n**8)*delta_t2
     x/del_t2**2*exp(-delta_t1**2/del_t1**2-delta_t2**2/del_t2**2)

      Trac_Jacob(2,2) = 
     x -2*(-a/delta_n**2*del_n**2+b/delta_n**8*del_n**8)/del_t1**2*e
     xxp(-delta_t1**2/del_t1**2-delta_t2**2/del_t2**2)+4*(-a/delta_n**2*
     xdel_n**2+b/delta_n**8*del_n**8)*delta_t1**2/del_t1**4*exp(-delta_t
     x1**2/del_t1**2-delta_t2**2/del_t2**2)

      Trac_Jacob(2,3) = 
     x 4*(-a/delta_n**2*del_n**2+b/delta_n**8*del_n**8)*delta_t1/del
     x_t1**2*delta_t2/del_t2**2*exp(-delta_t1**2/del_t1**2-delta_t2**2/d
     xel_t2**2)

      Trac_Jacob(3,3) = 
     x -2*(-a/delta_n**2*del_n**2+b/delta_n**8*del_n**8)/del_t2**2*e
     xxp(-delta_t1**2/del_t1**2-delta_t2**2/del_t2**2)+4*(-a/delta_n**2*
     xdel_n**2+b/delta_n**8*del_n**8)*delta_t2**2/del_t2**4*exp(-delta_t
     x1**2/del_t1**2-delta_t2**2/del_t2**2)
c
      Trac_Jacob(2,1) = Trac_Jacob(1,2)
      Trac_Jacob(3,1) = Trac_Jacob(1,3)
      Trac_Jacob(3,2) = Trac_Jacob(2,3)
c
      delta_n = delta_n - delta_star  
c
c     --- Penalty stiffness augmentation in compression
 900  continue
c
c     --- Add the penalty stiffness in compression if required
c         and requested
      if (delta_n .lt. 0.d0 .and. kcomp .eq. 1) then
         const = youngsmod/dlength0**2
         Trac(1) = Trac(1) - const*3.d0*delta_n**2
         Trac_Jacob(1,1) = Trac_Jacob(1,1) - const*6.d0*delta_n
      endif
c
 999  continue
c
      return
      end
c================================================================
      subroutine kvzero (vec,n)
c
c Funtion:
c -------
c  Zero out an input double precision vector, vec(*)
c
      include 'ABA_PARAM.INC'
c
      dimension vec(*)
c
      do i = 1, n
         vec(i)=0.d0
      enddo
c
      return
      end
c==================================================================
      subroutine kmaxb (amat,bmat,cmat, nrowa,nrowb,nrowc, n1,n2,n3)
c
      include 'ABA_PARAM.INC'
c
c Function:
c --------
c    Perform matrix multiply of the form:
c
c              [C] = [A]*[B]
c
c     n1       - number of rows in matrix c to fill
c     n2       - number of columns in matrix c to fill
c     n3       - number of columns in matrix a and
c                number of rows in matrix b
c
      dimension amat(nrowa,*),bmat(nrowb,*),cmat(nrowc,*)
c
      do i = 1, n1
         do j = 1, n2
            cmat(i,j) = 0.d0
            do k = 1, n3
               cmat(i,j) =   cmat(i,j) + amat(i,k)*bmat(k,j)
            enddo
         enddo
      enddo
c
      return
      end
c====================================================================
      subroutine kmatxb (amat,bmat,cmat, nrowa,nrowb,nrowc, n1,n2,n3)
c
      include 'ABA_PARAM.INC'
c
c Function:
c --------
c    Perform matrix multiply of the form:
c
c                       t
c              [C] = [A] * [B]
c
c     n1       - number of rows in matrix c to fill
c     n2       - number of columns in matrix c to fill
c     n3       - number of columns in matrix a and
c                number of rows in matrix b
c
      dimension amat(nrowa,*),bmat(nrowb,*),cmat(nrowc,*)
c
      do i = 1, n1
         do j = 1, n2
            cmat(i,j) = 0.d0
            do k = 1, n3
               cmat(i,j) = cmat(i,j) + amat(k,i)*bmat(k,j)
            enddo
         enddo
      enddo
c
      return
      end
c=====================================================================
      subroutine kvmult1 (vec,n,const)
c
      include 'ABA_PARAM.INC'
c
c Function:
c --------
c   Multiply a vector by a constant and return
c
      dimension vec(*)
c
      do i = 1, n
         vec(i) = vec(i)*const
      enddo
c
      return
      end
c======================================================================
      subroutine kmatxv (amat,vecin,vecout,nrowa,n1)
c
      include 'ABA_PARAM.INC'
c
c Function:
c --------
c   To form the product of a matrix transpose with vector
c
c                      t
c        {vecout} = [A]  * {vecin}
c
      dimension amat(nrowa,*),vecout(*),vecin(*)
c
      do i = 1, n1
         vecout(i) = 0.d0
         do j = 1, nrowa
            vecout(i) = vecout(i) + amat(j,i)*vecin(j)
         enddo
      enddo
c
      return
      end
c======================================================================
      subroutine kvapb1 (vecin1,vecin2,n)
c
      include 'ABA_PARAM.INC'
c
c Function:
c --------
c    To add a vector (vecin2) to another vector (vecin1)
c    and store in the same vector (vecin1)
c
c       {vecin1} = {vecin1} + {vecin2}
c
      dimension  vecin1(*),vecin2(*)
c
      do i = 1, n
         vecin1(i) = vecin1(i) + vecin2(i)
      enddo
c
      return
      end
c
c=================================================================
      subroutine kvmove (vecin1,vecin2,n)
c
      include 'ABA_PARAM.INC'
c
c Function:
c --------
c    To move a vector (vecin1) to another vector (vecin2)
c
c       {vecin2} = {vecin1}
c
      dimension  vecin1(*),vecin2(*)
c
      do i = 1, n
         vecin2(i) = vecin1(i)
      enddo
c
      return
      end
c
c======================================================
      subroutine kvimove (v1,inc1,v2,inc2,n)
c
      include 'ABA_PARAM.INC'
c
c
c     v1        - vector v1
c     inc1      - increment on vector v1
c     inc2      - increment on vector v2
c     n         - number of items to be moved
c
c     v2        - vector v2
c
      dimension  v1(*),v2(*)
c
      i1 = 1
      i2 = 1
      do  10  i = 1,n
         v2(i2) = v1(i1)
         i1 = i1 + inc1
         i2 = i2 + inc2
  10  continue
c
      return
      end
c=============================================================
      subroutine kQuadIntPts (intg,nini,GaussPts,GaussWts)
c
c       intg          - integration rule
c       nini          - number of integration pts 
c       GaussPts      - integration locations
c       GaussWts      - integration weights

      include 'ABA_PARAM.INC'
c                                                        
      dimension GaussPts(2,9),GaussWts(9)
c                                
      dimension 
     x x2cor(2),wt2(2),
     x x3cor(3),wt3(3),
     x x4cor(4),wt4(4),
     x x5cor(5),wt5(5),
     x x6cor(6),wt6(6),
     x x7cor(7),wt7(7),
     x x8cor(8),wt8(8),
     x x9cor(9),wt9(9)
                      
      data x2cor /-0.577350269189626d0, 0.577350269189626d0/
                                                            
      data wt2   / 1.000000000000000d0, 1.000000000000000d0/
                                                            
      data x3cor /-0.774596669241483d0, 0.000000000000000d0,
     x             0.774596669241483d0/
                                       
      data wt3   / 0.555555555555556d0, 0.888888888888889d0,
     x             0.555555555555556d0/
                                       
      data x4cor /-0.861136311594953d0,-0.339981043584856d0,
     x             0.339981043584856d0, 0.861136311594953d0/
                                                            
      data wt4   / 0.347854845137454d0, 0.652145154862546d0,
     x             0.652145154862546d0, 0.347854845137454d0/
                                                            
      data x5cor /-0.906179845938664d0,-0.538469310105683d0,
     x             0.000000000000000d0, 0.538469310105683d0,
     x             0.906179845938664d0/
                                       
      data wt5  /  0.236926885056189d0, 0.478628670499366d0,
     x             0.568888888888889d0, 0.478628670499366d0,
     x             0.236926885056189d0/
                                       
      data x6cor /-0.932469514203152d0,-0.661209386466265d0,
     x            -0.238619186083197d0, 0.238619186083197d0,
     x             0.661209386466265d0, 0.932469514203152d0/
                                                            
      data wt6   / 0.171324492379170d0, 0.360761573048139d0,
     x             0.467913934572691d0, 0.467913934572691d0,
     x             0.360761573048139d0, 0.171324492379170d0/
                                                            
      data x7cor /-0.949107912342759d0,-0.741531185599394d0,
     x            -0.405845151377397d0, 0.000000000000000d0,
     x             0.405845151377397d0, 0.741531185599394d0,
     x             0.949107912342759d0/
                                       
      data wt7   / 0.129484966168870d0, 0.279705391489277d0,
     x             0.381830050505119d0, 0.417959183673469d0,
     x             0.381830050505119d0, 0.279705391489277d0,
     x             0.129484966168870d0/
                                       
      data x8cor /-0.960289856497536d0,-0.796666477413627d0,
     x            -0.525532409916329d0,-0.183434642495650d0,
     x             0.183434642495650d0, 0.525532409916329d0,
     x             0.796666477413627d0, 0.960289856497536d0/
                                                            
      data wt8   / 0.101228536290376d0, 0.222381034453374d0,
     x             0.313706645877887d0, 0.362683783378362d0,
     x             0.362683783378362d0, 0.313706645877887d0,
     x             0.222381034453374d0, 0.101228536290376d0/
                                                            
      data x9cor /-0.968160239507626d0,-0.836031107326636d0,
     x            -0.613371432700590d0,-0.324253423403809d0,
     x             0.000000000000000d0, 0.324253423403809d0,
     x             0.613371432700590d0, 0.836031107326636d0,
     x             0.968160239507626d0/
                                       
      data wt9   / 0.081274388361574d0, 0.180648160694857d0,
     x             0.260610696402935d0, 0.312347077040003d0,
     x             0.330239355001260d0, 0.312347077040003d0,
     x             0.260610696402935d0, 0.180648160694857d0,
     x             0.081274388361574d0/
c
      if (intg.eq.2) then
         nini = 2
         call kvimove (x2cor(1),1,GaussPts(1,1),2,nini)
         call kvimove (x2cor(1),1,GaussPts(2,1),2,nini)
         call kvmove (wt2(1),GaussWts(1),nini)
      elseif (intg.eq.3) then
         nini = 3
         call kvimove (x3cor(1),1,GaussPts(1,1),2,nini)
         call kvimove (x3cor(1),1,GaussPts(2,1),2,nini)
         call kvmove (wt3(1),GaussWts(1),nini)
      elseif (intg.eq.4) then
         nini = 4
         call kvimove (x4cor(1),1,GaussPts(1,1),2,nini)
         call kvimove (x4cor(1),1,GaussPts(2,1),2,nini)
         call kvmove (wt4(1),GaussWts(1),nini)
      elseif (intg.eq.5) then
         nini = 5
         call kvimove (x5cor(1),1,GaussPts(1,1),2,nini)
         call kvimove (x5cor(1),1,GaussPts(2,1),2,nini)
         call kvmove (wt5(1),GaussWts(1),nini)
      elseif (intg.eq.6) then
         nini = 6
         call kvimove (x6cor(1),1,GaussPts(1,1),2,nini)
         call kvimove (x6cor(1),1,GaussPts(2,1),2,nini)
         call kvmove (wt6(1),GaussWts(1),nini)
      elseif (intg.eq.7) then
         nini = 7
         call kvimove (x7cor(1),1,GaussPts(1,1),2,nini)
         call kvimove (x7cor(1),1,GaussPts(2,1),2,nini)
         call kvmove (wt7(1),GaussWts(1),nini)
      elseif (intg.eq.8) then
         nini = 8
         call kvimove (x8cor(1),1,GaussPts(1,1),2,nini)
         call kvimove (x8cor(1),1,GaussPts(2,1),2,nini)
         call kvmove (wt8(1),GaussWts(1),nini)
      else
         nini = 9
         call kvimove (x9cor(1),1,GaussPts(1,1),2,nini)
         call kvimove (x9cor(1),1,GaussPts(2,1),2,nini)
         call kvmove (wt9(1),GaussWts(1),nini)
      endif
           
      return
      end
c=================================================================
