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

c     Purpose: Plane stress/strain Pian Sumihara Element

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)
      integer   ix(1)
      real*8 xl(ndm,*),tl(*),dm,ap
      real*8 d(*),ul(ndf,*),s(nst,*),p(*),sig(4),p1,p2,p3
      character*4 wd(2),yyy*80

      integer        maxa
      parameter     (maxa=32000)
      common /adata/ ap(maxa)
      common /cdata/  numnp,numel,nummat,nen,neq
      common /eldata/ dm,n,ma,mct,iel,nel
      common /elcom3/ xs,xt,xh,ys,yt,yh,xj0,xj1,xj2,a1(3),a2(3),beta(5)
      common /iofile/ ior,iow
      data wd/'ress','rain'/
c.... go to correct array processor
      go to(1,2,3,3,3,3,7,3), isw
c.... input/output material properties
   1  if(ior.lt.0) write(*,3000)
      call pintio(yyy,10)
      read(yyy,1000,err=11) (d(i),i=8,11),is
      is = max(1,min(is,2))
      d(12)= is
      d(2) = d(9)*d(8)/(1.+d(9))/(1.-is * d(9))
      d(4) = d(8)/(1.+d(9))
      d(3) = d(4)/2.
      d(1) = d(2) + d(4)
c.... set parameters for plane stress (is = 1)
      if(is.eq.1) then
        if(d(11).le.0.0d0) d(11) = 1.0
        d(5) =  1./d(8)
        d(6) = -d(9)/d(8)
        d(13)=  0.0
c.... set parameters for plane strain (is = 2)
      else
        d(11)=  1.0
        d(5) =  (1.-d(9))/d(4)
        d(6) = -d(9)/d(4)
        d(13)=  d(9)
      endif
      write(iow,2000) wd(is),(d(i),i=8,11)
      if(ior.lt.0) write(*,2000) wd(is),(d(i),i=8,11)
      d(10) = d(10)*d(11)
      return
 11   call perror('PCEL03',yyy)
      go to 1
c.... check mesh
   2  call ckisop(ix,xl,ap,ndm)
      return
c.... Compute the Pian-Sumihara arrays for elastic: compute jacobian
   3  xs = (-xl(1,1)+xl(1,2)+xl(1,3)-xl(1,4))/4.
      ys = (-xl(2,1)+xl(2,2)+xl(2,3)-xl(2,4))/4.
      xt = (-xl(1,1)-xl(1,2)+xl(1,3)+xl(1,4))/4.
      yt = (-xl(2,1)-xl(2,2)+xl(2,3)+xl(2,4))/4.
      xh = ( xl(1,1)-xl(1,2)+xl(1,3)-xl(1,4))/4.
      yh = ( xl(2,1)-xl(2,2)+xl(2,3)-xl(2,4))/4.
      xj0 = xs*yt - xt*ys
      xj1 = xs*yh - xh*ys
      xj2 = xh*yt - xt*yh
      if(isw.eq.5) go to 5
      ssa = xj1/xj0/3.
      tta = xj2/xj0/3.
c.... form stiffness for elastic part and compute the beta parameters
      call pian03(d,ul,s,p,nst,ndf,isw)
      if(isw.eq.4) go to 4
      if(isw.eq.8) go to 8
c.... compute symetric part of s
        do 334 i = 1,nst
          do 332 j = i,nst
            s(j,i) = s(i,j)
 332      continue
 334    continue
      return
c.... compute the stresses
   4  is = d(12)
c.... compute the stresses at the center and the specified points
      ssg    = -ssa*beta(5)
      ttg    = -tta*beta(4)
      sig(1) =  beta(1) + a1(1)*ttg + a2(1)*ssg
      sig(2) =  beta(2) + a1(2)*ttg + a2(2)*ssg
      sig(3) =  beta(3) + a1(3)*ttg + a2(3)*ssg
      sig(4) =  d(13)*(sig(1)+sig(2))
 330  continue
      call pstres(sig,p1,p2,p3)
      xx = (xl(1,1)+xl(1,2)+xl(1,3)+xl(1,4))/4.0
      yy = (xl(2,1)+xl(2,2)+xl(2,3)+xl(2,4))/4.0
      mct = mct - 1
      if(mct.gt.0) go to 450
      call prthed(iow)
      write(iow,2001) wd(is)
      if(ior.lt.0) write(*,2001) wd(is)
      mct = 25
450   write(iow,2002) n,ma,xx,yy,p1,p2,p3,(sig(i),i=1,4)
      if(ior.lt.0) write(*,2002) n,ma,xx,yy,p1,p2,p3,
     1             (sig(i),i=1,4)
      return
c.... compute a lumped mass matrix
   5  p(      1) = (xj0-(xj1+xj2)/3.)*d(10)
      p(  ndf+1) = (xj0+(xj1-xj2)/3.)*d(10)
      p(2*ndf+1) = (xj0+(xj1+xj2)/3.)*d(10)
      p(3*ndf+1) = (xj0-(xj1-xj2)/3.)*d(10)
      do 510 i = 2,nst,ndf
        p(i) = p(i-1)
510   continue
      return
c.... error estimator goes here!
   7  continue
      return
c.... compute the nodal stress values
   8  call stcn03(ix,d,ssa,tta,ap,ap(numnp+1),numnp)
      return
c.... format statements
1000  format(4f10.0,i10)
2000  format(5x,'Plane St',a4,' Element'//10x,'modulus      =',e12.5/
     1  10x,'poisson ratio=', f8.5/10x,'mass density =',e12.5/
     2  10x,'thickness    =', e12.5)
2001  format(5x,'Plane St',a4,' Stresses'//'  element material',
     1  5x,'1-coord',5x,'2-coord',8x,'sig1',8x,'sig2',3x,'angle'/
     2  38x,'s-11',8x,'s-12',8x,'s-22',8x,'s-33'/1x)
2002  format(2i9,2f12.4,2e12.4,f8.2/30x,4e12.4)
3000  format(' Input: e, nu, rho, th, is (1=pl.stress,2=pl.strain)'/
     1    3x,'mate>',$)
      end
c
      subroutine pian03(d,ul,s,p,nst,ndf,isw)

c     Pian-Sumihara stiffness matrix developed explicitly

      implicit real*8 (a-h,o-z)
      real*8 ul(ndf,1),s(nst,1),p(1),d(1),rr(5),x1(4),x2(4),y1(4),y2(4)
     1      ,ax(4),bx(4),cx(4),ay(4),by(4),cy(4)
      common /elcom3/ xs,xt,xh,ys,yt,yh,xj0,xj1,xj2,a1(3),a2(3),beta(5)
c..1.) set up stress interpolants for the 4-5 term
      r1 = xj1/xj0
      r2 = xj2/xj0
      a1(1) = xs*xs
      a1(2) = ys*ys
      a1(3) = xs*ys
      a2(1) = xt*xt
      a2(2) = yt*yt
      a2(3) = xt*yt
c..2.) set up shape function coefficients - jacobian weighted
      ax(1) = -yt + ys
      ax(2) =  yt + ys
      ax(3) = -ax(1)
      ax(4) = -ax(2)
      bx(1) = -yh - ys
      bx(2) = -bx(1)
      bx(3) =  yh - ys
      bx(4) = -bx(3)
      cx(1) =  yt + yh
      cx(2) = -yt + yh
      cx(3) = -cx(2)
      cx(4) = -cx(1)
      ay(1) =  xt - xs
      ay(2) = -xt - xs
      ay(3) = -ay(1)
      ay(4) = -ay(2)
      by(1) =  xh + xs
      by(2) = -by(1)
      by(3) = -xh + xs
      by(4) = -by(3)
      cy(1) = -xt - xh
      cy(2) =  xt - xh
      cy(3) = -cy(2)
      cy(4) = -cy(1)
c..3.) compute volume and stabilization h-array
      vol = 4.*xj0
      d11 = d(1)/vol
      d12 = d(2)/vol
      d33 = d(3)/vol
      hy  = vol*3.
      hx  = hy*d(5)
      h44 = hx*(1.-r2*r2/3.)*(a1(1)+a1(2))**2
      h55 = hx*(1.-r1*r1/3.)*(a2(1)+a2(2))**2
      h45 =-(r1*r2/3.)*(hx*(xs*xt+ys*yt)**2+d(6)*hy*(ys*xt-xs*yt)**2)
c..4.) Invert stabilization h-array
      hx  = h44*h55 - h45*h45
      hy  = h55/hx
      h55 = h44/hx
      h45 =-h45/hx
      h44 = hy
c..5.) Compute the current stress parameters
      call pconsd(rr,5,0.0d0)
      do 50 j = 1,4
        hx = cx(j) - r2*ax(j)
        hy = cy(j) - r2*ay(j)
        x1(j) = a1(1)*hx + a1(3)*hy
        x2(j) = a1(2)*hy + a1(3)*hx
        hx = bx(j) - r1*ax(j)
        hy = by(j) - r1*ay(j)
        y1(j) = a2(1)*hx + a2(3)*hy
        y2(j) = a2(2)*hy + a2(3)*hx
        rr(1) = rr(1) + ax(j)*ul(1,j)
        rr(2) = rr(2) + ay(j)*ul(2,j)
        rr(3) = rr(3) + ay(j)*ul(1,j) + ax(j)*ul(2,j)
c...... (stabilization terms)
        rr(4) = rr(4) + x1(j)*ul(1,j) + x2(j)*ul(2,j)
        rr(5) = rr(5) + y1(j)*ul(1,j) + y2(j)*ul(2,j)
50    continue
      beta(1) = d11*rr(1) + d12*rr(2)
      beta(2) = d12*rr(1) + d11*rr(2)
      beta(3) = d33*rr(3)
      beta(4) = (h44*rr(4) + h45*rr(5))*3.
      beta(5) = (h45*rr(4) + h55*rr(5))*3.
c..6.) Form stiffness matrix for 1-pt (constant) terms
      if(isw.eq.3) then
        d11 = d11*d(11)
        d12 = d12*d(11)
        d33 = d33*d(11)
        i1 = 1
        do 110 i = 1,2
          bd11 = ax(i)*d11
          bd12 = ax(i)*d12
          bd13 = ay(i)*d33
          bd21 = ay(i)*d12
          bd22 = ay(i)*d11
          bd23 = ax(i)*d33
          j1   = i1
          do 100 j = i,2
            s(i1  ,j1  ) = bd11*ax(j) + bd13*ay(j)
            s(i1  ,j1+1) = bd12*ay(j) + bd13*ax(j)
            s(i1+1,j1  ) = bd21*ax(j) + bd23*ay(j)
            s(i1+1,j1+1) = bd22*ay(j) + bd23*ax(j)
            j1 = j1 + ndf
100       continue
          i1 = i1 + ndf
110     continue
c..7.) Copy other parts from computed terms
        i1 = ndf + ndf
        do 120 i = 1,i1
        do 120 j = i,i1
          s(i   ,j+i1) =-s(i,j)
          s(i+i1,j+i1) = s(i,j)
120     continue
        j1 = i1 + ndf
        s(2    ,i1+1) =-s(2,1)
        s(ndf+2,j1+1) =-s(ndf+2,ndf+1)
        s(ndf+1,i1+1) =-s(1,ndf+1)
        s(ndf+1,i1+2) =-s(2,ndf+1)
        s(ndf+2,i1+1) =-s(1,ndf+2)
        s(ndf+2,i1+2) =-s(2,ndf+2)
c..8.) Add stabilization matrix
        h44 = h44*d(11)
        h45 = h45*d(11)
        h55 = h55*d(11)
        j1 = 1
        do 210 j = 1,4
          bd11 = h44*x1(j) + h45*y1(j)
          bd12 = h44*x2(j) + h45*y2(j)
          bd21 = h45*x1(j) + h55*y1(j)
          bd22 = h45*x2(j) + h55*y2(j)
          i1 = 1
          do 200 i = 1,j
            s(i1  ,j1  ) = s(i1  ,j1  ) + x1(i)*bd11 + y1(i)*bd21
            s(i1  ,j1+1) = s(i1  ,j1+1) + x1(i)*bd12 + y1(i)*bd22
            s(i1+1,j1  ) = s(i1+1,j1  ) + x2(i)*bd11 + y2(i)*bd21
            s(i1+1,j1+1) = s(i1+1,j1+1) + x2(i)*bd12 + y2(i)*bd22
            i1 = i1 + ndf
200       continue
          j1 = j1 + ndf
210     continue
      endif
c..9.) compute the residual force vector
      if(mod(isw,3).eq.0) then
c...a.) compute the constant part
        do 300 i = 1,5
          beta(i) = beta(i)*d(11)
300     continue
        do 310 i = 1,2
          p(2*i-1) = -(ax(i)*beta(1) + ay(i)*beta(3))
          p(2*i+3) = -p(2*i-1)
          p(2*i  ) = -(ay(i)*beta(2) + ax(i)*beta(3))
          p(2*i+4) = -p(2*i  )
310     continue
c...b.) compute the stabilization part
        do 320 i = 1,4
          p(2*i-1) = p(2*i-1) - (x1(i)*beta(4) + y1(i)*beta(5))/3.0
          p(2*i  ) = p(2*i  ) - (x2(i)*beta(4) + y2(i)*beta(5))/3.0
320     continue
      endif
      return
      end
c
      subroutine stcn03(ix,d,ssa,tta,dt,st,nnp)
      implicit real*8 (a-h,o-z)
c     integer*2 ix(1)
      integer   ix(1)
      real*8 dt(nnp),st(nnp,1),ss(4),tt(4)
      real*8 d(1)
      common /elcom3/ xs,xt,xh,ys,yt,yh,xj0,xj1,xj2,a1(3),a2(3),beta(5)
      data ss/-1.0,1.0,1.0,-1.0/,tt/-1.0,-1.0,1.0,1.0/
c.... compute stress projections
      do 200 jj = 1,4
        ll = abs(ix(jj))
        if(ll.gt.0) then
c.... compute weighted stresses at nodes
          xsj = xj0 + ss(jj)*xj1 + tt(jj)*xj2
          ssg = (ss(jj) - ssa)*beta(5)
          ttg = (tt(jj) - tta)*beta(4)
          dt(ll)   = dt(ll)   + xsj
          sig1     = (beta(1) + a1(1)*ttg + a2(1)*ssg)*xsj
          sig2     = (beta(2) + a1(2)*ttg + a2(2)*ssg)*xsj
          sig3     = (beta(3) + a1(3)*ttg + a2(3)*ssg)*xsj
          st(ll,1) = st(ll,1) + sig1
          st(ll,2) = st(ll,2) + sig3
          st(ll,3) = st(ll,3) + sig2
          st(ll,4) = st(ll,4) + d(13)*(sig1+sig2)
        endif
200   continue
      return
      end
