PROGRAM FEM8Q
C=======================================================================
C  2-DIM FEM PROGRAM FOR SOLVING LAPLACE EQUATION: DP/DX + DP/DY = 0.
C     USING UPPER HALF BANDED MATRIX AND 8-NODED ISO-PARAMETERIC
C SEQUENCE: (-1,-1),(0,-1),(+1,-1),(+1,0),(+1,+1),(0,+1),(-1,+1),(-1,0)
C                  7-----6-----5
C                  |     |     |
C                  |     |     |
C                  8     +     4
C                  |     |     |
C                  |     |     |
C                  1-----2-----3
C   ORIGINAL:1984 EIJI FUKUMORI BUFFALO NY & REVISED: 1994 ACHI
C   FEB. 15, 2013
C=======================================================================
      INCLUDE 'PARAM.DAT'
      IMPLICIT REAL*8 ( A-H , O-Z )
CCCC      PARAMETER (ND=8,MXE=3200,MXN=3400,MXB=2200,INTEPT=3,MXW=1040)
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN), IBND(MXB),
     * ITYPE(MXB), BVALUE(MXB),RK(MXN,MXW), RHS(MXN),SK(ND,ND), 
     * X(ND), Y(ND),BX(2,ND),SAI(INTEPT), W(INTEPT),F1(ND),F2(ND),
     * BPP(2,ND,INTEPT,INTEPT)
      CHARACTER INPFILE*14
      WRITE (*,*)' FEM8Q.FOR SOLVER'
C=======================================================================
      CALL GRULE ( INTEPT, SAI, W   )
      CALL DERIV ( ND, INTEPT, F1, F2, SAI, BPP )
C=======================================================================
      CALL INPUT ( INPFILE,ND,MXE,MXN,MXB, NE,NNODE,NB, 
     * EXX,EXY,EYY,NODEX, XCOORD, YCOORD, IBND, ITYPE, BVALUE )
C=======================================================================
      CALL BANDWID ( MXE, ND, NE,NODEX, NBW )
      IF ( NBW .GT. MXW ) STOP 'Error ... NBW > MXW'
C=======================================================================
      CALL GSM ( INTEPT,ND,MXE,MXN,MXW,NE,NNODE,NBW,
     * EXX,EXY,EYY,BPP,W,BX,SK,NODEX,XCOORD,YCOORD,RK )
C=======================================================================
      CALL FORMQ (MXN,MXB,MXW,NNODE,NB,NBW,RK,RHS,ITYPE,BVALUE,IBND)
C=======================================================================
      CALL SYSTEMQ ( MXN, MXW, NNODE, NBW, RK, RHS )
C=======================================================================
      CALL RESULT ( IW, INPFILE, ND,MXE,MXN,MXB,NE,NNODE,NB,
     * EXX,EXY,EYY,NODEX, XCOORD, YCOORD, IBND, ITYPE, BVALUE, RHS )
C=======================================================================
      CALL RESULT0 (IW,ND,MXE,MXN,NE,F1, NODEX, XCOORD, YCOORD,RHS )
C=======================================================================
      CLOSE (IW)
C=======================================================================
      CALL SOL ( MXN,NNODE,RHS )
C=======================================================================
      STOP 'NORMAL TERMINATION'
      END
C
C
      SUBROUTINE SOL ( MXN,NNODE,U )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION U(MXN)
      CHARACTER OUTFILE*12, EXFILE*3
      LOGICAL YES
C========> FILENAME SOLUTION.BIN
      OUTFILE = 'SOLUTION.BIN'
      OPEN ( 1, FILE=OUTFILE,STATUS='UNKNOWN',FORM='UNFORMATTED' )
      WRITE (1) ( U(I) , I = 1 , NNODE )
      CLOSE (1)
      RETURN
      END
C
C
      SUBROUTINE GSM ( INTEPT,ND,MXE,MXN,MXW,NE,NNODE,NBW,
     * EXX,EXY,EYY,BPP,W,BX,SK,NODEX,XCOORD,YCOORD,RK )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),RK(MXN,MXW),
     * SK(ND,ND), BX(2,ND),W(INTEPT),BPP(2,ND,INTEPT,INTEPT)
COMPUTATION OF LOCAL STIFFNESS MATRICES AND THE ASSEMBLY OF 
C  THE GLOBAL STIFFNESS MATRIX,RK.
C---- RESET
      DO I = 1 , NNODE
      DO J = 1 , NBW
      RK (I,J) = 0.D0
      END DO
      END DO
C------ INTEGRATION
      DO IEL = 1 ,NE
C---------RESET SK(I,J)
      DO I = 1 , ND
      DO J = 1 , ND
      SK (I,J) = 0.D0
      END DO
      END DO
      DO K = 1 , INTEPT
      DO L = 1 , INTEPT
      WEIGHT = W(K) * W(L)
C---------COMPUTATION OF B(1,J) = DN(J) / DX
C---------COMPUTATION OF B(2,J) = DN(J) / DY
C---------COMPUTATION OF JACOBIAN MATRIX [YAC].
      YAC11 = 0.D0
      YAC12 = 0.D0
      YAC21 = 0.D0
      YAC22 = 0.D0
      DO I = 1 , ND
      NODE = NODEX(IEL,I)
      YAC11 = YAC11 + BPP(1,I,K,L) * XCOORD(NODE)
      YAC12 = YAC12 + BPP(1,I,K,L) * YCOORD(NODE)
      YAC21 = YAC21 + BPP(2,I,K,L) * XCOORD(NODE)
      YAC22 = YAC22 + BPP(2,I,K,L) * YCOORD(NODE)
      END DO
C---------COMPUTATION OF THE INVERSE OF THE JACOBIAN MATRIX.
      DETJ = YAC11 * YAC22 - YAC12 * YAC21
      Y11 =   YAC22 / DETJ
      Y12 = - YAC12 / DETJ
      Y21 = - YAC21 / DETJ
      Y22 =   YAC11 / DETJ
C---------COMPUTATION OF [B] MATRIX.
      DO J = 1 , ND
      BX(1,J) = Y11 * BPP(1,J,K,L) + Y12 * BPP(2,J,K,L)
      BX(2,J) = Y21 * BPP(1,J,K,L) + Y22 * BPP(2,J,K,L)
      END DO
      BETA = WEIGHT * DETJ
      DO I = 1 ,  ND
      DO J = 1 ,  ND
      SK(I,J) = SK(I,J) + ((BX(1,I)*EXX+BX(2,I)*EXY)*BX(1,J)
     *                  +  (BX(1,I)*EXY+BX(2,I)*EYY)*BX(2,J))*BETA
      END DO
      END DO
C--------
      END DO
      END DO
C----ASSEMBLY
      DO K = 1 , ND
      I = NODEX(IEL,K)
      DO L = 1 , ND
      J = NODEX(IEL,L) - I + 1
      IF ( J .GT. 0 ) RK (I,J) = RK(I,J) + SK(K,L)
      END DO
      END DO
      END DO
      RETURN
      END
C
C
      SUBROUTINE BANDWID ( MXE, ND, NE, NODEX , NBW  )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND)
      NBW = 0
      DO 30 I = 1 , NE
      DO 20 J = 1 , ND-1
      DO 10 K = J+1 , ND
      NBW = MAX0(NBW,IABS(NODEX(I,J)-NODEX(I,K)))
   10 CONTINUE
   20 CONTINUE
   30 CONTINUE
      NBW = NBW + 1
      write(*,*) ' bandwidth =', nbw
      RETURN
      END
C
C
      SUBROUTINE SYSTEMQ ( MXN, MXW, NUMNP, MBAND, A, B )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION A(MXN,MXW) , B(MXN)
C---------- ELIMINATION ------------------
      DO N = 1 , NUMNP 
      DO L = 2 , MBAND 
      C = A(N,L) / A(N,1) 
      I = N + L - 1 
      IF ( I .LE. NUMNP ) THEN
      J = 0 
      DO K = L , MBAND 
      J = J + 1 
      A(I,J) = A(I,J) - C * A(N,K)
      END DO
      A(N,L) = C
      B(I) = B(I) - A(N,L) * B(N)
      ENDIF
      END DO
      B(N) = B(N) / A(N,1)
      END DO
C---------- BACKSUBSTITUTION -------------
      N = NUMNP
      DO WHILE ( N .GT. 0 )
      DO K = 2 , MBAND
      L = N + K - 1 
      IF ( L .LE. NUMNP ) B(N) = B(N)-A(N,K)*B(L) 
      END DO
      N = N - 1 
      END DO
      RETURN
      END 
C
C
      SUBROUTINE FORMQ (MXN,MXB,MXW,NNODE,NB,NBW,A,RHS,IBC,BV,IBND)
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION RHS(MXN),A(MXN,MXW),IBC(MXB),BV(MXB),IBND(MXB)
C  BOUNDARY CONDITIONS ARE COUPLED HERE.
C  ITYPE(I) = 1 ---> DIRICHLET, ITYPE(I) = 2 ---> NEUMANN
C----- REFORMATION OF VECTOR RHS DUE TO FIRST KIND B.C. 
C------ CLEAR RIGHT HAND SIDE.
      DO I = 1 , NNODE
      RHS (I) = 0.D0
      END DO
C------ DIRICHLET B.C.'S DO TO RHS.
      DO K = 1 , NB
        IF ( IBC(K) .EQ. 1 ) THEN
          I = IBND(K) 
          DO J = 2 , NBW 
              I = I - 1 
              IF ( I.GT. 0 ) THEN
              RHS(I) = RHS(I) - BV(K) * A(I,J)
              END IF
          END DO
            I = IBND(K) 
          DO J = 2 , NBW 
              I = I + 1 
              IF ( I .LE. NNODE ) THEN
              RHS(I) = RHS(I) - BV(K) * A(IBND(K),J)
              END IF
          END DO
        END IF
      END DO
C-----REFORMATION OF MATRIX A.
      DO K = 1 , NB
      I = IBND (K)
      IF ( IBC(K) .EQ. 1 ) THEN
      RHS(I) = BV(K) 
      A(I,1) = 1.D0
      DO J = 2 , NBW 
      L = I - J + 1 
      A(I,J) = 0.D0
      IF ( L  .GT. 0 ) THEN
      A( L ,J) = 0.D0
      END IF
      END DO
      ELSE
      RHS (I) = RHS(I) - BV(K)
      END IF
      END DO
      RETURN
      END
C
C
      SUBROUTINE GRULE ( N , SAI , W )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(N) , W(N)
      IF ( N .LT. 2 ) STOP'N<2'
      IF ( N .GT. 6 ) STOP'N>6'
      IF ( N .EQ. 2 ) THEN
      SAI(1) = DSQRT(3.D0)/3.D0
      W(1) = 1.D0
      SAI(2) = - SAI(1)
      W(2) = W(1)
      RETURN
      END IF
      IF ( N .EQ. 3 ) THEN
      SAI(1) = DSQRT(15.D0)/5.D0
      SAI(2) = 0.D0
      W(1) = 5.D0/ 9.D0
      W(2) = 8.D0/ 9.D0
      SAI(3) = - SAI(1)
      W(3) = W(1)
      RETURN
      END IF
      IF ( N .EQ. 4 ) THEN
      SAI(1) = 0.33998104358485D0
      SAI(2) = 0.86113631159405D0
        W(1) = 0.65214515486254D0
        W(2) = 0.34785484513745D0
      SAI(3) = -SAI(2)
      SAI(4) = -SAI(1)
      W(3) = W(2)
      W(4) = W(1)
      RETURN
      END IF
      IF ( N .EQ. 5 ) THEN
      SAI(1) = 0.90617984593866D0
      SAI(2) = 0.53846931010568D0
      SAI(3) = 0.D0
        W(1) = 0.23692688505619D0
        W(2) = 0.47862867049937D0
        W(3) = 5.12D0 / 9.D0
      SAI(4) = -SAI(2)
      SAI(5) = -SAI(1)
      W(4) = W(2)
      W(5) = W(1)
      RETURN
      END IF
      IF ( N .EQ. 6 ) THEN
      SAI(1) = 0.23861918608320D0
      SAI(2) = 0.66120938646626D0
      SAI(3) = 0.93246951420315D0
        W(1) = 0.46791393457269D0
        W(2) = 0.36076157304814D0
        W(3) = 0.17132449237917D0
      SAI(4) = -SAI(1)
      SAI(5) = -SAI(2)
      SAI(6) = -SAI(3)
        W(4) = W(1)
        W(5) = W(2)
        W(6) = W(3)
      END IF
      RETURN
      END
C
C
      SUBROUTINE DERIV ( ND, INTEPT, F1, F2, SAI, BPP )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(INTEPT),BPP(2,ND,INTEPT,INTEPT),F1(ND),F2(ND)
      DX = 0.5D0
      DO 40 K = 1 , INTEPT
      E1 = SAI (K)
      DO 30 L = 1 , INTEPT
      E2 = SAI (L)
COMPUTATION OF BP(J) = DN(J) / DETA1
      CALL ISOPARA ( ND, E1-DX , E2 , F1 )
      CALL ISOPARA ( ND, E1+DX , E2 , F2 )
      DO I = 1 , ND
      BPP(1,I,K,L) = ( F2(I) - F1(I) ) / ( 2.D0* DX )
      END DO
COMPUTATION OF BP(J)= DN(J)/DETA2
      CALL ISOPARA ( ND, E1 , E2-DX , F1 )
      CALL ISOPARA ( ND, E1 , E2+DX , F2 )
      DO I = 1 , ND
      BPP(2,I,K,L) = ( F2(I) - F1(I) ) / ( 2.D0* DX )
      END DO
   30 CONTINUE
   40 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE ISOPARA ( ND, E1 , E2 , F )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION F(ND)
      F(1) = -0.25D0*(1.D0- E1   )*(1.D0- E2   )*(E1+E2+1.D0)
      F(2) =  0.50D0*(1.D0- E1*E1)*(1.D0- E2   )
      F(3) =  0.25D0*(1.D0+ E1   )*(1.D0- E2   )*(E1-E2-1.D0)
      F(4) =  0.50D0*(1.D0+ E1   )*(1.D0- E2*E2)
      F(5) =  0.25D0*(1.D0+ E1   )*(1.D0+ E2   )*(E1+E2-1.D0)
      F(6) =  0.50D0*(1.D0- E1*E1)*(1.D0+ E2   )
      F(7) = -0.25D0*(1.D0- E1   )*(1.D0+ E2   )*(E1-E2+1.D0)
      F(8) =  0.50D0*(1.D0- E1   )*(1.D0- E2*E2)
      RETURN
      END
C
C
      SUBROUTINE INPUT ( INPFILE,ND,MXE,MXN,MXB, NE,NNODE,NB, 
     * EXX,EXY,EYY,NODEX, XCOORD, YCOORD, IBND, ITYPE, BVALUE )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBND(MXB),
     * ITYPE(MXB), BVALUE(MXB)
      LOGICAL YES
      CHARACTER INPFILE*14
      IR = 1
      IF ( ND .LE.  2 ) STOP 'NOT ACCEPTABLE ND'
      IF ( ND .EQ.  3 ) INPFILE = 'FEM03INPUT.DAT'
      IF ( ND .EQ.  4 ) INPFILE = 'FEM04INPUT.DAT'
      IF ( ND .EQ.  8 ) INPFILE = 'FEM08INPUT.DAT'
      IF ( ND .EQ.  9 ) INPFILE = 'FEM09INPUT.DAT'
      IF ( ND .EQ. 12 ) INPFILE = 'FEM12INPUT.DAT'
      INQUIRE ( FILE=INPFILE, EXIST=YES )
      IF ( .NOT.YES ) STOP' INPUT FILE DOES NOT EXIST'
      OPEN ( IR, FILE = INPFILE, STATUS = 'OLD' )
      READ (IR,*) EXX, EXY, EYY
      READ (IR,*) NE
      IF ( NE .GT. MXE ) STOP 'NE > MXE'
      READ (IR,*) (IEL,(NODEX(IEL,J),J=1,ND), I=1,NE)
      READ (IR,*) NNODE
      IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
      READ (IR,*) (NODE,XCOORD(NODE),YCOORD(NODE),J=1,NNODE)
      READ (IR,*) NB
      IF ( NB .GT. MXB ) STOP 'NB > MXB'
      READ (IR,*) (IBND(I), ITYPE(I), BVALUE(I),I=1,NB)
      RETURN
      END
C
C
      SUBROUTINE RESULT ( IW, INPFILE, ND,MXE,MXN,MXB,NE,NNODE,NB,
     * EXX,EXY,EYY,NODEX, XCOORD, YCOORD, IBND, ITYPE, BVALUE, RHS )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND), XCOORD(MXN), YCOORD(MXN), RHS(MXN)
      DIMENSION ITYPE (MXB), IBND(MXB), BVALUE(MXB)
      CHARACTER OUTFILE*12, EXFILE*3, INPFILE*14
      LOGICAL YES
C=================== ECHO AND RESULT PRINTS ======================
C--------- FILE INQUIRERY --------
      OUTFILE = 'SOLUTION.FEM'
      IW = 1
      INQUIRE ( FILE=OUTFILE, EXIST=YES )
      IF ( YES ) EXFILE = 'OLD'
      IF ( .NOT. YES ) EXFILE = 'NEW'
      OPEN ( IW, FILE = OUTFILE, STATUS = EXFILE )
C--------- CRT DUMP ------------
      WRITE (*,*)  ' INPUT  FILE: ', INPFILE
      WRITE (*,*)  ' OUTPUT FILE: ', OUTFILE
C--------- ECHO PRINT ------------
      WRITE (IW,*)  ' NAME OF INPUT FILE: ', INPFILE
      WRITE (IW,'(1X,78(1H-))')
      WRITE (IW,*) ' PROPERTY OF DOMAIN:'
      WRITE (IW,*) '   EXX =',EXX
      WRITE (IW,*) '   EXY =',EXY
      WRITE (IW,*) '   EYY =',EYY
C--------- ELEMENT ---------------
      WRITE (IW,'(1X,78(1H-))')
      WRITE (IW,*) ' TYPE OF ELEMENT:'
      WRITE (IW,*) '   NUMBER OF NODES AT EACH ELEMENT =', ND
C--------- DISCRETIZATION ---------------
      WRITE (IW,'(1X,78(1H-))')
      WRITE (IW,*) ' DISCRETIZATION OF DOMAIN INTO ELEMENTS:'
      WRITE (IW,*) '   NUMBER OF ELEMENTS(NE) =', NE
      WRITE (IW,'(10X,11HELEMENT NO ,9(2X,1H(,I1,1H)))') (I,I=1,ND)
      DO I = 1 , NE
      WRITE (IW, '(10X ,I10, 9I5)') I,(NODEX(I,J),J=1,ND)
      END DO
C--------- NODAL POINT COORDINATS -----
      WRITE (IW,'(1X,78(1H-))')
      WRITE (IW,*) ' COORDINATES OF NODAL POINTS:'
      WRITE (IW,*) '   NUMBER OF NODAL POINTS(NNODE) =', NNODE
      WRITE (IW,*) '   (I=NODAL POINT, X & Y = X- & Y-COORDINATES)'
      DO I = 1 , NNODE
      WRITE (IW,*) '   (I,X,Y)=',I,XCOORD(I), YCOORD(I)
      END DO
C--------- BOUNDARY CONDITIONS AND VALUES -------
      WRITE (IW,'(1X,78(1H-))')
      WRITE (IW,*) ' BOUNDARY CONDITIONS:'
      WRITE (IW,*) '  NUMBER OF BOUNDARY NODES =',NB
      WRITE (IW,*) '  (I=NODAL POINT, T=B.C. B=B.V.)'
      WRITE (IW,*) '  (DIRICHLET IF B.C.=1, NUEMANN IF B.C.=2)'
      DO I = 1 , NB
      WRITE (IW,*) '  (I,T,B)=',IBND(I),ITYPE(I),BVALUE(I)
      END DO
C-------- PRINT RESULT ---------
      WRITE (IW,'(1X,78(1H-))')
      WRITE (IW,*) ' RESULTS:'
      WRITE (IW,*) '   UNKNOWN VALUES AT NODEL POINTS'
      WRITE (IW,*) '   (I=NODAL POINTS, U=UNKNOWN VALUE)'
      DO I = 1 , NNODE
      WRITE (IW,*) '   (I,U)=', I,RHS(I)
      END DO
      RETURN
      END
C
C
      SUBROUTINE RESULT0 (IW,ND,MXE,MXN,NE,F,NODEX,XCOORD,YCOORD,RHS)
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND), XCOORD(MXN), YCOORD(MXN), RHS(MXN)
      DIMENSION F(ND)
C===EVALUATION AND PRINT OF  POTENTIAL AT CENTER OF ELEMENT ===
      E1 = 0.D0
      E2 = 0.D0
      CALL ISOPARA ( ND, E1 , E2 , F )
      WRITE (IW,'(1X,78(1H-))')
      WRITE (IW,*) ' POTENTIAL AT CENTER OF ELEMENT'
      DO I = 1 , NE
      X = 0.D0
      Y = 0.D0
      P = 0.D0
      DO J = 1 , ND
      X = X + F(J)*XCOORD(NODEX(I,J))
      Y = Y + F(J)*YCOORD(NODEX(I,J))
      P = P + F(J)*RHS(NODEX(I,J))
      END DO
      WRITE (IW,*) ' X,Y,P=', X, Y, P
      END DO
      RETURN
      END