PROGRAM FEM4Q
C=======================================================================
C 2-DIM FEM PROGRAM FOR SOLVING LAPLACE EQUATION
C USING UPPER HALF BANDED MATRIX
C EQUATION: DP/DX + DP/DY = 0.
C ELEMENT : 4-NODED ISO-PARAMETERIC
C NUMBERING ORDER: (-1,-1),(+1,-1),(+1,+1),(-1,+1)
C ORIGINAL:1984 EIJI FUKUMORI BUFFALO NY & REVISED: 1994 ACHI
C FEB. 13, 2013
C=======================================================================
INCLUDE 'PARAM.DAT'
IMPLICIT REAL*8 ( A-H , O-Z )
CCCCCC PARAMETER (ND=4,MXE=3400,MXN=3500,MXB=2300,INTEPT=2,MXW=520)
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBND(MXB),
* ITYPE(MXB),BVALUE(MXB),RK(MXN,MXW),RHS(MXN),SK(ND,ND),DUM1(ND),
*DUM2(ND),BX(2,ND),SAI(INTEPT),W(INTEPT),BPP(2,ND,INTEPT,INTEPT)
CHARACTER INPFILE*14
WRITE (*,*)' FEM4Q.FOR SOLVER'
C=======================================================================
CALL GRULE ( INTEPT, SAI, W )
CALL DERIV ( ND, INTEPT, DUM1, DUM2, 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 ( INPFILE, ND,MXE,MXN,MXB,NE,NNODE,NB,
* EXX,EXY,EYY,NODEX, XCOORD, YCOORD, IBND, ITYPE, BVALUE, RHS )
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 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
CCCCCCCCCCCCCCCCCCCCC SAI(1) = 0.8165D0
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 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)
F(2) = 0.25D0*(1.D0+ E1)*(1.D0- E2)
F(3) = 0.25D0*(1.D0+ E1)*(1.D0+ E2)
F(4) = 0.25D0*(1.D0- E1)*(1.D0+ E2)
RETURN
END
C
C
SUBROUTINE DERIV ( ND, INTEPT, F0, F1, SAI, BPP )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION SAI(INTEPT),BPP(2,ND,INTEPT,INTEPT),F0(ND),F1(ND)
DO K = 1 , INTEPT
E1 = SAI (K)
DO L = 1 , INTEPT
E2 = SAI (L)
COMPUTATION OF BP(J) = DN(J) / DETA1
CALL ISOPARA ( ND, E1+0.5D0 , E2 , F1 )
CALL ISOPARA ( ND, E1-0.5D0 , E2 , F0 )
DO I = 1 , ND
BPP(1,I,K,L) = F1(I) - F0(I)
END DO
COMPUTATION OF BP(J)= DN(J)/DETA2
CALL ISOPARA ( ND, E1 , E2+0.5D0 , F1 )
CALL ISOPARA ( ND, E1 , E2-0.5D0 , F0 )
DO I = 1 , ND
BPP(2,I,K,L) = F1(I) - F0(I)
END DO
END DO
END DO
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 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 SYSTEMQ ( MXN, MXW, NNODE, NBW, A, B )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION A(MXN,MXW) , B(MXN)
C---------- GAUSS ELIMINATION --------------
DO N = 1 , NNODE-1
DO L = 2 , MIN((NNODE-N+1),NBW)
C = A(N,L) / A(N,1)
IF ( C .NE. 0. ) THEN
I = N + L - 1
DO K = L , NBW
J = K-L+1
A(I,J) = A(I,J) - C * A(N,K)
END DO
A(N,L) = C
B(I) = B(I) - C * B(N)
END IF
END DO
B(N) = B(N) / A(N,1)
END DO
C------ SOLUTION OF LAST UNKNOWN ---------
B(NNODE) = B(NNODE) / A(NNODE,1)
C---------- BACK SUBSTITUTION -------------
DO I = NNODE-1 , 1 , -1
DO J = 2 , MIN((NNODE-I+1),NBW)
B(I) = B(I)-A(I,J)*B(I+J-1)
END DO
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.
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 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'
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*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
C-------- FORMATS --------------
RETURN
END