戻る

C ******************************************************************** C
C C
C F 2 D - F E M C
C C
C 2D Fluid Analysis Software by FEM C
C C
C ( uvp / Quadrilateral Elements ) C
C ============== C
C C
C ( Version 5.1 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ******************************************************************** C
C
PROGRAM MAIN
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
C ***** DATA ( 1 ) ***************************************************
C
C NNP : No. of Nodes NEL : No. of Elements
C
C ----- For Unstructured meshes --------------------------------------
C =====================
C
PARAMETER ( NNP = 969, NEL = 888, N = 3632,
C
C ----- For Structured meshes ( Max.case )----------------------------
C ===================
C
*** PARAMETER ( NNP = 2601, NEL = 2500, N = 10100,
C
C ********************************************************************
C
1 NNZ = NNP + N, NZ2 = NNP + 2* N)
C --------------------------------------------------------------------
C
C NNP, NEL, N NNP, NEL, N
C
C 10F 121, 100, 420 20F 441, 400, 1640
C 30F 961, 900, 3660 40F 1681, 1600, 6480
C 50F 2601, 2500, 10100
C
C ********************************************************************
C
DIMENSION U ( NNP), V ( NNP), P ( NNP), DLP( NNP),
1 SU ( NNP), SV ( NNP), U1 ( NNP), V1 ( NNP),
2 X ( NNP), Y ( NNP), A ( NNZ), B ( NNP),
3 F ( NNP), DI ( NNP), DLX( NEL,4), DLY( NEL,4),
4 NC1( NNP), NC2( NNP,10), NC3( NNP,10)
C
DIMENSION UVE( NEL), ARS( NEL), RLE( NEL), FF ( NNP),
1 GG ( NNP), HHX( NNP), HHY( NNP), DD ( NNZ),
2 PP ( NNZ), D ( NNZ), AA ( 10,NNP), NDE( NEL,4)
C
DIMENSION KCC( NNP), NCC( 2,NNZ), NCW( 2,NZ2), NBU( NNP),
1 NBV( NNP), NBP( NNP)
C
DIMENSION XLX( 51), YLY( 51), XI ( 4,2,2), W ( 4,2),
1 EPP( 4,4), EDD( 4,4)
C
C ----- Attention --- LOP = < 3000 ---
C
COMMON /CNST/ DT, COF, PHI, RE
COMMON /SDTA/ LOP, EPS, VRV( 3000)
COMMON /CST2/ ROP, ERO, ITR, NIT
COMMON /OTRS/ LPE, EPM, LYX, LYY, UMX, VMX, NPR
C
C ***** DATA ( 2 ) ***************************************************
C
C ----- For Unstructured meshes ( 30 X 30 ) --------------------------
C
OPEN ( 10,FILE ='F2D_FEM.DAT' )
C
CALL USMS ( X, Y, NDE, NNP, NEL, XMIN )
C
C ----- For Structured meshes ( 50 X 50) -----------------------------
C
*** OPEN ( 10, FILE ='50S_F2D_FEM.DAT' )
C
*** CALL STMS ( X, Y, NDE, NNP, NEL )
C
C ----------------------------------------------------------------------
C
*** CALL PRM10F ( XLX, YLY, LYX, LYY, N, NNP, NEL, NNZ, NZ2 )
C
*** CALL PRM20F ( XLX, YLY, LYX, LYY, N, NNP, NEL, NNZ, NZ2 )
C
*** CALL PRM30F ( XLX, YLY, LYX, LYY, N, NNP, NEL, NNZ, NZ2 )
C
*** CALL PRM40F ( XLX, YLY, LYX, LYY, N, NNP, NEL, NNZ, NZ2 )
C
*** CALL PRM50F ( XLX, YLY, LYX, LYY, N, NNP, NEL, NNZ, NZ2 )
C
C ***** DATA - File ( 3 )***********************************************
C
* OPEN ( 7, FILE ='UM30_R100_F2D_FEM_RSLT.OUT' )
* OPEN ( 11, FILE ='UM30_R100_F2D_FEM_PLT.OUT' )
C
* OPEN ( 7, FILE ='UM30_R400_F2D_FEM_RSLT.OUT' ) ! Re = 400, DT = 0.03 -- 715
* OPEN ( 11, FILE ='UM30_R400_F2D_FEM_PLT.OUT' )
C
OPEN ( 7, FILE ='UM30_R1000_F2D_FEM_RSLT.OUT' )
C ! Re = 1000, DT = 0.02 -- 1310
OPEN ( 11, FILE ='UM30_R1000_F2D_FEM_PLT.OUT' )
C
* OPEN ( 7, FILE ='M50_R1000_F2D_FEM_RSLT.OUT' )
* OPEN ( 11, FILE ='M50_R1000_F2D_FEM_PLT.OUT' )
C
C **********************************************************************
C
CALL INPT ( X, Y, U, V, P, NDE, PP, NBU, NBV,
1 NBP, NCC, NCW, KCC, ISL, EPP, EDD, DD, XLX,
2 YLY, XI, W, ARS, RLE, NC1, NC2, NC3, NNP,
3 NEL, NNZ, NZ2, XMIN )
C
LOP = 0
C
ITF = 0
C
C ----------------------------------------------------------
999 LOP = LOP + 1
C
NIT = 0
C
CALL FGEM ( FF, GG, X, Y, DT, COF, U, V, NC1, NC2,
1 NC3, NDE, XI, W, AA, HHX, HHY, NNP, NEL, LOP )
C
C ----- Pseudo-Velocities ---
C
CALL CMPA ( X, Y, U, V, U, NDE, P, PP, A, B,
1 D, F, U1, NBU, NCC, NCW, KCC, NNZ, NZ2, NNP,
2 NEL )
C
CALL CMPA ( X, Y, U, V, V, NDE, P, PP, A, B,
1 D, F, V1, NBV, NCC, NCW, KCC, NNZ, NZ2, NNP,
2 NEL )
C
CALL PSSN ( X, Y, U1, V1, P, NDE, A, B, D, F,
1 NCC, NCW, KCC, NBP, DD, NNZ, NZ2, NNP, NEL )
C
C ----- Velocity "U" ---
C
CALL CMPB ( X, Y, U, V, U, NDE, P, PP, A, B,
1 D, F, RLE, 1, SU, NBU, NCC, NCW, KCC, ISL,
2 FF, NNZ, NZ2, NNP, NEL )
C
C ----- Velocity "V" ---
C
CALL CMPB ( X, Y, U, V, V, NDE, P, PP, A, B,
1 D, F, RLE, 2, SV, NBV, NCC, NCW, KCC, ISL,
2 FF, NNZ, NZ2, NNP, NEL )
C
CALL PSSN ( X, Y, SU, SV, DLP, NDE, A, B, D, F,
1 NCC, NCW, KCC, NBP, DD, NNZ, NZ2, NNP, NEL )
C
CALL VLPR ( X, Y, P, DLP, DI, NDE, NC3, NC2, NC1, SU,
1 SV, NBU, NBV, DLX, DLY, NNP, NEL )
C
C ----- Check ------
C
CALL CKWR ( U, V, SU, SV, P, X, Y, NDE, ISL, KYK,
1 NNZ, NZ2, NNP, NEL, ITF, *333 )
C
IF ( LOP.LT.LPE ) GO TO 999
C ----------------------------------------------------------
333 CONTINUE
C
CALL FSBR ( FF, U, V, UVE, RLE, P, X, Y, NDE, ITF,
1 NNP, NEL )
C
CLOSE ( 7)
CLOSE (10)
CLOSE (11)
C
STOP
END
C **********************************************************************
C
SUBROUTINE STMS ( X, Y, NDE, NNP, NEL )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION X( NNP), Y( NNP), NDE( NEL,4)
C
READ(10,1000) NNPE, NELM
1000 FORMAT(2I5)
READ(10,1100) ( MM, ( NDE( I,J), J = 1, 4), I = 1, NELM )
1100 FORMAT(10X,5I5)
C
C --------------------------------------
DO 100 I = 1, NEL
DO 100 J = 1, 4
C
JJ = J + 1
IF ( J.EQ.4) THEN
JJ = 1
END IF
C
100 CONTINUE
C --------------------------------------
XMIN = 1.D0
DO 200 I = 2, NNP
C
XX = X( I)
IF ( XX.LE.0.0001D0) GO TO 200
XMIN = DMIN1( XMIN, XX)
C
200 CONTINUE
C --------------------------------------
* WRITE(6,2000) XMIN
* 2000 FORMAT(/,' * Xmin.=', F7.4)
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CKWR ( U, V, SU, SV, P, X, Y, NDE,
1 ISL, KYK, NNZ, NZ2, NNP, NEL, ITF, * )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION U ( NNP), V ( NNP), P ( NNP), X( NNP), Y( NNP),
1 SU( NNP), SV( NNP), NDE( NEL,4)
C
COMMON /CNST/ DT, COF, PHI, RE
COMMON /CST2/ ROP, ERO, ITR, NIT
COMMON /SDTA/ LOP, EPS, VRV( 3000)
COMMON /OTRS/ LPE, EPM, LYX, LYY, UMX, VMX, NPR
C
C --------------------------------------------------
C
CALL RCHK ( U, V, SU, SV, P, PMX, PMN, NNP )
C
C --------------------------------------------------
C
KYK = 0
C
ITF = ITF + NIT
C -----------------------------------
C
IF ( UMX.GT.1.D0) GO TO 100
IF ( VMX.GT.1.D0) GO TO 100
IF ( EPS.LT.EPM ) GO TO 200
IF ( MOD( LOP,NPR).EQ.0)
1 WRITE(6,2000) LOP, NIT, PMN, PMX, VRV( LOP)
2000 FORMAT(' ** LOP =', I5,' TIT=',I4,' Pmin.=',F7.3,' Pmax.=',F7.3,
1 ' EPS =',F7.3,' (%)')
GO TO 990
C
100 WRITE(6,*) ' ** Error ** Velocity is over 1 '
STOP
C
200 CONTINUE
WRITE(6,2100) EPM
2100 FORMAT( ' *** Converged ! ',1P,D11.5,' * End * ' )
WRITE(6,*) ' '
WRITE(6,2000) LOP, NIT, PMN, PMX, VRV( LOP)
C
WRITE(7,*) '*****************************************************'
WRITE(7,*) '* *'
WRITE(7,*) '* F 2 D - F E M *'
WRITE(7,*) '* *'
WRITE(7,*) '* ( V.5.1 ) *'
WRITE(7,*) '* *'
WRITE(7,*) '* 2D Fluid Analysis by FEM *'
WRITE(7,*) '* *'
WRITE(7,*) '* ( uvp / Quadrilateral Elements ) *'
WRITE(7,*) '* *'
WRITE(7,*) '* Copyright 2011 : Yasuhiro MATSUDA *'
WRITE(7,*) '* *'
WRITE(7,*) '*****************************************************'
WRITE(7,*) ' '
WRITE(7,*) ' --- Ctrl Data ---'
WRITE(7,*) ' '
WRITE(7,2200) RE, DT, LPE,ISL
2200 FORMAT( 1X,' RE =',F10.0,2X,' DT =',F10.5,2X,' LOP =', I5,
1 2X,' ISL =', I5)
WRITE(7,*) ' '
WRITE(7,*) ' --- Size Information ---'
WRITE(7,*) ' '
WRITE(7,*) ' NNP =',NNP,' NEL =',NEL,' NNZ =',NNZ,' NZ2=',NZ2
WRITE(7,*) ' '
WRITE(7,*) ' --- NDE Information ---'
WRITE(7,*) ' '
WRITE(7,*) 'NDE NO. X Y NDE NO. X Y '
WRITE(7,2300) ( I, X( I), Y( I), I = 1, NNP )
2300 FORMAT( 5( 2X, I5,2F8.3))
WRITE(7,*) ' '
WRITE(7,*) ' --- Element Information ---'
WRITE(7,*)
1' ELMT NO. 1 2 3 4'
WRITE(7,2400) ( I, ( NDE( I,J), J = 1, 4), I = 1, NEL )
2400 FORMAT(3X, I5,3X,4I6,3X, I5,3X,4I6,3X, I5,3X,4I6)
WRITE(7,*) ' '
WRITE(7,*) ' --- Velocity --- LOP =',LOP
WRITE(7,*) ' NDE NO. U V NDE NO. U V'
WRITE(7,2500) ( I, U( I), V( I), I= 1, NNP )
2500 FORMAT(5(3X, I5,2F8.4))
WRITE(7,*) ' '
WRITE(7,*) ' --- Pressure --- LOP =',LOP
WRITE(7,*) ' NDE NO. P NDE NO. P NDE NO. P'
WRITE(7,2600) ( I, P( I), I= 1, NNP )
2600 FORMAT(6(3X, I6,F9.4))
WRITE(7,*) ' '
WRITE(7,2700) PMN, PMX
2700 FORMAT(5X,' Pmin.=',F15.4,3X,' Pmax.=',F15.4)
C
RETURN 1
C
990 RETURN
C
END
C **********************************************************************
C
SUBROUTINE CMPB ( X, Y, U, V, C, NDE, P, PP, A,
1 B, D, F, RLE, ID, CNW, NBC, NCC, NCW,
2 KCC, ISL, FF, NNZ, NZ2, NNP, NEL )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION X ( NNP), Y ( NNP), U ( NNP), V ( NNP),
1 C ( NNP), FF ( NNP), P ( NNP), PP ( NNZ),
2 A ( NNZ), B ( NNP), D ( NNZ), F ( NNP),
3 CNW( NNP), RLE( NEL), NBC( NNP), NCC( 2,NNZ),
4 NCW( 2,NZ2), KCC( NNP), NDE( NEL,4)
C
DIMENSION XX ( 4), YY ( 4), PE( 4), EXX( 4,4), EYY( 4,4),
1 EFX( 4,4), EFY( 4,4)
C
COMMON /CNST/ DT, COF, PHI, RE
C
C ----- MGM - Formulation ( Matrix-A & Vector-B ) ---
C
C ----------------------------
DO 100 I = 1, NNP
C
F ( I) = 0.D0
B ( I) = 0.D0
CNW( I) = C( I)
C
100 CONTINUE
C ----------------------------
DO 200 I = 1, NNZ
C
D( I) = 0.D0
200 CONTINUE
C ----------------------------
C
FSX = 0.D0
FSY = 0.D0
C
IF ( ID.EQ.1) FSX = 1.D0
IF ( ID.EQ.2) FSY = 1.D0
C
C ----------------------------------------------------------
DO 300 INL = 1, NEL
C
UE = 0.D0
VE = 0.D0
FF2 = 0.D0
C
C ------------------------------------------------
DO 400 I = 1, 4
C
N = NDE( INL,I)
XX( I) = X( N)
YY( I) = Y( N)
PE( I) = P( N)
UE = UE + U( N)/ 4.D0
VE = VE + V( N)/ 4.D0
C
400 CONTINUE
C ------------------------------------------------
DO 450 I = 1, 4
C
FF2 = FF2 + FF( NDE( INL,I))/ 4.D0
C
450 CONTINUE
C ------------------------------------------------
C
CALL JCDF ( XX, YY, EXX, EYY, EFX, EFY )
C
C ------------------------------------------------
DO 300 I = 1, 4
DO 300 J = 1, 4
C
II = NDE( INL,I)
JJ = NDE( INL,J)
FK = COF* FF2
IF ( ISL.EQ.1) DS = FK* ( EXX( I,J) + EYY( I,J))
C
C ------------------------------------------------
IF ( ISL.EQ.0) THEN
C
C ----- Attention ---
C
DX = RLE( INL)
DY = RLE( INL)
C
C ---------------------------------
FKX = FK - DX* DX/ DT /6.D0
FKY = FK - DY* DY/ DT /6.D0
C ---------------------------------
C
DS = FKX* EXX( I,J) + FKY* EYY( I,J)
C
ENDIF
C ------------------------------------------------
C
CALL DCON ( DS, II, JJ, D, NCC, NNZ, L )
C
C ------------------------------------------------
C
F( II) = F( II) + ( UE* EFX( I,J) + VE* EFY( I,J))* C( JJ)
1 + FSX* EFX( I,J)* PE( J) + FSY* EFY( I,J)* PE( J)
C
300 CONTINUE
C ----------------------------------------------------------
C
C ----- Simultaneous Linear Equations by MGM ---
C
JS = 1
C ----------------------------------------------------------
IF ( ISL.EQ.1) THEN
C
C ----------------------------------------------------------
DO 500 I = 1, NNP
C
JL = KCC( I)
C
C ------------------------------------------------
DO 600 J = JS, JL
C
KPT = NCW( 1,J)
KCL = NCW( 2,J)
A( KPT) = PP( KPT)/ DT + PHI* D( KPT)
B( I) = B ( I) + ( PP( KPT)/ DT
1 - ( 1.D0 - PHI)* D( KPT))* C( KCL)
C
600 CONTINUE
C ------------------------------------------------
C
B( I) = B( I) - F( I)
JS = JL + 1
C
500 CONTINUE
C ------------------------------------------------------------
C
CALL SOLV ( B, CNW, A, NCW, KCC, NBC, NNP, NNZ, NZ2 )
C
C ------------------------------------------------------------
ELSE
C --------------------------------------
DO 700 I = 1, NNP
C
JL = KCC( I)
FW = 0.D0
C
C ----------------------------
DO 800 J = JS, JL
C
KPT = NCW( 1,J)
KCL = NCW( 2,J)
FW = FW + PP( KPT)
B( I) = B( I) + ( PP( KPT)/ DT - D( KPT))* C( KCL)
C
800 CONTINUE
C ----------------------------
IF ( NBC( I).NE.9999) CNW( I) = ( B( I) - F( I))/ FW* DT
JS = JL + 1
C
700 CONTINUE
C --------------------------------------
C
END IF
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CMPA ( X, Y, U, V, C, NDE, P, PP,
1 A, B, D, F, CNW, NBC, NCC, NCW,
2 KCC, NNZ, NZ2, NNP, NEL )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION X ( NNP), Y ( NNP), U ( NNP), V ( NNP),
1 P ( NNP), F ( NNP), B ( NNP), C ( NNP),
2 CNW( NNP), PP ( NNZ), D ( NNZ), A ( NNZ),
3 KCC( NNP), NBC( NNP), NDE( NEL,4), NCC( 2,NNZ),
4 NCW( 2,NZ2)
C
DIMENSION XX ( 4), YY ( 4), EXX( 4,4), EYY( 4,4),
1 EFX( 4,4), EFY( 4,4)
C
COMMON /CNST/ DT, COF, PHI, RE
C
C ----- MGM Formulation ( Matrix-A, Vector-B ) ---
C
C --------------------------------------
DO 100 I = 1, NNP
C
F ( I) = 0.D0
B ( I) = 0.D0
CNW( I) = C( I)
C
100 CONTINUE
C --------------------------------------
DO 200 I = 1, NNZ
C
D( I) = 0.D0
C
200 CONTINUE
C --------------------------------------
C
DO 300 INL = 1, NEL
C
UE = 0.D0
VE = 0.D0
C
C ----------------------------
DO 400 I = 1, 4
C
N = NDE( INL,I)
XX( I) = X( N)
YY( I) = Y( N)
UE = UE + U( N)/ 4.D0
VE = VE + V( N)/ 4.D0
C
400 CONTINUE
C ----------------------------
C
C ----- Matrix generation ( D, F ) -------------------------------
C
CALL JCDF ( XX, YY, EXX, EYY, EFX, EFY )
C
C ----------------------------------------------------------------
DO 300 I = 1, 4
DO 300 J = 1, 4
C
II = NDE( INL,I)
JJ = NDE( INL,J)
DS = COF* ( EXX( I,J) + EYY( I,J))
C
C ------------------------------------------------
C
CALL DCON ( DS, II, JJ, D, NCC, NNZ, L )
C
C ------------------------------------------------
C
F( II) = F( II) + ( UE* EFX( I,J) + VE* EFY( I,J))* C( JJ)
C
300 CONTINUE
C ----------------------------------------------------------------
C
JS = 1
C ----------------------------------------------------------
DO 500 I = 1, NNP
C
JL = KCC( I)
C
C ------------------------------------------------
DO 600 J = JS, JL
C
KPT = NCW( 1,J)
KCL = NCW( 2,J)
A( KPT) = PP( KPT)/ DT
B( I) = B( I) + ( PP( KPT)/ DT - D( KPT))* C( KCL)
C
600 CONTINUE
C ------------------------------------------------
B( I) = B( I) - F( I)
JS = JL + 1
C
500 CONTINUE
C ----------------------------------------------------------
C
CALL SOLV( B, CNW, A, NCW, KCC, NBC, NNP, NNZ, NZ2 )
C
C ----------------------------------------------------------
C
RETURN
C
END
C **********************************************************************
C
SUBROUTINE DCON ( VAL, I, J, RMT, NCC, NNZ, IC )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION RMT( NNZ), NCC( 2,NNZ)
C
IC = 0
IF ( I.GT.J) RETURN
C
N1 = 1
N2 = NNZ
333 IF ( ( N1+ N2)/ 2.EQ.IC) N2 = N2 + 1
IC = ( N1+ N2)/ 2
C
C ----------------------------
IF ( NCC( 1,IC).LT.I.OR.
1 ( NCC( 1,IC) .EQ. I .AND. NCC( 2,IC).LT.J)) THEN
N1 = IC
GO TO 333
C
ELSE IF ( NCC( 1,IC).GT.I.OR.
1 ( NCC( 1,IC).EQ.I.AND.NCC( 2,IC).GT.J)) THEN
N2 = IC
GO TO 333
C
END IF
C ----------------------------
C
RMT( IC) = RMT( IC) + VAL
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DCNL ( PDV, I, J, PD, NCC, K, NNZ )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION PD( NNZ), NCC( 2,NNZ)
C
IF ( I.GT.J) RETURN
C
IE = K
N1 = 1
IF ( K.EQ.0) N2 = K + 1
IF ( K.EQ.0) GO TO 10
N2 = K
C
10 CONTINUE
N = N2
IF ( IE.EQ.0) IE = 1
C
C ----------------------------------------------------------
DO 100 II = 1, IE
C
IF ( K.EQ.0) GO TO 43
IF ( I.LT.NCC( 1,N).OR.( I.EQ.NCC( 1,N).AND.J.LT.NCC( 2,N)))
1 GO TO 21
IF ( I .EQ.NCC( 1,N).AND. J.EQ.NCC( 2,N)) GO TO 22
GO TO 23
C
21 CONTINUE
N2 = N
N = N2 - ( N2 - N1)/ 2
IF ( N.NE.N2) GO TO 100
IF ( I.LT.NCC( 1,N1).OR.( I.EQ.NCC( 1,N1).AND.J.LT.NCC( 2, N1)))
1 GO TO 31
IF ( I.EQ.NCC( 1,N1).AND.J.EQ.NCC( 2,N1)) GO TO 32
GO TO 33
C
31 N = N1
33 CONTINUE
C ---------------------------------------------------
C
CALL SUBS ( PDV, I, J, N, PD, NCC, K, NNZ )
C
C ---------------------------------------------------
RETURN
C
32 CONTINUE
N = N1
GO TO 22
C
42 N = N2
22 CONTINUE
PD( N) = PD( N) + PDV
RETURN
C
23 CONTINUE
N1 = N
N = N1 + ( N2 - N1)/ 2
IF ( N.NE.N1) GO TO 100
IF ( I.LT.NCC( 1,N2).OR.( I.EQ.NCC( 1,N2).AND.J.LT.NCC( 2,N2)))
1 GO TO 41
IF ( I.EQ.NCC( 1,N2).AND.J.EQ.NCC( 2,N2)) GO TO 42
GO TO 43
C
41 N = N2
GO TO 33
C
43 K = K + 1
IF ( K.GT.NNZ) WRITE(6,2000)
2000 FORMAT(' ** Error in DCNL ** ''NNZ''in''SOR''is too small **')
C
PD ( K) = PDV
NCC( 1,K) = I
NCC( 2,K) = J
RETURN
C
100 CONTINUE
C ----------------------------------------------------------
C
WRITE(6,2100) K
2100 FORMAT(' * Error in DCNL -- DO LOP Strange K =', I8)
C
STOP
END
C **********************************************************************
C
SUBROUTINE DCIJ ( KCC, NCC, NZ, NNP, NZ2 )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION KCC( NNP), NCC( 2,NZ2)
C
C NP, NZ : Current length of KCC, NCC
C
C NNP, NZ2 : Max. length of KCC, NCC
C
C -------------------------------------------
NP = 1
C
DO 100 I = 1, NZ
C
IF ( NP.EQ.NCC( 1,I)) GO TO 10
C
C ----- To set KCC ( NP ) ---
C
KCC( NP) = I - 1
NP = NCC( 1,I)
IF ( NP.GT.NNP ) GO TO 30
10 NCC( 1,I) = I
C
100 CONTINUE
C -------------------------------------------
C ------ To set last row's counter ---
C
KCC( NP) = NZ
C
C ------ To set ( I,J) in NCC for I > J ---
C
C ----------------------------------------------------------
DO 200 I = 2, NP
C
I1 = I - 1
IS = KCC( I1)
C
C ------------------------------------------------
DO 300 J = 1,I1
C
JS = 1
C
IF ( J.GT.1) JS = KCC( J-1) + 1
C
JL = KCC( J)
C
C --------------------------------------
DO 400 K = JS, JL
C
IF ( NCC( 2,K).NE.I) GO TO 400
K1 = K
GO TO 20
C
400 CONTINUE
C --------------------------------------
C
C ------ ( I,J) Not exist ---
GO TO 300
C
20 CONTINUE
C
C ------ To shift NCC ---
C
NZIS = NZ - IS
C
C --------------------------------------
DO 500 L = 1, NZIS
C
NCC( 1, NZ-L+2) = NCC( 1, NZ-L+1)
NCC( 2, NZ-L+2) = NCC( 2, NZ-L+1)
C
500 CONTINUE
C --------------------------------------
C
NZ = NZ + 1
IF ( NZ.GT.NZ2) GO TO 40
C
C --------------------------------------
DO 600 L = I, NP
C
KCC( L) = KCC( L) + 1
C
600 CONTINUE
C --------------------------------------
C
C ------ To insert ( I,J) element ---
C
IS = IS + 1
NCC( 1,IS) = NCC( 1,K1)
NCC( 2,IS) = J
C
300 CONTINUE
C ------------------------------------------------
C
200 CONTINUE
C ----------------------------------------------------------
RETURN
C
30 WRITE(6,2000) NNP, NP
STOP
C
40 WRITE(6,2000) NZ2, NZ
2000 FORMAT(' *** Error ( DCIJ)', 2I5)
STOP
C
END
C **********************************************************************
C
SUBROUTINE TTLE
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
WRITE(6,*) ' '
WRITE(6,*) '*****************************************************'
WRITE(6,*) '* *'
WRITE(6,*) '* F 2 D - F E M *'
WRITE(6,*) '* *'
WRITE(6,*) '* ( V.5.1 ) *'
WRITE(6,*) '* *'
WRITE(6,*) '* 2D Fluid Analysis by FEM *'
WRITE(6,*) '* *'
WRITE(6,*) '* ( uvp / Quadrilateral Elements ) *'
WRITE(6,*) '* *'
WRITE(6,*) '* Copyright 2011 : Yasuhiro MATSUDA *'
WRITE(6,*) '* *'
WRITE(6,*) '*****************************************************'
C
* WRITE(6,2000) XMIN
* 2000 FORMAT(/,' * Xmin.=', F7.4)
C ----- 1999 -
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DLXY ( X, Y, C, NDE, DLX, DLY, NNP, NEL )
C
C ----- Differentiation ---
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION X ( NNP), Y ( NNP), C( NNP), NDE( NEL,4),
1 DLX( NEL,4), DLY( NEL,4)
C
DIMENSION XX ( 4), YY ( 4), ENN( 4,4), ENZ( 2,4,4),
1 EXY( 2,4,4), YAC( 2,2,4), YIN( 2,2,4), CC ( 4)
C
COMMON /XE/ XI( 4), ETA( 4)
C
C ----------------------------------------------------------
DO 100 INL= 1, NEL
DO 100 M = 1, 4
C ----------------------------------------------------------
C
C ------------------------------------------------
DO 200 I = 1, 4
C
N = NDE( INL,I)
CC( I) = C( N)
XX( I) = X( N)
YY( I) = Y( N)
C
200 CONTINUE
C ------------------------------------------------
C
C ----- Shape functions / Derivatives for shape function ---
C
ENN( 1,M) = ( 1.D0 - XI( M))* ( 1.D0 - ETA( M))/ 4.D0
ENN( 2,M) = ( 1.D0 + XI( M))* ( 1.D0 - ETA( M))/ 4.D0
ENN( 3,M) = ( 1.D0 + XI( M))* ( 1.D0 + ETA( M))/ 4.D0
ENN( 4,M) = ( 1.D0 - XI( M))* ( 1.D0 + ETA( M))/ 4.D0
C
ENZ( 1,1,M) = - ( 1.D0 - ETA( M))/ 4.D0
ENZ( 1,2,M) = ( 1.D0 - ETA( M))/ 4.D0
ENZ( 1,3,M) = ( 1.D0 + ETA( M))/ 4.D0
ENZ( 1,4,M) = - ( 1.D0 + ETA( M))/ 4.D0
C
ENZ( 2,1,M) = - ( 1.D0 - XI ( M))/ 4.D0
ENZ( 2,2,M) = - ( 1.D0 + XI ( M))/ 4.D0
ENZ( 2,3,M) = ( 1.D0 + XI ( M))/ 4.D0
ENZ( 2,4,M) = ( 1.D0 - XI ( M))/ 4.D0
C
C ----- Jacobian matrix ---
C
C ------------------------------------------------
DO 300 I = 1, 2
DO 300 J = 1, 2
C
YAC( I,J, M) = 0.D0
C
300 CONTINUE
C ------------------------------------------------
DO 400 I = 1, 4
DO 400 K = 1, 2
C
YAC( K,1,M) = YAC( K,1,M) + ENZ( K,I,M)* XX( I)
YAC( K,2,M) = YAC( K,2,M) + ENZ( K,I,M)* YY( I)
C
400 CONTINUE
C ------------------------------------------------
C
C ----- Inverse Jacobian ---
C
DTJ = YAC( 1,1,M)* YAC( 2,2,M) - YAC( 1,2,M)* YAC( 2,1,M)
C
YIN( 1,1,M) = YAC( 2,2,M)/ DTJ
YIN( 1,2,M) = - YAC( 1,2,M)/ DTJ
YIN( 2,1,M) = - YAC( 2,1,M)/ DTJ
YIN( 2,2,M) = YAC( 1,1,M)/ DTJ
C
C ------ Transformation from "ξη" to "ｘｙ" co-ordinate
C
C --- ( Derivatives of shape functionS) ---
C
C ------------------------------------------------
DO 500 I = 1, 2
DO 500 J = 1, 4
C
EXY( I,J,M) = 0.D0
500 CONTINUE
C ------------------------------------------------
C
C ----------------------------------------------------------
DO 100 I = 1, 2
DO 100 J = 1, 4
DO 100 K = 1, 2
C ----------------------------------------------------------
C
EXY( I,J,M) = EXY( I,J,M) + YIN( I,K,M)* ENZ( K,J,M)
C
DLX( INL,M) = EXY( 1,1,M)* CC( 1) + EXY( 1,2,M)* CC( 2)
1 + EXY( 1,3,M)* CC( 3) + EXY( 1,4,M)* CC( 4)
C
DLY( INL,M) = EXY( 2,1,M)* CC( 1) + EXY( 2,2,M)* CC( 2)
1 + EXY( 2,3,M)* CC( 3) + EXY( 2,4,M)* CC( 4)
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
FUNCTION DVALX ( X, C, NDE, NNP, NEL, I )
C
C ----- Differentiation ---
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION X( NNP), C( NNP), NDE( NEL,4)
C
HX = DABS( X( NDE( I,3)) - X( NDE( I,1)))
C
C1 = C( NDE( I,1))
C2 = C( NDE( I,2))
C3 = C( NDE( I,3))
C4 = C( NDE( I,4))
C
DVALX = ( - C1 + C2 + C3 - C4)/ ( 2.D0* HX)
C
RETURN
END
C **********************************************************************
C
FUNCTION DVALY ( Y, C, NDE, NNP, NEL, I)
C
C ----- Differentiation ---
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION Y( NNP), C( NNP), NDE( NEL,4)
C
HY = DABS( Y( NDE( I,3)) - Y( NDE( I,1)))
C
C1 = C( NDE( I,1))
C2 = C( NDE( I,2))
C3 = C( NDE( I,3))
C4 = C( NDE( I,4))
C
DVALY = ( - C1 - C2 + C3 + C4 )/ ( 2.D0* HY )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE RCHK ( U, V, SU, SV, P, PMX, PMN, NNP )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION U( NNP), V( NNP), SU( NNP), SV( NNP), P( NNP)
C
COMMON /SDTA/ LOP, EPS, VRV( 3000)
COMMON /OTRS/ LPE, EPM, LYX, LYY, UMX, VMX, NPR
C
C ----- Attention --- LOP = < 3000 ---
C
EPS = 0.D0
UMX = 0.D0
VMX = 0.D0
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
UN = SU( I)
VN = SV( I)
C
UMX = DMAX1( UMX,UN)
VMX = DMAX1( VMX,VN)
C
IF ( DABS( UN).GT.0.01D0)
1 EPS = DMAX1( EPS, DABS( ( U( I) - UN)/ UN))
C
IF ( DABS( VN).GT.0.01D0)
1 EPS = DMAX1( EPS, DABS( ( V( I) - VN)/ VN))
C
U( I) = UN
V( I) = VN
C
100 CONTINUE
C ------------------------------------------------
C
C ----- EPS: Max.var.of velocity ---
C
EPS = EPS* 100.D0
C
VRV( LOP) = EPS
C
PMN = 1000.D0
PMX = - 1000.D0
C --------------------------------------
DO 200 I = 1, NNP
C
PMN = DMIN1( PMN,P( I))
PMX = DMAX1( PMX,P( I))
C
200 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FCHK ( FF, NNP )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION FF( NNP)
C
FMAX = 0.D0
FMIN = 100.D0
C
DO 100 I = 1, NNP
C
IF ( FF( I).GE.FMAX) FMAX = FF( I)
IF ( FF( I).LE.FMIN) FMIN = FF( I)
C
100 CONTINUE
C
WRITE(6,2000) FMAX, FMIN
2000 FORMAT(/,' ** f-max.=',F7.3,' ** f-min.=',F7.3)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FGEM ( FF, GG, X, Y, DT, COF, U, V,
1 NC1, NC2, NC3, NDE, XI, W, AA, HHX,
2 HHY, NNP, NEL, LOP )
C
C ----- f & g ---
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION FF ( NNP), GG ( NNP), X ( NNP),
1 Y ( NNP), U ( NNP), V ( NNP),
2 NC1( NNP), NC2( NNP,10), NC3( NNP,10),
3 NDE( NEL,4), XI ( 4,2,2), W ( 4,2),
4 AA ( 10,NNP ), HHX( NNP), HHY( NNP)
C ***************
C
DIMENSION ELX ( 4,4), ELY ( 4,4), EKX( 4,4),
1 EKY ( 4,4), DPX ( 4), DPY( 4),
2 DXDS( 2,2), DSDX( 2,2), PSI( 4),
3 DPSI( 4,2), XX ( 2, 4), EC ( 4,4),
4 PP ( 20,4), DX ( 20,4), DY ( 20,4),
5 FX ( 20,4), FY ( 20,4), WW ( 20,4),
6 WWW ( 20,4), WX ( 20,4), WY ( 20,4)
C
DIMENSION NC11( 20000), X1( 100), Y1( 100), X2( 100), Y2( 100),
C ==============
1 FF1( 10000), FF2( 10000), GG1( 10000)
C
C -------------------
C
PHI = 1.D0
C
C ----------------------------------------------------------
DO 200 ID = 1, NNP
C
KEY = NC1( ID)
C
C ----- Boundary ---
C
IF ( KEY.LT.3) GO TO 910
IF ( LOP.NE.1) GO TO 920
C
C ----- Important ---
C
C ------------------------------------------------
DO 210 I = 1, 20000
C
NC11( I) = 0
C
210 CONTINUE
C ------------------------------------------------
DO 220 NLE = 1, KEY
C
X1( NLE) = 0.D0
Y1( NLE) = 0.D0
C
X2( NLE) = 0.D0
Y2( NLE) = 0.D0
C
C ---------------------------------
DO 225 I = 1, 4
C
NC11( NDE( NC2( ID,NLE ),I))
1 = NC11( NDE( NC2( ID,NLE),I)) + 1
C
225 CONTINUE
C ---------------------------------
C
220 CONTINUE
C -----------------------------------------------
DO 230 NLE = 1, KEY
DO 230 I = 1, 4
C
IF ( NC11( NDE( NC2( ID,NLE ), I)).EQ.1) THEN
C
X1( NLE) = X( NDE( NC2( ID,NLE ), I))
Y1( NLE) = Y( NDE( NC2( ID,NLE ), I))
C
ENDIF
C
230 CONTINUE
C ------------------------------------------------
C
X2MN = 10.D0
C
C ------------------------------------------------
DO 240 NLE = 1, KEY
C
IF ( Y1( NLE).LT.Y( ID)) THEN
X2( NLE) = X1( NLE)
Y2( NLE) = Y1( NLE)
X2MN = DMIN1( X2MN, X2( NLE))
ENDIF
C
240 CONTINUE
C ------------------------------------------------
DO 250 NLE = 1, KEY
C
ABSDL = DABS( X2MN - X2( NLE))
C
IF ( ABSDL.LE.0.0001D0) THEN
C
HHX( ID) = DABS( X( ID) - X2( NLE))
HHY( ID) = DABS( Y( ID) - Y2( NLE))
C
END IF
C
250 CONTINUE
C -----------------------------------------------
C
C ----------------------------------------------------------
DO 100 NLE = 1, KEY
C
C ---------------------------
DO 110 I = 1,4
C
XX( 1,I) = X( NDE( NC2( ID,NLE),I)) - X( ID)
XX( 2,I) = Y( NDE( NC2( ID,NLE),I)) - Y( ID)
C
110 CONTINUE
C ---------------------------
C
NGY = NC3( ID, NLE)
C
C ----- To initialize element arrays ---
C
C ----------------------------
DO 120 I = 1, 4
DO 120 J = 1, 4
C ----------------------------
C
EC ( I,J) = 0.D0
EKX( I,J) = 0.D0
EKY( I,J) = 0.D0
ElX( I,J) = 0.D0
ElY( I,J) = 0.D0
C
120 CONTINUE
C ----------------------------
NL = 4
C
C ----------------------------------------------------------
DO 50 L = 1, NL
C
IF ( NL.EQ.4) NN = 2
C
CALL FGS4 ( XI( L,1,NN), XI( L,2,NN), PSI, DPSI )
C
C ----- To compute DXDS ---
C
C -----------------------------------------------------
DO 60 I = 1, 2
DO 60 J = 1, 2
C -----------------------------------------------------
C
DXDS( I,J) = 0.D0
C
C -----------------------------------------------------
DO 60 K = 1, 4
C -----------------------------------------------------
C
DXDS( I,J) = DXDS( I,J) + DPSI( K,J)* XX( I,K)
C
60 CONTINUE
C -----------------------------------------------------
C
C ----- To calculate DSDX ---
C
DTJ = DXDS( 1,1)* DXDS( 2,2) - DXDS( 1,2)* DXDS( 2,1)
C
DSDX( 1,1) = DXDS( 2,2)/ DTJ ! A11
DSDX( 2,2) = DXDS( 1,1)/ DTJ ! A22
DSDX( 1,2) = - DXDS( 1,2)/ DTJ ! A21
DSDX( 2,1) = - DXDS( 2,1)/ DTJ ! A12
C
C ----- To calculate D (PSI)/ DX ---
C
C -----------------------------------------------------
DO 70 I = 1, 4
C
DPX( I) = 0.D0
DPY( I) = 0.D0
DPX( I) = DPSI( I,1)* DSDX( 1,1) + DPSI( I,2)* DSDX( 2,1)
DPY( I) = DPSI( I,1)* DSDX( 1,2) + DPSI( I,2)* DSDX( 2,2)
C
70 CONTINUE
C -----------------------------------------------------
C
C ----- To accumulate integration point value of integarals ---
C
FAC = DTJ* W( L,NN)
C
C ----------------------------
DO 80 I = 1, 4
DO 80 J = 1, 4
C
EC ( I,J) = EC( I,J) + PSI( I)* PSI( J)* FAC
C
EKX( I,J) = EKX( I,J) + DPX( I)* DPX( J)* HHX( ID)* HHX( ID)* FAC
EKY( I,J) = EKY( I,J) + DPY( I)* DPY( J)* HHY( ID)* HHY( ID)* FAC
C
ELX( I,J) = ELX( I,J) + PSI( I)* DPX( J)* HHX( ID)* FAC
ELY( I,J) = ELY( I,J) + PSI( I)* DPY( J)* HHY( ID)* FAC
C
80 CONTINUE
C -----------------------------
C
50 CONTINUE
C ----------------------------------------------------------
C
IF ( NLE.EQ.1) THEN
C
C -------------------------------------
DO 130 I = 1, 10
C
AA( I,ID) = 0.D0
C
130 CONTINUE
C --------------------------------------
C
ENDIF
C
C ------------------------------------------------
DO 140 J = 1, 4
C
PP( NLE,J) = EC ( NGY,J)
DX( NLE,J) = EKX( NGY,J)
DY( NLE,J) = EKY( NGY,J)
FX( NLE,J) = ELX( NGY,J)
FY( NLE,J) = ELY( NGY,J)
C
140 CONTINUE
C ------------------------------------------------
DO 150 I = 1, 4
C
WX( NLE,I) = XX( 1,I)/ HHX( ID)
WY( NLE,I) = XX( 2,I)/ HHY( ID)
C
150 CONTINUE
C ------------------------------------------------
DO 160 I = 1, 4
C
WW ( NLE,I) = WX( NLE,I) + WY( NLE,I)
WWW( NLE,I) = WW( NLE,I)** 2
C
160 CONTINUE
C ------------------------------------------------
C
C ------------------------------------------------
DO 170 I = 1, 4
C
AA( 1,ID) = AA( 1,ID) + PP( NLE,I)
AA( 2,ID) = AA( 2,ID) + PP( NLE,I)* WW ( NLE,I)
C
AA( 3,ID) = AA( 3,ID) + FX( NLE,I)* WW ( NLE,I)
AA( 4,ID) = AA( 4,ID) + FY( NLE,I)* WW ( NLE,I)
AA( 5,ID) = AA( 5,ID) + FX( NLE,I)* WWW( NLE,I)
AA( 6,ID) = AA( 6,ID) + FY( NLE,I)* WWW( NLE,I)
C
AA( 7,ID) = AA( 7,ID) + DX( NLE,I)* WW ( NLE,I)
AA( 8,ID) = AA( 8,ID) + DY( NLE,I)* WW ( NLE,I)
AA( 9,ID) = AA( 9,ID) + DX( NLE,I)* WWW( NLE,I)
AA(10,ID) = AA(10,ID) + DY( NLE,I)* WWW( NLE,I)
C
170 CONTINUE
C ------------------------------------------------
C
100 CONTINUE
C ----------------------------------------------------------
C
920 CONTINUE
C
BBX = DT* U( ID)/ HHX( ID)
BBY = DT* V( ID)/ HHY( ID)
C
RRX = DT* COF / HHX( ID)/ HHX( ID)
RRY = DT* COF / HHY( ID)/ HHY( ID)
C
ABBX = DABS( BBX)
ABBY = DABS( BBY)
C
C ------ In case of small velocity ---
C
IF ( ABBX.LE.0.0001D0 .AND. ABBY.LE. 0.0001D0 ) GO TO 910
C
C ------
C
FF1( ID) = AA(1,ID)* ( BBX + BBY )
1 * ( AA(5,ID)* BBX + AA(6,ID)* BBY )
C
2 + ( AA(3,ID)* BBX + AA(4,ID)* BBY)* ( AA( 1,ID)
3 * ( 2.D0* ( RRX + RRY ) + ( BBX + BBY )** 2 )
C
4 - 2.D0* AA( 2,ID)* ( BBX + BBY ))
C
C -----
C
FF2( ID) = - ( AA( 3,ID)* BBX + AA( 4,ID)* BBY )
1 *( AA( 9,ID)* RRX + AA(10,ID)* RRY )
C
2 + ( AA(7,ID)* RRX + AA(8,ID)* RRY )
3 *( ( AA(5,ID)* BBX + AA(6,ID)* BBY )
C
4 + 2.D0* PHI* ( BBX + BBY )
5 * ( AA( 3,ID)* BBX + AA( 4,ID)* BBY ))
C
FF( ID) = FF1( ID)/ FF2( ID)
C
C ************************************************
C
GG1( ID) = - AA(1,ID)* ( BBX + BBY )
1 * ( AA(9,ID)* RRX + AA(10,ID)* RRY )
2 - ( AA(7,ID)* RRX + AA(8,ID)* RRY )
3 * ( AA(1,ID)* (2.D0* ( RRX + RRY )
C
4 + ( 1.D0 - 2.D0* PHI)* ( BBX + BBY )** 2 )
5 - 2.D0* AA(2,ID)* ( BBX + BBY ))
C
GG( ID) = GG1( ID)/ FF2( ID)
C
GO TO 930
C
910 CONTINUE
C
FF( ID) = 1.D0
GG( ID) = 1.D0
C
930 CONTINUE
C
200 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FGS4 ( XI, YI, PSI, DPSI )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION PSI( 4), DPSI( 4,2)
C
PSI( 1) = 0.25D0* ( 1.D0 - XI )* ( 1.D0 - YI )
PSI( 2) = 0.25D0* ( 1.D0 + XI )* ( 1.D0 - YI )
PSI( 3) = 0.25D0* ( 1.D0 + XI )* ( 1.D0 + YI )
PSI( 4) = 0.25D0* ( 1.D0 - XI )* ( 1.D0 + YI )
C
C ----- To compute Derivatibes of Shape Functions ---
C
DPSI( 1,1) = 0.25D0* ( YI - 1.D0 )
DPSI( 1,2) = 0.25D0* ( XI - 1.D0 )
DPSI( 2,1) = 0.25D0* ( 1.D0 - YI )
DPSI( 2,2) = 0.25D0* (-1.D0 - XI )
C
DPSI( 3,1) = 0.25D0* ( 1.D0 + YI )
DPSI( 3,2) = 0.25D0* ( 1.D0 + XI )
DPSI( 4,1) = 0.25D0* (-1.D0 - YI )
DPSI( 4,2) = 0.25D0* ( 1.D0 - XI )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE IPGN ( LX, LY, X, Y, NDE, U, V, P, NBU,
1 NBV, XLX, YLY, NNP, NEL )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION U ( NNP), V ( NNP), X( NNP), Y( NNP),
1 NBU( NNP), NBV( NNP), P( NNP), NDE( NEL,4),
2 XLX( 51), YLY( 51)
C
COMMON /CNST/ DT, COF, PHI, RE
COMMON /CST2/ ROP, ERO, ITR, NIT
C
PHI = 0.9D0
ERO = 0.0001D0
ITR = 200
C
C ----- Data generation ( X, Y, Z, NDE etc.) ( DGEN) ---
C
C ----------------------------
J = 1
C
DO 100 IY = 1, LY
DO 100 IX = 1, LX
C
X( J) = XLX( IX)
Y( J) = YLY( IY)
J = J + 1
C
100 CONTINUE
C ----------------------------
C
J = 1
C
DO 200 IY = 1, LY-1
DO 200 IX = 1, LX-1
C
ST = IX + ( IY-1)* LX
NDE( J,1) = ST
NDE( J,2) = ST + 1
NDE( J,3) = ST + LX + 1
NDE( J,4) = ST + LX
J = J + 1
C
200 CONTINUE
C ----------------------------
C
C ----- Fixed boundary ( SETINI) --- Nodes for Boundary ---
C
J = 1
C
C *******************************************
DO 300 I = 1, NNP
C
C ----- Velocities for corners ---
C
C Moving Wall : U = 1, V = 0
C
P( J) = 1.D0
C ******************
C
DXX0 = DABS( X( I))
DXX1 = DABS( X( I) - 1.D0)
C
DYY0 = DABS( Y( I))
DYY1 = DABS( Y( I) - 1.D0)
C
IF ( DYY1.LE.0.0001D0 .AND.
1 ( DXX0.LE.0.0001D0 .OR. DXX1.LE.0.0001D0)) THEN
C ***********************************************************
C
C ----- Velocity direction ( 1 ) ---
C
U( J) = 1.D0
NBU( J) = 9999
C
V( J) = 0.D0
NBV( J) = 9999
C
ELSE IF ( DYY1.LE.0.0001D0) THEN
C *****************************************
C
C ----- Velocity direction ( 2 ) ---
C
U( J) = 1.D0
C
NBU( J) = 9999
C
V ( J) = 0.D0
NBV( J) = 9999
C
C---- U = 0, V = 0 -----------
C
ELSE IF ( DXX0.LE.0.0001D0 .OR.
1 DXX1.LE.0.0001D0 .OR. DYY0.LE.0.0001D0) THEN
C -------------------------------------------------------------
C
U ( J) = 0.D0
NBU( J) = 9999
C
V ( J) = 0.D0
NBV( J) = 9999
C
ENDIF
C
J = J+1
C
300 CONTINUE
C ***************************************
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INPT ( X, Y, U, V, P, NDE, PP,
1 NBU, NBV, NBP, NCC, NCW, KCC, ISL,
2 EPP, EDD, DD, XLX, YLY, XI, W,
3 ARS, RLE, NC1, NC2, NC3, NNP, NEL,
4 NNZ, NZ2, XMIN )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION X ( NNP), Y ( NNP), U ( NNP), V( NNP),
1 P ( NNP), NDE( NEL,4), PP( NNZ), NBU( NNP),
2 NBV( NNP), NBP ( NNP), NCC( 2,NNZ), NCW( 2,NZ2),
3 KCC( NNP)
C
DIMENSION EPP( 4,4), EDD( 4,4), DD( NNZ), XLX( 51), YLY( 51),
1 XX( 4), YY ( 4)
C
DIMENSION XI ( 4,2,2), W( 4,2), NC1( NNP), NC2( NNP,10),
1 NC3( NNP,10)
C
DIMENSION ARS( NEL), RLE( NEL)
C
DIMENSION X0( 3), Y0( 3)
C
COMMON /CNST/ DT, COF, PHI, RE
COMMON /CST2/ ROP, ERO, ITR, NIT
COMMON /OTRS/ LPE, EPM, LYX, LYY, UMX, VMX, NPR
C
C ----- Title ---
C
CALL TTLE
C
C ----- Initialization ---
C
DO 100 I = 1, NNP
C
NBU( I) = 0
NBV( I) = 0
NBP( I) = 0
C
100 CONTINUE
C
DO 200 I = 1, NNP
C
U( I) = 0.D0
V( I) = 0.D0
C
200 CONTINUE
C
333 WRITE(6,2000)
2000 FORMAT(/,' *** Method (?) : * Implicit = 1 * Explicit = 0 ')
*** READ(5,*) ISL
ISL = 1
C
WRITE(6,2100) ISL
2100 FORMAT(' ** Method =',I2 )
C
IF ( ISL.NE.0.AND.ISL.NE.1) GO TO 333
C
720 WRITE(6,2200)
2200 FORMAT(' ** Reynolse Number = ?')
C
*** READ(5,*) RE
C
RE = 100.D0
C
RE = 400.D0
C
RE = 1000.D0
C =================
C
IF ( RE.LT.1) GOTO 720
COF = 1.D0/ RE
C
*** READ(5,*) DT
C
DT = 0.04D0 ! RE = 100
DT = 0.03D0 ! RE = 400
DT = 0.03D0
C
DT = 0.02D0 ! RE = 1000
C ==============================
C
*** DT = 0.01D0 ! For Structured meshes ( 50 X 50) ---
C
*** DT = 1.D0 * OK Stability point :
C For Non-structured meshes ( 30 X 30) ---
C
** DT = 0.04D0 ! For Non-structured meshes ( 30 X 30) ---
C
IF ( ISL.EQ.0) DT = 0.025D0
C
730 CONTINUE
WRITE(6,2400)
2400 FORMAT(' ** No. of Loop = ? & NPR = ?')
C
*** READ(5,*) LPE, NPR
C
LPE = 800 ! For Re= 400
NPR = 50
C
LPE = 1500 ! For Re=1000
NPR = 100
C
* LPE = 20
* NPR = 10
C -----------------
C
ROP = 1.9D0
C *****************
C
IF ( LPE.LT.1) GO TO 730
C
WRITE(6,2600) NNP, NEL, RE, DT, LPE, ROP
2600 FORMAT(/,1X,' NNP =', I5,' NEL =', I5,' * Re =',F8.1,' DT =',F6.3,
1 ' LOP =', I5,' * ROP =',F7.3)
C
WRITE(6,*) ' '
C
C ----- Data ( IPGN ) ---
C
EPM = 0.1D0
C ****************
C
C -------------------------------------------------------------------
C
CALL IPGN ( LYX, LYY, X, Y, NDE, U, V, P, NBU, NBV,
1 XLX, YLY, NNP, NEL )
C
C -------------------------------------------------------------------
C
C ----- Weighting factor ----------------------------
C
CALL SJAC
C
C ---------------------------------------------------
C
CALL STNC ( NCC, NDE, 4, PP, NNP, NEL, NNZ )
C
C ---------------------------------------------------
DO 300 I = 1, NNZ
C
NCW( 1,I) = NCC( 1,I)
NCW( 2,I) = NCC( 2,I)
C
300 CONTINUE
C ------------------------------------------------
NZ = NNZ
C
CALL DCIJ ( KCC, NCW, NZ, NNP, NZ2 )
C
C ------------------------------------------------
DO 400 I = 1, NNZ
C
DD( I) = 0.D0
PP( I) = 0.D0
C
400 CONTINUE
C ------------------------------------------------
K = 0
C
C ------------------------------------------------
DO 500 INL = 1, NEL
C
C -------------------------------------
DO 600 I = 1, 4
C
N = NDE( INL,I)
XX( I) = X( N)
YY( I) = Y( N)
C
600 CONTINUE
C -------------------------------------
C
CALL JACP ( XX, YY, EPP, EDD )
C
C --------------------------------------
DO 700 I = 1, 4
DO 700 J = 1, 4
C
II = NDE( INL,I)
JJ = NDE( INL,J)
C
C --------------------------------------------------------
C
CALL DCON ( EDD( I,J), II, JJ, DD, NCC, NNZ, L )
C
C --------------------------------------------------------
C
IF ( L.NE.0) PP( L) = PP( L) + EPP( I,J)
C
700 CONTINUE
C -------------------------------------
C
500 CONTINUE
C -------------------------------------------------------
C
CALL SNCW ( XI, W, NC1, NC2, NC3, NDE, NNP, NEL )
C
C -------------------------------------------------------
C
C ----- To compute "Area" ---
C
C --------------------------------------
DO 800 I = 1, NEL
C
X0( 1) = X( NDE(I,1))
X0( 2) = X( NDE(I,2))
X0( 3) = X( NDE(I,3))
Y0( 1) = Y ( NDE(I,1))
Y0( 2) = Y ( NDE(I,2))
Y0( 3) = Y ( NDE(I,3))
C
A1 = X0( 1)* ( Y0( 2) - Y0( 3))
A2 = X0( 2)* ( Y0( 3) - Y0( 1))
A3 = X0( 3)* ( Y0( 1) - Y0( 2))
D1 = ( A1 + A2 + A3 )/ 2.D0
C
X0( 1) = X( NDE(I,3))
X0( 2) = X( NDE(I,4))
X0( 3) = X( NDE(I,1))
Y0( 1) = Y( NDE(I,3))
Y0( 2) = Y( NDE(I,4))
Y0( 3) = Y( NDE(I,1))
C
A1 = X0( 1)* ( Y0( 2) - Y0( 3))
A2 = X0( 2)* ( Y0( 3) - Y0( 1))
A3 = X0( 3)* ( Y0( 1) - Y0( 2))
D2 = ( A1 + A2 + A3 )/ 2.D0
C
ARS( I) = D1 + D2
RLE( I) = DSQRT( ARS( I))
C
800 CONTINUE
C --------------------------------------
C
C ----- To check computation area ---
C
ARSS = 0.D0
C --------------------------------------
DO 900 I = 1, NEL
C
ARSS = ARSS + ARS( I)
C
900 CONTINUE
C --------------------------------------
C
*** WRITE(6,*) ' *** ARSS =', ARSS ! Exact : 1
*** WRITE(6,*) ' '
C
RETURN
END
C **********************************************************************
C
SUBROUTINE JCDF ( XX, YY, EXX, EYY, EFX, EFY )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION XX ( 4), YY ( 4), EXX( 4,4), EYY( 4,4),
1 EFX( 4,4), EFY( 4,4), ENN( 4), ENZ( 2,4),
2 EXY( 2,4), YAC( 2,2), YIN( 2,2)
C
COMMON /XE/ XI( 4), ETA( 4)
COMMON /PW/ PG( 8), WG ( 8), NPW
C
C ----- Construction for D and F ---
C
C --------------------------------
DO 100 I = 1, 4
DO 100 J = 1, 4
C
EXX( I,J) = 0.D0
EYY( I,J) = 0.D0
EFX( I,J) = 0.D0
EFY( I,J) = 0.D0
C
100 CONTINUE
C --------------------------------
C
C ----------------------------------------------------------
DO 200 II = 1, NPW
DO 200 JI = 1, NPW
C ----------------------------------------------------------
C
C ----- Shape function etc. ---
C
C ------------------------------------------------
DO 300 I = 1, 4
C
XO = PG( II)* XI ( I)
YO = PG( JI)* ETA( I)
ENN( I) = ( 1.D0 + XO)* ( 1.D0 + YO)/ 4.D0
ENZ( 1,I) = XI ( I)* ( 1.D0 + YO)/ 4.D0
ENZ( 2,I) = ETA( I)* ( 1.D0 + XO)/ 4.D0
C
300 CONTINUE
C ------------------------------------------------
C
C ----- Jacobian matrix ---
C
DO 400 I = 1,2
DO 400 J = 1,2
C
YAC( I,J) = 0.D0
C
400 CONTINUE
C ------------------------------------------------
DO 500 I = 1, 4
DO 500 K = 1, 2
C
YAC( K,1) = YAC( K,1) + ENZ( K,I)* XX( I)
YAC( K,2) = YAC( K,2) + ENZ( K,I)* YY( I)
C
500 CONTINUE
C ------------------------------------------------
C
C ----- Inversed Jacobian matrix ---
C
DTJ = YAC( 1,1)* YAC( 2,2) - YAC( 1,2)* YAC( 2,1)
YIN( 1,1) = YAC( 2,2)/ DTJ
YIN( 1,2) = - YAC( 1,2)/ DTJ
YIN( 2,1) = - YAC( 2,1)/ DTJ
YIN( 2,2) = YAC( 1,1)/ DTJ
C
C ----- Transformation fROP "ξηζ" to "ｘｙｚ" co-ordinate
C
C --- ( Derivatives of shape functionS) ---
C
C ------------------------------------------------
DO 600 I = 1, 2
DO 600 J = 1, 4
C
EXY( I,J) = 0.D0
C
600 CONTINUE
C ------------------------------------------------
DO 700 I = 1, 2
DO 700 J = 1, 4
DO 700 K = 1, 2
C
EXY( I,J) = EXY( I,J) + YIN( I,K)* ENZ( K,J)
C
700 CONTINUE
C ------------------------------------------------
C
C ----- Numerical integration ( Element matriX) ---
C
WGT = WG( II)* WG( JI)* DTJ
C
C ----------------------------------------------------------
DO 200 I = 1, 4
DO 200 J = 1, 4
C ----------------------------------------------------------
C
EXX( I,J) = EXX( I,J) + EXY( 1,I)* EXY( 1,J)* WGT
EYY( I,J) = EYY( I,J) + EXY( 2,I)* EXY( 2,J)* WGT
C
EFX( I,J) = EFX( I,J) + ENN( I) * EXY( 1,J)* WGT
EFY( I,J) = EFY( I,J) + ENN( I) * EXY( 2,J)* WGT
C
200 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE JACP ( XX, YY, EPP, EDD)
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION XX ( 4), YY ( 4), EPP( 4,4), EDD( 4,4), ENN( 4),
1 ENZ( 2,4), EXY( 2,4), YAC( 2,2), YIN( 2,2)
C
COMMON /XE/ XI( 4), ETA( 4)
COMMON /PW/ PG( 8), WG ( 8), NPW
C
C ----- PD-Matrixes by Jacobian formulation ---
C
C ----------------------------------------------------------
DO 100 I = 1, 4
DO 100 J = 1, 4
C
EPP( I,J) = 0.D0
EDD( I,J) = 0.D0
C
100 CONTINUE
C ----------------------------------------------------------
DO 200 II = 1, NPW
DO 200 JI = 1, NPW
C ----------------------------------------------------------
C
C ----- Shape functions, Derivatives of shape functions ---
C
C ------------------------------------------------
DO 300 I = 1, 4
C
XO = PG( II)* XI ( I)
YO = PG( JI)* ETA( I)
ENN( I) = ( 1.D0 + XO )* ( 1.D0 + YO )/ 4.D0
C
ENZ( 1,I) = XI ( I)* ( 1.D0 + YO )/ 4.D0
ENZ( 2,I) = ETA( I)* ( 1.D0 + XO )/ 4.D0
C
300 CONTINUE
C ------------------------------------------------
C
C ----- Jacobian matrix ---
C
DO 400 I = 1, 2
DO 400 J = 1, 2
C
YAC( I,J) = 0.D0
C
400 CONTINUE
C -----------------------------------------------
DO 500 I = 1, 4
DO 500 K = 1, 2
C
YAC( K,1) = YAC( K,1) + ENZ( K,I)* XX( I)
YAC( K,2) = YAC( K,2) + ENZ( K,I)* YY( I)
C
500 CONTINUE
C -----------------------------------------------
C
C ----- Inversed Jacobian Matrix ---
C
DTJ = YAC( 1,1)* YAC( 2,2) - YAC( 1,2)* YAC( 2,1)
YIN( 1,1) = YAC( 2,2)/ DTJ
YIN( 1,2) = - YAC( 1,2)/ DTJ
YIN( 2,1) = - YAC( 2,1)/ DTJ
YIN( 2,2) = YAC( 1,1)/ DTJ
C
C ----- Transformation fROP "ξηζ" to "ｘｙｚ" co-ordinate
C
C --- ( Derivatives of shape functionS) ---
C
C ------------------------------------------------
DO 600 I = 1, 2
DO 600 J = 1, 4
C
EXY( I,J) = 0.D0
C
600 CONTINUE
C ------------------------------------------------
DO 700 I = 1, 2
DO 700 J = 1, 4
DO 700 K = 1, 2
C
EXY( I,J) = EXY( I,J) + YIN( I,K)* ENZ( K,J)
C
700 CONTINUE
C -----------------------------------------------
C
C ----- Numerical Integration ---
C
WGT = WG( II)* WG( JI)* DTJ
C
C ----------------------------------------------------------
DO 200 I = 1, 4
DO 200 J = 1, 4
C ----------------------------------------------------------
C
EPP(I,J) = EPP(I,J) + ENN( I)* ENN( J)* WGT
C
EDD(I,J) = EDD(I,J)
1 + ( EXY(1,I)* EXY(1,J) + EXY(2,I)* EXY(2,J))* WGT
C
200 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE USMS ( X, Y, NDE, NNP, NEL, XMIN )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION X( NNP), Y( NNP), NDE( NEL,4)
C
READ(10,1000) NNPE, NELM
1000 FORMAT(2I5)
READ(10,1100) ( MM, ( NDE( I,J), J = 1, 4), I = 1, NELM)
1100 FORMAT(5I5)
READ(10,1200) ( NN, X( I), Y( I), I = 1, NNPE )
1200 FORMAT(I5,2E15.7)
C
C --------------------------------------
DO 100 I = 1, NEL
DO 100 J = 1, 4
C
JJ = J + 1
IF ( J.EQ.4) THEN
JJ = 1
END IF
C
100 CONTINUE
C --------------------------------------
C
XMIN = 1.D0
DO 200 I = 2, NNP
C
XX = X( I)
IF ( XX.LE.0.0001D0) GO TO 200
XMIN = DMIN1( XMIN, XX)
C
200 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRM10F ( XLX, YLY, LYX, LYY, N, NNP, NEL, NNZ, NZ2 )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION XLX( 51), YLY( 51)
C
NEX = 10
NEY = 10
C
LYX = NEX + 1
LYY = NEY + 1
C
NNP = LYX* LYY
NEL = NEX* NEY
C
N = 4* ( NEX* NEY) + NEX + NEY
NNZ = NNP + N
NZ2 = NNP + 2* N
C
XLX( 1) = 0.D0
C --------------------------------------
DO 100 I = 2, 5
C
XLX( I) = ( DEXP( 2.D0* DFLOAT( I-1)* 2.D0/DFLOAT( NEX)) - 1.D0)
1 / ( 2.D0* ( DEXP( 2.D0)-1.D0))
C
100 CONTINUE
C --------------------------------------
XLX( 6) = 0.5D0
C
C ---------------------------------------
DO 200 I = 7, 10
C
XLX( I) = 1.D0 - XLX( 12-I)
200 CONTINUE
C --------------------------------------
XLX( 11) = 1.D0
C
C --------------------------------------
DO 300 I = 1, 11
C
YLY( I) = XLX( I)
300 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRM20F ( XLX, YLY, LYX, LYY, N, NNP, NEL, NNZ, NZ2 )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION XLX( 51), YLY( 51)
C
NEX = 20
NEY = 20
LYX = NEX + 1
LYY = NEY + 1
NNP = LYX* LYY
NEL = NEX* NEY
C
N = 4* ( NEX* NEY ) + NEX + NEY
NNZ = NNP + N
NZ2 = NNP + 2* N
C
C --------------------------------------
XLX( 1) = 0.D0
DO 100 I = 1,10
C
XLX( I) = ( DEXP( 2.D0* DFLOAT( I-1)*2.D0 / DFLOAT( NEX))- 1.D0 )
1 / ( 2.D0* ( DEXP( 2.D0) - 1.D0))
C
100 CONTINUE
C --------------------------------------
XLX( 1) = 0.D0
XLX(11) = 0.5D0
DO 200 I = 12, 20
C
XLX( I) = 1.D0 - XLX( 22-I)
200 CONTINUE
C --------------------------------------
XLX( 21) = 1.D0
DO 300 I = 1, 21
C
YLY( I) = XLX( I)
300 CONTINUE
C --------------------------------------
RETURN
END
C **********************************************************************
C
SUBROUTINE PRM30F ( XLX, YLY, LYX, LYY, N, NNP, NEL, NNZ, NZ2 )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION XLX( 51), YLY( 51)
C
NEX = 30
NEY = 30
LYX = NEX + 1
LYY = NEY + 1
NNP = LYX* LYY
NEL = NEX* NEY
C
N = 4* ( NEX* NEY ) + NEX + NEY
NNZ = NNP + N
NZ2 = NNP + 2* N
C
C --------------------------------------
XLX( 1) = 0.D0
C
DO 100 I = 1, 15
C
XLX( I) = ( DEXP( 2.D0* DFLOAT( I-1)* 2.D0 / DFLOAT( NEX))
1 - 1.D0)/ ( 2.D0* ( DEXP ( 2.D0) - 1.D0))
C
100 CONTINUE
C --------------------------------------
XLX( 1) = 0.D0
XLX(16) = 0.5D0
DO 200 I = 17, 30
C
XLX( I) = 1.D0 - XLX( 32-I)
200 CONTINUE
C --------------------------------------
XLX(31) = 1.D0
C
DO 300 I = 1, 31
C
YLY( I) = XLX( I)
300 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRM40F ( XLX, YLY, LYX, LYY, N, NNP, NEL, NNZ, NZ2 )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION XLX( 51), YLY( 51)
C
NEX = 40
NEY = 40
LYX = NEX + 1
LYY = NEY + 1
NNP = LYX* LYY
NEL = NEX* NEY
C
N = 4* ( NEX* NEY ) + NEX + NEY
NNZ = NNP + N
NZ2 = NNP + 2* N
C
C --------------------------------------
XLX( 1) = 0.D0
C
DO 100 I = 1,20
C
XLX( I) =( DEXP( 2.D0* DFLOAT( I-1)* 2.D0 / DFLOAT( NEX)) - 1.D0)
1 / ( 2.D0* ( DEXP( 2.D0) - 1.D0 ))
C
100 CONTINUE
C --------------------------------------
XLX( 1) = 0.D0
XLX(21) = 0.5D0
DO 200 I = 22, 40
C
XLX( I) = 1.D0 - XLX( 42-I)
200 CONTINUE
C --------------------------------------
XLX(41) = 1.D0
DO 300 I = 1, 41
C
YLY( I) = XLX( I)
300 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRM50F ( XLX, YLY, LYX, LYY, N, NNP, NEL, NNZ, NZ2 )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION XLX( 51), YLY( 51)
C
NEX = 50
NEY = 50
LYX = NEX + 1
LYY = NEY + 1
NNP = LYX* LYY
NEL = NEX* NEY
C
N = 4* ( NEX* NEY ) + NEX + NEY
NNZ = NNP + N
NZ2 = NNP + 2* N
C
C --------------------------------------
XLX( 1) = 0.D0
DO 100 I = 1, 25
C
XLX( I) =( DEXP( 2.D0* DFLOAT( I-1)* 2.D0/ DFLOAT( NEX)) - 1.D0 )
1 / ( 2.D0* ( DEXP( 2.D0) - 1.D0 ))
C
100 CONTINUE
C --------------------------------------
XLX( 1) = 0.D0
XLX(26) = 0.5D0
DO 200 I = 27, 50
C
XLX( I) = 1.D0 - XLX( 52-I)
200 CONTINUE
C --------------------------------------
C
XLX(51) = 1.D0
DO 300 I = 1, 51
C
YLY( I) = XLX( I)
300 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FSBR ( FF, U, V, UVE, RLE, P, X, Y, NDE,
1 ITF, NNP, NEL )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION FF( NNP), U ( NNP), V( NNP), P( NNP), X( NNP),
1 Y ( NNP), NDE( NEL,4)
C
DIMENSION UVE( NEL), RLE( NEL)
C
C ----- Attention --- LOP = < 3000 ---
C
COMMON /CNST/ DT, COF, PHI, RE
COMMON /CST2/ ROP, ERO, ITR, NIT
COMMON /SDTA/ LOP, EPS, VRV( 3000)
C
WRITE(6,2000) ROP, ITF
2000 FORMAT(/,' ** ROP =',F7.3,' ** Total Iteration Number =',I8)
WRITE(6,*) ' '
WRITE(6,*) '*** END *** '
C
C ------------------------------------------------
XMIN = 1.D0
DO 100 I = 2, NNP
C
XX = X( I)
AXX = DABS( X( I))
IF ( AXX.GE.0.0001D0) XMIN = DMIN1( XMIN, XX)
C
100 CONTINUE
C ------------------------------------------------
UMX = 0.D0
DO 200 I = 1, NNP
C
IF ( DABS( U(I)) .GE. UMX) UMX = DABS( U(I))
200 CONTINUE
C
*** WRITE(6,*) ' ** Abs_Umax =', UMX
C ------------------------------------------------
C
C ----- Estimation ---
C
C ----------------------------------------------------------
DO 300 I = 1, NEL
C
UU = ( U( NDE(I,1)) + U( NDE(I,2)) + U( NDE(I,3)) + U( NDE(I,4)))
1 / 4.D0
C
VV = ( V( NDE(I,1)) + V( NDE(I,2)) + V( NDE(I,3)) + V( NDE(I,4)))
1 / 4.D0
C
UVE( I) = DSQRT( UU* UU + VV* VV )
C
300 CONTINUE
C ----------------------------------------------------------
C
BMAX = 0.D0
DO 400 I = 1, NEL
C
BVR = UVE( I)* DT/ RLE( I)
IF ( BVR.GE.BMAX) THEN
BMAX = BVR
NSR = I
END IF
C
400 CONTINUE
C ----------------------------------------------------------
WRITE(6,2100) XMIN, BMAX, NSR
2100 FORMAT(/,' * Xmin.=',F7.4,' * b_max.=',F7.3,' at Node =',I5)
IF ( BMAX.GE.1.5D0)
1 WRITE(6,*) ' *** Attention : b_max. >= 1.5 ***'
C
C ----------- Correction Coefficient of "f" ---
C
CALL FCHK ( FF, NNP )
C
C ------------------------------------------------
WRITE(11,2200) LOP, DT, RE, NNP, NEL
2200 FORMAT(I4, F6.3, F8.1,2I4)
WRITE(11,2300) ( X( I), Y( I), I = 1, NNP )
2300 FORMAT(10F9.4)
WRITE(11,2400) ( ( NDE( I,J), J = 1, 4), I = 1, NEL )
2400 FORMAT(20I4)
WRITE(11,2300) ( VRV( I), I = 1, LOP )
WRITE(11,2300) ( U( I), V( I), I = 1, NNP )
WRITE(11,2300) ( P( I), I = 1, NNP )
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PSSN ( X, Y, U, V, P, NDE, A, B,
1 D, F, NCC, NCW, KCC, NBP, DD, NNZ,
2 NZ2, NNP, NEL )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION X( NNP), Y( NNP), U( NNP), V( NNP), P( NNP),
1 F( NNP), A( NNZ), B( NNP), D( NNZ), NBP( NNP),
2 DD( NNZ), NDE( NEL,4), KCC( NNP), NCC( 2,NNZ),
3 NCW( 2,NZ2)
C
DIMENSION XX ( 4), YY( 4), EXX( 4,4), EYY( 4,4), EFX( 4,4),
1 EFY( 4,4)
C
COMMON /CNST/ DT, COF, PHI, RE
COMMON /CST2/ ROP, ERO, ITR, NIT
C
C ----------------------------
DO 100 I = 1, NNP
C
F( I) = 0.D0
B( I) = 0.D0
C
100 CONTINUE
C ----------------------------
DO 200 I = 1, NNZ
C
D( I) = 0.D0
C
200 CONTINUE
C -------------------------------------------------
DO 300 INL= 1, NEL
C
UE = 0.D0
VE = 0.D0
C
C ----------------------------
DO 400 I = 1, 4
C
N = NDE( INL,I)
XX( I) = X( N)
YY( I) = Y( N)
C
400 CONTINUE
C ------------------------------------------------
C
CALL JCDF ( XX, YY, EXX, EYY, EFX, EFY )
C
C ------------------------------------------------
DO 300 I = 1, 4
DO 300 J = 1, 4
C ------------------------------------------------
C
II = NDE( INL,I)
JJ = NDE( INL,J)
DS = EXX( I,J) + EYY( I,J)
C
CALL DCON ( DS, II, JJ, D, NCC, NNZ, L )
C
F( II) = F( II)
1 + ( U( JJ)* EFX( I,J) + V( JJ)* EFY( I,J))/ DT
C
300 CONTINUE
C ------------------------------------------------
C
JS = 1
C ------------------------------------------------
DO 500 I = 1, NNP
C
JL = KCC( I)
C --------------------------------------
DO 600 J = JS, JL
C
KPT = NCW( 1,J)
KCL = NCW( 2,J)
A( KPT) = D( KPT)
C
600 CONTINUE
C --------------------------------------
B( I) = B( I) - F( I)
JS = JL + 1
C
500 CONTINUE
C ------------------------------------------------
C
C ------------------------------------------------------------
C
CALL SOLV ( B, P, DD, NCW, KCC, NBP, NNP, NNZ, NZ2 )
C
C ------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SJAC
C
C ----- Gaussian Weghting Factors ---
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
COMMON /XE/ XI( 4), ETA( 4)
COMMON /PW/ PG( 8), WG ( 8), NPW
C
C ----- Quadrilateral element for ξ-η coordinates ---
C
XI ( 1) = - 1.D0
ETA( 1) = - 1.D0
XI ( 2) = 1.D0
ETA( 2) = - 1.D0
C
XI ( 3) = 1.D0
ETA( 3) = 1.D0
XI ( 4) = - 1.D0
ETA( 4) = 1.D0
C
C ----- Weighting coefficients for Gaussian integlation ---
C
NPW = 2
C ============
C
IF ( NPW.EQ.2) GO TO 10
IF ( NPW.EQ.3) GO TO 20
IF ( NPW.EQ.4) GO TO 30
IF ( NPW.EQ.8) GO TO 40
C
10 CONTINUE
C
PG( 1) = - 0.577350269189626D0
PG( 2) = - PG( 1)
WG( 1) = 1.D0
WG( 2) = 1.D0
C
RETURN
C
20 CONTINUE
C
PG( 1) = - 0.774596669241483D0
PG( 2) = 0.D0
PG( 3) = - PG( 1)
C
WG( 1) = 5.D0 / 9.D0
WG( 2) = 8.D0 / 9.D0
WG( 3) = WG( 1)
C
RETURN
C
30 CONTINUE
C
PG( 1) = - 0.861136311594053D0
PG( 2) = - 0.339981043584856D0
PG( 3) = - PG( 2)
PG( 4) = - PG( 1)
C
WG( 1) = 0.347854845137454D0
WG( 2) = 0.652145154862546D0
WG( 3) = WG( 2)
WG( 4) = WG( 1)
RETURN
C
40 CONTINUE
C
PG( 1) = - 0.9602898564D0
PG( 2) = - 0.7966664774D0
PG( 3) = - 0.5255324099D0
PG( 4) = - 0.1834346424D0
C
PG( 5) = - PG( 4)
PG( 6) = - PG( 3)
PG( 7) = - PG( 2)
PG( 8) = - PG( 1)
C
WG( 1) = 0.1012285362D0
WG( 2) = 0.2223810344D0
WG( 3) = 0.3137066458D0
WG( 4) = 0.3626837833D0
C
WG( 5) = WG( 4)
WG( 6) = WG( 3)
WG( 7) = WG( 2)
WG( 8) = WG( 1)
RETURN
END
C **********************************************************************
C
SUBROUTINE STNC ( NCC, NDE, NP, P, NNP, NEL, NNZ )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION P( NNZ), NCC( 2,NNZ), NDE( NEL,4)
C
K = 0
C --------------------------------------------------
DO 100 IN = 1, NEL
DO 100 I = 1, NP
C --------------------------------------------------
C
II = NDE( IN,I)
C --------------------------------------------------
DO 100 J = 1, NP
C
JJ = NDE( IN,J)
CALL DCNL ( 0.D0, II, JJ, P, NCC, K, NNZ )
C
100 CONTINUE
C --------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SNCW ( XI, W, NC1, NC2, NC3, NDE, NNP, NEL )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION XI ( 4,2,2), W ( 4,2), NDE( NEL,4), NC1( NNP),
1 NC2( NNP,10), NC3( NNP,10)
C
C --------------------------------------
DO 100 I = 1, NNP
C
NC1( I) = 0
C
100 CONTINUE
C --------------------------------------
DO 200 I = 1, NNP
DO 200 J = 1, 10
C
NC2( I,J) = 0
NC3( I,J) = 0
C
200 CONTINUE
C --------------------------------------
DO 300 I = 1, NEL
DO 300 J = 1, 4
C
N = NDE( I,J)
NC1( N) = NC1( N) + 1
NC2( N,NC1( N)) = I
NC3( N,NC1( N)) = J
C
300 CONTINUE
C --------------------------------------
C
C ----- Two-point quadrature --- Attention
C
XI( 1,1,2) = - 0.577350269189626D0
XI( 2,1,2) = - XI( 1,1,2)
XI( 3,1,2) = XI( 1,1,2)
XI( 4,1,2) = XI( 2,1,2)
C
XI( 1,2,2) = - 0.577350269189626D0
XI( 2,2,2) = XI( 1,2,2)
XI( 3,2,2) = - XI( 1,2,2)
XI( 4,2,2) = - XI( 1,2,2)
C
W( 1,2) = 1.D0
W( 2,2) = 1.D0
W( 3,2) = 1.D0
W( 4,2) = 1.D0
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SOLV ( FVC, C, TSM, NCC, KCC, NBC, NNP, NNZ, NZ2 )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
COMMON /CST2/ ROP, ERO, ITR, NIT
C
DIMENSION FVC( NNP), C ( NNP), KCC( NNP), NBC( NNP),
1 TSM( NNZ), NCC( 2,NZ2)
C
C ----- S.O.R. Method ---
C
CMX = 0.D0
C --------------------------------------
DO 100 I = 1, NNP
C
AC = DABS( C( I))
IF ( AC.GE.CMX) CMX = AC
C
100 CONTINUE
C --------------------------------------
C
EPS = CMX* 0.01D0
LMR = 3* NNP
KTR = 2* NNP
ITR = 0
C
C ----------------------------------------------------------
900 ITR = ITR + 1
ERX = 0.D0
C
C ------------------------------------------------
DO 200 I = 1, NNP
C
IF ( NBC( I).EQ.9999) GO TO 200
IF ( I.EQ.1) GO TO 250
MS = KCC( I-1) + 1
GO TO 270
C
250 MS = 1
270 ML = KCC( I)
IF ( MS.GT.ML) GO TO 1000
DK = FVC( I)
C
C --------------------------------------
DO 300 M = MS, ML
C
M1 = NCC( 1,M)
M2 = NCC( 2,M)
IF ( M2.EQ.I) DSM = 1.D0/ TSM( M1)
IF ( M2.EQ.I) GO TO 300
DK = DK - TSM( M1)* C( M2)
C
300 CONTINUE
C --------------------------------------
C
DC = - C( I) + DSM* DK
C( I) = C( I) + ROP* DC
C ==============================
C
TEMP = DABS( C( I))
IF ( TEMP.LE.EPS) GO TO 200
TEMP = DABS( DC/ C( I))
C
IF ( TEMP.LE.ERX) GO TO 200
ERX = TEMP
NMX = I
C
200 CONTINUE
C ------------------------------------------------
IF ( ERX.LE.ERO ) GO TO 710
IF ( ITR.GE.KTR ) GO TO 701
IF ( ITR.GE.LMR ) GO TO 800
GO TO 900
C ----------------------------------------------------------
C
800 CONTINUE
WRITE(6,2000) EPS, ITR, NMX, ERX
2000 FORMAT('SOLV'' EPS=',E11.4,' Iter=', I6,2X,' Node No.=', I5,2X,
1 ' Error =',E11.4)
WRITE(6,2100) ( C( I), I = 1, NNP )
2100 FORMAT(' ',6E18.6)
GO TO 900
C
701 CONTINUE
WRITE(6,2300) ITR
2300 FORMAT('0',' *** No. of iterations exceeded, Iter =',I4)
GO TO 710
C
1000 CONTINUE
WRITE(6,2400) I, MS, ML
2400 FORMAT(' * SOLV * Error * ',3I4)
710 CONTINUE
C ----------------------
C
NIT = NIT + ITR
C
C ----------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SUBS ( PDV, I, J, N, PD, NCC, K, NNZ )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION PD( NNZ), NCC( 2,NNZ)
C
K = K + 1
IF ( K.GT.NNZ) WRITE(6,2000)
2000 FORMAT(' * Error K.GT.NNZ *')
C
C --------------------------------------
JE = K - N
C
DO 100 JJ = 1, JE
C
J1 = K - JJ + 1
PD ( J1) = PD ( J1-1)
NCC( 1,J1) = NCC( 1,J1-1)
NCC( 2,J1) = NCC( 2,J1-1)
C
100 CONTINUE
C --------------------------------------
C
PD ( N) = PDV
NCC( 1,N) = I
NCC( 2,N) = J
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VLPR ( X, Y, P, DLP, DI, NDE, NC3, NC2,
1 NC1, SU, SV, NBU, NBV, DLX, DLY, NNP,
2 NEL )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION X ( NNP), Y ( NNP), P ( NNP), DLP( NNP),
1 DI ( NNP), NDE( NEL,4), NC3( NNP,10), NC2( NNP,10),
2 NC1( NNP), SU ( NNP), SV ( NNP), NBU( NNP),
3 NBV( NNP), DLX( NEL,4), DLY( NEL,4)
C
COMMON /CNST/ DT, COF, PHI, RE
C
C ----------------------------------------------------------
C
CALL DLXY ( X, Y, DLP, NDE, DLX, DLY, NNP, NEL )
C
C ---------------------------------------------------------------
C
CALL VLND ( DI, DLX, NC1, NC2, NC3, X, Y, NDE, NNP, NEL )
C
C ---------------------------------------------------------------
DO 100 I = 1, NNP
C
IF ( NBU( I).NE.9999) SU( I) = SU( I) - DT* DI( I)
100 CONTINUE
C ---------------------------------------------------------------
C
CALL VLND ( DI, DLY, NC1, NC2, NC3, X, Y, NDE, NNP, NEL )
C
C ---------------------------------------------------------------
DO 200 I = 1, NNP
C
IF ( NBV( I).NE.9999) SV( I)= SV( I) - DT* DI( I)
200 CONTINUE
C
C ----- Pressure -----------------
C
DO 300 I = 1, NNP
C
P( I) = P( I) + DLP( I)
C
300 CONTINUE
C --------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VLND ( CI, CE, NC1, NC2, NC3, X, Y, NDE, NNP, NEL )
C
IMPLICIT REAL*8 ( A-H,O-Z )
C
DIMENSION CI ( NNP), CE( NEL,4), NC1( NNP), NC2( NNP,10),
1 NC3( NNP,10), X ( NNP), Y ( NNP), NDE( NEL,4)
C
DIMENSION SNNP( 10)
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
NQ = NC1( I)
CT = 0.D0
SPL = 0.D0
C
C --------------------------------------
DO 150 J = 1, NQ
C
JJ = NC2( I,J)
JI = NC3( I,J)
C
SNNP( J)
C
1 = ( X( NDE( JJ,2))* Y( NDE( JJ,4))
2 + X( NDE( JJ,1))* Y( NDE( JJ,2))
3 + X( NDE( JJ,4))* Y( NDE( JJ,1))
C
4 - X( NDE( JJ,1))* Y( NDE( JJ,4))
5 - X( NDE( JJ,2))* Y( NDE( JJ,1))
6 - X( NDE( JJ,4))* Y( NDE( JJ,2))
C
7 + X( NDE( JJ,2))* Y( NDE( JJ,3))
8 + X( NDE( JJ,3))* Y( NDE( JJ,4))
9 + X( NDE( JJ,4))* Y( NDE( JJ,2))
C
A - X( NDE( JJ,2))* Y( NDE( JJ,4))
B - X( NDE( JJ,3))* Y( NDE( JJ,2))
C - X( NDE( JJ,4))* Y( NDE( JJ,3)))/ 2.D0
C
CT = CT + CE( JJ,JI)* SNNP( J)
SPL = SPL + SNNP( J)
C
150 CONTINUE
C --------------------------------------
C
CI( I) = CT/ SPL
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************