C
C PROGRAM 7.1 SOLUTION OF LAPLACE'S EQUATION FOR
C PLANE FREE-SURFACE FLOW USING 8-NODED QUADRILATERALS
C
C ALTER NEXT LINE TO CHANGE PROBLEM SIZE
C
PARAMETER (IKV=416000,ILOADS=3500,INF=3000,INO=100,INX=60)
C
REAL PERMX
REAL PERMY
REAL DET
REAL QUOT
REAL REACT
REAL JAC(2,2),JAC1(2,2),KAY(2,2),SAMP(3,2),DTKD(4,8),KP(4,8),
+COORD(8,2),DER(2,8),DERIV(2,8),DERIVT(8,2),KDERIV(2,8),
+FUN(8),VAL(INO),KVH(IKV),KV(IKV),LOADS(ILOADS),DISPS(ILOADS),
+OLDPOT(ILOADS),WIDTH(INX),SURF(INX)
INTEGER G(4),NO(INO),NF(INF,1)
DATA IT,IJAC,IJAC1,IKAY,IDER,IDERIV,IKDERV/7*2/,ISAMP/3/
DATA IDTKD,IKP,ICOORD,IDERVT,NOD/5*4/,NODOF/1/
C
C INPUT AND INITIALISATION
C
READ (5,FMT=*) NXE,NYE,N,IW,NN,NR,NGP,ITS,PERMX,PERMY
READ (5,FMT=*) (WIDTH(I),I=1,NXE+1)
READ (5,FMT=*) (SURF(I),I=1,NXE+1)
CALL READNF(NF,INF,NN,NODOF,NR)
IR = N* (IW+1)
CALL NULVEC(OLDPOT,N)
CALL NULL(KAY,IKAY,IT,IT)
KAY(1,1) = PERMX
KAY(2,2) = PERMY
CALL GAUSS(SAMP,ISAMP,NGP)
C
C ITERATE FOR POSITION OF FREESURFACE
C
ITERS = 0
10 ITERS = ITERS + 1
CALL NULVEC(KV,IR)
C
C ELEMENT INTEGRATION AND ASSEMBLY
C
DO 20 IP = 1,NXE
DO 20 IQ = 1,NYE
CALL WELGEO(IP,IQ,NXE,NYE,WIDTH,SURF,COORD,ICOORD,G,NF,
+INF)
CALL NULL(KP,IKP,NOD,NOD)
DO 30 I = 1,NGP
DO 30 J = 1,NGP
CALL FORMLN(DER,IDER,FUN,SAMP,ISAMP,I,J)
CALL MATMUL(DER,IDER,COORD,ICOORD,JAC,IJAC,IT,NOD,
+IT)
CALL TWOBY2(JAC,IJAC,JAC1,IJAC1,DET)
CALL MATMUL(JAC1,IJAC1,DER,IDER,DERIV,IDERIV,IT,
+IT,NOD)
CALL MATMUL(KAY,IKAY,DERIV,IDERIV,KDERIV,IKDERV,
+IT,IT,NOD)
CALL MATRAN(DERIVT,IDERVT,DERIV,IDERIV,IT,NOD)
CALL MATMUL(DERIVT,IDERVT,KDERIV,IKDERV,DTKD,
+IDTKD,NOD,IT,NOD)
QUOT = DET*SAMP(I,2)*SAMP(J,2)
CALL MSMULT(DTKD,IDTKD,QUOT,NOD,NOD)
30 CALL MATADD(KP,IKP,DTKD,IDTKD,NOD,NOD)
20 CALL FORMKV(KV,KP,IKP,G,N,NOD)
CALL VECCOP(KV,KVH,IR)
C
C SPECIFY FIXED POTENTIALS AND REDUCE EQUATIONS
C
IF (ITERS.EQ.1) READ (5,FMT=*) IFIX, (NO(I),I=1,IFIX)
CALL NULVEC(LOADS,N)
DO 40 I = 1,IFIX
KV(NO(I)) = KV(NO(I)) + 1.E20
40 LOADS(NO(I)) = KV(NO(I))*SURF(NXE+1)
DO 50 IQ = 1,NYE - 1
J = IQ* (NXE+1) + 1
50 LOADS(J) = KV(J)* (NYE-IQ)*SURF(1)/NYE
CALL BANRED(KV,N,IW)
C
C SOLVE EQUATIONS
C
CALL BACSUB(KV,LOADS,N,IW)
CALL VECCOP(LOADS,SURF,NXE)
C
C CHECK CONVERGENCE
C
CALL CHECON(LOADS,OLDPOT,N,0.001,ICON)
IF (ITERS.NE.ITS .AND. ICON.EQ.0) GO TO 10
CALL LINMUL(KVH,LOADS,DISPS,N,IW)
REACT = 0.
DO 60 I = 1,NYE
60 REACT = REACT + DISPS(I* (NXE+1))
REACT = REACT + DISPS(N)
CALL PRINTV(LOADS,N)
WRITE (6,FMT='(E12.4)') REACT
WRITE (6,FMT='(I10)') ITERS
STOP
END
SUBROUTINE READNF(NF,INF,NN,NODOF,NR)
C
C THIS SUBROUTINE READS THE NODAL FREEDOM DATA
C
INTEGER NF(INF,*)
DO 1 I = 1,NN
DO 1 J = 1,NODOF
1 NF(I,J) = 1
IF (NR.GT.0) READ (5,FMT=*) (K, (NF(K,J),J=1,NODOF),I=1,NR)
N = 0
DO 2 I = 1,NN
DO 2 J = 1,NODOF
IF (NF(I,J).NE.0) THEN
N = N + 1
NF(I,J) = N
END IF
2 CONTINUE
RETURN
END
SUBROUTINE GAUSS(SAMP,ISAMP,NGP)
C
C THIS SUBROUTINE PROVIDES THE WEIGHTS AND SAMPLING POINTS
C FOR GAUSS-LEGENDRE QUADRATURE
C
REAL SAMP(ISAMP,*)
GO TO (1,2,3,4,5,6,7) NGP
1 SAMP(1,1) = 0.
SAMP(1,2) = 2.
GO TO 100
2 SAMP(1,1) = 1./SQRT(3.)
SAMP(2,1) = -SAMP(1,1)
SAMP(1,2) = 1.
SAMP(2,2) = 1.
GO TO 100
3 SAMP(1,1) = .2*SQRT(15.)
SAMP(2,1) = .0
SAMP(3,1) = -SAMP(1,1)
SAMP(1,2) = 5./9.
SAMP(2,2) = 8./9.
SAMP(3,2) = SAMP(1,2)
GO TO 100
4 SAMP(1,1) = .86113631
SAMP(2,1) = .33998104
SAMP(3,1) = -SAMP(2,1)
SAMP(4,1) = -SAMP(1,1)
SAMP(1,2) = .34785484
SAMP(2,2) = .65214515
SAMP(3,2) = SAMP(2,2)
SAMP(4,2) = SAMP(1,2)
GO TO 100
5 SAMP(1,1) = .90617984
SAMP(2,1) = .53846931
SAMP(3,1) = .0
SAMP(4,1) = -SAMP(2,1)
SAMP(5,1) = -SAMP(1,1)
SAMP(1,2) = .23692688
SAMP(2,2) = .47862867
SAMP(3,2) = .56888888
SAMP(4,2) = SAMP(2,2)
SAMP(5,2) = SAMP(1,2)
GO TO 100
6 SAMP(1,1) = .93246951
SAMP(2,1) = .66120938
SAMP(3,1) = .23861918
SAMP(4,1) = -SAMP(3,1)
SAMP(5,1) = -SAMP(2,1)
SAMP(6,1) = -SAMP(1,1)
SAMP(1,2) = .17132449
SAMP(2,2) = .36076157
SAMP(3,2) = .46791393
SAMP(4,2) = SAMP(3,2)
SAMP(5,2) = SAMP(2,2)
SAMP(6,2) = SAMP(1,2)
GO TO 100
7 SAMP(1,1) = .94910791
SAMP(2,1) = .74153118
SAMP(3,1) = .40584515
SAMP(4,1) = .0
SAMP(5,1) = -SAMP(3,1)
SAMP(6,1) = -SAMP(2,1)
SAMP(7,1) = -SAMP(1,1)
SAMP(1,2) = .12948496
SAMP(2,2) = .27970539
SAMP(3,2) = .38183005
SAMP(4,2) = .41795918
SAMP(5,2) = SAMP(3,2)
SAMP(6,2) = SAMP(2,2)
SAMP(7,2) = SAMP(1,2)
100 CONTINUE
RETURN
END
SUBROUTINE WELGEO(IP,IQ,NXE,NYE,WIDTH,SURF,COORD,ICOORD,G,NF,INF)
C
C THIS SUBROUTINE FORMS THE COORDINATES AND STEERING VECTOR
C FOR 4-NODE QUADS NUMBERING IN THE X-DIRECTION
C LAPLACE'S EQUATION,VARIABLE MESH, 1-FREEDOM PER NODE
C
REAL BB
REAL BVAR
REAL COORD(ICOORD,*),WIDTH(*),SURF(*)
INTEGER NUM(4),G(*),NF(INF,*)
NUM(1) = IQ* (NXE+1) + IP
NUM(2) = (IQ-1)* (NXE+1) + IP
NUM(3) = NUM(2) + 1
NUM(4) = NUM(1) + 1
DO 1 I = 1,4
1 G(I) = NF(NUM(I),1)
BB = SURF(IP)/NYE
BVAR = SURF(IP+1)/NYE
COORD(1,1) = WIDTH(IP)
COORD(2,1) = WIDTH(IP)
COORD(3,1) = WIDTH(IP+1)
COORD(4,1) = WIDTH(IP+1)
COORD(1,2) = (NYE-IQ)*BB
COORD(2,2) = (NYE-IQ+1)*BB
COORD(3,2) = (NYE-IQ+1)*BVAR
COORD(4,2) = (NYE-IQ)*BVAR
RETURN
END
SUBROUTINE NULL(A,IA,M,N)
C
C THIS SUBROUTINE NULLS A 2-D ARRAY
C
REAL A(IA,*)
DO 1 I = 1,M
DO 1 J = 1,N
1 A(I,J) = 0.0
RETURN
END
SUBROUTINE FORMLN(DER,IDER,FUN,SAMP,ISAMP,I,J)
C
C THIS SUBROUTINE FORMS THE SHAPE FUNCTIONS AND
C THEIR DERIVATIVES FOR 4-NODED QUADRILATERAL ELEMENTS
C
REAL ETA
REAL XI
REAL ETAM
REAL ETAP
REAL XIM
REAL XIP
REAL DER(IDER,*),FUN(*),SAMP(ISAMP,*)
ETA = SAMP(I,1)
XI = SAMP(J,1)
ETAM = .25* (1.-ETA)
ETAP = .25* (1.+ETA)
XIM = .25* (1.-XI)
XIP = .25* (1.+XI)
FUN(1) = 4.*XIM*ETAM
FUN(2) = 4.*XIM*ETAP
FUN(3) = 4.*XIP*ETAP
FUN(4) = 4.*XIP*ETAM
DER(1,1) = -ETAM
DER(1,2) = -ETAP
DER(1,3) = ETAP
DER(1,4) = ETAM
DER(2,1) = -XIM
DER(2,2) = XIM
DER(2,3) = XIP
DER(2,4) = -XIP
RETURN
END
SUBROUTINE TWOBY2(JAC,IJAC,JAC1,IJAC1,DET)
C
C THIS SUBROUTINE FORMS THE INVERSE OF A 2 BY 2 MATRIX
C
REAL DET
REAL JAC(IJAC,*),JAC1(IJAC1,*)
DET = JAC(1,1)*JAC(2,2) - JAC(1,2)*JAC(2,1)
JAC1(1,1) = JAC(2,2)
JAC1(1,2) = -JAC(1,2)
JAC1(2,1) = -JAC(2,1)
JAC1(2,2) = JAC(1,1)
DO 1 K = 1,2
DO 1 L = 1,2
1 JAC1(K,L) = JAC1(K,L)/DET
RETURN
END
SUBROUTINE MATRAN(A,IA,B,IB,M,N)
C
C THIS SUBROUTINE FORMS THE TRANSPOSE OF A MATRIX
C
REAL A(IA,*),B(IB,*)
DO 1 I = 1,M
DO 1 J = 1,N
1 A(J,I) = B(I,J)
RETURN
END
SUBROUTINE MATMUL(A,IA,B,IB,C,IC,L,M,N)
C
C THIS SUBROUTINE FORMS THE PRODUCT OF TWO MATRICES
C
REAL X
REAL A(IA,*),B(IB,*),C(IC,*)
DO 1 I = 1,L
DO 1 J = 1,N
X = 0.0
DO 2 K = 1,M
2 X = X + A(I,K)*B(K,J)
C(I,J) = X
1 CONTINUE
RETURN
END
SUBROUTINE MSMULT(A,IA,C,M,N)
C
C THIS SUBROUTINE MULTIPLIES A MATRIX BY A SCALAR
C
REAL C
REAL A(IA,*)
DO 1 I = 1,M
DO 1 J = 1,N
1 A(I,J) = A(I,J)*C
RETURN
END
SUBROUTINE MATADD(A,IA,B,IB,M,N)
C
C THIS SUBROUTINE ADDS TWO EQUAL SIZED ARRAYS
C
REAL A(IA,*),B(IB,*)
DO 1 I = 1,M
DO 1 J = 1,N
1 A(I,J) = A(I,J) + B(I,J)
RETURN
END
SUBROUTINE FORMKV(BK,KM,IKM,G,N,IDOF)
C
C THIS SUBROUTINE FORMS THE GLOBAL STIFFNESS MATRIX
C STORING THE UPPER TRIANGLE AS A VECTOR BK(N*(IW+1))
C
REAL BK(*),KM(IKM,*)
INTEGER G(*)
DO 1 I = 1,IDOF
IF (G(I).EQ.0) GO TO 1
DO 5 J = 1,IDOF
IF (G(J).EQ.0) GO TO 5
ICD = G(J) - G(I) + 1
IF (ICD-1) 5,4,4
4 IVAL = N* (ICD-1) + G(I)
BK(IVAL) = BK(IVAL) + KM(I,J)
5 CONTINUE
1 CONTINUE
RETURN
END
SUBROUTINE NULVEC(VEC,N)
C
C THIS SUBROUTINE NULLS A COLUMN VECTOR
C
REAL VEC(*)
DO 1 I = 1,N
1 VEC(I) = 0.
RETURN
END
SUBROUTINE BANRED(BK,N,IW)
C
C THIS SUBROUTINE PERFORMS GAUSSIAN REDUCTION OF
C THE STIFFNESS MATRIX STORED AS A VECTOR BK(N*(IW+1))
C
REAL SUM
REAL BK(*)
DO 1 I = 2,N
IL1 = I - 1
KBL = IL1 + IW + 1
IF (KBL-N) 3,3,2
2 KBL = N
3 DO 1 J = I,KBL
IJ = (J-I)*N + I
SUM = BK(IJ)
NKB = J - IW
IF (NKB) 4,4,5
4 NKB = 1
5 IF (NKB-IL1) 6,6,8
6 DO 7 M = NKB,IL1
NI = (I-M)*N + M
NJ = (J-M)*N + M
7 SUM = SUM - BK(NI)*BK(NJ)/BK(M)
8 BK(IJ) = SUM
1 CONTINUE
RETURN
END
SUBROUTINE BACSUB(BK,LOADS,N,IW)
C
C THIS SUBROUTINE PERFORMS THE GAUSSIAN BACK-SUBSTITUTION
C
REAL SUM
REAL BK(*),LOADS(*)
LOADS(1) = LOADS(1)/BK(1)
DO 1 I = 2,N
SUM = LOADS(I)
I1 = I - 1
NKB = I - IW
IF (NKB) 2,2,3
2 NKB = 1
3 DO 4 K = NKB,I1
JN = (I-K)*N + K
SUM = SUM - BK(JN)*LOADS(K)
4 CONTINUE
LOADS(I) = SUM/BK(I)
1 CONTINUE
DO 5 JJ = 2,N
I = N - JJ + 1
SUM = 0.
I1 = I + 1
NKB = I + IW
IF (NKB-N) 7,7,6
6 NKB = N
7 DO 8 K = I1,NKB
JN = (K-I)*N + I
8 SUM = SUM + BK(JN)*LOADS(K)
LOADS(I) = LOADS(I) - SUM/BK(I)
5 CONTINUE
RETURN
END
SUBROUTINE VECCOP(A,B,N)
C
C THIS SUBROUTINE COPIES VECTOR A INTO VECTOR B
C
REAL A(*),B(*)
DO 1 I = 1,N
1 B(I) = A(I)
RETURN
END
SUBROUTINE CHECON(LOADS,OLDLDS,N,TOL,ICON)
C
C THIS SUBROUTINE SETS ICON TO ZERO IF THE RELATIVE CHANGE
C IN VECTORS 'LOADS' AND 'OLDLDS' IS GREATER THAN 'TOL'
C
REAL TOL
REAL BIG
REAL LOADS(*),OLDLDS(*)
ICON = 1
BIG = 0.
DO 1 I = 1,N
1 IF (ABS(LOADS(I)).GT.BIG) BIG = ABS(LOADS(I))
DO 2 I = 1,N
IF (ABS(LOADS(I)-OLDLDS(I))/BIG.GT.TOL) ICON = 0
2 OLDLDS(I) = LOADS(I)
RETURN
END
SUBROUTINE LINMUL(BK,DISPS,LOADS,N,IW)
C
C THIS SUBROUTINE MULTIPLIES A MATRIX BY A VECTOR
C THE MATRIX IS SYMMETRICAL WITH ITS UPPER TRIANGLE
C STORED AS A VECTOR
C
REAL X
REAL BK(*),DISPS(*),LOADS(*)
DO 1 I = 1,N
X = 0.
DO 2 J = 1,IW + 1
IF (I+J.LE.N+1) X = X + BK(N* (J-1)+I)*DISPS(I+J-1)
2 CONTINUE
DO 3 J = 2,IW + 1
IF (I-J+1.GE.1) X = X + BK((N-1)* (J-1)+I)*DISPS(I-J+1)
3 CONTINUE
LOADS(I) = X
1 CONTINUE
RETURN
END
Write, Run & Share Fortran code online using OneCompiler's Fortran online compiler for free. It's one of the robust, feature-rich online compilers for Fortran language, running on the latest version 7. Getting started with the OneCompiler's Fortran compiler is simple and pretty fast. The editor shows sample boilerplate code when you choose language as Fortran and start coding.
OneCompiler's Fortran online editor supports stdin and users can give inputs to programs using the STDIN textbox under the I/O tab. Following is a sample Fortran program which takes name as input and prints hello message with your name.
program hello
character :: name*30
read *, name
print *, "Hello ", name
end program hello
Fortran language was initially developed for scientific calculations by IBM in 1957. It has a number of in-built functions to perform mathematical calculations and is ideal for applications which has more mathematical calculations.
| Data type | Description | Usage |
|---|---|---|
| Integer | To store integer variables | integer :: x |
| Real | To store float values | real :: x |
| Complex | To store complex numbers | complex :: x,y |
| Logical | To store boolean values True or false | logical :: x=.True. , logical :: x = .FALSE. |
| Character | To store characters and strings | character :: x |
Variable is a name given to the storage area in order to manipulate them in our programs.
data type :: variable_name
Array is a collection of similar data which is stored in continuous memory addresses.
data-type, dimension (x,y) :: array-name
integer, dimension(3,3) :: cube
Do is used to execute a set of statement(s) iteratively when a given condition is true and the loop variable must be an integer.
do i = start, stop [,step]
! code
end do
Do-While is used to execute a set of statement(s) iteratively when a given condition is true.
do while (condition)
!Code
end do
If is used to execute a set of statements based on a condition.
if (logical-expression) then
!Code
end if
If is used to execute a set of statements based on a condition and execute another set of statements present in else block, if condition specified in If block fails.
if (logical-expression) then
!code when the condition is true
else
!code when the condition fails
end if
Case is similar to switch in C language.
[name:] select case (regular-expression)
case (value1)
! code for value 1
... case (value2)
! code for value 2
...
case default
! default code
...
end select [name]