戻る

C ******************************************************************** C
C C
C H F 2 D - F E M C
C C
C ( V.5.0 ) C
C C
C Two-Dimensional Thermal Fluid Analysis by F.E.M. C
C C
C ( U-V-P-C / For Quardrilateral Element ) C
C ======================== 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 ----- ( A ) Unstructured meshes ( 30 x 30 ) --------------------------
C
PARAMETER ( NNP= 969, NEL=888, N=3632, NNZ= NNP+N, NZ2= NNP+2*N )
C
C ----- ( B ) Structured meshes ----------------------------------------
C
** PARAMETER ( NEX = 50, NEY = 50, ! Attention " CMNU"
C ===================
C
** 1 LYX = NEX + 1, LYY = NEY + 1,
** 2 NNP = ( NEX+1)* ( NEY+1), NEL = NEX* NEY,
** 3 N = 4* NEX* NEY + NEX + NEY,
** 4 NNZ = NNP + N, NZ2 = NNP + 2* N )
C
** DIMENSION XLX( LYX), YLY( LYY) ! Attention
C
C **********************************************************************
C
DIMENSION U ( NNP), V ( NNP), T ( NNP), P ( NNP),
1 X ( NNP), Y ( NNP),
2 B ( NNP), F ( NNP), SU ( NNP), SV ( NNP),
3 ST ( NNP), GU ( NNP), GV ( NNP), DI ( NNP),
4 U1 ( NNP), V1 ( NNP), DLP( NNP), DPX( NEL,4),
5 DPY( NEL,4), NC1( NNP), NC2( NNP,10), NC3( NNP,10)
C
DIMENSION FF ( NNP), GG ( NNP), AA( 10,NNP), DD ( NNZ),
1 PP ( NNZ), D ( NNZ), A ( NNZ), HHX( NNP),
2 HHY( NNP), NDE( NEL,4)
C
DIMENSION KCC( NNP), NCC( 2,NNZ), NBU( NNP), NBV( NNP),
1 NBP( NNP), NBT( NNP), NCW( 2,NZ2), NBA( NNP)
C
COMMON /DEF/ XI( 4,2), WS( 4), EPP( 4,4), EDD( 4,4)
C
COMMON /CS2/ XMN, NPR, KSR
COMMON /CS1/ OMG, ERR, NTR, ITT
COMMON /CS3/ UVX, UMX, VMX, RA, EPS, EPM
COMMON /ABC/ HGN, RLW, ANU, HNX, HNI, RNX, RNI
COMMON /CST/ DT, VCF, PHI, TDF, FT, ALF, TB, LPE, LOP
C
C ----- Attention ( LOP < = 100,000 ) ---
C
COMMON /VRP/ VRV( 100000), VNA( 100000), VNB( 100000),
1 VDN( 100000)
C
CALL TTLE
C
C ***** Data ( 2 ) *****************************************************
C
C ----- ( A ) Unstructured meshes ( 30 X 30 ) --------------------------
C
OPEN ( 10, FILE='HF2D_FEM.DAT' )
C
KSR = 2 ! Unstructured meshes
C
CALL MSHG ( X, Y, NDE, NNP, NEL ) ! KSR = 2
C
* OPEN ( 7, FILE='T_HF2D_FEM_UE30_RA7100_R' ) ! DT = 0.001
* OPEN ( 8, FILE='T_HF2D_FEM_UE30_RA7100_Nu' )
* OPEN ( 9, FILE='T_HF2D_FEM_UE30_RA7100_OUT')
C
* OPEN ( 7, FILE='T_HF2D_FEM_UE30_RA10_5_R' ) ! DT = 0.0001
* OPEN ( 8, FILE='T_HF2D_FEM_UE30_RA10_5_Nu' )
* OPEN ( 9, FILE='T_HF2D_FEM_UE30_RA10_5_OUT')
C
OPEN ( 7, FILE='HF2D_FEM_UE30_RA10_5_R' ) ! DT = 0.0001
OPEN ( 8, FILE='HF2D_FEM_UE30_RA10_5_Nu' )
OPEN ( 9, FILE='HF2D_FEM_UE30_RA10_5_OUT')
C
C ----- ( B ) Structured meshes ----------------------------------------
C
** KSR = 1 ! Structured meshes
C
C ----- ( 1 ) Mesh with unequaled length -------------------
C

C ----- Comp. of Nu : OK ---
C
** CALL PRM50F ( XLX, YLY, LYX, LYY, N, NNZ, NZ2 )
C
C ----------------------------------------------------------
C
* CALL PRM4F ( XLX, YLY, LYX, LYY, N, NNZ, NZ2 )
C
* CALL PRM10F ( XLX, YLY, LYX, LYY, N, NNZ, NZ2 )
C
* CALL PRM20F ( XLX, YLY, LYX, LYY, N, NNZ, NZ2 )
C
* CALL PRM30F ( XLX, YLY, LYX, LYY, N, NNZ, NZ2 )
C
* CALL PRM40F ( XLX, YLY, LYX, LYY, N, NNZ, NZ2 )
C
* CALL PRM100F ( XLX, YLY, LYX, LYY, N, NNZ, NZ2 )
C
* CALL PRM200F ( XLX, YLY, LYX, LYY, N, NNZ, NZ2 )
*
* NNP , NEL , N NNP , NEL , N
* 10F 121 , 100 , 420 40F 1681 , 1600 , 6480
* 20F 441 , 400 , 1640 50F 2601 , 2500 , 10100
* 30F 961 , 900 , 3660
*
* ----- Mesh with equaled length ---------------------------
*
* CALL PRM10H ( XLX, YLY )
*
* CALL PRM20H ( XLX, YLY )
C
C ----- ( 2 ) ----------------------------------------------
C
** CALL SMSH ( X, Y, NDE, XLX, YLY, LYX, LYY, NNP, NEL )
C
C ----- ( 3 ) ----------------------------------------------
C
* OPEN ( 7, FILE='HF2D_FEM_M4_RA710_R' )
* OPEN ( 8, FILE='HF2D_FEM_M4_RA710_Nu' )
* OPEN ( 9, FILE='HF2D_FEM_M4_RA710_OUT')
C
* OPEN ( 7, FILE='T_HF2D_FEM_M4_RA710_R' )
* OPEN ( 8, FILE='T_HF2D_FEM_M4_RA710_Nu' )
* OPEN ( 9, FILE='T_HF2D_FEM_M4_RA710_OUT')
C
C ----------------------------------------------------------
C
* OPEN ( 7, FILE='HF2D_FEM_M50_RA710_R' ) ! DT = 0.001
* OPEN ( 8, FILE='HF2D_FEM_M50_RA710_Nu')
* OPEN ( 9, FILE='HF2D_FEM_M50_RA710_OUT')
C
* OPEN ( 7, FILE='HF2D_FEM_M50_RA7100_R' ) ! DT = 0.001
* OPEN ( 8, FILE='HF2D_FEM_M50_RA7100_Nu')
* OPEN ( 9, FILE='HF2D_FEM_M4_RA710_OUT' )
C
** OPEN ( 7, FILE='HF2D_FEM_M50_RA10_5_R' ) ! DT = 0.0001
** OPEN ( 8, FILE='HF2D_FEM_M50_RA10_5_Nu')
** OPEN ( 9, FILE='HF2D_FEM_M50_RA10_5_OUT')
C
* OPEN ( 7, FILE='HF2D_FEM_M50_RA10_6_R' ) ! DT = 0.000 01'
* OPEN ( 8, FILE='HF2D_FEM_M50_RA10_6_Nu')
* OPEN ( 9, FILE='HF2D_FEM_M50_RA10_6_OUT')
C
* OPEN ( 7, FILE='HF2D_FEM_M50_RA10_7_R' ) ! DT = 0.000 001
* OPEN ( 8, FILE='HF2D_FEM_M50_RA10_7_Nu')
C
C **********************************************************************
C
CALL INPT ( X, Y, U, V, T, P, DD, NDE, PP,
1 NBU, NBV, NBT, NBP, NBA, NCC, NCW, KCC, ISL,
2 NNP, NEL, NNZ, NZ2 )
C
CALL INIT ( SU, SV, ST, DPX, DPY, DLP, NNP, NEL )
C
CALL SNCN ( NDE, NC1, NC2, NC3, NNP, NEL )
C
WRITE(6,*) '*** Computation started ***'
WRITE(6,*) ' '
C
LOP = 0
IYF = 0
C --------------------------------------------------------------------
999 CONTINUE
LOP = LOP + 1
ITT = 0
C
CALL FGEM ( FF, GG, X, Y, TDF, U, V, NC1, NC2, NC3,
1 NDE, AA, HHX, HHY, NNP, NEL )
C
CALL STUV ( U, V, GU, GV, NNP )
C
C ----- Temperature ---
C
CALL UVTC ( X, Y, U, V, T, T, NDE, P, PP, A,
1 B, D, F, 0, ST, NBT, NCC, NCW, KCC, ISL,
2 FF, GU, GV, NBA, NNP, NEL, NNZ, NZ2 )
C
CALL FGEM ( FF, GG, X, Y, VCF, U, V, NC1, NC2, NC3,
1 NDE, AA, HHX, HHY, NNP, NEL )
C
C ----- Pseudo-Velocities ---
C
CALL CNC1 ( X, Y, U, V, ST, U, NDE, P, PP, A,
1 B, D, F, 1, U1, NBU, NCC, NCW, KCC, NNP,
2 NEL, NNZ, NZ2 )
C
CALL CNC1 ( X, Y, U, V, ST, V, NDE, P, PP, A,
1 B, D, F, 2, V1, NBV, NCC, NCW, KCC, NNP,
2 NEL, NNZ, NZ2 )
C
C ----- Predicted Pressure ---
C
CALL POIS ( X, Y, U1, V1, P, NDE, A, B, D, F,
1 DD, NCC, NCW, KCC, NBP, NNP, NEL, NNZ, NZ2 )
C
C ----- Predicted Velocity "u" ---
C
CALL UVTC ( X, Y, U, V, ST, U, NDE, P, PP, A,
1 B, D, F, 1, SU, NBU, NCC, NCW, KCC, ISL,
2 FF, GU, GV, NBA, NNP, NEL, NNZ, NZ2 )
C
C ----- Predicted Velocity "v" ---
C
CALL UVTC ( X, Y, U, V, ST, V, NDE, P, PP, A,
1 B, D, F, 2, SV, NBV, NCC, NCW, KCC, ISL,
2 FF, GU, GV, NBA, NNP, NEL, NNZ, NZ2 )
C
C ----- For "DLP" ---
C
CALL POIS ( X, Y, SU, SV, DLP, NDE, A, B, D, F,
1 DD, NCC, NCW, KCC, NBP, NNP, NEL, NNZ, NZ2 )
C
C ----- For Pressure ---
C
CALL PRSS ( P, DLP, NNP )
C
C ----- For Velocities ---
C
CALL VPRS ( X, Y, P, DLP, DI, NDE, NC1, NC2, NC3,
1 SU, SV, NBU, NBV, DPX, DPY, NNP, NEL )
C
C ----- Final routine ---
C
CALL FSBR ( X, Y, U, V, T, P, DLP, SU, SV,
1 ST, DI, NC1, NC2, NC3, NBU, NBV, DPX, DPY,
2 NDE, ISL, LYX, LYY, NNP, NEL, IYF, *999 )
C
CALL PLTD ( U, V, T, FF, NNP )
C
CLOSE ( 7)
CLOSE ( 8)
CLOSE ( 9)
C
STOP
END
C **********************************************************************
C
SUBROUTINE CMNU ( X, Y, T, LYX, LYY, NNP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CS2/ XMN, NPR, KSR
COMMON /ABC/ HGN, RLW, ANU, HNX, HNI, RNX, RNI
COMMON /VRP/ VRV( 100000), VNA( 100000), VNB( 100000),
1 VDN( 100000)
C
DIMENSION X( NNP), Y( NNP), T( NNP)
C
DIMENSION DDT( 100), DT1( 100), SSY( 100), SYY( 100),
1 XX ( 6,100), YY ( 6,100), TT ( 6,100)
C
C ----- Computation for Nu- Value ---
C
C ( Attention : Str. Meshes of 50 x 50 or Unstr.Meshes of 30 x 30 ONLY )
C ==================================================================
C
C === > Attention : Definition of computation nodes ===
C
C ----- For Structured meshes ( KSR = 1 ) of 50 X 50 -------
C
IF ( KSR.EQ.1 .AND. LYX.NE.51) THEN
ANU = 0.D0
WRITE(6,*) ' *** Comp. for Nu is omitted ***'
GO TO 990
END IF
C
IF ( KSR.EQ.1) CALL HDSB ( X, Y, T, XX, YY, TT, NNP )
C
C ----- For Unstructured meshes ( KSR = 2 ) of 30 X 30 -----
C
IF ( KSR.EQ.2) THEN
C
CALL HDSR ( X, Y, T, XX, YY, TT, NNP )
C
HNX = - 10000.D0
HNI = 10000.D0
RNX = - 10000.D0
RNI = 10000.D0
C
NEX = 30 ! Un-structured meshes
NEY = 30 ! Un-structured meshes
LYX = NEX + 1
LYY = NEY + 1
C
END IF
C ----------------------------------------------------------
JJ = 1
C
DO 100 I = 1, LYX
C
H1 = XX( 2,I) - XX( 1,I) ! ( Left:3points )
H2 = XX( 3,I) - XX( 2,I)
H3 = XX( 6,I) - XX( 5,I) ! ( Right: 3 points )
H4 = XX( 5,I) - XX( 4,I)
C
A = - ( 2.D0* H1 + H2)/( H1* ( H1 + H2))
B = ( H1 + H2)/( H1* H2)
C = - ( H1/( H2* ( H1 + H2)))
A1 = ( 2.D0* H3 + H4)/( H3* ( H3 + H4))
B1 = - ( H3 + H4)/( H3* H4)
C1 = ( H3/( H4* ( H3 + H4)))
C
C ----- Second Derivatives ---
C
DDT( I) = A * TT( 1,I) + B * TT( 2,I) + C * TT( 3,I)
DT1( I) = A1* TT( 6,I) + B1* TT( 5,I) + C1* TT( 4,I)
C
C ---- First Derivatives ---
C
JJ = JJ + LYX
HNX = DMAX1( HNX,- DDT( I))
HNI = DMIN1( HNI,- DDT( I))
RNX = DMAX1( RNX,- DT1( I))
RNI = DMIN1( RNI,- DT1( I))
C
100 CONTINUE
C ----------------------------------------------------------
C
J = 1
HGN = 0.D0
RLW = 0.D0
C
C ----- Second-orderd Integration ---
C
DO 200 I = 1, LYY-2, 2
C
Y0 = 0.D0
Y1 = YY( 1,I+1) - YY( 1,I) ! + Value
Y2 = YY( 1,I+2) - YY( 1,I) ! + Value
C
C ----- Attention : Important --- CHECK ---
C
*** Y1 = - Y1
*** Y2 = - Y2
C
SSY1 = DDT( I)* ( - Y2** 3 + 3.D0* Y1* Y2** 2 - 2.D0* Y0** 3
1 + 3.D0* Y0** 2* Y2 + 3.D0* Y0** 2* Y1 - 6.D0* Y0* Y1* Y2 )
2 /( 6.D0* ( Y0-Y1)* ( Y0-Y2))
C
SSY2 = DDT( I+1)* ( - Y2** 3 + Y0** 3 + 3.D0* Y0* Y2** 2
1 - 3.D0* Y0** 2* Y2)/( 6.D0* ( Y1-Y0)* ( Y1-Y2))
C
SSY3 = DDT( I+2)* ( 2.D0* Y2** 3 + Y0** 3 - 3.D0* Y1* Y2** 2
1 - 3.D0* Y0* Y2** 2 + 6.D0* Y0* Y1* Y2 - 6.D0* Y0** 2* Y1 )
2 /( 6.D0* ( Y2-Y0)* ( Y2-Y1))
C
C ----------------------------------------
C
SY11 = DT1( I)* ( - Y2** 3 + 3.D0* Y1* Y2** 2 - 2.D0* Y0** 3
1 + 3.D0* Y0** 2* Y2 + 3.D0* Y0** 2* Y1 - 6.D0* Y0* Y1* Y2 )
2 /( 6.D0*( Y0-Y1)* ( Y0-Y2))
C
SY22 = DT1( I+1)* ( - Y2** 3 + Y0** 3 + 3.D0* Y0* Y2** 2
1 - 3.D0* Y0** 2* Y2 )/( 6.D0* ( Y1-Y0)* ( Y1-Y2))
C
SY33 = DT1( I+2)* ( 2.D0* Y2** 3 + Y0** 3 - 3.D0* Y1* Y2** 2
1 - 3.D0* Y0* Y2** 2 + 6.D0* Y0* Y1* Y2 - 6.D0* Y0** 2* Y1 )
2 /( 6.D0* ( Y2-Y0)* ( Y2-Y1))
C
SSY( J) = SSY1 + SSY2 + SSY3
SYY( J) = SY11 + SY22 + SY33
HGN = HGN + SSY( J)
RLW = RLW + SYY( J)
J = J + 1
C
200 CONTINUE
C
ANU = DABS( HGN + RLW)/ 2.D0
C
990 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CNC1 ( X, Y, U, V, T, C, NDE, P,
1 PP, A, B, D, F, ID, CNW, NBC,
2 NCC, NCW, KCC, NNP, NEL, NNZ, NZ2 )
C
C ----- MGM ( Matrix A & Vector B ) ---
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X ( NNP), Y ( NNP), U ( NNP), V ( NNP),
1 T ( NNP), C ( NNP), NDE( NEL,4), P ( NNP),
2 PP( NNZ), A ( NNZ), B ( NNP), D ( NNZ),
3 F ( NNP), CNW( NNP), NBC( NNP), NCC( 2,NNZ),
4 NCW( 2,NZ2), KCC( NNP)
C
DIMENSION EDY( 4,4), EFX( 4,4), EFY( 4,4), EDX( 4,4),
1 EPP( 4,4), PE ( 4), XX ( 4), YY ( 4),
2 TT ( 4)
C
COMMON /CST/ DT, VCF, PHI, TDF, FT, ALF, TB, LPE, LOP
C
IF ( ID.EQ.1) FH = 0.D0
IF ( ID.EQ.2) FH = 1.D0
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 ------------------------------------------------
DO 300 IEL = 1, NEL
C
UE = 0.D0
VE = 0.D0
C
C --------------------------------------
DO 350 I = 1, 4
C
N = NDE( IEL,I)
XX( I) = X( N)
YY( I) = Y( N)
PE( I) = P( N)
TT( I) = T( N)
UE = UE + U( N)/ 4.D0
VE = VE + V( N)/ 4.D0
C
350 CONTINUE
C --------------------------------------
C
C ----- D, F ------------------------------------------
C
CALL JPDF ( XX, YY, EDX, EDY, EFX, EFY, EPP )
C
C -----------------------------------------------------
DO 300 I = 1, 4
DO 300 J = 1, 4
C ------------------------------------------------
C
II = NDE( IEL,I)
JJ = NDE( IEL,J)
DS = VCF* ( EDX( I,J) + EDY( 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)
1 - FH* EPP( I,J)* TT( J)* FT
C
300 CONTINUE
C -----------------------------------------------
C
JS = 1
C ------------------------------------------------
DO 400 I = 1, NNP
C
JL = KCC( I)
C -------------------------------------
DO 450 J = JS, JL
C
KPR = NCW( 1,J)
KCL = NCW( 2,J)
A( KPR)= PP( KPR)/ DT
B( I) = B( I) + ( PP( KPR)/ DT - D( KPR))* C( KCL)
C
450 CONTINUE
C -------------------------------------
B( I) = B( I) - F( I)
JS = JL + 1
C
400 CONTINUE
C -------------------------------------------------------------
C
CALL SLVE( B, CNW, A, NCW, KCC, NBC, NNP, NNZ, NZ2 )
C
C -------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CPNU ( X, Y, T, HGN, RLW, ANU, HNX,
1 HNI, RNX, RNI, NDE, NC1, NC2, NC3,
2 NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X ( NNP), Y ( NNP), T ( NNP), DDT( 100),
1 DT1( 100), SSY( 100), SYY( 100), NDE( NEL,4),
2 NC1( NNP), NC2( NNP,10), NC3( NNP,10)
C
DIMENSION YA( 100,3), NA( 100,3), YA2( 100,3), NA2( 100,3)
C
HNX = -10000.D0
HNI = 10000.D0
RNX = -10000.D0
RNI = 10000.D0
C
JJ = 1
C ------------------------------------------
DO 100 I = 1, NNP
C
IF ( DABS( X( I)).LE.0.001D0) THEN
NA( JJ,1) = I
YA( JJ,1) = Y( I)
JJ = JJ + 1
ENDIF
C
100 CONTINUE
C ------------------------------------------
J1 = JJ - 1
DO 150 I = 1, J1-1
DO 150 J = I+1, J1
C ------------------------------------------
C
IF ( YA( I,1).GT.YA( J,1)) THEN
YA1 = YA( I,1)
NA1 = NA( I,1)
YA( I,1) = YA( J,1)
NA( I,1) = NA( J,1)
YA( J,1) = YA1
NA( J,1) = NA1
END IF
C
150 CONTINUE
C ------------------------------------------
DO 200 I = 1, J1-1
KEY = NC1( NA( I,1))
C -------------------------------------------
DO 200 J = 1, KEY
C
IF ( NDE( NC2( NA( I,1),J),1).EQ.NA( I,1))
1 NA( I,2) = NDE( NC2( NA( I,1),J),2)
C
200 CONTINUE
C ------------------------------------------
NA( J1,2) = NDE( NC2( NA( J1,1),1), 3)
DO 300 I = 1, J1-1
KEY = NC1( NA( I,2))
DO 300 J = 1, KEY
C -------------------------------------------
IF ( NDE( NC2( NA( I,2),J), 1).EQ.NA( I,2))
1 NA( I,3) = NDE( NC2( NA( I,2),J), 2)
C
300 CONTINUE
C -------------------------------------------
KEY = NC1( NA( J1,2))
DO 350 I = 1, KEY
C
IF ( NDE( NC2( NA( J1,2),I), 4).EQ.NA( J1,2))
1 NA( J1,3) = NDE( NC2( NA( J1,2),I), 3)
350 CONTINUE
C -------------------------------------------
JJ = 1
DO 400 I = 1, NNP
C
IF ( DABS( X( I)-1.D0).LE.0.001D0) THEN
NA2( JJ,1) = I
YA2( JJ,1) = Y( I)
JJ = JJ + 1
ENDIF
C
400 CONTINUE
C -------------------------------------------
J2 = JJ - 1
DO 450 I = 1, J2-1
DO 450 J = I+1, J2
C
IF ( YA2( I,1).GT.YA2( J,1))THEN
YA3 = YA2( I,1)
NA3 = NA2( I,1)
YA2( I,1) = YA2( J,1)
NA2( I,1) = NA2( J,1)
YA2( J,1) = YA3
NA2( J,1) = NA3
END IF
C
450 CONTINUE
C -------------------------------------------
DO 500 I = 1, J2-1
KEY = NC1( NA2( I,1))
DO 500 J = 1, KEY
C --------------------------------------------
C
IF ( NDE( NC2( NA2( I,1),J), 2).EQ.NA2( I,1))
1 NA2( I,2) = NDE( NC2( NA2( I,1),J), 1)
C
500 CONTINUE
C --------------------------------------------
NA2( J2,2) = NDE( NC2( NA2( J2,1),1), 4)
C
DO 600 I = 1, J2-1
KEY = NC1( NA2( I,2))
C
DO 600 J = 1, KEY
C -------------------------------------------
C
IF ( NDE( NC2( NA2( I,2),J),2).EQ.NA2( I,2))
1 NA2( I,3) = NDE( NC2( NA2( I,2),J), 1)
C
600 CONTINUE
C -------------------------------------------
KEY = NC1( NA2( J2,2))
DO 650 I = 1, KEY
C
IF ( NDE( NC2( NA2( J2,2),I),3).EQ. NA2( J2,2))
1 NA2( J2,3) = NDE( NC2( NA2( J2,2),I), 4 )
650 CONTINUE
C --------------------------------------------
DO 700 I = 1, J1
C
H1 = X( NA( I,2)) - X( NA( I,1))
H2 = X( NA( I,3)) - X( NA( I,2))
A = - ( 2.D0* H1 + H2)/( H1* ( H1 + H2))
B = ( H1 + H2)/( H1* H2 )
C = - ( H1/( H2* ( H1 + H2)))
C
C ------ Second Derivatives ---
C
DDT( I) = A* T( NA( I,1)) + B* T( NA( I,2)) + C* T( NA( I,3))
HNX = DMAX1( HNX, - DDT( I))
HNI = DMIN1( HNI, - DDT( I))
C
700 CONTINUE
C --------------------------------------------
DO 750 I = 1, J2
C
H3 = X( NA2( I,1)) - X( NA2( I,2))
H4 = X( NA2( I,2)) - X( NA2( I,3))
A1 = ( 2.D0* H3 + H4)/( H3* ( H3 + H4))
B1 = - ( H3 + H4)/( H3* H4)
C1 = ( H3/( H4* ( H3 + H4)))
C
C ------ Second Derivatives ---
C
DT1( I)= A1* T( NA2( I,1)) + B1* T( NA2( I,2)) + C1* T( NA2( I,3))
RNX = DMAX1( RNX, - DT1( I))
RNI = DMIN1( RNI, - DT1( I))
C
750 CONTINUE
C -------------------------------------------
C
C ----- Second-orderd Integration ---
C
J = 1
HGN = 0.D0
C -------------------------------------------
DO 800 I = 1, J1-2, 2
C
Y0 = 0.D0
Y1 = Y( NA( I+1,1)) - Y( NA( I,1))
Y2 = Y( NA( I+2,1)) - Y( NA( I,1))
C
SSY1 = DDT( I)* ( - Y2** 3 + 3.D0* Y1* Y2** 2 - 2.D0* Y0** 3
1 + 3.D0* Y0** 2* Y2 + 3.D0* Y0** 2* Y1 - 6.D0* Y0* Y1* Y2 )
2 /( 6.D0* ( Y0-Y1)* ( Y0-Y2))
C
SSY2 = DDT( I+1)* ( - Y2** 3 + Y0** 3 + 3.D0* Y0* Y2** 2
1 - 3.D0* Y0** 2* Y2 )/( 6.D0* ( Y1-Y0)* ( Y1-Y2))
C
SSY3 = DDT( I+2)* ( 2.D0* Y2** 3 + Y0** 3 - 3.D0* Y1* Y2** 2
1 - 3.D0* Y0* Y2** 2 + 6.D0* Y0* Y1* Y2 - 6.D0* Y0** 2* Y1)
2 /( 6.D0* ( Y2-Y0)* ( Y2-Y1))
C
SSY( J) = SSY1 + SSY2 + SSY3
HGN = HGN + SSY( J)
J = J + 1
C
800 CONTINUE
C -------------------------------------------
C
J = 1
RLW = 0.D0
C -------------------------------------------
DO 850 I = 1, J2-2, 2
C
Y0 = 0.D0
Y1 = Y( NA2( I+1,1)) - Y( NA2( I,1))
Y2 = Y( NA2( I+2,1)) - Y( NA2( I,1))
C
SY11= DT1( I)* ( - Y2** 3 + 3.D0* Y1* Y2** 2 - 2.D0* Y0** 3
1 + 3.D0* Y0** 2* Y2 + 3.D0* Y0** 2* Y1 - 6.D0* Y0* Y1* Y2 )
2 /( 6.D0* ( Y0-Y1 )* ( Y0-Y2 ))
C
SY22 = DT1( I+1)* ( - Y2** 3 + Y0** 3 + 3.D0* Y0* Y2** 2
1 - 3.D0* Y0** 2* Y2 )/( 6.D0* ( Y1-Y0 )* ( Y1-Y2 ))
C
SY33 = DT1( I+2)* ( 2.D0* Y2** 3 + Y0** 3 - 3.D0* Y1* Y2** 2
1 - 3.D0* Y0* Y2** 2 + 6.D0* Y0* Y1* Y2 - 6.D0* Y0** 2* Y1 )
2 /( 6.D0* ( Y2-Y0 )* ( Y2-Y1 ))
C
SYY( J) = SY11 + SY22 + SY33
RLW = RLW + SYY( J)
J = J + 1
C
850 CONTINUE
C -------------------------------------------
C
ANU = DABS( HGN + RLW)/ 2.D0
C -------------------------------------------
C
RETURN
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 : Currnt length of KCC, NCC
C NNP, NZ2 : Max. length of KCC, NCC
C
CHARACTER*4 ANNP, ANZ2
C
DATA ANNP/'NNP '/, ANZ2/'NZ2'/
C
NP = 1
C
C --------------------------------------
DO 100 I = 1, NZ
C
IF ( NP.EQ.NCC( 1,I)) GO TO 10
C ----- To set KCC ( NP) ---
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
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
N = I - 1
IS = KCC( N)
C
C --------------------------------------
DO 300 J = 1, N
C
JS = 1
IF ( J.GT.1) JS = KCC( J-1) + 1
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 ----- ( I,J) not exist ---
C
GO TO 300
C
20 CONTINUE
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 ----------------------------
NZ = NZ + 1
IF ( NZ.GT.NZ2) GO TO 40
C
C ----- To increase KCC ---
C
C ----------------------------
DO 600 L = I, NP
C
KCC( L) = KCC( L) + 1
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 ------------------------------------------------
C
RETURN
C
30 WRITE(6,2000) ANNP, NNP, ANNP, NP
STOP
C
40 WRITE(6,2000) ANZ2,NZ2, ANZ2, NZ
STOP
C
2000 FORMAT('** Error ( DCIJ) ',A3,'(',I5,') too small',A2,'=',I5 )
C
END
C **********************************************************************
C
SUBROUTINE DCNL ( PDV, I, J, PP, NCC, K, NNZ )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION PP( 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
10 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)
1 .AND.J.LT.NCC( 2,N))) 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)
1 .AND.J.LT.NCC( 2,N1))) 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 -----------------------------------------------------
C
CALL SUBS ( PDV, I, J, N, PP, NCC, K, NNZ )
C
C -----------------------------------------------------
RETURN
C
32 N = N1
GO TO 22
C
42 N = N2
22 CONTINUE
PP( N) = PP( N) + PDV
RETURN
C
23 CONTINUE
N1 = N
N = N1 + ( N2-N1)/ 2
C
IF ( N.NE.N1) GO TO 100
IF ( I.LT.NCC( 1,N2).OR.( I.EQ.NCC( 1,N2)
1 .AND.J.LT.NCC( 2,N2))) GO TO 41
IF ( I.EQ.NCC( 1,N2).AND. J.EQ.NCC( 2,N2)) GO TO 42
GO TO 43
41 N = N2
GO TO 33
C
43 K = K + 1
IF ( K.GT.NNZ) WRITE(6,2000)
2000 FORMAT('** Error DCNL **''NNZ'' in ''SOR'' is too small **')
C
PP ( K) = PDV
NCC( 1,K) = I
NCC( 2,K) = J
RETURN
C
100 CONTINUE
C ----------------------------------------------------------
WRITE(6,2100) K
2100 FORMAT(' * Error DCNL --DO LOP strange K =',I8 )
STOP
END
C **********************************************************************
C
SUBROUTINE DCON ( VAL, I, J, RMAT, NCC, NNZ, IC )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION RMAT( NNZ), NCC( 2,NNZ)
C
IC = 0
IF ( I.GT.J) RETURN
C
N1 = 1
N2 = NNZ
100 CONTINUE
IF (( N1+N2)/2 .EQ.IC) N2 = N2 + 1
IC = ( N1+N2)/ 2
IF ( NCC( 1,IC).LT.I.OR.
1 ( NCC( 1,IC).EQ.I.AND.NCC( 2,IC).LT.J)) THEN
N1 = IC
GO TO 100
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 100
C
END IF
C
RMAT( IC) = RMAT( IC) + VAL
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DPXY ( X, Y, C, NDE, DPX, DPY, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X ( NNP), Y ( NNP), C( NNP), NDE( NEL,4),
1 DPX( NEL,4), DPY( NEL,4)
C
DIMENSION CC ( 4), XX ( 4), YY ( 4),
1 ENN2( 4,4), EZG2( 2,4,4), EXY2( 2,4,4),
2 YAC2( 2,2,4), YIN2( 2,2,4)
C
COMMON /XE/ XIR( 4), ETA( 4)
C
C ---------------------------------------
DO 100 IEL = 1, NEL
DO 100 M = 1, 4
C ---------------------------------------
C
C ----------------------------
DO 200 I = 1, 4
C
N = NDE( IEL,I)
CC( I) = C( N)
XX( I) = X( N)
YY( I) = Y( N)
C
200 CONTINUE
C ----------------------------
C
C ----- Shape functions and its derivatives ---
C
ENN2( 1,M) = ( 1.D0 - XIR( M))* ( 1.D0 - ETA( M))/ 4.D0
ENN2( 2,M) = ( 1.D0 + XIR( M))* ( 1.D0 - ETA( M))/ 4.D0
ENN2( 3,M) = ( 1.D0 + XIR( M))* ( 1.D0 + ETA( M))/ 4.D0
ENN2( 4,M) = ( 1.D0 - XIR( M))* ( 1.D0 + ETA( M))/ 4.D0
C
EZG2( 1,1,M) = - ( 1.D0 - ETA( M))/ 4.D0
EZG2( 1,2,M) = ( 1.D0 - ETA( M))/ 4.D0
EZG2( 1,3,M) = ( 1.D0 + ETA( M))/ 4.D0
EZG2( 1,4,M) = - ( 1.D0 + ETA( M))/ 4.D0
C
EZG2( 2,1,M) = - ( 1.D0 - XIR ( M))/ 4.D0
EZG2( 2,2,M) = - ( 1.D0 + XIR ( M))/ 4.D0
EZG2( 2,3,M) = ( 1.D0 + XIR ( M))/ 4.D0
EZG2( 2,4,M) = ( 1.D0 - XIR ( M))/ 4.D0
C
C ----- Jacobians ------------
C
DO 400 I = 1, 2
DO 400 J = 1, 2
C
YAC2( I,J,M) = 0.D0
400 CONTINUE
C ----------------------------
DO 500 I = 1, 4
DO 500 K = 1, 2
C
YAC2( K,1,M) = YAC2( K,1,M) + EZG2( K,I,M)* XX( I)
YAC2( K,2,M) = YAC2( K,2,M) + EZG2( K,I,M)* YY( I)
C
500 CONTINUE
C ----------------------------
C
C ----- Inverse Jacobian ---
C
DTJ = YAC2( 1,1,M)*YAC2( 2,2,M) - YAC2( 1,2,M)*YAC2( 2,1,M)
C
YIN2( 1,1,M) = YAC2( 2,2,M)/ DTJ
YIN2( 1,2,M) = - YAC2( 1,2,M)/ DTJ
YIN2( 2,1,M) = - YAC2( 2,1,M)/ DTJ
YIN2( 2,2,M) = YAC2( 1,1,M)/ DTJ
C
C ----- Conversion from ξηζ- coordinates to ｘｙｚ coordinates ---
C
C ( Derivatives of a shape function )
C
C -----------------------------
DO 600 I = 1, 2
DO 600 J = 1, 4
C
EXY2( I,J,M) = 0.D0
600 CONTINUE
C -----------------------------
C
C -------------------------------------
DO 100 I = 1, 2
DO 100 J = 1, 4
DO 100 K = 1, 2
C --------------------------------------
C
EXY2( I,J,M) = EXY2( I,J,M) + YIN2( I,K,M)* EZG2( K,J,M)
C
DPX( IEL,M) = EXY2( 1,1,M)* CC( 1) + EXY2( 1,2,M)* CC( 2)
1 + EXY2( 1,3,M)* CC( 3) + EXY2( 1,4,M)* CC( 4)
C
DPY( IEL,M) = EXY2( 2,1,M)* CC( 1) + EXY2( 2,2,M)* CC( 2)
1 + EXY2( 2,3,M)* CC( 3) + EXY2( 2,4,M)* CC( 4)
C
100 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
FUNCTION DVLX ( X, C, NDE, NNP, NEL, I )
C
C ----Derivatives ( 2-Dimensional ) ----
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)))
C1 = C( NDE( I,1))
C2 = C( NDE( I,2))
C3 = C( NDE( I,3))
C4 = C( NDE( I,4))
C
DVLX = ( - C1 + C2 + C3 - C4 )/( 2.D0* HX)
C
RETURN
END
C **********************************************************************
C
FUNCTION DVLY ( Y, C, NDE, NNP, NEL, I )
C
C ----- Derivatives ( 2-Dimensional ) ---
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
DVLY = ( - C1 - C2 + C3 + C4)/( 2.D0* HY)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EPSL ( U, V, T, SU, SV, ST, NNP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- Attention --- < = 100000 ---
C
COMMON /CS3/ UVX, UMX, VMX, RA, EPS, EPM
COMMON /CST/ DT, VCF, PHI, TDF, FT, ALF, TB, LPE, LOP
COMMON /VRP/ VRV( 100000), VNA( 100000), VNB( 100000),
1 VDN( 100000)
C
DIMENSION U ( NNP), V( NNP), T( NNP), SU( NNP), SV( NNP),
1 ST( NNP)
C
EPS = 0.D0
UVX = - 10000.D0
C ------------------------
DO 100 I = 1, NNP
C
UVX = DMAX1( UVX, DABS( SU( I)))
UVX = DMAX1( UVX, DABS( SV( I)))
C
100 CONTINUE
C --------------------------------------
DO 200 I = 1, NNP
C
UN = SU( I)
VN = SV( I)
TN = ST( I)
IF ( DABS(UN).GT.UVX*0.01D0) EPS = DMAX1(EPS, DABS((U(I)-UN)/UN))
IF ( DABS(VN).GT.UVX*0.01D0) EPS = DMAX1(EPS, DABS((V(I)-VN)/VN))
U( I) = UN
V( I) = VN
T( I) = TN
C
200 CONTINUE
C --------------------------------------
EPS = EPS* 100.D0 ! (%) Expression
VRV( LOP) = EPS
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FGEM ( FF, GG, X, Y, ADF, U, V, NC1,
1 NC2, NC3, NDE, AA, HHX, HHY, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION FF ( NNP), GG ( NNP), X ( NNP), Y ( NNP),
1 U ( NNP), V ( NNP), AA ( 10,NNP), HHX( NNP),
2 HHY( NNP), NC1( NNP), NC2( NNP,10), NC3( NNP,10),
3 NDE( NEL,4)
C
DIMENSION ELX ( 4,4), ELY ( 4,4), EKX( 4,4), EKY( 4,4),
1 DPSX( 4), DPSY( 4), DXS( 2,2), DSX( 2,2)
C
COMMON /DEF / XI( 4,2), WS( 4), EPP( 4,4), EDD( 4,4)
C
DIMENSION PSI( 4), DPSI( 4,2), XX ( 2, 4), EC( 4,4),
1 PP ( 20,4), DX ( 20,4), DY ( 20,4), FX( 20,4),
2 FY ( 20,4), WW ( 20,4), WWW( 20,4), WX( 20,4),
3 WY ( 20,4)
C
C ----- Important --- Attention: ( NNP < = 30000 ) ---
C
DIMENSION NC11( 20000), X1 ( 100), Y1 ( 100), X2 ( 100),
1 Y2 ( 100), FF1( 30000), FF2( 30000), GG1( 30000)
C
COMMON /CST/ DT, VCF, PHI, TDF, FT, ALF, TB, LPE, LOP
C
C ----- Check ---
C
IF ( NNP.GT.30000) WRITE(6,*) ' *** Check in FGEM *** STOP ***'
IF ( NNP.GT.30000) STOP
C
PHI = 1.D0
C
C ----- Attention ---
C
NN = 2
C ---------------------------------------------------------------
DO 200 ID = 1, NNP
C
KEY = NC1( ID)
IF ( KEY.LT.3) GO TO 910
C
C ----------------------------
DO 210 I = 1, NNP
C
NC11( I) = 0
210 CONTINUE
C ----------------------------
C
IF ( LOP.NE.1) GO TO 920
C
C --------------------------------------
DO 220 NLE = 1, KEY
C
X1( NLE) = 0.D0
Y1( NLE) = 0.D0
X2( NLE) = 0.D0
Y2( NLE) = 0.D0
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 ----------------------------
220 CONTINUE
C --------------------------------------
CV
C --------------------------------------
DO 230 NLE = 1, KEY
DO 230 I = 1, 4
C
IF ( NC11( NDE( NC2( ID,NLE),I)).EQ.1 ) THEN
X1( NLE) = X( NDE( NC2( ID,NLE),I))
Y1( NLE) = Y( NDE( NC2( ID,NLE),I))
ENDIF
C
230 CONTINUE
C --------------------------------------
C
C ----- X2MIN = 10.---
C
X2MIN = 10.D0
C --------------------------------------
DO 240 NLE = 1, KEY
C
IF ( Y1( NLE).LT.Y( ID)) THEN
X2( NLE) = X1( NLE)
Y2( NLE) = Y1( NLE)
X2MIN = DMIN1( X2MIN, X2( NLE))
ENDIF
C
240 CONTINUE
C --------------------------------------
DO 250 NLE = 1, KEY
C
DIFF = DABS( X2MIN - X2( NLE))
IF ( DIFF.LE.0.0001D0) THEN
HHX( ID) = DABS ( X( ID)-X2( NLE))
HHY( ID) = DABS ( Y( ID)-Y2( NLE))
END IF
C
250 CONTINUE
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 --------------------------------------
NGL = NC3( ID,NLE)
C
C ----- To initialize element arrays ---
C
C --------------------------------------
DO 120 I = 1, 4
DO 120 J = 1, 4
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 --------------------------------------
C
NL = 4
C ------------------------------------------------
DO 130 L = 1, NL
C
CALL FGS4 ( XI( L,1), XI( L,2), PSI, DPSI)
C
C ----- To calculate DXS ---
C
C --------------------------------------
DO 20 I = 1, 2
DO 20 J = 1, 2
C
DXS( I,J) = 0.D0
C --------------------------------------
DO 20 K = 1, 4
C
DXS( I,J) = DXS( I,J) + DPSI( K,J)* XX( I,K)
20 CONTINUE
C --------------------------------------
C
C ------ To compute DSX ---
C
DTJ = DXS( 1,1)* DXS( 2,2) - DXS( 1,2)* DXS( 2,1)
DSX( 1,1) = DXS( 2,2)/ DTJ ! A11
DSX( 2,2) = DXS( 1,1)/ DTJ ! A22
DSX( 1,2) = - DXS( 1,2)/ DTJ ! A21
DSX( 2,1) = - DXS( 2,1)/ DTJ ! A12
C
C ----- To calculate D( PSI)/ DX
C
C --------------------------------------
DO 40 I = 1, 4
C
DPSX( I) = 0.D0
DPSY( I) = 0.D0
DPSX( I) = DPSI( I,1)* DSX( 1,1) + DPSI( I,2)* DSX( 2,1)
DPSY( I) = DPSI( I,1)* DSX( 1,2) + DPSI( I,2)* DSX( 2,2)
C
40 CONTINUE
C --------------------------------------
C
C ----- To accumlate integration point value of integrals ---
C
** FAC = DTJ * W( L, NN)
C
FAC = DTJ* WS( L)
C --------------------------------------
DO 60 I = 1, 4
DO 60 J = 1, 4
C
EC( I,J) = EC( I,J) + PSI( I)* PSI( J)* FAC
C
EKX( I,J) = EKX( I,J) + DPSX( I)* DPSX( J)* HHX( ID)** 2* FAC
EKY( I,J) = EKY( I,J) + DPSY( I)* DPSY( J)* HHY( ID)** 2* FAC
ELX( I,J) = ELX( I,J) + PSI( I)* DPSX( J)* HHX( ID)* FAC
ELY( I,J) = ELY( I,J) + PSI( I)* DPSY( J)* HHY( ID)* FAC
C
60 CONTINUE
C --------------------------------------
C
130 CONTINUE
C ------------------------------------------------
C
IF ( NLE.EQ.1) THEN
C
C --------------------------------------
DO 135 I = 1, 10
C
AA ( I,ID) = 0.D0
135 CONTINUE
C --------------------------------------
ENDIF
C
C --------------------------------------
DO 140 J = 1, 4
C
PP( NLE,J) = EC ( NGL,J)
DX( NLE,J) = EKX( NGL,J)
DY( NLE,J) = EKY( NGL,J)
FX( NLE,J) = ELX( NGL,J)
FY( NLE,J) = ELY( NGL,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 --------------------------------------
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)
C
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)
C
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 ----- In case of LOP GE.2 ---
C
BBX = DT* U( ID)/ HHX( ID)
BBY = DT* V( ID)/ HHY( ID)
RRX = DT* ADF/ HHX( ID)/ HHX( ID)
RRY = DT* ADF/ HHY( ID)/ HHY( ID)
C
IF ( DABS( BBX).LE.0.0001D0 .AND. DABS( BBY ).LE.0.0001D0)
1 GO TO 910
C
C ------------------------------------------------
C
FF1( ID) = AA( 1,ID)* ( BBX + BBY)
1 * ( AA( 5,ID)* BBX + AA( 6,ID)* BBY )
2 + ( AA( 3,ID)* BBX + AA( 4,ID)* BBY )
C
3 * ( AA( 1,ID)* ( 2.D0* ( RRX + RRY ) + ( BBX + BBY )** 2 )
4 - 2.D0* AA( 2,ID)* ( BBX + BBY ))
C
FF2( ID) = - ( AA( 3,ID)* BBX + AA( 4,ID)* BBY )
1 * ( AA( 9,ID)* RRX + AA(10,ID)* RRY )
2 + ( AA( 7,ID)* RRX + AA( 8,ID)* RRY )
3 * (( AA( 5,ID)* BBX + AA( 6,ID)* BBY )
4 + 2.D0* PHI* ( BBX + BBY )* ( 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 )
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
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 FGST ( XI, WS )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XI( 4,2), WS( 4)
C
C ----- Two - point Quadrature --- Gauss-Legendre formulae ---
C
XI( 1,1) = - 0.577350269189626D0
XI( 2,1) = 0.577350269189626D0
XI( 3,1) = - 0.577350269189626D0
XI( 4,1) = 0.577350269189626D0
C
XI( 1,2) = - 0.577350269189626D0
XI( 2,2) = - 0.577350269189626D0
XI( 3,2) = 0.577350269189626D0
XI( 4,2) = 0.577350269189626D0
C
WS( 1) = 1.D0
WS( 2) = 1.D0
WS( 3) = 1.D0
WS( 4) = 1.D0
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 calculate derivatives 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 FSBR ( X, Y, U, V, T, P, DLP, SU,
1 SV, ST, DI, NC1, NC2, NC3, NBU, NBV,
2 DPX, DPY, NDE, ISL, LYX, LYY, NNP, NEL,
3 IYF, * )
C
C --- UMX. on a vertical centered axis
C
C VMX. on a horizontal centered axis ---
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X ( NNP), Y ( NNP), U ( NNP), V ( NNP),
1 T ( NNP), P ( NNP), DLP( NNP), SU( NNP),
2 SV( NNP), ST( NNP), NDE( NEL,4)
C
DIMENSION DI ( NNP), NC1( NNP), NC2( NNP,10), NC3( NNP,10),
1 NBU( NNP), NBV( NNP), DPX( NEL,4), DPY( NEL,4)
C
DIMENSION KNUM( 100), KNVM( 100)
C ----- Attention : Max.Size = < 100 ---
C
COMMON /CS2/ XMN, NPR, KSR
COMMON /CS1/ OMG, ERR, NTR, ITT
COMMON /CS3/ UVX, UMX, VMX, RA, EPS, EPM
COMMON /ABC/ HGN, RLW, ANU, HNX, HNI, RNX, RNI
COMMON /CST/ DT, VCF, PHI, TDF, FT, ALF, TB, LPE, LOP
COMMON /VRP/ VRV( 100000), VNA( 100000), VNB( 100000),
1 VDN( 100000)
C
IYF = IYF + ITT
C
C -----------------------------------------------------
C
CALL EPSL ( U, V, T, SU, SV, ST, NNP )
C
C -----------------------------------------------------
C
C ----- Comp. for Nu number ---------------------------
C
C ***** Attention : Structured meshes ( KSR = 1 ) of 50 X 50
C
C or Unstructured meshes ( KSR = 2 ) of 30 X 30 ONLY ***
C
C ----- To be changed for Other Meshes ---
C
IF ( KSR.EQ.1.AND.NEL.EQ.2500)
C
1 CALL CMNU ( X, Y, T, LYX, LYY, NNP )
C
IF ( KSR.EQ.2) CALL CMNU ( X, Y, T, LYX, LYY, NNP )
C
*** IF ( KSR.EQ.2) CALL CPNU ( X, Y, T, LYX, LYY, NNP )
C
C -----------------------------------------------------------------------
UMX = - 1000.D0
VMX = - 1000.D0
K1 = 0
K2 = 0
C
DO 100 I = 1, NNP
C
IF ( DABS( X( I)-0.5D0).LE.0.001D0) THEN
K1 = K1 + 1
UMX = DMAX1( UMX, U( I))
KNUM( K1) = I
END IF
C ----------------------------
IF ( DABS( Y( I)-0.5D0).LE.0.001D0) THEN
K2 = K2 + 1
VMX = DMAX1( VMX, V( I))
KNVM( K2) = I
END IF
C
100 CONTINUE
C ------------------------------------------------
UMX = - 1000.D0
C
DO 150 I = 1, K1
C
NN = KNUM( I)
C
C ----- X( NN), Y( NN), U( NN) ---
C
IF ( U( NN).GE.UMX) THEN
UMX = U( NN)
K_YV = NN
END IF
C
150 CONTINUE
C ------------------------------------------------
VMX = - 1000.D0
C
DO 170 I = 1, K2
C
NN = KNVM( I)
C
C ----- X( NN), Y( NN), U( NN) ---
C
IF ( V( NN).GE.VMX) THEN
VMX = V( NN)
K_XV = NN
END IF
C
170 CONTINUE
C ------------------------------------------------
C
DNU = DABS( HGN) - DABS( RLW)
C ----------------------------------
C
IF ( MOD( LOP,NPR).EQ.0) THEN
WRITE(6,2000) LOP, EPS, ITT, DABS( HGN), DABS( RLW), DNU
2000 FORMAT('* L=',I6,' EPS=', F8.3,'(%) TIT=',I5,' Nu_0=',F7.3,
1 ' Nu_1=',F7.3,' D_Nu=',F7.3 )
WRITE(9,2000) LOP, EPS, ITT, DABS( HGN), DABS( RLW), DNU
END IF
C
KKC = 0
C
C ----- Attention : LOP < = 100,000 --- Abs. Value ---
C
VNA( LOP) = DABS( HGN)
VNB( LOP) = DABS( RLW)
VDN( LOP) = DNU
C
RLN = 0.001D0* ( VNA( LOP) + VNB( LOP))/ 2.D0
C
! 0.1(%) Criterion
C
C ----- Attention ----------------------------------------------------
C
IF ( EPS.LE.EPM ) KKC = 1
C
C Check for Var. of Max. Velocity
C
IF ( DNU.LE.RLN ) KKC = 2 ! Check for Difference of Nu_Values
C
C ======================================================================
C
IF ( KKC.EQ.1.OR.KKC.EQ.2) GO TO 200
IF ( UVX.GT.2000.D0) GO TO 300
IF ( LOP.LT.LPE) RETURN 1
GO TO 400
C
200 CONTINUE
C
IF ( KKC.EQ.1) THEN
WRITE(6,2100)
WRITE(9,2100)
2100 FORMAT(/,' *** Converged from Stationarity Check : End ***')
WRITE(6,2000) LOP, EPS, ITT, DABS( HGN), DABS( RLW), DNU
WRITE(9,2000) LOP, EPS, ITT, DABS( HGN), DABS( RLW), DNU
END IF
C
IF ( KKC.EQ.2) THEN
WRITE(6,2220) LOP
WRITE(9,2220) LOP
2220 FORMAT(/,' *** Converged from Nu-Value Check at LOP =',I4 )
WRITE(6,2000) LOP, EPS, ITT, DABS( HGN), DABS( RLW), DNU
WRITE(9,2000) LOP, EPS, ITT, DABS( HGN), DABS( RLW), DNU
END IF
GO TO 400
C
300 WRITE(6,2200) LOP
2200 FORMAT(' *** Error : Velocity is over 1000 * LOP =', I6 )
WRITE(6,*) '** STOP **'
STOP
C
400 CONTINUE
WRITE(6,2300) IYF
WRITE(9,2300) IYF
2300 FORMAT(/,' *** IYF =',I8)
C
C ---------------------------
C
CRNX = UMX* DT/ XMN
CRNY = VMX* DT/ XMN
C
C ---------------------------
C
WRITE(6,*) ' '
WRITE(9,*) ' '
WRITE(6,2400) UMX, Y( K_YV), VMX, X( K_XV), CRNX, CRNY
WRITE(9,2400) UMX, Y( K_YV), VMX, X( K_XV), CRNX, CRNY
2400 FORMAT(' * At X = 0.5 : Umax.=', F8.3,' at Y =',F8.4,/,
1 ' * At Y = 0.5 : Vmax.=', F8.3,' at X =',F8.4,/,
2 ' * Cour_X =',F6.3,' * Cour_Y =',F6.3)
C
C --------------------------------------
PMIN = 100000.D0
PMAX = - 100000.D0
C
DO 500 I = 1, NNP
C
PMIN = DMIN1( PMIN, P( I))
PMAX = DMAX1( PMAX, P( I))
C
500 CONTINUE
C --------------------------------------
WRITE(6,2500) PMIN, PMAX
WRITE(9,2500) PMIN, PMAX
2500 FORMAT(' * Pmin.=',F12.3,' Pmax.=',F12.3)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE HDSB ( X, Y, T, XX, YY, TT, NNP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X ( NNP), Y ( NNP), T ( NNP), XX( 6,100),
1 YY( 6,100), TT( 6,100)
C
DO 100 I = 1, 6
DO 100 J = 1, 100
C
XX( I,J) = 0.D0
YY( I,J) = 0.D0
TT( I,J) = 0.D0
C
100 CONTINUE
C
C ===== Attention : For 50 X 50 Structured meshes ===
C =============================
K1 = 0
C
C --------------------------------------
DO 200 I = 1, NNP
C
IF ( DABS( X( I)).LE.0.00001D0) THEN
K1 = K1 + 1
XX( 1,K1) = X( I)
YY( 1,K1) = Y( I)
TT( 1,K1) = T( I)
END IF
C
200 CONTINUE
C --------------------------------------
K2 = 0
DO 300 I = 1, NNP
C
IF ( DABS( X( I) - 0.0065D0 ).LE.0.0001D0 ) THEN
C ===== Attention ===============================
C
K2 = K2 + 1
C
XX( 2,K2) = X( I)
YY( 2,K2) = Y( I)
TT( 2,K2) = T( I)
C
END IF
C
300 CONTINUE
C --------------------------------------
K3 = 0
C
DO 400 I = 1, NNP
C
IF ( DABS( X( I) - 0.01358D0 ).LE.0.0001D0 ) THEN
C ====== Attention ===============================
C
K3 = K3 + 1
XX( 3,K3) = X( I)
YY( 3,K3) = Y( I)
TT( 3,K3) = T( I)
END IF
C
400 CONTINUE
C
K6 = 0
C
DO 500 I = 1, NNP
C
IF ( DABS( X( I) - 1.D0).LE. 0.00001D0) THEN
C ===========================================
C
K6 = K6 + 1
XX( 6,K6) = X( I)
YY( 6,K6) = Y( I)
TT( 6,K6) = T( I)
C
END IF
C
500 CONTINUE
C
K5 = 0
C
DO 600 I = 1, NNP
C
IF ( DABS( X( I) - 0.9935D0 ).LE.0.0001D0 ) THEN
C ====== Attention ==============================
K5 = K5 + 1
XX( 5,K5) = X( I)
YY( 5,K5) = Y( I)
TT( 5,K5) = T( I)
END IF
C
600 CONTINUE
C
K4 = 0
DO 700 I = 1, NNP
C
IF ( DABS( X( I) - 0.9864D0 ).LE.0.0001D0 ) THEN
C ====== Attention ==============================
K4 = K4 + 1
XX( 4,K4) = X( I)
YY( 4,K4) = Y( I)
TT( 4,K4) = T( I)
END IF
C
700 CONTINUE
C
C ----- Check ---
C
IF ( K1.NE.K2) THEN
WRITE(6,*) ' *** K1=',K1,' K2=', K2,' K3=', K3,' K4=',K4,' K5=',
1 K5,' K6=', K6
WRITE(6,*) ' '
WRITE(6,*) ' * STOP : Error: K1.NE.K2 from ''HDSB'' '
STOP
END IF
C
C ----- Change order ---
C
C ----- Left hand side ---
C
DO 800 I = 1, K1-1
DO 800 J = I+1, K1
C
IF ( YY( 1,I).GT.YY( 1,J)) THEN
C
A = YY( 1,I)
YY( 1,I) = YY( 1,J)
YY( 1,J) = A
C
C = TT( 1,I)
D = TT( 1,J)
TT( 1,I) = C
TT( 1,J) = D
C
E = XX( 1,I)
F = XX( 1,J)
XX( 1,I) = F
XX( 1,J) = E
C
END IF
C
800 CONTINUE
C
C ----- 1 from left ---
C
DO 900 I = 1, K1-1
DO 900 J = I+1, K1
C
IF ( YY( 2,I).GT.YY( 2,J)) THEN
C
A = YY( 2,I)
YY( 2,I) = YY( 2,J)
YY( 2,J) = A
C
C = TT( 2,I)
D = TT( 2,J)
TT( 2,I) = C
TT( 2,J) = D
C
E = XX( 2,I)
F = XX( 2,J)
XX( 2,I) = F
XX( 2,J) = E
C
END IF
C
900 CONTINUE
C
C ----- 2 from left ---
C
DO 920 I = 1, K1-1
DO 920 J = I+1, K1
C
IF ( YY( 3,I).GT.YY( 3,J)) THEN
C
A = YY( 3,I)
YY( 3,I) = YY( 3,J)
YY( 3,J) = A
C
C = TT( 3,I)
D = TT( 3,J)
TT( 3,I) = C
TT( 3,J) = D
C
E = XX( 3,I)
F = XX( 3,J)
XX( 3,I) = F
XX( 3,J) = E
C
END IF
C
920 CONTINUE
C ------------------------------------------------
C
C ----- Right edge ---
C
DO 940 I = 1, K1-1
C
DO 940 J = I+1, K1
C
IF ( YY( 6,I).GT.YY( 6,J)) THEN
C
A = YY( 6,I)
YY( 6,I) = YY( 6,J)
YY( 6,J) = A
C
C = TT( 6,I)
D = TT( 6,J)
TT( 6,I) = C
TT( 6,J) = D
C
E = XX( 6,I)
F = XX( 6,J)
XX( 6,I) = F
XX( 6,J) = E
C
END IF
C
940 CONTINUE
C
C ------ 1 from right ---
C
DO 950 I = 1, K1-1
DO 950 J = I+1, K1
C
IF ( YY( 5,I).GT.YY( 5,J)) THEN
C
A = YY( 5,I)
YY( 5,I) = YY( 5,J)
YY( 5,J) = A
C
C = TT( 5,I)
D = TT( 5,J)
TT( 5,I) = C
TT( 5,J) = D
C
E = XX( 5,I)
F = XX( 5,J)
XX( 5,I) = F
XX( 5,J) = E
C
END IF
C
950 CONTINUE
C
C ----- 2 from right ---
C
DO 960 I = 1, K1-1
C
DO 960 J = I+1, K1
C
IF ( YY( 4,I).GT.YY( 4,J)) THEN
C
A = YY( 4,I)
YY( 4,I) = YY( 4,J)
YY( 4,J) = A
C
C = TT( 4,I)
D = TT( 4,J)
TT( 4,I) = C
TT( 4,J) = D
C
E = XX( 4,I)
F = XX( 4,J)
XX( 4,I) = F
XX( 4,J) = E
C
END IF
C
960 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE HDSR ( X, Y, T, XX, YY, TT, NNP )

IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X ( NNP), Y( NNP), T( NNP), XX( 6,100),
1 YY( 6,100), TT( 6,100)
C
C ----- Attention : For Unstructured meshes ( 30X30) ---
C
DO 100 I = 1, 6
DO 100 J = 1, 100
C
XX( I,J) = 0.D0
YY( I,J) = 0.D0
TT( I,J) = 0.D0
C
100 CONTINUE
C
K = 0
C ------------------------------
DO 200 I = 1, NNP
C
IF ( DABS( X( I)).LE.0.0001D0) THEN
C
K = K + 1
XX( 1,K) = X( I)
YY( 1,K) = Y( I)
TT( 1,K) = T( I)
C
END IF
C
200 CONTINUE
C -----------------------------
C
K = 0
DO 300 I = 1, NNP
C
IF ( 0.0248D0.LT.X( I) .AND. X( I).LT.0.0316D0) THEN
C ----- Attention ------------------------------------
C
K = K + 1
XX( 2,K) = X( I)
YY( 2,K) = Y( I)
TT( 2,K) = T( I)
END IF
C
300 CONTINUE
C ----------------------------------------------------------
C
K = 0
C
DO 400 I = 1, NNP
C
IF ( 0.0498D0.LT.X( I) .AND. X( I).LT.0.06377D0) THEN
C ----- Attention ------------------------------------
K = K + 1
XX( 3,K) = X( I)
YY( 3,K) = Y( I)
TT( 3,K) = T( I)
END IF
C
400 CONTINUE
C
K = 0
DO 500 I = 1, NNP
C
IF ( DABS ( X( I) - 1.D0).LE. 0.0001D0) THEN
C
K = K + 1
XX( 6,K) = X( I)
YY( 6,K) = Y( I)
TT( 6,K) = T( I)
C
END IF
C
500 CONTINUE
C
K = 0
C
DO 600 I = 1, NNP
C
IF ( 0.9679D0.LT.X( I) .AND. X( I).LT.0.976D0) THEN
C ===== Attention ==================================
C
K = K + 1
XX( 5,K) = X( I)
YY( 5,K) = Y( I)
TT( 5,K) = T( I)
C
END IF
C
600 CONTINUE
C
K = 0
C
DO 700 I = 1, NNP
C
IF ( 0.9353D0.LT.X( I) .AND. X( I).LT.0.951D0) THEN
C ===== Attention ==================================
C
K = K + 1
XX( 4,K) = X( I)
YY( 4,K) = Y( I)
TT( 4,K) = T( I)
C
END IF
C
700 CONTINUE
C
C ----- Change order --- ************* ( ???? ) ****
C
C ----- Left hand side ---
C
C ------------------------------------------------
DO 800 I = 1, 41
C
C --------------------------------------
DO 850 J = I+1, 41
C
IF ( YY( 1,I).GT.YY( 1,J)) THEN
C
A = YY( 1,I)
YY( 1,I) = YY( 1,J)
YY( 1,J) = A
C
C = TT( 1,I)
D = TT( 1,J)
TT( 1,I) = C
TT( 1,J) = D
C
E = XX( 1,I)
F = XX( 1,J)
XX( 1,I) = F
XX( 1,J) = E
C
END IF
C
850 CONTINUE
C --------------------------------------
C
800 CONTINUE
C ------------------------------------------------
C
C ----- 1 from left ---
C
C ------------------------------------------------
DO 900 I = 1, 41
C
C ------------------------------------------------
DO 900 J = I+1, 41
C
IF ( YY( 2,I).GT.YY( 2,J)) THEN
C
A = YY( 2,I)
YY( 2,I) = YY( 2,J)
YY( 2,J) = A
C
C = TT( 2,I)
D = TT( 2,J)
TT( 2,I) = C
TT( 2,J) = D
C
E = XX( 2,I)
F = XX( 2,J)
XX( 2,I) = F
XX( 2,J) = E
C
END IF
C
900 CONTINUE
C ------------------------------------------------
C
C ----- 2 from left ---
C
C ------------------------------------------------
DO 950 I = 1, 41
C
DO 950 J = I+1, 41
C ------------------------------------------------
C
IF ( YY( 3,I).GT.YY( 3,J)) THEN
C
A = YY( 3,I)
YY( 3,I) = YY( 3,J)
YY( 3,J) = A
C
C = TT( 3,I)
D = TT( 3,J)
TT( 3,I) = C
TT( 3,J) = D
C
E = XX( 3,I)
F = XX( 3,J)
XX( 3,I) = F
XX( 3,J) = E
C
END IF
C
950 CONTINUE
C ------------------------------------------------
C
C ----- Right edge ---
C
DO 960 I = 1, 41
C
DO 960 J = I+1, 41
C
IF ( YY( 6,I).GT.YY( 6,J)) THEN
C
A = YY( 6,I)
YY( 6,I) = YY( 6,J)
YY( 6,J) = A
C
C = TT( 6,I)
D = TT( 6,J)
TT( 6,I) = C
TT( 6,J) = D
C
E = XX( 6,I)
F = XX( 6,J)
XX( 6,I) = F
XX( 6,J) = E
C
END IF
C
960 CONTINUE
C
C ------ 1 from right ---
C
C ------------------------------------------------
DO 970 I = 1, 41
DO 970 J = I+1, 41
C ------------------------------------------------
C
IF ( YY( 5,I).GT.YY( 5,J)) THEN
C
A = YY( 5,I)
YY( 5,I) = YY( 5,J)
YY( 5,J) = A
C
C = TT( 5,I)
D = TT( 5,J)
TT( 5,I) = C
TT( 5,J) = D
C
E = XX( 5,I)
F = XX( 5,J)
XX( 5,I) = F
XX( 5,J) = E
C
END IF
C
970 CONTINUE
C ------------------------------------------------
C
C ----- 2 from right ---
C
C ------------------------------------------------
DO 980 I = 1, 41
DO 980 J = I+1, 41
C
IF ( YY( 4,I).GT.YY( 4,J)) THEN
C
A = YY( 4,I)
YY( 4,I) = YY( 4,J)
YY( 4,J) = A
C
C = TT( 4,I)
D = TT( 4,J)
TT( 4,I) = C
TT( 4,J) = D
C
E = XX( 4,I)
F = XX( 4,J)
XX( 4,I) = F
XX( 4,J) = E
C
END IF
C
980 CONTINUE
C ------------------------------------------------
C
RETURN
END
C *********************************************************************
C
SUBROUTINE IDGN ( U, V, P, T, X, Y, NBU, NBV, NBP,
1 NBT, NBA, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION U ( NNP), V ( NNP), P ( NNP), T ( NNP),
1 X ( NNP), Y ( NNP), NBU( NNP), NBV( NNP),
2 NBP( NNP), NBT( NNP), NBA( NNP)
C
COMMON /CS1/ OMG, ERR, NTR, ITT
COMMON /CST/ DT, VCF, PHI, TDF, FT, ALF, TB, LPE, LOP
C
C ------ Attention ---------------
C
DO 100 I = 1, NNP
C
U( I) = 0.D0
V( I) = 0.D0
P( I) = 0.D0
T( I) = 0.D0
C
100 CONTINUE
C --------------------------------
DO 200 I = 1, NNP
C
NBU( I) = 0
NBV( I) = 0
NBP( I) = 0
NBT( I) = 0
NBA( I) = 0
C
200 CONTINUE
C --------------------------------
C
C ----- Fixed value at boundary ---
C
J = 1
C ------------------------------------------------
DO 300 I = 1, NNP
C
C ----- Velocity at each corner ---
C
IF ( DABS( Y( I)-1.D0).LE.0.0001D0 .AND.
1 ( DABS( X( I)) .LE. 0.0001D0).OR.
2 ( DABS( X( I)-1.D0).LE.0.0001D0 )) THEN
C
U( J) = 0.D0
NBU( J) = 9999
V( J) = 0.D0
NBV( J) = 9999
C
ELSE IF ( DABS( Y( I)-1).LE.0.0001D0 ) THEN
U( J) = 0.D0
NBU( J) = 9999
C
V( J) = 0.D0
NBV( J) = 9999
C
ELSE
1 IF( DABS( X( I)).LE.0.0001D0 .OR.
2 DABS( X( I) - 1.D0).LE.0.0001D0 .OR.
3 DABS( Y( I)).LE.0.0001D0) THEN
C
U( J) = 0.D0
NBU( J) = 9999
V( J) = 0.D0
NBV( J) = 9999
END IF
C
IF ( DABS(X( I)) .LE. 0.0001D0) THEN
T( J) = 1.D0
NBT( J) = 9999
C
ELSE IF
1 ( DABS( X( I)-1.D0).LE.0.0001D0) THEN
T ( J) = 0.D0
NBT( J) = 9999
ENDIF
C
C ----- Heat convection at boundary ---
C
C ----- Upper Wall ---
C
IF ( DABS( Y( I)-1.D0).LE.0.0001D0 .AND.
1 ( DABS( X( I)).LE.0.0001D0 .OR.
2 DABS( X( I)-1.D0) .LE.0.0001D0)) THEN
C
NBA( J) = 9999
C
ELSE
1IF( DABS( Y( I) - 1.D0).LE.0.0001D0) THEN
NBA( J) = 9999
ENDIF
C
J = J + 1
C
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C *********************************************************************
C
SUBROUTINE INIT ( SU, SV, ST, DPX, DPY, DLP, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION SU ( NNP), SV ( NNP), ST( NNP), DPX( NEL,4),
1 DPY( NEL,4), DLP( NNP)
C
C ----------------------------
DO 100 I = 1, NNP
C
SU( I) = 0.D0
SV( I) = 0.D0
ST( I) = 0.D0
C
DLP( I) = 0.D0
C
100 CONTINUE
C ----------------------------
DO 200 I = 1, NEL
DO 200 J = 1, 4
C
DPX( I,J) = 0.D0
DPY( I,J) = 0.D0
C
200 CONTINUE
C ----------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INPT ( X, Y, U, V, T, P, DD, NDE,
1 PP, NBU, NBV, NBT, NBP, NBA, NCC, NCW,
2 KCC, ISL, NNP, NEL, NNZ, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X ( NNP), Y ( NNP), U ( NNP), V ( NNP),
1 T ( NNP), P ( NNP), DD ( NNZ), NDE( NEL,4),
2 PP ( NNZ), NBU( NNP), NBV( NNP), NBT( NNP),
3 NBP( NNP), NBA( NNP), NCC( 2,NNZ), NCW( 2,NZ2),
4 KCC( NNP)
C
COMMON /DEF/ XI( 4,2), WS( 4), EPP( 4,4), EDD( 4,4)
C
DIMENSION XX( 4), YY( 4)
C
COMMON /CS2/ XMN, NPR, KSR
COMMON /CS1/ OMG, ERR, NTR, ITT
COMMON /CS3/ UVX, UMX, VMX, RA, EPS, EPM
COMMON /CST/ DT, VCF, PHI, TDF, FT, ALF, TB, LPE, LOP
C
UVX = 0.D0
C
WRITE(6,2000)
2000 FORMAT(' *** Ra = ?')
*** READ( 5,*) RA
C
RA = 710.D0
RA = 1.D+5
** RA = 7100
C
PR = 0.71D0
GR = RA/ PR
RE = 1.D0/ PR
** RE = 400.D0
C
TDF = 1.D0/( RE* PR)
FT = GR/( RE* RE)
VCF = PR
C
C ----- ALF : Coeff. of Heat Convection / TB : Ref. temperature ---
C
WRITE(6,2100)
2100 FORMAT(' ** ALF = ?')
**** READ( 5,*) ALF
ALF = 0.D0
C
** WRITE(6,2200)
** 2200 FORMAT(' ** Tb = ?')
**** READ( 5,*) TB
C
TB = 0.D0
C ===== Attention ===
C
ISL = 1
C
WRITE(6,2400)
2400 FORMAT(' * Time Step = ? & Loop = ? NPR = ? ')
*** READ( 5,*) DT, LPE, NPR
C
DT = 0.0001D0
LPE = 3000
NPR = 100
C
IF ( KSR.EQ.2) THEN ! KSR = 2 Unstructured Meshes
DT = 0.01D0
LPE = 20
C
DT = 0.0001D0
LPE = 3000
NPR = 100
C
* DT = 0.001D0
* LPE = 2000
* NPR = 50
C
END IF
C ===================
C
EPM = 0.01D0 ! EPM: Expression by (%)
C ===========================================
C
PHI = 1.D0
ERR = 0.0001D0
NTR = 200
C
C -----------------
OMG = 1.7D0
C -----------------
WRITE(6,2700) NNP, NEL, RA, PR, ALF, DT, LPE, NPR, OMG
2700 FORMAT(/,2X,' NNP =',I5,' NEL =',I5,' * Ra =',F8.1,' Pr =',
1 F5.2,1X,' ALF =',F5.2,/,' DT =',F9.6,' Max.Loop =',I6,
2 ' NPR=',I4,' ** OMG =',F5.2/ )
C
WRITE(7,2040) DT, RA, NNP, NEL, NPR
2040 FORMAT(F12.9, F12.1, 3I5)
WRITE(7,2050) ( X( I), Y( I), I = 1, NNP)
2050 FORMAT( 10F8.4)
WRITE(7,2060) (( NDE( I,J), J = 1, 4), I= 1, NEL)
2060 FORMAT( 20I4)
C
C ----------------------------------------------------------------------
C
CALL IDGN ( U, V, P, T, X, Y, NBU, NBV, NBP, NBT,
1 NBA, NNP, NEL )
C
C ----------------------------------------------------------------------
C
C ----- To set Gaussian "Weighting" factors ---
C
CALL STJC
C
C ---------------------------------------------
C
CALL STNC ( NCC, NNZ, NDE, NEL, 4, PP )
C
C ---------------------------------------------
DO 100 J = 1, NNZ
C
NCW( 1,J) = NCC( 1,J)
NCW( 2,J) = NCC( 2,J)
C
100 CONTINUE
C ------------------------------------------------
C
NZ = NNZ
C
CALL DCIJ ( KCC, NCW, NZ, NNP, NZ2 )
C
C ------------------------------------------------
DO 200 I = 1,NNZ
C
DD( I) = 0.D0
PP( I) = 0.D0
C
200 CONTINUE
C ------------------------------------------------
C
K = 0
C ------------------------------------------------
DO 300 IEL = 1, NEL
C
C ------------------------------------
DO 400 I = 1, 4
C
N = NDE( IEL,I)
XX( I) = X( N)
YY( I) = Y( N)
C
400 CONTINUE
C ------------------------------------
C
CALL JACP ( XX, YY, EPP, EDD)
C
C ------------------------------------
DO 500 I = 1, 4
DO 500 J = 1, 4
C
II = NDE( IEL,I)
JJ = NDE( IEL,J)
C
CALL DCON ( EDD( I,J), II, JJ, DD, NCC, NNZ, L )
C
IF ( L.NE.0) PP( L) = PP( L) + EPP( I,J)
C
500 CONTINUE
C --------------------------------------
C
300 CONTINUE
C ------------------------------------------------
C
CALL FGST ( XI, WS )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE JACP ( XX, YY, EPP, EDD)
C
C ----- PD - Matrixes by Jacobian ---
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XX ( 4), YY ( 4), EPP( 4,4), EDD( 4,4), ENN( 4),
1 ENG( 2,4), EXY( 2,4), YAC( 2,2), YIN( 2,2)
C
COMMON /XE/ XIR( 4), ETA( 4)
COMMON /PW/ PG ( 8), WG ( 8), NPW
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
X0 = PG( II)* XIR( I)
Y0 = PG( JI)* ETA( I)
ENN( I) = ( 1.D0 + X0)* ( 1.D0 + Y0)/ 4.D0
ENG( 1,I) = XIR( I)* ( 1.D0 + Y0)/ 4.D0
ENG( 2,I) = ETA( I)* ( 1.D0 + X0)/ 4.D0
C
300 CONTINUE
C --------------------------------------
C
C ----- Jacobian ---
C
C --------------------------------------
DO 400 I = 1, 2
DO 400 J = 1, 2
C
YAC( I,J) = 0.D0
400 CONTINUE
C --------------------------------------
DO 500 I = 1, 4
DO 500 K = 1, 2
C
YAC( K,1) = YAC( K,1) + ENG( K,I)* XX( I)
YAC( K,2) = YAC( K,2) + ENG( 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)
C
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 from ξη- coordinate
C
C ( Derivatives of shape functions ) into ｘｙ- coordinate ---
C
C --------------------------------------
DO 600 I = 1, 2
DO 600 J = 1, 4
C
EXY( I,J) = 0.D0
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)* ENG( K,J)
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
EPP( I,J) = EPP( I,J) + ENN( I)* ENN( J)* WGT
C
EDD( I,J) = EDD( I,J) + ( EXY( 1,I)* EXY( 1,J)
1 + EXY( 2,I)* EXY( 2,J))* WGT
C
200 CONTINUE
C -----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE JPDF ( XX, YY, EDX, EDY, EFX, EFY, EPP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ------DF-Matrixes by Jacobian ---
C
DIMENSION XX ( 4), YY ( 4), EDX( 4,4), EDY( 4,4),
1 EFX( 4,4), EFY( 4,4), EPP( 4,4)
C
DIMENSION ENN( 4), ENG( 2,4), EXY( 2,4), YAC( 2,2), YIN( 2,2)
C
COMMON /XE/ XIR( 4), ETA( 4)
COMMON /PW/ PG ( 8), WG ( 8), NPW
C
C ----------------------------
DO 100 I = 1, 4
DO 100 J = 1, 4
C
EPP( I,J) = 0.D0
EDX( I,J) = 0.D0
EDY( I,J) = 0.D0
EFX( I,J) = 0.D0
EFY( I,J) = 0.D0
C
100 CONTINUE
C -------------------------------------------------------
DO 200 II = 1, NPW
DO 200 JI = 1, NPW
C -------------------------------------------------------
C
C ----- Shape functions & its Derivatives ---
C
C ------------------------------------------------
DO 300 I = 1, 4
C
X0 = PG( II)* XIR( I)
Y0 = PG( JI)* ETA( I)
ENN( I) = ( 1.D0 + X0)* ( 1.D0 + Y0)/ 4.D0
ENG( 1,I) = XIR( I)* ( 1.D0 + Y0)/ 4.D0
ENG( 2,I) = ETA( I)* ( 1.D0 + X0)/ 4.D0
C
300 CONTINUE
C ------------------------------------------------
C
C ----- Jacobian matrix ---
C
C ----------------------------
DO 400 I = 1, 2
DO 400 J = 1, 2
C
YAC( I,J) = 0.D0
400 CONTINUE
C ----------------------------
DO 500 I = 1, 4
DO 500 K = 1, 2
C
YAC( K,1) = YAC( K,1) + ENG( K,I)* XX( I)
YAC( K,2) = YAC( K,2) + ENG( K,I)* YY( I)
C
500 CONTINUE
C ----------------------------
C
C ----- Inversed Jacobian ---
C
DTJ = YAC( 1,1)* YAC( 2,2) - YAC( 1,2)* YAC( 2,1)
C
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 from ξ-η to x-y coordinates
C
C & Derivatives of shape functions ---
C
C --------------------------------------
DO 600 I = 1, 2
DO 600 J = 1, 4
C
EXY( I,J) = 0.D0
600 CONTINUE
C --------------------------------------
DO 700 I = 1, 2
DO 700 J = 1, 4
DO 700 K = 1, 2
C --------------------------------------
C
EXY( I,J) = EXY( I,J) + YIN( I,K)* ENG( K,J)
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
EPP( I,J) = EPP( I,J) + ENN( I) * ENN( J) * WGT
EDX( I,J) = EDX( I,J) + EXY( 1,I)* EXY( 1,J)* WGT
EDY( I,J) = EDY( I,J) + EXY( 2,I)* EXY( 2,J)* WGT
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 MSHG ( X, Y, NDE, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CS2/ XMN, NPR, KSR
C
DIMENSION X( NNP), Y( NNP), NDE( NEL,4)
C
C ----- Key for Un-structured meshes ---
C
KSR = 2
C =============
C
READ(10,1000) NNP, NEL
1000 FORMAT(2I5)
READ(10,1100) (( NDE( I,J),J = 1, 4), I = 1, NEL)
1100 FORMAT( 4I5)
READ(10,1200) ( X( I), Y( I), NN, I = 1, NNP)
1200 FORMAT( 2E15.7,I5)
C
XMN = 1.D0
DO 100 I = 2, NNP
C
XX = X( I)
IF ( DABS( XX).GT.0.0001D0) XMN = DMIN1( XMN, XX)
C
100 CONTINUE
C
WRITE(6,2000) XMN
2000 FORMAT(' * Xmin.=', F7.4)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLTD ( U, V, T, FF, NNP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION U( NNP), V( NNP), T( NNP), FF( NNP)
C
COMMON /CST/ DT, VCF, PHI, TDF, FT, ALF, TB, LPE, LOP
COMMON /VRP/ VRV( 100000), VNA( 100000), VNB( 100000),
1 VDN( 100000)
C
C ------------------------------------------------
IF ( LOP.GE.30000) THEN
C
C ----- For Plot. Routines N_Plot = < 30000 ---
C
NL = 0
DO 100 I = 1, LOP, 3
C
NL = NL + 1
VRV( NL) = VRV( I)
VDN( NL) = VDN( I)
VNA( NL) = VNA( I)
VNB( NL) = VNB( I)
C
100 CONTINUE
C ------------------------------------------------
WRITE(6,*) ' *** Attention: No_Plot_LOP =', NL
WRITE(9,*) ' *** Attention: No_Plot_LOP =', NL
LOP = NL
C
END IF
C ------------------------------------------------
WRITE(7,2000) LOP
2000 FORMAT(I10)
C
WRITE(8,2100) ( VNA( I), VNB( I), I=1, LOP )
2100 FORMAT(10F8.4)
C
WRITE(9,*) ' '
WRITE(7,2200) ( VRV( I), I= 1, LOP)
** WRITE(6,*) ( VRV( I), I= 1, LOP)
2200 FORMAT(10F11.4)
C
WRITE(7,2200) ( VDN( I), I= 1, LOP)
WRITE(7,2200) ( U( I), V( I), I= 1, NNP)
WRITE(7,2200) ( T( I), I= 1, NNP)
C
C ------------------------------------------------
RMIN = 100.D0
RMAX = 0.D0
DO 200 I = 1, NNP
C
IF ( FF( I).GE.RMAX) RMAX = FF( I)
IF ( FF( I).LE.RMIN) RMIN = FF( I)
C
200 CONTINUE
C ------------------------------------------------
C
WRITE(6,2300) RMAX, RMIN
WRITE(9,2300) RMAX, RMIN
2300 FORMAT(' *** f-Max.=',F7.3,' f-Min.=',F7.3)
WRITE(6,2400)
WRITE(9,2400)
2400 FORMAT(/,' *** End ***')
C
RETURN
END
C **********************************************************************
C
SUBROUTINE POIS ( X, Y, U, V, P, NDE, A, B,
1 D, F, DD, NCC, NCW, KCC, NBP, NNP,
2 NEL, NNZ, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X ( NNP), Y ( NNP), U ( NNP), V ( NNP),
1 P ( NNP), NDE( NEL,4), A ( NNZ), B ( NNP),
2 D ( NNZ), F ( NNP), DD ( NNZ),
3 NCC( 2,NNZ), NCW( 2,NZ2), KCC( NNP), NBP( NNP)
C
DIMENSION XX ( 4), YY ( 4), EDX( 4,4), EDY( 4,4), PE( 4),
1 EPP( 4,4), EFX( 4,4), EFY( 4,4)
C
COMMON /CS1/ OMG, ERR, NTR, ITT
COMMON /CST/ DT, VCF, PHI, TDF, FT, ALF, TB, LPE, LOP
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
200 CONTINUE
C --------------------------------------
DO 300 IEL = 1, NEL
C
UE = 0.D0
VE = 0.D0
C ----------------------------
DO 350 I = 1, 4
C
N = NDE( IEL,I)
XX( I) = X( N)
YY( I) = Y( N)
PE( I) = P( N)
C
350 CONTINUE
C ----------------------------
C
C ----- D, F ------------------------------------------
C
CALL JPDF ( XX, YY, EDX, EDY, EFX, EFY, EPP)
C
C -----------------------------------------------------
DO 300 I = 1, 4
DO 300 J = 1, 4
C --------------------------------------
C
II = NDE( IEL,I)
JJ = NDE( IEL,J)
DS = EDX( I,J) + EDY( I,J)
C
C ----------------------------
C
CALL DCON ( DS, II,JJ, D, NCC, NNZ, L )
C
C ----------------------------
C
F( II) = F( II)
1 + ( U( JJ)* EFX( I,J) + V( JJ)* EFY( I,J))/ DT
C
C ----- Simultaneous Equations by MGM ---
C
300 CONTINUE
C --------------------------------------
C
JS = 1
C ---------------------------------
DO 400 I = 1, NNP
C
JL = KCC( I)
C ----------------------------
DO 450 J = JS, JL
C
KPR = NCW( 1,J)
KCL = NCW( 2,J)
A( KPR) = D( KPR)
C
450 CONTINUE
C ----------------------------
B( I) = B( I) - F( I)
JS = JL + 1
C
400 CONTINUE
C -----------------------------------------------------------
C
CALL SLVE ( B, P, DD, NCW, KCC, NBP, NNP, NNZ, NZ2 )
C
C -----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRM4F ( XLX, YLY, LYX, LYY, N, NNZ, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XLX( LYX), YLY( LYY)
C
NEX = LYX - 1
C
C ------------------------------------------------
XLX( 1) = 0.D0
DO 100 I = 1, 2
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( 3) = 0.5D0
DO 200 I = 4, 5
C
XLX( I) = 1.D0 - XLX( 6-I)
200 CONTINUE
C ------------------------------------------------
XLX( 5) = 1.D0
WRITE(6,*) ' * XLX =', ( XLX( I), I= 1, 5 )
C ------------------------------------------------
DO 300 I = 1, 5
C
YLY( I) = XLX( I)
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRM10F ( XLX, YLY, LYX, LYY, N, NNZ, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XLX( 11), YLY( 11)
C
NEX = 10
NEY = 10
LYX = NEX + 1
LYY = NEY + 1
C
NNP = LYX* LYY
NEL = NEX* NEY
N = 4* NEX* NEY + NEX + NEY
C
NNZ = NNP + N
NZ2 = NNP + 2* N
C
XLX( 1) = 0.D0
C ------------------------------------------------
DO 100 I = 2, 5
C
IM = I - 1
XLX( I) = ( DEXP( 2.D0* DFLOAT( IM)* 2.D0/ DFLOAT( NEX)) - 1.D0)
1 /( 2.D0* ( DEXP( 2.D0) - 1.D0 ))
C
100 CONTINUE
C ------------------------------------------------
XLX( 6) = 0.5D0
C
DO 200 I = 7, 10
C
XLX( I) = 1.D0 - XLX( 12-I)
200 CONTINUE
C ------------------------------------------------
C
XLX( 11) = 1.D0
DO 300 I = 1, 11
C
YLY( I) = XLX( I)
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRM10H ( XLX, YLY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XLX( 11), YLY( 11)
C
NEX = 10
NEY = 10
C
XLX( 1) = 0.D0
XLX( 2) = 0.1D0
XLX( 3) = 0.2D0
XLX( 4) = 0.3D0
XLX( 5) = 0.4D0
C
XLX( 6) = 0.5D0
XLX( 7) = 0.6D0
XLX( 8) = 0.7D0
XLX( 9) = 0.8D0
XLX(10) = 0.9D0
XLX(11) = 1.D0
C
YLY( 1) = 0.D0
YLY( 2) = 0.1D0
YLY( 3) = 0.2D0
YLY( 4) = 0.3D0
YLY( 5) = 0.4D0
C
YLY( 6) = 0.5D0
YLY( 7) = 0.6D0
YLY( 8) = 0.7D0
YLY( 9) = 0.8D0
YLY(10) = 0.9D0
YLY(11) = 1.D0
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRM20F ( XLX, YLY, LYX, LYY, N, NNZ, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XLX( 21), YLY( 21)
C
NEX = 20
NEY = 20
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
C ------------------------------------------------
XLX( 1) = 0.D0
DO 100 I = 1, 10
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( 11) = 0.5D0
C
DO 200 I = 12, 20
C
XLX( I) = 1.D0 - XLX( 22-I)
200 CONTINUE
C ------------------------------------------------
XLX( 21) = 1.D0
C
DO 300 I = 1, 21
C
YLY( I) = XLX( I)
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRM20H ( XLX, YLY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XLX( 21), YLY( 21)
C
XLX( 1) = 0.D0
XLX( 2) = 0.05D0
XLX( 3) = 0.1D0
XLX( 4) = 0.15D0
C
XLX( 5) = 0.2D0
XLX( 6) = 0.25D0
XLX( 7) = 0.3D0
XLX( 8) = 0.35D0
XLX( 9) = 0.4D0
XLX( 10) = 0.45D0
XLX( 11) = 0.5D0
C
XLX( 12) = 1.D0 - XLX(10)
XLX( 13) = 1.D0 - XLX( 9)
XLX( 14) = 1.D0 - XLX( 8)
XLX( 15) = 1.D0 - XLX( 7)
XLX( 16) = 1.D0 - XLX( 6)
XLX( 17) = 1.D0 - XLX( 5)
XLX( 18) = 1.D0 - XLX( 4)
XLX( 19) = 1.D0 - XLX( 3)
XLX( 20) = 1.D0 - XLX( 2)
XLX( 21) = 1.D0 - XLX( 1)
C
C ----------------------------
DO 100 I = 1, 21
C
YLY( I) = XLX( I)
100 CONTINUE
C ----------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRM30F ( XLX, YLY, LYX, LYY, N, NNZ, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XLX( 31), YLY( 31)
C
NEX = 30
NEY = 30
LYX = NEX + 1
LYY = NEY + 1
NNP = LYX * LYY
C
C ------------------------------------------------
XLX( 1) = 0.D0
DO 100 I = 1, 15
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( 16) = 0.5D0
C
DO 200 I = 17, 30
C
XLX( I) = 1.D0 - XLX( 32-I)
200 CONTINUE
C ------------------------------------------------
C
XLX( 31) = 1.D0
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, NNZ, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XLX( 41), YLY( 41)
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
DO 100 I = 1, 20
C
IM = I-1
XLX( I)=( DEXP( 2.D0* DFLOAT( IM)* 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
C
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, NNZ, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XLX( LYX), YLY( LYY)
C
NEX = LYX - 1
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
XLX( 51) = 1.D0
C ------------------------------------------------
DO 300 I = 1, 51
C
YLY( I) = XLX( I)
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRM100F ( XLX, YLY, LYX, LYY, N, NNZ, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XLX( 101), YLY( 101)
C
NEX = 100
NEY = 100
C
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, 50
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( 51) = 0.5D0
C
DO 200 I = 52, 100
C
XLX( I) = 1.D0 - XLX ( 102-I)
200 CONTINUE
C
XLX( 101) = 1.D0
C ------------------------------------------------
DO 300 I = 1, 101
C
YLY( I) = XLX ( I)
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRM200F ( XLX, YLY, LYX, LYY, N, NNZ, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XLX( 201), YLY( 201)
C
NEX = 200
NEY = 200
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, 100
C
XLX( I) = ( DEXP ( 2.D0 * DFLOAT( I-1)* 2.D0
1 / DFLOAT( NEX)) - 1.D0)
2 /( 2.D0* ( DEXP ( 2.D0 ) - 1.D0 ))
C
100 CONTINUE
C --------------------------------------
XLX( 1) = 0.D0
XLX( 101) = 0.5D0
C
DO 200 I = 102, 200
C
XLX( I) = 1.D0 - XLX ( 202-I)
200 CONTINUE
C
XLX ( 201) = 1.D0
C --------------------------------------
DO 300 I = 1, 201
C
YLY( I) = XLX( I)
300 CONTINUE
C --------------------------------------
C
RETURN
END
C *********************************************************************
C
SUBROUTINE PRSS ( P, DLP, NNP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION P( NNP), DLP( NNP)
C
DO 100 I = 1, NNP
C
P( I) = P( I) + DLP( I)
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SLVE ( FVEC, C, TSM, NCC, KCC, NBC, NNP, NNZ, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CS1/ OMG, ERR, NTR, ITT
COMMON /CST/ DT, VCF, PHI, TDF, FT, ALF, TB, LPE, LOP
C
DIMENSION FVEC( NNP), C ( NNP), KCC( NNP), NBC( NNP),
1 TSM ( NNZ), NCC( 2,NZ2)
C
C ----- S.O.R. Method ---
C
CMAX = 0.D0
C
DO 100 I = 1, NNP
C
AC = DABS( C( I))
IF ( AC.GE.CMAX) CMAX = AC
C
100 CONTINUE
C -----------------------------------
EPS = CMAX* 0.01D0
LMT = 3* NNP
KIT = 2* NNP
C
NIR = 0
900 NIR = NIR + 1
EMX = 0.D0
C ----------------------------------------------------
DO 200 I = 1, NNP
C
IF ( NBC( I).EQ.9999) GO TO 200
IF ( I.EQ.1) GO TO 604
MS = KCC( I-1) + 1
GO TO 605
C
604 MS = 1
605 ML = KCC( I)
IF ( MS.GT.ML ) GO TO 1000
DK = FVEC( I)
C
C --------------------------------------
DO 220 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 220
DK = DK - TSM( M1)* C( M2)
C
220 CONTINUE
C --------------------------------------
C
DC = - C( I) + DSM* DK
C( I)= C( I) + OMG* DC
TMP = DABS( C( I))
C
IF ( TMP.LE.EPS) GO TO 200
TMP = DABS( DC/ C( I))
C
IF ( TMP.LE.EMX) GO TO 200
EMX = TMP
NRX = I
C
200 CONTINUE
C ---------------------------------------------------
IF ( EMX.LE.ERR) GO TO 710
IF ( NIR.GE.KIT) GO TO 701
IF ( NIR.GE.LMT) GO TO 800
GO TO 900
C
800 CONTINUE
WRITE(6,2000) EPS, NIR, NRX, EMX
2000 FORMAT(' * EPS=',D11.4,' Iter =',I6,2X,' NDE=',I5,2X,' Error=',
1 E11.4)
WRITE(6,2100) ( C( I),I = 1,NNP)
2100 FORMAT(' ',6E18.6)
GO TO 900
C
701 CONTINUE
WRITE(6,2200) NIR
2200 FORMAT('0',' **** No. of Iterations exceeded, NIR=',I5 )
GO TO 710
C
1000 CONTINUE
WRITE(6,2300) I, MS, ML
2300 FORMAT('1','* Error ',3I6)
710 CONTINUE
C
C ----- Total Iteration Numbers ---
C
ITT = ITT + NIR
C
C ---------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SMSH ( X, Y, NDE, XLX, YLY, LX, LY, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X( NNP), Y( NNP), NDE( NEL,4), XLX( LX), YLY( LY )
C
COMMON /CS2/ XMN, NPR, KSR
COMMON /CS1/ OMG, ERR, NTR, ITT
COMMON /CST/ DT, VCF, PHI, TDF, FT, ALF, TB, LPE, LOP
C
C ----- For structured meshes ----------
C
XMN = XLX( 2) - XLX( 1)
C --------------------------------------
J = 1
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 --------------------------------------
J = 1
DO 200 IY = 1, LY-1
DO 200 IX = 1, LX-1
C
IST = IX + ( IY-1)* LX
NDE( J,1) = IST
NDE( J,2) = IST + 1
NDE( J,3) = IST + LX + 1
NDE( J,4) = IST + LX
J = J + 1
C
200 CONTINUE
C -------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SNCN ( NDE, NC1, NC2, NC3, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION NDE( NEL,4), NC1( NNP), NC2( NNP,10), NC3( NNP,10)
C
C --------------------------------------
DO 100 I = 1, NNP
C
NC1( I) = 0
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
RETURN
END
C **********************************************************************
C
SUBROUTINE STJC
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- To set Gaussian weighting parameters ---
C
COMMON /XE/ XIR( 4), ETA( 4)
COMMON /PW/ PG ( 8), WG ( 8), NPW
C
C ----- Quadrilateral element for ξ−η coordinates ---
C
XIR( 1) = -1.D0
ETA( 1) = -1.D0
C
XIR( 2) = 1.D0
ETA( 2) = -1.D0
C
XIR( 3) = 1.D0
ETA( 3) = 1.D0
C
XIR( 4) = -1.D0
ETA( 4) = 1.D0
C
C ---- Positions & Weighting parameters
C for Gaussian integration formulae ---
C
NPW = 2
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
PG( 1) = - 0.577350269189626 D0
PG( 2) = - PG( 1)
WG( 1) = 1.D0
WG( 2) = 1.D0
RETURN
C
20 CONTINUE
PG( 1) = - 0.774596669241483 D0
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)
RETURN
C
30 CONTINUE
PG( 1) = - 0.8611 3631 1594 053D0
PG( 2) = - 0.3399 8104 3584 856D0
PG( 3) = - PG( 2)
PG( 4) = - PG( 1)
C
WG( 1) = 0.3478 5484 5137 454D0
WG( 2) = 0.6521 4515 4862 546D0
WG( 3) = WG( 2)
WG( 4) = WG( 1)
RETURN
C
40 CONTINUE
PG( 1) = - 0.9602 8985 64D0
PG( 2) = - 0.7966 6647 74D0
PG( 3) = - 0.5255 3240 99D0
PG( 4) = - 0.1834 3464 24D0
C
PG( 5) = - PG( 4)
PG( 6) = - PG( 3)
PG( 7) = - PG( 2)
PG( 8) = - PG( 1)
C
WG( 1) = 0.1012 2853 62D0
WG( 2) = 0.2223 8103 44D0
WG( 3) = 0.3137 0664 58D0
WG( 4) = 0.3626 8378 33D0
C
WG( 5) = WG( 4)
WG( 6) = WG( 3)
WG( 7) = WG( 2)
WG( 8) = WG( 1)
RETURN
END
C **********************************************************************
C
SUBROUTINE STNC ( NCC, NNZ, NDE, NEL, NP, PP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION NDE( NEL,NP), PP( NNZ), NCC( 2,NNZ)
C
K = 0
C -----------------------------------------------------
DO 100 INE = 1, NEL
DO 100 I = 1, NP
C -----------------------------------------------------
II = NDE( INE,I)
C -----------------------------------------------------
DO 100 J = 1, NP
C
JJ = NDE( INE,J)
C
CALL DCNL ( 0.D0, II, JJ, PP, NCC, K, NNZ )
C
C ------------------------------------------------
C
100 CONTINUE
C -----------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE STUV ( U, V, GU, GV, NNP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION U( NNP), V( NNP), GU( NNP), GV( NNP)
C
DO 100 I = 1, NNP
C
GU( I) = U( I)
GV( I) = V( I)
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SUBS ( PDV, I, J, N, PP, NCC, K, NNZ )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION PP( NNZ), NCC( 2,NNZ)
C
K = K + 1
IF ( K.GT.NNZ) WRITE(6,2000)
2000 FORMAT(' * Error K.GT.NNZ * ')
JE = K - N
C
C --------------------------------------
DO 100 JJ = 1,JE
C
J1 = K - JJ + 1
PP ( J1) = PP ( J1-1)
NCC( 1,J1) = NCC( 1,J1-1)
NCC( 2,J1) = NCC( 2,J1-1)
C
100 CONTINUE
C --------------------------------------
PP ( N) = PDV
NCC( 1,N) = I
NCC( 2,N) = J
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
IMPLICIT REAL*8 (A-H,O-Z)
C
WRITE(6,*) ' '
WRITE(6,*) 'C ************************************************* C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C H F 2 D - F E M C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( Version 5.0 ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C 2D Thermal Fluid Analysis by FEM C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( Quardrilateral Element ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C Copyright 2011 : Yasuhiro MATSUDA C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ************************************************* C'
WRITE(6,*) ' '
C
RETURN
END
C **********************************************************************
C
SUBROUTINE UVTC ( X, Y, U, V, T, C, NDE, P,
1 PP, A, B, D, F, ID, CNW, NBC,
2 NCC, NCW, KCC, ISL, FF, GU, GV, NBA,
3 NNP, NEL, NNZ, NZ2 )
C
C ----- Modified Galerkin Formulation ( Matrix A & Vector B ) ---
C
IMPLICIT REAL*8 (A-H,O-Z)
c
DIMENSION X ( NNP), Y ( NNP), U ( NNP), V ( NNP),
1 T ( NNP), C ( NNP), NDE( NEL,4), P ( NNP),
2 PP( NNZ), A ( NNZ), B ( NNP), D ( NNZ),
3 F ( NNP), CNW( NNP), NBC( NNP), NCC( 2,NNZ),
4 NCW( 2,NZ2), KCC( NNP), FF ( NNP), GU ( NNP),
5 GV( NNP), NBA( NNP)
C
DIMENSION XX ( 4), YY ( 4), EDX( 4,4), EDY( 4,4),
1 TT ( 4), EFX( 4,4), EFY( 4,4), PE ( 4),
2 EPP( 4,4), NA ( 3)
C
COMMON/CST/ DT, VCF, PHI, TDF, FT, ALF, TB, LPE, LOP
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
FH = 0.D0
CVF = VCF
C
IF ( ID.EQ.0) CVF = TDF
IF ( ID.EQ.1) FSX = 1.D0
IF ( ID.EQ.2) FSY = 1.D0
IF ( ID.EQ.2) FH = 1.D0
C
C ------------------------------------------------
DO 300 IEL = 1, NEL
C
UE = 0.D0
VE = 0.D0
GUE = 0.D0
GVE = 0.D0
FF2 = 0.D0
C
NA( 1) = 0
NA( 2) = 0
NA( 3) = 0
JN = 1
C
C ----------------------------------
DO 400 I = 1, 4
C
N = NDE( IEL,I)
XX( I) = X( N)
YY( I) = Y( N)
TT( I) = T( N)
PE( I) = P( N)
C
UE = UE + U ( N)/ 4.D0
VE = VE + V ( N)/ 4.D0
GUE = GUE + GU( N)/ 4.D0
GVE = GVE + GV( N)/ 4.D0
FF2 = FF2 + FF( N)/ 4.D0
C
IF ( NBA( N).EQ.9999) THEN
NA( JN) = I
JN = JN + 1
END IF
C
400 CONTINUE
C ------------------------------------
C
C ----- DF Matrix -------------------------------------
C
CALL JPDF ( XX, YY, EDX, EDY, EFX, EFY, EPP )
C
C -----------------------------------------------------
DO 360 I = 1, 4
C
C --------------------------------------
DO 350 J = 1, 4
C
II = NDE( IEL,I)
JJ = NDE( IEL,J)
DS = CVF* FF2* ( EDX( I,J) + EDY( I,J))
C
C ---------------------------------
C
CALL DCON ( DS, II, JJ, D, NCC, NNZ, L )
C
C ---------------------------------
F( II) = F( II)
1 + ( GUE* EFX( I,J) + GVE* EFY( I,J))* C( JJ)
2 + FSX* EFX( I,J)* PE( J) + FSY* EFY( I,J)* PE( J)
3 - FH* EPP( I,J)* TT( J)* FT
C
350 CONTINUE
C --------------------------------------
C
360 CONTINUE
C ------------------------------------------------
C
C ----- Heat Convection ---
C
IF ( ID.NE.0) GO TO 300
IF ( DABS( ALF).LE.0.0001D0) GO TO 300
C
C ------------------------------------------------
IF ( NA( 1).EQ.1.AND.NA( 2).EQ.2) THEN
C
X12 = ( XX( 2) - XX( 1))** 2
Y12 = ( YY( 2) - YY( 1))** 2
DET = DSQRT ( X12 + Y12)/ 2.D0
C
T1 = ALF* ( 2.D0* TT( 1)/ 3.D0 + TT( 2)/ 3.D0 - TB )* DET
T2 = ALF* ( TT( 1)/ 3.D0 + 2.D0* TT( 2)/ 3.D0 - TB )* DET
C
F( NDE( IEL,1)) = F( NDE( IEL,1)) + T1
F( NDE( IEL,2)) = F( NDE( IEL,2)) + T2
C
C ------------------------------------------------
IF ( NA( 3).EQ.3 ) THEN
C
X23 = ( XX( 3)-XX( 2))** 2
Y23 = ( YY( 3)-YY( 2))** 2
DET = DSQRT( X23 + Y23)/ 2.D0
C
T2 = ALF* ( 2.D0* TT( 2)/ 3.D0 + TT( 3)/ 3.D0 - TB )* DET
T3 = ALF* ( TT( 2)/ 3.D0 + 2.D0* TT( 3)/ 3.D0 - TB )* DET
C
F( NDE( IEL,2))= F( NDE( IEL,2)) + T2
F( NDE( IEL,3))= F( NDE( IEL,3)) + T3
C
ENDIF
C ------------------------------------------------
IF ( NA( 3).EQ.4 ) THEN
C
X41 = ( XX( 4)- XX( 1))** 2
Y41 = ( YY( 4)- YY( 1))** 2
DET = DSQRT( X41+Y41)/ 2.D0
C
T1 = ALF* ( 2.D0* TT( 1)/ 3.D0 + TT( 4)/ 3.D0 - TB)* DET
T4 = ALF* ( TT( 1)/ 3.D0 + 2.D0* TT( 4)/ 3.D0 - TB)* DET
F( NDE( IEL,1)) = F( NDE( IEL,1)) + T1
F( NDE( IEL,4)) = F( NDE( IEL,4)) + T4
C
ENDIF
C ------------------------------------------------
C
ENDIF
C ------------------------------------------------
IF ( NA( 1).EQ.2.AND.NA( 2).EQ.3) THEN
C
X23 = ( XX( 3) - XX( 2))** 2
Y23 = ( YY( 3) - YY( 2))** 2
DET = DSQRT( X23 + Y23)/ 2.D0
C
T2 = ALF* ( 2.D0* TT( 2)/ 3.D0 + TT( 3)/ 3.D0 - TB)* DET
T3 = ALF* ( TT( 2)/ 3.D0 + 2.D0* TT( 3)/ 3.D0 - TB)* DET
F( NDE( IEL,2)) = F( NDE( IEL,2)) + T2
F( NDE( IEL,3)) = F( NDE( IEL,3)) + T3
C
C --------------------------------------
IF ( NA( 3).EQ.4) THEN
C
X34 = ( XX( 3) - XX( 4))** 2
Y34 = ( YY( 3) - YY( 4))** 2
DET = DSQRT ( X34 + Y34)/ 2.D0
C
T3 = ALF* ( 2.D0* TT( 3)/ 3.D0 + TT( 4)/ 3.D0 - TB)* DET
T4 = ALF* ( TT( 3)/ 3.D0 + 2.D0* TT( 4)/ 3.D0 - TB)* DET
C
F( NDE( IEL,3)) = F( NDE( IEL,3)) + T3
F( NDE( IEL,4)) = F( NDE( IEL,4)) + T4
C
END IF
C --------------------------------------
C
END IF
C ------------------------------------------------
IF (( NA( 1).EQ.3.AND.NA( 2).EQ.4)
1 .OR.( NA( 1).EQ.1.AND.NA( 2).EQ.3))
C
1THEN
C
X34 = ( XX( 3) - XX( 4))** 2
Y34 = ( YY( 3) - YY( 4))** 2
DET = DSQRT( X34 + Y34)/ 2.D0
C
T3 = ALF* ( 2.D0* TT( 3)/ 3.D0 + TT( 4)/ 3.D0 - TB )* DET
T4 = ALF* ( TT( 3)/ 3.D0 + 2.D0* TT( 4)/ 3.D0 - TB )* DET
C
F( NDE( IEL,3)) = F( NDE( IEL,3)) + T3
F( NDE( IEL,4)) = F( NDE( IEL,4)) + T4
C
IF ( NA( 3).EQ.4.AND.NA( 1).EQ.1) THEN
C
X41 = ( XX( 4) - XX( 1))** 2
Y41 = ( YY( 4) - YY( 1))** 2
DET = DSQRT( X41 + Y41)/ 2.D0
C
T1 = ALF* ( 2.D0* TT( 1)/ 3.D0 + TT( 4)/ 3.D0 - TB )* DET
T4 = ALF* ( TT( 1)/ 3.D0 + 2.D0* TT( 4)/ 3.D0 - TB )* DET
F( NDE( IEL,1)) = F( NDE( IEL,1)) + T1
F( NDE( IEL,4)) = F( NDE( IEL,4)) + T4
C
ENDIF
C ------------------------------------------------
C
ENDIF
C ----------------------------------------------------------
IF ( NA( 1).EQ.1.AND.NA( 2).EQ.4) THEN
C
X41 = ( XX( 4) - XX( 1))** 2
Y41 = ( YY( 4) - YY( 1))** 2
DET = DSQRT( X41 + Y41)/ 2.D0
C
T1 = ALF* ( 2.D0* TT( 1)/ 3.D0 + TT( 4)/ 3.D0 - TB )* DET
T4 = ALF* ( TT( 1)/ 3.D0 + 2.D0* TT( 4)/ 3.D0 - TB )* DET
F( NDE( IEL,1)) = F( NDE( IEL,1)) + T1
F( NDE( IEL,4)) = F( NDE( IEL,4)) + T4
C
ENDIF
C ----------------------------------------------------------
C
300 CONTINUE
C ----------------------------------------------------------
C
C ----- Simultaneous Eq. by Modified Galerkin Formulation ---
C
JS = 1
IF ( ISL.EQ.1) THEN
C
C ------------------------------------------------
DO 500 I = 1, NNP
C
JL = KCC( I)
C --------------------------------------
DO 600 J = JS, JL
C
KPR = NCW( 1,J)
KCL = NCW( 2,J)
A( KPR) = PP( KPR)/ DT + PHI* D( KPR)
B( I) = B( I) + ( PP( KPR)/ DT - ( 1.D0 - PHI)* D( KPR))* C( KCL)
C
600 CONTINUE
C --------------------------------------
B( I) = B( I) - F( I)
JS = JL+1
C
500 CONTINUE
C ---------------------------------------------------------------
C
CALL SLVE ( B, CNW, A, NCW, KCC, NBC, NNP, NNZ, NZ2 )
C
ELSE
C ---------------------------------------------------------------
DO 700 I = 1, NNP
C
JL = KCC( I)
FW = 0.D0
C
C ------------------------------------------------
DO 800 J = JS, JL
C
KPR = NCW( 1,J)
KCL = NCW( 2,J)
FW = FW + PP( KPR)
B( I)= B( I) + ( PP( KPR)/ DT - D( KPR))* C( KCL )
C
800 CONTINUE
C ------------------------------------------------
C
IF ( NBC( I).NE.9999) CNW( I)=( B( I)- F( I))/ FW* DT
JS = JL + 1
C
700 CONTINUE
C ----------------------------------------------------------
C
ENDIF
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) ! Attention OK
C
C ------------------------------------------------
DO 100 I = 1, 10
C
SNNP( I) = 0.D0
100 CONTINUE
C ------------------------------------------------
DO 200 I = 1, NNP
C
NQ = NC1( I)
CT = 0.D0
SPAL = 0.D0
C --------------------------------------
DO 250 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))
4 - X( NDE( JJ,1))* Y( NDE( JJ,4))
C
5 - X( NDE( JJ,2))* Y( NDE( JJ,1))
6 - X( NDE( JJ,4))* Y( NDE( JJ,2))
7 + X( NDE( JJ,2))* Y( NDE( JJ,3))
8 + X( NDE( JJ,3))* Y( NDE( JJ,4))
C
9 + X( NDE( JJ,4))* Y( NDE( JJ,2))
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)))
C
D / 2.D0
C
CT = CT + CE( JJ,JI)* SNNP( J)
SPAL = SPAL + SNNP( J)
C
250 CONTINUE
C --------------------------------------
C
CI( I) = CT/ SPAL
C
200 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VPRS ( X, Y, P, DLP, DI, NDE, NC1, NC2,
1 NC3, SU, SV, NBU, NBV, DPX, DPY, 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), NC1( NNP), NC2( NNP,10),
2 NC3( NNP,10), SU ( NNP), SV ( NNP), NBU( NNP),
3 NBV( NNP), DPX( NEL,4), DPY( NEL,4)
C
COMMON/CST/ DT, VCF, PHI, TDF, FT, ALF, TB, LPE, LOP
C
C ------------------------------------------------------------------
C
CALL DPXY ( X, Y, DLP, NDE, DPX, DPY, NNP, NEL )
C
C ------------------------------------------------------------------
C
CALL VLND ( DI, DPX, 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, DPY, 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
RETURN
END
C **********************************************************************