subroutine elmt05(d,ul,xl,ix,tl,s,p,ndf,ndm,nst,isw)
      use shapefunctions

      implicit none
      include 'cdata.h'
      include 'hdata.h'
      include 'eldata.h'
      include 'iofile.h'
      include 'eltran.h'
      include 'comblk.h'
      include 'elplot.h'
      
      real*8  ul(ndf,nel),xl(ndm,nel),s(nst,nst),p(nst),d(2)
      real*8  cel(6,6), elas, tl(*) 
      real*8  sig(6), eps(6),dsig(6),deps(6)
      real*8  gp(3), u(24)
      real*8  A(6,24),B(6,24), detJ, Bt(24,6)
      real*8  K(24,24),shp(4,8) 
      integer ix(*), isw, i, j,ndf,ndm,nst

      save

      if(isw.eq.1) then
      call inpt05(d)
      write(*,3000) d(1),d(2)
      write(iow,3000) d(1),d(2)
           
      else if (isw .eq. 2) then
      call ckbrk8 ( n, ix, xl, ndm, 8, shp )

      else if((isw.eq.3).or.(isw.eq.6).or.(isw.eq.4).or.(isw.eq.8))then

      elas=d(1)/((1.00D0+d(2))*(1.00D0-2*d(2)))
      cel=0.00d0
      do i=1,3
      cel(i,i)=1.00D0-d(2)
      cel(i+3,i+3)=(1.00D0-2*d(2))/2.00D0
      enddo
      cel(1,2)=d(2)
      cel(1,3)=d(2)
      cel(2,1)=d(2)
      cel(2,3)=d(2)
      cel(3,1)=d(2)
      cel(3,2)=d(2)
      cel=elas*cel

      do i=0,7
      do j=1,3
      u(3*i+j)=ul(j,i+1)
      end do 
      end do
c--------------------------------------------------------      
      open(unit=101,file="displacements.txt",status="old")
      write(101,*) "Displacements of element :",n
      write(101,'(8E12.4)') ul(1,:)
      write(101,'(8E12.4)') ul(2,:)
      write(101,'(8E12.4)') ul(3,:)
c--------------------------------------------------------   
      open(unit=102,file="constitutive.txt",status="old")
      write(102,*) "Constitutive matrix of element :",n
      do i=1,6
      write(102,'(6F12.4)') cel(i,:)
      end do
c--------------------------------------------------------
      gp=(/0.00D0,0.00D0,0.00D0/)
      call B3D8N(xl,gp,B,detJ)   
c--------------------------------------------------------
      open(unit=103,file="bmat.txt",status="old")
      write(103,*) "B matrix of element :",n
      do i=1,6
      write(103,'(24F12.4)') B(i,:)
      end do
c--------------------------------------------------------      
      Bt=transpose(B)
      A=matmul(cel,B)*detJ
      K=matmul(Bt,A)
c--------------------------------------------------------      
      open(unit=104,file="kmat.txt",status="old")
      write(104,*) "K matrix of element :",n
      do i=1,24
      write(104,'(24F12.4)') K(i,:)
      end do
c-------------------------------------------------------- 
      ! STIFFNESS
      do i=1,24
      do j=1,24
            s(i,j)=K(i,j)
      end do
      end do
      ! STRAIN
      do i=1,6
      eps(i)=0.0000D0
      do j=1,24
      eps(i)=eps(i)+B(i,j)*u(j)
      end do
      tt(i)=eps(i)
      end do
      ! STRESS
      do i=1,6
      sig(i)=0.0000D0
      do j=1,6
      sig(i)=sig(i)+cel(i,j)*eps(j)
      end do
      tt(i+6)=sig(i)
      end do
      ! RESIDUAL
      do i=1,24
      p(i)=0.0000D0
      do j=1,6
      p(i)=p(i)+Bt(i,j)*sig(j)*detJ
      end do
c      p(i)=p(i)-q(n,i)
c      q(n,i)=p(i)
      end do

c--------------------------------------------------------      
      open(unit=105,file="element.txt",status="old")
      write(105,*) "Strain at the centroid of element :",n
      write(105,'(6F12.4)') eps
      write(105,*) "Stress at the centroid of element :",n
      write(105,'(6F12.4)') sig
      write(105,*) "Residual at the centroid of element :",n
      write(105,'(24F12.4)') p
c--------------------------------------------------------       
      elseif (isw.eq.12) then  !update element history data
      do i=0,23
      hr(nh1+i) = u(i+1)
      enddo
      do i=24,29
      hr(nh1+i) = eps(i-23)
      hr(nh1+i+6) = sig(i-23)
      enddo

      
      elseif(isw.eq.14) then
      do i=0,23
      hr(nh1+i)=0.00D0
      end do
      endif      

         
3000  format(/5x,'Elmt 4: Isotropic Linear Elastic Analysis'
     & ' with one-point numerical integration - email@mubeen.info'/// 
     &       10x,'Elastic Modulus   = ',e14.5/
     &       10x,'Poisson ratio     = ',e14.5/x)
      end

      subroutine inpt05(d)
      implicit none
       include 'iofile.h'
       real*8   td, d(2)
       logical pcomp, errck, tinput       
       character*9 text
       save
       text = 'start'
       do while (.not.pcomp(text,'    ',4))
          errck = tinput(text,1,td,1)          
          if (pcomp(text,'Emod',4)) then
            d(1) = td
          elseif (pcomp(text,'neu',3)) then
            d(2) = td
          endif
      end do       

      end SUBROUTINE