C-----==================================================================
C-----ONE-DIMENSIONAL HEAT EQAUTION ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
C-----GALERKIN FINITE ELEMENT METHOD WITH EULER IMPLICIT ]]]]]]]]]]]]]]]
C-----GAUSS QUADRATURE IS USED FOR INTEGRATION ]]]]]]]]]]]]]]]]]]]]]]]]]
C-----==================================================================
      IMPLICIT NONE
      INTEGER NP,NE,ITE
      PARAMETER (NP=21,NE=NP-1,ITE=1000)
      INTEGER NEC(NE,2)
      REAL*8 X1(NP)
      REAL*8 A(NP,NP),B(NP),T(NP),E(NP)
      REAL*8 TIME,DELX,DELT,PI,BETA,ERROR,SUM
      REAL*8 N1,N2,DN1,DN2
      REAL*8 XL(2),XW(2)
      INTEGER I,L,M,N,IT,IPIV(2*NP),INFO


C-----TIME STEP ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      DELT=0.1D00/DBLE(ITE)
      PI=4D00*ATAN(1D00)


C-----GAUSS WEIGHTS AND LOCATIONS ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      XL(1)=-1D00/SQRT(3D00)
      XL(2)= 1D00/SQRT(3D00)
      XW(1)= 1D00
      XW(2)= 1D00


C-----READ MESH DATA ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      CALL MESH(X1,NEC)


C-----INITIAL CONDITIONS ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      T(:)=0.0D00
      T(NP)=1.D00



C-----==================================================================
C-----TIME INTEGRATION ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
C-----==================================================================
      TIME=0.0D00
      DO 100 IT=1,ITE
      TIME=TIME+DELT
      WRITE(6,*) "TIME=",TIME

C-----CONSTRUCT DISCRETIZED EQUATIONS ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      A(:,:)=0.0D00
      B(:)=0.000D00
      DO I=1,NE
      DO L=1,2

C-----SHAPE FUNCTIONS N1 AND N2 AND THEIR DERIVATIONS ]]]]]]]]]]]]]]]]]]
      N1=0.5D00*(1D00-XL(L)) !SHAPE FUNCTION N1
      N2=0.5D00*(1D00+XL(L)) !SHAPE FUNCTION N2
      DN1=-0.5D00            !DERIVATION OF SHAPE FUNCTION N1
      DN2= 0.5D00            !DERIVATION OF SHAPE FUNCTION N2
      DELX=X1(NEC(I,2))-X1(NEC(I,1))

C-----MASS MATRIX PART ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      M=NEC(I,1)
      N=NEC(I,2)
      A(M,M)=A(M,M)+N1*N1*(DELX/2D00)*XW(L)/DELT
      A(M,N)=A(M,N)+N1*N2*(DELX/2D00)*XW(L)/DELT
      B(M)=B(M)+N1*N1*(DELX/2D00)*XW(L)/DELT*T(M)
     *         +N1*N2*(DELX/2D00)*XW(L)/DELT*T(N)

      A(N,M)=A(N,M)+N2*N1*(DELX/2D00)*XW(L)/DELT
      A(N,N)=A(N,N)+N2*N2*(DELX/2D00)*XW(L)/DELT
      B(N)=B(N)+N2*N1*(DELX/2D00)*XW(L)/DELT*T(M)
     *         +N2*N2*(DELX/2D00)*XW(L)/DELT*T(N)

C-----DIFFUSION PART ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      A(M,M)=A(M,M)+DN1*DN1*(2D00/DELX)*XW(L)
      A(M,N)=A(M,N)+DN1*DN2*(2D00/DELX)*XW(L)
      
      A(N,M)=A(N,M)+DN2*DN1*(2D00/DELX)*XW(L)
      A(N,N)=A(N,N)+DN2*DN2*(2D00/DELX)*XW(L)
      END DO
      END DO


C-----IMPOSE DRICHLET BOUNDARY CONDITIONS ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      A(1,:)=0.D00
      A(1,1)=1.D00
      B(1)=0.0D0
      A(NP,:)=0.0D00
      A(NP,NP)=1.D00
      B(NP)=1.0D0


C-----SOLVE SQUARE SYSTEM (LAPACK SUBROUTINE) ]]]]]]]]]]]]]]]]]]]]]]]]]]
C     CALL DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )
      CALL DGESV( NP, 1, A, NP, IPIV, B, NP, INFO ) 
      IF (INFO.NE.0) WRITE(6,*) "INFO=",INFO
      T(:)=B(:)
  100 CONTINUE
      
      
C-----==================================================================
C-----ANALYTICAL SOLUTION ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
C-----==================================================================
      WRITE(6,*)
      WRITE(6,*) "NUMERICAL AND ANALYTICAL SOLUTION..."
      DO I=1,NP
      SUM=0.0D00
      DO N=1,20
      SUM=SUM+2D00*(-1D00)**N/(DBLE(N)*PI) * 
     *    EXP(-DBLE(N)**2*PI**2*TIME)*SIN(DBLE(N)*PI*X1(I))
      END DO
      E(I)=X1(I)+SUM
      WRITE(6,11) X1(I),T(I),E(I)
      END DO

      
C-----COMPUTE ERROR ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      WRITE(6,*)
      ERROR=0.0D00
      DO I=1,NP
      ERROR=ERROR+(T(I)-E(I))**2
      END DO
      ERROR=ERROR/DSQRT(DBLE(NP))
      WRITE(6,*) "DELX=",DELX,"ERROR=",ERROR


C-----==================================================================
C-----TECPLOT FILE FORMAT ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
C-----==================================================================
      OPEN(10,FILE='ANIME.DAT',FORM='FORMATTED')
      WRITE(10,*) 'TITLE="HEAT EQUATION"'
      WRITE(10,*) 'VARIABLES="X","NUMERIC","ANALYTIC"'
      WRITE(10,*) 'ZONE  F=FEPOINT N=',NP,'E=',NE,'ET=LINESEG'
      DO I=1,NP
      WRITE(10,*) X1(I),T(I),E(I)
      END DO
      DO I=1,NE
      WRITE(10,*) NEC(I,1),NEC(I,2)
      END DO
      CLOSE(10)

   11 FORMAT(D18.12,2X,D18.12,2X,D18.12)
      END





C     ==================================================================
C     [][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]
C     ==================================================================
      SUBROUTINE MESH(X1,NEC)
      IMPLICIT NONE
      INTEGER NP,NE
      PARAMETER (NP=21,NE=NP-1)
      INTEGER NEC(NE,2)
      INTEGER I
      REAL*8  X1(NP)

C-----MESH NODES ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      DO I=1,NP
      X1(I)=(DBLE(I-1)/DBLE(NP-1))
      END DO

C-----NODE ELEMENT CONNECTIVITY ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      DO I=1,NE
      NEC(I,1)=I
      NEC(I,2)=I+1
      END DO
      RETURN
      END