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

c     Purpose: Plane linear elastic displacement model element routine

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  none
      integer   ndf,ndm,nst,isw
      integer   ityp, j,j1, k,k1, l,lint
      
      real*8    a11,a21,a31,a41,a12,a32,a42, w11,w12,w22
      real*8    e,xnu,alp,t0, xx,yy,zz, dv, xsj, sigr4

      integer   ix(*)
      real*8    xl(ndm,*),tl(*)
      real*8    d(*),ul(ndf,*),s(nst,*),p(*)
      real*8    eps(4),sigr(6),shp(3,9),sg(16),tg(16),wg(16),ang
      character wd(3)*12,yyy*80

      integer    maxa
      parameter (maxa=32000)

      real*8         aa
      common /adata/ aa(maxa)

      integer        numnp,numel,nummat,nen,neq
      common /cdata/ numnp,numel,nummat,nen,neq

      real*8          dm
      integer            n,ma,mct,iel,nel
      common /eldata/ dm,n,ma,mct,iel,nel

      integer         ior,iow
      common /iofile/ ior,iow

      data wd/'Plane Stress','Plane Strain','Axisymmetric'/

c     Go to correct array processor

      l    = d(5)
      k    = d(6)
      ityp = d(15)
      go to(1,2,3,3,3,3,7,8), isw

c     Input material properties

1     if(ior.lt.0) write(*,5000)
      call pintio(yyy,10)
      read(yyy,1000,err=110) e,xnu,d(4),l,k,ityp
11    if(ior.lt.0) write(*,5001)
      call pintio(yyy,10)
      read(yyy,1001,err=111) d(14),d(11),d(12),alp,t0

c     Set material parameter type and flags

      ityp = max(1,min(ityp,3))
      j = min(ityp,2)
      d(1) = e*(1.+(1-j)*xnu)/(1.+xnu)/(1.-j*xnu)
      d(2) = xnu*d(1)/(1.+(1-j)*xnu)
      d(3) = e/2./(1.+xnu)
      d(13)= d(2)*(j-1)
      if(d(14).le.0.0d0 .or. ityp.ge.2) d(14) = 1.0
      d(15) = ityp
      d(16) = e
      d(17) = xnu
      d(18) = -xnu/e
      l = min(4,max(1,l))
      k = min(4,max(1,k))
      d(5) = l
      d(6) = k
      d(9) = t0
      d(10)= e*alp/(1.-j*xnu)
      lint = 0
      write(iow,2000) wd(ityp),d(16),d(17),d(4),l,k,d(14),d(11),d(12),
     &             alp,t0
      if(ior.lt.0) write(*,2000) wd(ityp),d(16),d(17),d(4),l,k,
     &             d(14),d(11),d(12),alp,t0
      d(4) = d(4)*d(14)
      return

c     Read error messages

110   call perror('PCELM1',yyy)
      go to 1
111   call perror('PCELM1',yyy)
      go to 11

2     call ckisop(ix,xl,shp,ndm)
      return

c     Stiffness/residual computation

3     if(isw.eq.4) l = k
      if(l*l.ne.lint) call pgauss(l,lint,sg,tg,wg)

c     Compute integrals of shape functions

      do 340 l = 1,lint
        call shape(sg(l),tg(l),xl,shp,xsj,ndm,nel,ix,.false.)
        call stre01(d,xl,ul,tl,shp,eps,sigr,xx,yy,ndm,ndf,nel,ityp)
        xsj = xsj*wg(l)*d(14)

c       Compute jacobian correction

        if(ityp.le.2) then
          dv  = xsj
          xsj = 0.0
          zz  = 0.0
          sigr4 = -d(11)*dv
        else
          dv  = xsj*xx
          zz  = 1./xx
          sigr4 = sigr(4)*xsj - d(11)*dv
        endif
        j1 = 1

c       Compute the mass term

        if(isw.eq.5) then
          dv = dv*d(4)
          do 315 j = 1,nel
            p(j1  ) = p(j1) + shp(3,j)*dv
            p(j1+1) = p(j1)
            j1 = j1 + ndf
315       continue
        elseif(isw.eq.4) then
          call pstres(sigr,sigr(5),sigr(6),ang)

c         Output stresses and strains

          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,xx,sigr,ma,yy,eps,ang
          if(ior.lt.0) write(*,2002) n,xx,sigr,ma,yy,eps,ang
        else

c         Loop over rows

          do 330 j = 1,nel
            w11 = shp(1,j)*dv
            w12 = shp(2,j)*dv
            w22 = shp(3,j)*xsj

c           Compute the internal forces

            p(j1  ) = p(j1  ) - (shp(1,j)*sigr(1)+shp(2,j)*sigr(2))*dv
     &                        -  shp(3,j)*sigr4
            p(j1+1) = p(j1+1) - (shp(1,j)*sigr(2)+shp(2,j)*sigr(3))*dv
     &                        + d(12)*shp(3,j)*dv

c           Loop over columns (symmetry noted)

            if(isw.eq.3) then
              k1 = j1
              a11 = d(1)*w11 + d(2)*w22
              a21 = d(2)*w11 + d(1)*w22
              a31 = d(2)*(w11+w22)
              a41 = d(3)*w12
              a12 = d(2)*w12
              a32 = d(1)*w12
              a42 = d(3)*w11
              do 320 k = j,nel
                w11 = shp(1,k)
                w12 = shp(2,k)
                w22 = shp(3,k)*zz
                s(j1  ,k1  ) = s(j1  ,k1  ) + w11*a11+w22*a21+w12*a41
                s(j1+1,k1  ) = s(j1+1,k1  ) + (w11 + w22)*a12+w12*a42
                s(j1  ,k1+1) = s(j1  ,k1+1) + w12*a31 + w11*a41
                s(j1+1,k1+1) = s(j1+1,k1+1) + w12*a32 + w11*a42
                k1 = k1 + ndf
320           continue
            endif
            j1 = j1 + ndf
330       continue
        endif
340   continue

c     Make stiffness symmetric and compute a residual

      if(isw.eq.3) then
        do 360 j = 1,nst
        do 360 k = j,nst
          s(k,j) = s(j,k)
360     continue
      endif
      return

c     Compute the stress errors

7     if(l*l.ne.lint) call pgauss(l,lint,sg,tg,wg)
      call ster01(ix,d,xl,ul,tl,shp,aa(numnp+1),ndf,ndm,
     &   numnp,numel,sg,tg,wg,sigr,eps,lint,ityp)
      return

c     Compute the nodal stress values

8     if(l*l.ne.lint) call pgauss(l,lint,sg,tg,wg)
      call stcn01(ix,d,xl,ul,tl,shp,aa,aa(numnp+1),ndf,ndm,nel,
     &   numnp,sg,tg,sigr,eps,lint,ityp)
      return

c     Formats for input-output

1000  format(3f10.0,3i10)
1001  format(8f10.0)
2000  format(/5x,a12,' Linear Elastic Element'//
     & 10x,'Modulus',e18.5/10x,'Poisson ratio',f8.5/10x,'Density',e18.5/
     & 10x,'Gauss pts/dir',i3/10x,'Stress pts',i6/10x,'Thickness',e16.5/
     & 10x,'1-gravity',e16.5/10x,'2-gravity',e16.5/10x,'Alpha',e20.5/
     & 10x,'Base temp',e16.5/)
2001  format(5x,'Element Stresses'//' elmt  1-coord',2x,'11-stress',2x,
     & '12-stress',2x,'22-stress',2x,'33-stress',3x,'1-stress',3x,
     & '2-stress'/' matl  2-coord',2x,'11-strain',2x,'12-strain',2x,
     & '22-strain',2x,'33-strain',6x,'angle'/39(' -'))
2002  format(i5,0p1f9.3,1p6e11.3/i5,0p1f9.3,1p4e11.3,0p1f11.2/)
5000  format(' Input: E, nu, rho, pts/stiff, pts/stre',
     & ', type(1=stress,2=strain,3=axism)',/3x,'>',$)
5001  format(' Input: Thickness, 1-body force, 2-body force, alpha,'
     &       ,' Temp-base'/3x,'>',$)

      end

      subroutine stre01(d,xl,ul,tl,shp,eps,sig,xx,yy,ndm,ndf,nel,ityp)

c     Purpose: Compute stress and strain for point

c     Inputs:
c        d(*)     - Material set parameters
c        xl(ndm,*)- Nodal coordinates for element
c        ul(ndf,*) - Nodal solution parameters for element
c        tl(*)     - Nodal temperatures for element
c        shp(3,*)  - SHape functions and derivatives for point
c        ndm       - Spatial dimension of FEM mesh
c        ndf       - Number dof/node
c        nel       - Number of nodes on element
c        ityp      - Problem type

c     Outputs:
c        eps(*)    - Strains at point
c        sig(*)    - Stresses at point
c        xx        - x-coordinate at point
c        yy        - y-coordinate at point

      implicit none

      integer  j,ndm,ndf,nel,ityp
      real*8   xx,yy,ta

      real*8 xl(ndm,*),tl(*)
      real*8 d(*),ul(ndf,*),shp(3,*),eps(4),sig(4)

c     Compute strains and coordinates

      call pconsd(eps,4,0.0d0)
      xx = 0.0
      yy = 0.0
      ta = -d(9)
      do 100 j = 1,nel
        xx = xx + shp(3,j)*xl(1,j)
        yy = yy + shp(3,j)*xl(2,j)
        ta = ta + shp(3,j)*tl(j)
        eps(1) = eps(1) + shp(1,j)*ul(1,j)
        eps(2) = eps(2) + shp(1,j)*ul(2,j) + shp(2,j)*ul(1,j)
        eps(3) = eps(3) + shp(2,j)*ul(2,j)
        if(ityp.eq.3) eps(4) = eps(4) + shp(3,j)*ul(1,j)
100   continue
      ta = ta*d(10)

c     Compute stresses

      if(ityp.gt.2) then
        if(xx.ne.0.0d0) then
          eps(4) = eps(4)/xx
        else
          eps(4) = eps(1)
        endif
        sig(4) = d(1)*eps(4) + d(2)*(eps(1) + eps(3)) - ta
      else
        sig(4) = d(13)*(eps(1) + eps(3)) - ta
      endif
      sig(1) = d(1)*eps(1) + d(2)*(eps(3) + eps(4)) - ta
      sig(2) = d(3)*eps(2)
      sig(3) = d(1)*eps(3) + d(2)*(eps(1) + eps(4)) - ta
      if(ityp.eq.1) eps(4) = d(18)*(sig(1) + sig(3))

      end

      subroutine stcn01(ix,d,xl,ul,tl,shp,dt,st,ndf,ndm,nel,numnp,
     &                  sg,tg,sig,eps,lint,ityp)

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        ix(*)     - Nodal connections for element
c        tl(*)     - Nodal temperatures for element
c        shp(3,*)  - Shape function and derivatives
c        ndf       - Number dof/node
c        ndm       - Spatial dimension of FEM mesh
c        nel       - Number of nodes on element
c        numnp     - 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        ityp      - Problem type

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

      implicit none

      integer  ndf,ndm,nel,numnp,lint,ityp
      integer  i, l,ll
      real*8   gr,xsj,xsji, xx,yy

      integer  ix(*)
      real*8   dt(numnp),st(numnp,*),xl(ndm,*),tl(*)
      real*8   eps(4),sig(4),ul(ndf,*),shp(3,4),d(*),sg(9),tg(9)

      real*8          eerror,elproj,ecproj,efem,enerr,ebar
      common /errind/ eerror,elproj,ecproj,efem,enerr,ebar

      gr = (1.+d(17))/d(16)
      do 130 l = 1,lint

        call shape(sg(l),tg(l),xl,shp,xsj,ndm,nel,ix,.false.)
        call stre01(d,xl,ul,tl,shp,eps,sig,xx,yy,ndm,ndf,nel,ityp)
        enerr = enerr+((sig(1)+sig(3)+sig(4))**2)*d(18)*xsj 
     &     + gr*xsj*sig(2)**2

        do 110 i = 1,4
          enerr  = enerr  + gr*sig(i)**2*xsj
110     continue

        do 120 i = 1,nel
          ll = abs(ix(i))
          if(ll.gt.0) then
            xsji = xsj*shp(3,i)
            dt(ll) = dt(ll) + xsji
            st(ll,1) = st(ll,1) + sig(1)*xsji
            st(ll,2) = st(ll,2) + sig(2)*xsji
            st(ll,3) = st(ll,3) + sig(3)*xsji
            st(ll,4) = st(ll,4) + sig(4)*xsji
          endif
120     continue
130   continue

      end

      subroutine ster01(ix,d,xl,ul,tl,shp,st,ndf,ndm,
     &                  numnp,numel,sg,tg,wg,sig,eps,lint,ityp)

c     Purpose: Compute error estimator values

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        ix(*)     - Nodal connections for element
c        tl(*)     - Nodal temperatures for element
c        shp(3,*)  - Shape function and derivatives
c        st(*,*)   - Projected nodal stresses
c        ndf       - Number dof/node
c        ndm       - Spatial dimension of FEM mesh
c        nel       - Number of nodes on element
c        numnp     - Number of nodes in mesh
c        numel     - Number of elements in mesh
c        sg(*)     - Quadrature points
c        tg(*)     - Quadrature points
c        wg(*)     - quadrature weights
c        sig(4)    - Stresses at quadrature points
c        eps(4)    - Stresses at quadrature points
c        lint      - number of quadrature points
c        ityp      - Problem type

c     Outputs:
c        none      - done through common blocks

      implicit  none

      integer   ndf,ndm,numnp,numel,lint,ityp
      integer   i,ii, j, ll
      real*8    gr,xsj, psi,psis, xx,yy

      integer   ix(*)
      real*8    st(numnp,*),xl(ndm,*),tl(*)
      real*8    d(*),ul(ndf,*)
      real*8    shp(3,4),sig(6),sigp(4),dsig(4),eps(4)
      real*8    sg(16),tg(16),wg(16)

      integer         ior,iow
      common /iofile/ ior,iow

      real*8          dm
      integer            n,ma,mct,iel,nel
      common /eldata/ dm,n,ma,mct,iel,nel

      real*8          eerror,elproj,ecproj,efem,enerr,ebar
      common /errind/ eerror,elproj,ecproj,efem,enerr,ebar

c     Stress error computations

      psis = 0.0
      psi  = 0.0
      gr   = (1.+d(17))/d(16)
      do 200 ii = 1,lint
        call shape(sg(ii),tg(ii),xl,shp,xsj,ndm,nel,ix,.false.)
        call stre01(d,xl,ul,tl,shp,eps,sig,xx,yy,ndm,ndf,nel,ityp)
        xsj = xsj*wg(ii)
        do 100 i = 1,4
          sigp(i)= 0.0d0
100     continue
        do 110 i = 1,nel
          ll = abs(ix(i))
          if(ll.ne.0) then
            do 105 j = 1,4
              sigp(j) = sigp(j) + shp(3,i)*st(ll,j)
105         continue
          endif
110     continue

c       Compute the integral of the stress squares for error indicator use

        do 120 i = 1,4
          dsig(i) = sigp(i)-sig(i)
          efem    = efem   + sig(i)*sig(i)*xsj
          ecproj  = ecproj + sigp(i)*sigp(i)*xsj
          psis    = psis   + (dsig(i)**2)*xsj
          psi     = psi    +  gr*(dsig(i)**2)*xsj
120     continue
        psi = psi +  gr*dsig(2)**2*xsj + 
     &      ((dsig(1)+dsig(3)+dsig(4))**2)*d(18)*xsj
200   continue

      eerror = eerror + psis
      if(elproj.ne.0.0d0) then
        psi    = sqrt(abs(psi))/ebar
        psis   = 20.0*sqrt(abs(psis)/elproj*numel)
        if(mct.eq.0) then
          write(iow,2000)
          if(ior.lt.0) write(*,2000)
          mct = 50
        endif
        mct = mct - 1
        write(iow,2001) n,psis,psi
        if(ior.lt.0) write(*,2001) n,psis,psi
      endif

2000  format('   M e s h   R e f i n e m e n t s   f o r   5%',
     &       '   E r r o r'//'    elmt   h-sigma    h-energy'/)
2001  format(i8,1p2e12.4)

      end
