      PROGRAM FEM3
C=======================================================================
C     A BASIC FEM 2-DIM PROGRAM USING 3-NODED TRIANGULAR ELEMENTS
C     FOR LAPLACE EQUATION.             SPRING, 1980 EIJI FUKUMORI
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=3, MXE=32, MXN=25, MXB=20 )
      DIMENSION SK(ND,ND),X(ND),Y(ND),A(ND),B(ND),IBND(MXB),ITYPE(MXB),
     * BVALUE(MXB),NODEX(MXE,ND), XCOORD(MXN),YCOORD(MXN), RK(MXN,MXN),
     * RHS(MXN)
      CHARACTER INPFILE*12
C=======================================================================
      CALL INPUT ( INPFILE, ND, MXE,MXN,MXB, NE,NNODE,NB, EXX,EXY,EYY,
     * NODEX, XCOORD, YCOORD, IBND, ITYPE, BVALUE )
C=======================================================================
      CALL GLBSTIFF ( ND, MXE, MXN, NE, NNODE, EXX,EXY,EYY,
     * X,Y,A, B, SK,NODEX, XCOORD,YCOORD, RK )
C=======================================================================
      CALL FORM ( MXN, MXB, NNODE, NB,IBND,ITYPE,BVALUE, RK, RHS )
C=======================================================================
      CALL SYSTEM ( MXN, NNODE, RK, RHS )
C=======================================================================
      CALL RESULT ( INPFILE, ND,MXE,MXN,MXB,NE,NNODE,NB,EXX,EXY,EYY,
     * NODEX, XCOORD, YCOORD, IBND, ITYPE, BVALUE, RHS )
      STOP
      END
C
C
      SUBROUTINE GLBSTIFF ( ND, MXE, MXN, NE, NNODE, EXX,EXY,EYY,
     * X,Y,A, B, SK,NODEX, XCOORD,YCOORD, RK )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN), RK(MXN,MXN),
     * SK(ND,ND),  X(ND), Y(ND), A(ND), B(ND)
C  COMPUTATION OF LOCAL STIFFNESS MATRICES AND THE ASSEMBLY 
C  OF THE GLOBAL STIFFNESS MATRIX,RK.
C---- RESET
      DO I = 1 , NNODE
      DO J = 1 , NNODE
      RK(I,J) = 0.
      END DO
      END DO
C
      DO IEL = 1 ,NE
      DO I = 1 , ND
      X (I) = XCOORD (NODEX(IEL,I))
      Y (I) = YCOORD (NODEX(IEL,I))
      END DO
      CALL STIFMTX ( ND,  X, Y, A, B, EXX, EXY, EYY, SK  )
C----ASSEMBLY
      DO K = 1 , ND
      I = NODEX(IEL,K)
      DO L = 1 , ND
      J = NODEX(IEL,L)
      RK (I,J) = RK(I,J) + SK(K,L)
      END DO
      END DO
      END DO
      RETURN
      END
C
C
      SUBROUTINE STIFMTX (  ND, X, Y, A, B, EXX, EXY, EYY, STIFF)
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION X(ND), Y(ND), STIFF(ND,ND), A(ND), B(ND)
C   COMPUTATION OF DN/DX
      B(1) = Y(2) - Y(3)
      B(2) = Y(3) - Y(1)
      B(3) = Y(1) - Y(2)
C   COMPUTATION OF DN/DY
      A(1) = X(3) - X(2)
      A(2) = X(1) - X(3)
      A(3) = X(2) - X(1)
C   COMPUTATION OF AREA
      AREA = ( A(3)*B(2) - B(3)*A(2) ) / 2.
C   COMPUTATION OF LOCAL STIFFNESS MATRIX
      DO I = 1 , ND
      DO J = 1 , ND
      STIFF(I,J) = ( (B(I)*EXX+A(I)*EXY) * B(J)
     *           +   (B(I)*EXY+A(I)*EYY) * A(J) )/(4.*AREA)
      END DO
      END DO
      RETURN
      END
C
C
      SUBROUTINE FORM ( MXN, MXB, NNODE, NB,IBND,ITYPE,BVALUE, RK, RHS )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION RHS(MXN),RK(MXN,MXN),IBND(MXB),ITYPE(MXB),BVALUE(MXB)
C------ CLEAR RIGHT HAND SIDE --------
      DO I = 1 , NNODE
      RHS (I) = 0.
      END DO
C  BOUNDARY CONDITION IMPLEMENTATION
C  ITYPE(I) = 1 ---> DIRICHLET,    ITYPE(I) = 2 ---> NEUMANN
C
      DO K = 1 , NB
      L = IBND (K)
      IF ( ITYPE (K) .EQ. 1 ) THEN
      DO J = 1 , NNODE
      RK (L,J) = 0.
      END DO
      RK (L,L) = 1.
      RHS (L) =  BVALUE (K)
      ELSE
      RHS (L) = -BVALUE (K)
      ENDIF
      END DO
      RETURN
      END
C
C
      SUBROUTINE SYSTEM ( MXN, N , A , C )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION A (MXN,MXN) , C (MXN)
      N1 = N - 1
      DO K = 1, N1
      L = K + 1
      DO I = L , N
      P = A (I,K) / A (K,K)
      IF ( P .NE. 0. ) THEN
      DO J = L , N
      A (I,J) = A (I,J) - P * A ( K , J )
      END DO
      C ( I ) = C ( I) - P * C ( K )
      END IF
      END DO
      END DO
C---- BACK SUBSTITUTION
      C (N) = C (N) / A (N,N)
      DO K = 1 , N1
      I = N - K
      L = I + 1
      P = C ( I )
      DO J = L , N
      P = P - C (J) * A (I,J)
      END DO
      C ( I ) = P / A (I,I)
      END DO
      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*12
      IR = 1
      IF ( ND .LE. 2 ) STOP 'NOT ACCEPTABLE ND'
      IF ( ND .EQ. 3 ) INPFILE = 'FEM3DATA.QQQ'
      IF ( ND .EQ. 4 ) INPFILE = 'FEM4DATA.QQQ'
      IF ( ND .EQ. 8 ) INPFILE = 'FEM8DATA.QQQ'
      IF ( ND .EQ. 9 ) INPFILE = 'FEM9DATA.QQQ'
      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 ( 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*12
      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
C-------- FORMATS --------------
      RETURN
      END
