PROGRAM FEM3Q
C=======================================================================
C FEM 2-DIM PROGRAM FOR SOLVING LAPLACE EQUATION
C USING UPPER HALF BANDED MATRIX.
C ELEMENT: 3-NODED TRIANGULAR
C ORIGINAL: EIJI FUKUMORI 1983 BUFFALO NY, REVISED 1994 AICHI
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=3,MXE=600,MXN=290,MXB=100,MXW=20 )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),RK(MXN,MXW),
* RHS(MXN),IBND(MXB),ITYPE(MXB),BVALUE(MXB),
* SK(ND,ND),X(ND),Y(ND),A(ND),B(ND)
CHARACTER INPFILE*12
IF ( ND .NE. 3 ) STOP 'ND MUST BE 3.'
C=======================================================================
WRITE (*,*)' FEM3Q.FOR SOLVER'
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'NBW>MXW'
C=======================================================================
CALL GLBSTIFF ( ND,MXE,MXN,MXW,NE,NNODE,NBW,EXX,EXY,EYY,
* X,Y,A,B,SK,NODEX,XCOORD,YCOORD,RK )
C=======================================================================
CALL FORMQ ( MXN,MXB,MXW,NNODE,NB,NBW,RK,RHS,ITYPE,BVALUE,IBND )
C=======================================================================
CALL SYSTEM ( MXN, MXW, NNODE, NBW, RK, RHS )
C=======================================================================
CALL RESULT ( INPFILE, ND,MXE,MXN,MXB,NE,NNODE,NB,
* EXX,EXY,EYY,NODEX, XCOORD, YCOORD, IBND, ITYPE, BVALUE, RHS )
C=======================================================================
STOP
END
C
C
SUBROUTINE GLBSTIFF ( ND,MXE,MXN,MXW,NE,NNODE,NBW,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,MXW),
* SK(ND,ND),X(ND),Y(ND),A(ND),B(ND)
DO I = 1 , NNODE
DO J = 1 , NBW
RK ( I , J ) = 0.
END DO
END DO
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 , A , B, X , Y , EXX,EXY,EYY , SK )
DO K = 1 , ND
I = NODEX(IEL,K)
DO L = 1 , ND
J = NODEX(IEL,L) - I + 1
IF ( J .GE. 1 ) 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 I = 1 , NE
DO J = 1 , ND-1
DO K = J+1 , ND
NBW = MAX0(NBW,IABS(NODEX(I,J)-NODEX(I,K)))
END DO
END DO
END DO
NBW = NBW + 1
WRITE(*,*) 'HALF BANDWIDTH =', NBW
RETURN
END
C
C
SUBROUTINE SYSTEM ( MXN, MXW, NUMNP, MBAND, A, B )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION A(MXN,MXW) , B(MXN)
C---------- ELIMINATION ------------------
DO 30 N = 1 , NUMNP
DO 20 L = 2 , MBAND
C = A(N,L) / A(N,1)
I = N + L - 1
IF ( I .GT. NUMNP ) GO TO 20
J = 0
DO 10 K = L , MBAND
J = J + 1
A(I,J) = A(I,J) - C * A(N,K)
10 CONTINUE
A(N,L) = C
B(I) = B(I) - A(N,L) * B(N)
20 CONTINUE
B(N) = B(N) / A(N,1)
30 CONTINUE
C---------- BACKSUBSTITUTION -------------
40 DO 50 K = 2 , MBAND
L = N + K - 1
IF ( L .GT. NUMNP ) GO TO 60
B(N) = B(N) - A(N,K) * B(L)
50 CONTINUE
60 N = N - 1
IF ( N .GT. 0 ) GO TO 40
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)
DO I = 1 , NNODE
RHS (I) = 0.
END DO
DO 50 K = 1 , NB
IF ( IBC(K) .EQ. 2 ) GO TO 50
I = IBND(K)
DO 20 J = 2 , NBW
I = I - 1
IF ( I.LE. 0 ) GO TO 30
RHS(I) = RHS(I) - BV(K) * A(I,J)
20 CONTINUE
30 I = IBND(K)
DO 40 J = 2 , NBW
I = I + 1
IF ( I .GT. NNODE ) GO TO 50
RHS(I) = RHS(I) - BV(K) * A(IBND(K),J)
40 CONTINUE
50 CONTINUE
C-----REFORMATION OF MATRIX A.
DO 70 K = 1 , NB
I = IBND (K)
IF ( IBC(K) .EQ. 2 ) GO TO 65
RHS(I) = BV(K)
A(I,1) = 1.
DO 60 J = 2 , NBW
L = I - J + 1
A(I,J) = 0.
IF ( L .LE. 0 ) GO TO 60
A( L ,J) = 0.
60 CONTINUE
GO TO 70
65 RHS (I) = RHS(I) - BV(K)
70 CONTINUE
RETURN
END
C
C
SUBROUTINE STIFMTX ( ND, A, B, X, Y, EXX,EXY,EYY, STIFF )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION X(ND), Y(ND), STIFF(ND,ND), A(ND), B(ND)
A(1) = X(3) - X(2)
A(2) = X(1) - X(3)
A(3) = X(2) - X(1)
B(1) = Y(2) - Y(3)
B(2) = Y(3) - Y(1)
B(3) = Y(1) - Y(2)
AREA = ( A(3)*B(2) - B(3)*A(2) ) / 2.
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 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