      subroutine elmt02(d,ul,xl,ix,tl,s,p,ndf,ndm,nst,isw)

c     Purpose: Plane/Axisymmetric Linear Element -- bbar formulation

c     Inputs:
c        d(*)      - Material set parameters
c        ul(ndf,*) - Solution variables for element
c        xl(ndf,*) - Nodal coordinates for element
c        ix(*)     - Nodal connections for element
c        tl(*)     - Nodal temperatures for element
c        ndf       - Number dof/node
c        ndm       - Spatial dimension of FEM mesh
c        nst       - Size of element matrix
c        isw       - Switch parameter to control action taken by routine

c     Outputs:
c        s(nst,*)  - Element matrix
c        p(*)      - Element vector

      implicit real*8 (a-h,o-z)

      logical flg
      integer   ix(*)
      real*8 xl(ndm,*),tl(*),aa,dm
      real*8 d(*),ul(ndf,*),s(nst,*),p(*),eps(4),sig(6),bbar(4,2,4),
     1   bbd(4,2),sg(4),tg(4),wg(4),shp3(4,4),shp(3,4,4),xsj(4),ang,
     2   siga(6)

      integer        maxa
      parameter     (maxa=32000)
      common /adata/ aa(maxa)
      common /cdata/ numnp,numel,nummat,nen,neq
      common /eldata/ dm,n,ma,mct,iel,nel
      common /elcom2/ g(2,4),ad(4,4)
      common /iofile/ ior,iow
c.... go to correct array processor
      go to(1,2,3,3,5,3,7,3), isw
      return
c.... input material properties
   1  call matl02(d,ityp,ib)
      d(5) = ityp
      d(6) = ib
      l    = 2
      call pgauss(l,lint,sg,tg,wg)
      return
c.... check for element errors
   2  call ckisop(ix,xl,shp,ndm)
      return
c
c.... compute tangent stiffness and residual force vector
   3  type = d(5)
      ib   = d(6)
c.... compute volumetric integrals
      call pconsd(g,8,0.0d0)
      do 300 l = 1,lint
        call shape(sg(l),tg(l),xl,shp(1,1,l),xsj(l),ndm,4,ix,.false.)
        call gvc02(shp(1,1,l),shp3(1,l),xsj(l),wg(l),xl,type,ndm)
300   continue
      vol  = xsj(1) + xsj(2) + xsj(3) + xsj(4)
      do 310 i = 1,4
        g(1,i) = g(1,i)/vol
        g(2,i) = g(2,i)/vol
310   continue
      if(isw.eq.4) go to 4
      if(isw.eq.8) go to 8
      flg  = isw .eq. 3
      do 400 l = 1,lint
c.... compute stress, strain, and material moduli
        call strn02(shp(1,1,l),xl,ul,type,xr0,xz0,ndm,ndf,eps)
        call modl02(d,ul,eps,sig,xsj(l),ndf,ib)
        i1 = 0
        do 380 i = 1,4
c.... compute the internal stress divergence term
          p(i1+1) = p(i1+1) - shp(1,i,l)*sig(1) - shp(2,i,l)*sig(2)
     1                      - shp3(i,l)*sig(4)
          p(i1+2) = p(i1+2) - shp(2,i,l)*sig(3) - shp(1,i,l)*sig(2)
c.... compute stiffness
          if(flg) then
c.... compute b-bar matrix
            call bmat02(shp3(i,l),shp(1,i,l),g(1,i),bbar(1,1,i),ib)
            do 320 ii = 1,2
            do 320 jj = 1,4
              bbd(jj,ii) = dot(bbar(1,ii,i),ad(1,jj),4)
320         continue
            j1 = 0
            do 360 j  = 1,i
              do 340 ii = 1,2
              do 340 jj = 1,2
                s(ii+i1,jj+j1) = s(ii+i1,jj+j1)
     1                         + dot(bbd(1,ii),bbar(1,jj,j),4)
340           continue
360         j1 = j1 + ndf
          endif
380     i1 = i1 + ndf
400   continue
c.... form lower part by symmetry
      if(flg) then
        do 410 i = 1,nst
        do 410 j = 1,i
          s(j,i) = s(i,j)
410     continue
      endif
      return
c
c.... compute stresses at center of element
   4  call pconsd(sig,6,0.0d0)
      rr = 0.0
      zz = 0.0
      do 430 l = 1,lint
c.... compute stress, strain, and material moduli
        call strn02(shp(1,1,l),xl,ul,type,xr0,xz0,ndm,ndf,eps)
        call modl02(d,ul,eps,siga,1.d0,ndf,ib)
        do 420 i = 1,5
          sig(i) = sig(i) + 0.25d0*siga(i)
420     continue
        rr = rr + 0.25d0*xr0
        zz = zz + 0.25d0*xz0
430   continue
      yld = sig(5)
      call pstres(sig,sig(5),sig(6),ang)
c.... output stresses
      mct = mct - 2
      if(mct.le.0) then
        call prthed(iow)
        write(iow,2001)
        if(ior.lt.0) write(*,2001)
        mct = 50
      endif
      write(iow,2002) n,ma,sig,rr,zz,yld,ang
      if(ior.lt.0) write(*,2002) n,ma,sig,rr,zz,yld,ang
      return
c
c.... compute lumped mass matrix
   5  do 530 l = 1,lint
        call shape(sg(l),tg(l),xl,shp,xsj,ndm,4,ix,.false.)
c.... compute radius and multiply into jacobian for axisymmetry
        if(d(5).ne.0.0d0) then
          rr = 0.0d0
          do 500 j = 1,4
            rr = rr + shp(3,j,1)*xl(1,j)
500       continue
          xsj(1) = xsj(1)*rr
        endif
        xsj(1) = wg(l)*xsj(1)*d(4)
c.... for each node j compute db = rho*shape*dv
        j1 = 1
        do 520 j = 1,4
          p(j1  ) = p(j1) + shp(3,j,1)*xsj(1)
          p(j1+1) = p(j1)
          j1 = j1 + ndf
520     continue
530   continue
   7  return
c
c.... stress computations for nodes
   8  call stcn02(ix,d,xl,ul,shp,aa,aa(numnp+1),ndf,ndm,
     1   numnp,sg,tg,sig,eps,lint,type,ib)
c
c.... formats for input-output
1000  format(2i5)
1001  format(a5,1x,2i10,f10.0/9f10.0)
2001  format('  Element Stresses'//'  elmt  matl  11-stress  12-stress',
     1 '  22-stress  33-stress   1-stress   2-stress'/'  1-coord',
     2 '  2-coord  ',33x,'yield    angle')
2002  format(2i6,1p6e11.3/0p2f9.3,30x,1p1e11.3,0pf8.2/1x)
      end
c
      subroutine gvc02(shp,shp3,xsj,wg,xl,type,ndm)

c     Purpose: Compute volumetric integrals for b-bar

c     Inputs:
c        shp(3,*)    - Shape fucntions and derivatives
c        shp3(*)     - Shape function for axisymmetric, otherwise 1.0
c        xsj         - Jacobian
c        xl(ndm,*)   - Nodal coordinates for element
c        type        - Type of problem
c        ndm         - Spatial dimension of mesh

c     Outputs:
c        none        - stored in common block

      implicit real*8 (a-h,o-z)

      real*8 xl(ndm,4)
      real*8 shp(3,4),shp3(4)
      common /elcom2/ g(2,4),ad(4,4)
      if(type.ne.0.0d0) then
        rr = 0.0d0
        do 100 i = 1,4
          rr = rr + shp(3,i)*xl(1,i)
100     continue
        xsj = xsj*rr
      endif
      xsj = xsj*wg
      do 110 i = 1,4
        shp3(i) = 0.0d0
        if(type.ne.0.0d0) shp3(i) = shp(3,i)/rr
        g(1,i) = g(1,i) + (shp(1,i) + shp3(i))*xsj
        g(2,i) = g(2,i) +  shp(2,i)*xsj
110   continue
      end
c
      subroutine bmat02(sh3,shp,g,bbar,ib)

c     Purpose: Form B-bar matrix for plane and axisymmetric problems

c     Inputs:
c        sh3       - Shape functioin for axisymmetry, otherwise 1.0
c        shp(30    - Shape functions and derivatives for node
c        g(*)      - b-bar for volumetric terms
c        ib        - Switch:  =0 for B-bar; not 0 for displacement

c     Outputs:
c        bbar(4,2) - B-bar matrix for node

      implicit real*8 (a-h,o-z)

      real*8 shp(3),g(1),bbar(4,2)

      bbar(1,1) = shp(1)
      bbar(2,1) = 0.0d0
      bbar(3,1) = sh3
      bbar(4,1) = shp(2)
      bbar(1,2) = 0.0d0
      bbar(2,2) = shp(2)
      bbar(3,2) = 0.0d0
      bbar(4,2) = shp(1)

c.... correct for the b-bar effects

      if(ib.eq.0) then
        bb1 = (g(1) - shp(1) - sh3)/3.0d0
        bb2 = (g(2) - shp(2))/3.0d0
        do 100 i = 1,3
          bbar(i,1) = bbar(i,1) + bb1
          bbar(i,2) = bbar(i,2) + bb2
100     continue
      endif
      end
c
      subroutine strn02(shp,xl,ul,type,xr0,xz0,ndm,ndf,eps)

c     Purpose: Compute strain for point

c     Inputs:
c        shp(3,*)  - Shape functions and derivatives for point
c        xl(ndm,*) - Nodal coordinates for element
c        ul(ndf,*) - Nodal solution parameters for element
c        type      - Type of problem
c        ndm       - Spatial dimension of FEM mesh
c        ndf       - Number dof/node

c     Outputs:
c        xr0       - x-(r)-coordinate at point
c        xz0       - y-(z)-coordinate at point
c        eps(*)    - Strains at point

      implicit real*8 (a-h,o-z)
      real*8 xl(ndm,1)
      real*8 ul(ndf,1),eps(4),shp(3,1)
c.... compute strain and incremental tensors for constitutive equations
      xr0   = 0.0d0
      xz0   = 0.0d0
      call pconsd(eps,4,0.0d0)
      do 310 k = 1,4
        xr0   = xr0   + shp(3,k)*xl(1,k)
        xz0   = xz0   + shp(3,k)*xl(2,k)
        eps(1) = eps(1) + ul(1,k)*shp(1,k)
        eps(2) = eps(2) + ul(2,k)*shp(2,k)
        eps(3) = eps(3) + shp(3,k)*ul(1,k)
        eps(4) = eps(4) + (ul(2,k)*shp(1,k) + ul(1,k)*shp(2,k))/2.
310   continue
      eps(3) = type*eps(3)/xr0
      end
c
      subroutine stcn02(ix,d,xl,ul,shp,dt,st,ndf,ndm,
     &                  nnp,sg,tg,sig,eps,lint,type,ib)

c     Purpose: Project stresses onto nodes

c     Inputs
c        ix(*)     - Nodes connected to element
c        d(*)      - Material set parameters
c        xl(ndf,*) - Nodal coordinates for element
c        ul(ndf,*) - Solution variables for element
c        shp(3,*)  - Shape function and derivatives
c        ndf       - Number dof/node
c        ndm       - Spatial dimension of FEM mesh
c        nnp       - Number of nodes in mesh
c        sg(*)     - Quadrature points
c        tg(*)     - Quadrature points
c        sig(4)    - Stresses at quadrature points
c        eps(4)    - Stresses at quadrature points
c        lint      - number of quadrature points
c        type      - Problem type
c        ib        - Formulation: =0 B-bar; else displacement model

c     Outputs:
c        dt(*)     - Diagonal projection matrix
c        st(*,*)   - Integral of FEM stresses over element

      implicit real*8 (a-h,o-z)
c     integer*2 ix(1)
      integer   ix(1)
      real*8 dt(nnp),st(nnp,1),xl(ndm,1)
      real*8 d(1),ul(ndf,1),sig(1),eps(1),shp(3,4),sg(1),tg(1)
      do 200 jj = 1,lint
c.... compute stresses at nodes from history terms
        call shape(sg(jj),tg(jj),xl,shp,xsj,ndm,4,ix,.false.)
        call strn02(shp,xl,ul,type,rr,zz,ndm,ndf,eps)
        call modl02(d,ul,eps,sig,1.0d0,ndf,ib)
        do 100 ii = 1,4
          ll = abs(ix(ii))
          if(ll.gt.0) then
            shpj = shp(3,i)*xsj
            dt(ll) = dt(ll) + shpj
            do 50 i = 1,4
 50         st(ll,i) = st(ll,i) + sig(i)*shpj
          endif
100     continue
200   continue
      end
c
      subroutine matl02(d,it,ib)

c     Purpose:  Input material set parameters for elements

c     Outputs:
c           d(*)  -  Constants for constitutive equation.
c           it    -  Geometry type (0=plane, 1=axisymmetric)
c           ib    -  Formulation type (ib = 0 b-bar,
c                      ib = 1 normal b matrix).
c           nh1   -  number of history variables needed for each
c                      stress point (i.e., at gauss points in element).

      implicit real*8 (a-h,o-z)
      real*8 d(15)
      character*24 wa(2),yyy*80
c.... parameter specification for FEAP materials
      common /hdata/ nh1,nh2
      common /iofile/ ior,iow
      common /ydata/ yyy

      data wa/' P l a n e   S t r a i n',' A x i s y m m e t r i c'/
c
c....  input the material parameters
  1   if(ior.lt.0) write(*,3000)
      call pintio(yyy,10)
      read(yyy,1000,err=101) it,d(4)
      it  = max(0,min(1,it))
      ib  = 0
      nh  = 9
      nh1 = 36
      if(ior.lt.0) write(*,3001)
      call pintio(yyy,10)
      read(yyy,1001,err=101) ee,xnu,d(11),d(12),d(13)
      d(1)    = ee/(1. - 2.*xnu)/3.0d0
      d(2)    = ee/(1.+xnu)/2.
      if(d(11).ne.0.0d0) then
        write(iow,2000) wa(it+1),ee,xnu,d(4),d(11),d(12),d(13)
        if(ior.lt.0) write(*,2000) wa(it+1),ee,xnu,d(4),d(11),
     1                                             d(12),d(13)
      else
        write(iow,2001) wa(it+1),ee,xnu,d(4)
        if(ior.lt.0) write(*,2001) wa(it+1),ee,xnu
      endif
      return
101   call perror('PCELM2',yyy)
      go to 1
c.... formats
1000  format(i10,f10.0)
1001  format(8f10.0)
2000  format(2x,'E l a s t i c / P l a s t i c   M a t e r i a l'/
     +       2x,'v o n  M i s e s  Y i e l d - ',a24/
     +  10x,'Youngs Modulus (E)',e17.5/
     +  10x,'Poisson Ratio (nu)',e17.5/
     +  10x,'Mass Density (rho)',e17.5/
     +  10x,'Yield Stress      ',e17.5/
     +   8x,'Linear hardening moduli'/
     +  10x,'Isotropic Hardening',e16.5/
     +  10x,'Kinematic Hardening',e16.5/)
2001  format(2x,'E l a s t i c   M a t e r i a l'/
     +       2x,a24,'   A n a l y s i s'/
     +  10x,'Youngs Modulus (E)',e17.5/
     +  10x,'Poisson Ratio (nu)',e17.5/
     +  10x,'Mass Density (rho)',e17.5/)
3000  format(' Input: it, rho'/
     +  4x,'it = 0: Plane'/4x,'it = 1: Axisymmetric'/' >',$)
3001  format(' Input (Elastic/Plastic) : E, nu, Y, H-iso, H-kine'/
     +       '       (Linear  Elastic) : E, nu'/3x,'>',$)
      end
c
      subroutine modl02(d,ul,eps,sig,xsj,ndf,ib)

c     Purpose:  Form material tangent and stress state at point

c     Inputs:
c          d(*)    -  material parameters
c          eps(4)  -  current strains at point
c          h(*)    -  history terms at point
c          nh      -  number of history terms

c     Ouputs:
c          ad(4,4) -  current material tangent moduli
c          sig(4)  -  stresses at point.
c          sig(5)  -  yield state at point.

      implicit real*8 (a-h,o-z)
      real*8 d(*),ul(ndf,*),eps(4),sig(5)
      common /elcom2/ g(2,4),ad(4,4)
      common /hdata/  nh1,nh2
      common h(1)

c....  compute material moduli and stresses
      call pconsd(ad,16,0.0d0)
      call elpl02(d,eps,h(nh2),h(nh2+4),h(nh2+8),sig,ul,ndf,ib)
      nh2 = nh2 + 9
c
c....  multiply by jacobian
      do 110 i = 1,4
        sig(i) = sig(i)*xsj
        do 100 j = 1,4
           ad(i,j) = ad(i,j)*xsj
100     continue
110   continue
c
c.... reorder stresses for stress divergence calculations and prints
      temp   = sig(4)
      sig(4) = sig(3)
      sig(3) = sig(2)
      sig(2) = temp
      end
c
      subroutine elpl02(d,eps,epln,alph,ep,sig,ul,ndf,ib)
c
c     Purpose: Elasto-plastic model with isotropic / kinematic hardening
c
c     Inputs: - for t-n
c        d      - array of material constants
c        eps    - strains (at t-n+1)
c        sig    - stresses
c        alph   - back stress
c        ep     - effective plastic strain

c     Outputs: - at t-n+1
c        ep     - effective plastic strain
c        sig    - stress tensor
c        alph   - back stress tensor
c        ad     - "tangent" matrix

      implicit real*8 (a-h,o-z)
      real*8 dm
      real*8 d(*),ul(ndf,*),eps(4),sig(5),alph(4),epln(4),psi(4),en(4)
      common /eldata/ dm,n,ma,mct,iel,nel
      common /elcom2/ g(2,4),ad(4,4)

c.... set parameters
      data tt/.816496580927726D+00/
      oneg = d(2)
      twog = d(2) + d(2)
      bulk = d(1)
      elam = bulk - twog/3.d0
c
c.... compute the trial deviator stress
      treps  = (eps(1) + eps(2) + eps(3))/3.0d0
      do 100 i = 1,3
        eps(i) = eps(i) - treps
100   continue
      do 110 i = 1,4
        sig(i) = twog*(eps(i) - epln(i))
        psi(i) =       sig(i) - alph(i)
110   continue

c.... set up elastic tangent
      do 130 i = 1,3
        do 120 j = 1,3
          ad(i,j) = elam
120     continue
        ad(i,i) = elam + twog
130   continue
      ad(4,4) = oneg
c
c.... compute the yield state - J2d
      if(d(11).gt.0.0d0) then
        radius = tt*(d(11) + d(12)*ep)
        psitr  = sqrt(psi(1)**2+psi(2)**2+psi(3)**2+2.d0*psi(4)**2)
        sig(5) = psitr/tt/d(11)
c
c.... compute plasticity solution state
        if (psitr.gt.radius) then
          beta   = 1.d0/(1.d0 + (d(12) + d(13))/(3.d0*oneg))
c         gamn   = (beta*(psitr - radius) - 1.d-10*d(11))/twog
          gamn   = beta*(psitr - radius)/twog
		  gamn   = gamn*(1.d0 -1.d-10)
          ep     = ep + tt*gamn

c.... gam1 ensures stress is slightly outside yield surface.
          gam1   = gamn*twog
          gam2   = (d(13)+d(13))*gamn/3.d0
          sig(5) = (psitr - gam1 - gam2)/tt/d(11)
          do 140 i = 1,4
            en(i)   = psi(i)/psitr
            sig(i)  =  sig(i)  - gam1*en(i)
            alph(i) =  alph(i) + gam2*en(i)
            epln(i) =  epln(i) + gamn*en(i)
140       continue

c.... plastic modification for tangent
          gam1 = gam1/psitr*twog
          gam2 = gam1/3.0d0
c.... I-dev part
          do 160 i = 1,3
            do 150 j = 1,3
              ad(i,j) = ad(i,j) + gam2
 150        continue
            ad(i,i) = ad(i,i) - gam1
 160      continue
          ad(4,4) = ad(4,4) - 0.5*gam1
c....  n x n part
          gam1 = gam1 - twog*beta
          do 180 i = 1,4
            gam2    = gam1*en(i)
            do 170 j = 1,4
              ad(i,j) = ad(i,j) + gam2*en(j)
 170        continue
 180      continue
        endif
      endif
c
c.... compute trace of strain and add pressure term
      if(ib.eq.0) then
        treps = 0.0d0
        do 190 i = 1,4
          treps = treps + g(1,i)*ul(1,i) + g(2,i)*ul(2,i)
 190    continue
      else
        treps = treps*3.d0
      endif
      press = bulk*treps
      do 200 i = 1,3
        sig(i) = sig(i) + press
 200  continue
c
      end
