戻る

C *******************************************************************C
C C
C H F 2 D - F E R C
C C
C Thermal Fluid Analysis Software by Finite Element Method C
C C
C ( Rectangular Linear Element ) C
C C
C ( V.4.7 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C *******************************************************************C
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, ITB, IEB, KC1, MTR, IBP, IBV, IBT, LBC,
1 NTB, NCN, NEB
C
C ***** Problem : 121 Nodes / 100 Elements ***
C
DIMENSION X ( 121), Y ( 121), U ( 100), V ( 100),
1 UO ( 100), VO ( 100), UN ( 121), VN ( 121),
2 NDE( 100,4), ARS( 100), MTR( 100), LBC( 40),
3 P1 ( 541), P2 ( 541), DP ( 541), D ( 541),
4 FVC( 121), RKX( 100), RKY( 100),
5 NC1( 961), NC2( 2,961), NCW( 541)
C
DIMENSION PSI( 121), PSO( 121), IBP( 121), VOR( 121),
1 VRO( 121), IBV( 121), NEB( 121,10), TEP( 121),
2 TPO( 121), IBT( 121)
C
DIMENSION DSH( 40), SH ( 40), ITB( 121), IEB( 100),
1 NTB( 40), NCN( 121), KC1( 121), QWK( 100,4),
2 QE ( 100), ARW( 121), HXX( 12), HYY( 12),
3 FND( 121), FEL( 100)
C
DIMENSION WNU( 100,16), WNV( 100,16)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /TIM/ DT, LPM, LOP, TSP
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, PRN, RE
COMMON /FLW/ NVC, IVC( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 NTC, ITC( 50), TPC( 50), COF( 100)
C
OPEN ( 7, FILE='HFE_R_INP10.DAT' )
OPEN (10, FILE='HF2D_FER_PLT.DAT')
C
CALL TTLE
C
CALL INPT ( X, Y, U, V, DSH, SH, NDE, ITB, IEB,
1 NTB, MTR, LBC, NNZ, NZ2, ICN, NE2, HXX, HYY,
2 HKE, FND )
C
LOP = 0
T = 0.D0
C
CALL INSB ( X, Y, PSI, PSO, VOR, VRO, TEP, TPO, U,
1 V, UN, VN, MTR, RKX, RKY, QE, NC1, NDE,
2 IBP, IBV, IBT, LBC, NCN, NEB, NNZ, ICN, NE2,
3 HXX, HYY, UO, VO, FND, FEL, ARW )
C
C ----- P1 - Matrix ---
C
CALL PMAT ( P1, NC1, NCW, X, Y, NDE, KC1, NC2, NNZ, NZ2, ARS )
C
C ----- DP & D - Matrix ---
C
CALL DMAT ( X, Y, RKX, RKY, DP, NDE, NCW, NNZ, 1 )
C
CALL DMAT ( X, Y, RKX, RKY, D, NDE, NCW, NNZ, 2 )
C
C ----- VOR on a wall ---
C
CALL NBCN ( SH, DSH, PSI, VRO, LBC, NTB, NNP, LUV )
C
9999 LOP = LOP + 1
T = T + DT
C
C ----- Stream function ---
C
CALL LPSI ( U, V, UN, VN, PSI, PSO, IBP, VRO, DP,
1 RKX, RKY, NCN, FVC, NC2, KC1, NDE, MTR, X,
2 Y, IBV, NNZ, NZ2, NE2, HXX, HYY, FND, FEL,
3 ARW )
C
C ----- VOR on walls ---
C
CALL NBCN ( SH, DSH, PSI, VOR, LBC, NTB, NNP, LUV )
C
IF ( ICN.EQ.0) GO TO 100
C
C ----- To modify D & D2 ---
C
CALL MODD ( D, P1, P2, X, Y, RKX, RKY, NCW, NDE, NNZ, 3 )
C
CALL LPTP ( UN, VN, TEP, TPO, IBT, D, P2, QWK, FVC, QE,
1 NC2, KC1, NDE, X, Y, WNU, WNV, NNZ, NZ2 )
C
100 CONTINUE
C
C ----- To modify D & D2 ---
C
CALL MODD ( D, P1, P2, X, Y, RKX, RKY, NCW, NDE, NNZ, 2 )
C
CALL LPVR ( UN, VN, VOR, VRO, IBV, D, P2, FVC, NC2,
1 KC1, NDE, TEP, QWK, X, Y, NNZ, NZ2, ICN,
2 DMX, DN, WNU, WNV )
C
C ----- Stationarity Check ---
C
CALL CHCK ( PSI, PSO, VOR, VRO, TEP, TPO, U, V, UN, VN,
1 ICK, ICN, UO, VO, ARS, VRU, *999 )
C
IF ( LOP.LT.LPM) GO TO 9999
C
999 CONTINUE
C
CALL PLTD ( PSI, VOR, TEP, UN, VN, X, Y, NDE, ICN )
C
CLOSE ( 7)
CLOSE (10)
C
STOP
END
C **********************************************************************
C
SUBROUTINE BCS10 ( DSH, SH, NTB )
C
C ----- Attention : Boundary condition for Cavity problem of 121/100 ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NTB
C
DIMENSION DSH( LUV), SH( LUV), NTB( LUV)
C
C ----- Attention NTB ( LUV) ---
C
COMMON /NNL/ NNP, NEL, LUV
C
DO 100 I= 1, LUV
C
SH ( I) = 30.D0 ! ( = 3/l )
DSH( I) = 300.D0 ! ( = 3/l**2 )
C
100 CONTINUE
C
NTB( 1) = 13
NTB( 2) = 13
NTB( 3) = 24
NTB( 4) = 35
NTB( 5) = 46
NTB( 6) = 57
NTB( 7) = 68
NTB( 8) = 79
NTB( 9) = 90
NTB(10) =101
C
NTB(11) =101
NTB(12) =101
NTB(13) =102
NTB(14) =103
NTB(15) =104
NTB(16) =105
NTB(17) =106
NTB(18) =107
NTB(19) =108
NTB(20) =109
C
NTB(21) =109
NTB(22) =109
NTB(23) = 98
NTB(24) = 87
NTB(25) = 76
NTB(26) = 65
NTB(27) = 54
NTB(28) = 43
NTB(29) = 32
NTB(30) = 21
C
NTB(31) = 21
NTB(32) = 21
NTB(33) = 20
NTB(34) = 19
NTB(35) = 18
NTB(36) = 17
NTB(37) = 16
NTB(38) = 15
NTB(39) = 14
NTB(40) = 13
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CHCK ( PSI, PSO, VOR, VRO, TEP, TPO, U, V,
1 UN, VN, ICK, ICN, UO, VO, ARS, VRU,
2 * )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION PSI( NNP), PSO( NNP), VOR( NNP), VRO( NNP),
1 TEP( NNP), TPO( NNP), U ( NEL), V ( NEL),
2 UN ( NNP), VN ( NNP), UO ( NEL), VO ( NEL),
3 ARS( NEL)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /TIM/ DT, LPM, LOP, TSP
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, PRN, RE
COMMON /VPT/ VARV( 2000), VART( 2000), VRUV( 2000), VARP( 2000)
C
EPS = 0.1D0
C
ICK = 0
C ------------------------------------------------
C
VRMX = -100.D0
DO 100 I = 1, NNP
C
VRMX = DMAX1( VRMX, DABS( PSI( I)))
100 CONTINUE
VMX = VRMX* 0.001D0 ! 0.1(%) Up
C ------------------------------------------------
C
DMX = -100.D0
DO 200 I = 1, NNP
C
IF ( DABS( PSI( I)).LE.VMX) GO TO 200
DMX = DMAX1 ( DMX, DABS((PSI( I)-PSO( I))/ PSI( I)))
C
200 CONTINUE
VARP( LOP) = DMX* 100.D0
C ------------------------------------------------
C
VRMX = -100.D0
DO 300 I = 1, NNP
C
VRMX = DMAX1( VRMX, DABS( VOR( I)))
300 CONTINUE
VMX = VRMX* 0.001D0 ! 0.1(%) Up
C ------------------------------------------------
C
DMX = -100.D0
DO 400 I = 1, NNP
C
IF ( DABS( VOR( I)).LE.VMX) GO TO 400
DMX = DMAX1 ( DMX, DABS((VOR( I)-VRO( I))/ VOR( I)))
C
400 CONTINUE
VARV( LOP) = DMX* 100.D0
C ------------------------------------------------
C
IF ( ICN.EQ.1) THEN
C
TPMX = -100.D0
DO 500 I = 1, NNP
C
TPMX = DMAX1( TPMX, DABS( TEP( I)))
500 CONTINUE
C
TMX = TPMX* 0.001D0 ! 0.1(%) Up
C ------------------------------------------------
DO 600 I = 1 ,NNP
C
IF ( DABS( TEP( I)).LE.TMX) GO TO 600
TMX = DMAX1( TMX, DABS((TEP( I) - TPO( I))/ TEP( I)))
C
600 CONTINUE
VART( LOP) = TMX* 100.D0
C ------------------------------------------------
C
END IF
C
C ------------------------------------------------
VMIN = 0.01D0
VRU = -1000.D0
C
DO 700 I = 1, NEL
C
IF ( DABS( U( I)).LE.VMIN) GO TO 700
VRU = DMAX1 ( VRU, DABS(( U( I) - UO( I))/ U( I)))
C
700 CONTINUE
VRUV( LOP) = VRU* 100.D0
C ------------------------------------------------
C
IF ( ICN.EQ.1)
1WRITE(6,2000) LOP, VARV( LOP), VART( LOP), VRUV( LOP)
2000 FORMAT(' * LOP =', I4,' Var_Vr=', F8.3,' (%)',' Var_T =',F8.3,
1 ' (%)',' Var_Vel=',F8.3,' (%)')
C
IF ( ICN.EQ.0)
1WRITE(6,2100) LOP, VARP( LOP), VARV( LOP), VRUV( LOP)
2100 FORMAT(' * LOP =', I4,' Var_P=', F8.3,' (%)',' Var_Vr=', F8.3,
1 ' (%)',' Var_Vel=',F8.3,' (%)')
C
C ----- Stationarity Check ------
C
IF ( ICN.EQ.1) THEN
C
IF ( VART( LOP).LE.EPS .AND. VARV( LOP).LE.EPS) THEN
WRITE(6,*) ' '
WRITE(6,*) ' ** Stationary condition reached **'
IEND = 1
IF ( IEND.EQ.1) WRITE(6,2200) LOP, VARV( LOP), VART( LOP)
2200 FORMAT (' Stationary Condition Reached ',' LOP =',I5,' * VARV ='
1 ,F9.4,' * VART=', F9.4 )
RETURN 1
END IF
C
END IF
C ------------------------------------------------
IF ( ICN.EQ.0 ) THEN
C
IF ( VARV( LOP).LE.EPS) THEN
WRITE(6,*) ' '
WRITE(6,*) ' ** Stationary condition reached **'
IEND = 1
IF ( IEND.EQ.1) WRITE(6,2300) LOP, VARV( LOP)
2300 FORMAT (/' Stationary Condition Reached ',' LOP =',I5,' * VARV ='
1 ,F9.4 )
RETURN 1
END IF
C
END IF
C --------------------------------------
IF ( LOP.GE.LPM) THEN
WRITE(6,2400)
2400 FORMAT (/,' LOOP Limit Over')
RETURN 1
END IF
C
C ------------------------------------------------
DO 800 I = 1, NNP
C
IF ( DABS( UN( I)).GE.10.D0) GO TO 900
IF ( DABS( VN( I)).GE.10.D0) GO TO 900
C
800 CONTINUE
C ------------------------------------------------
C
CALL MOVR( PSI, 1.D0, PSO, NNP )
C
CALL MOVR( VOR, 1.D0, VRO, NNP )
C
IF ( ICN.EQ.1) CALL MOVR( TEP, 1.D0 , TPO, NNP )
C
CALL MOVR( U, 1.D0, UO, NEL )
C
CALL MOVR( V, 1.D0, VO, NEL )
C
RETURN
C
900 WRITE(6,2500)
2500 FORMAT(/,' *** U. or V.GE.10 * STOP ***')
STOP
C
END
C **********************************************************************
C
SUBROUTINE CKAK ( RKX, RKY, U, V, MTR, NE2, HXX, HYY,
1 UN, VN, NDE, FND, FEL )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, MTR
C
DIMENSION RKX( NEL), RKY( NEL), U ( NEL), V ( NEL),
1 MTR( NEL), HXX( NE2), HYY( NE2), UN ( NNP),
2 VN ( NNP), NDE( NEL,4 ), FND( NNP), FEL( NEL)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /TIM/ DT, LPM, LOP, TSP
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, PRN, RE
COMMON /FLW/ NVC, IVC( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 NTC, ITC( 50), TPC( 50), COF(100)
C
C ----- To estimate Corection Cofficients by Error Analysis ---
C
C ----- FND( I) : Nodal-"f" Coefficients ---
C
L = 0
C ------------------------------------------------
DO 100 J = 1, NE2 - 1
DO 100 I = 1, NE2 - 1
C
L = L + 1
C
AX = HXX( I)/ HXX( I+1)
AY = HYY( J+1)/ HYY( J)
C
BX = UN( L)* DT/ HXX(I+1) ! bx
BY = VN( L)* DT/ HYY( J) ! by
C
RX = COF( 1)* DT/( HXX( I+1)** 2) ! rx
RY = COF( 1)* DT/( HYY( J)** 2) ! ry
C
C ----- Attention : Correction Coefficient --------------------------
C
FND( L) = 1.D0 + 0.5D0* ( BX + BY)** 2/( RX + RY)
C
1 + (1.D0/ 6.D0)* ( ( AX-1.D0)* BX + ( AY-1.D0)* BY)/ ( RX + RY)
C
C -------------------------------------------------------------------
C
100 CONTINUE
C ------------------------------------------------
C
C ----- Correction Coefficients by Element -----------------
C
DO 200 I = 1, NEL
C
C ----- Averaged by element ---
C
FEL( I) = 0.25D0* ( FND( NDE( I,1)) + FND( NDE( I,2))
1 + FND( NDE( I,3)) + FND( NDE( I,4)))
C
200 CONTINUE
C ----------------------------------------------------------
C
C ----- Attention : Modified Diffusivity ---------
C
DO 300 I = 1, NEL
C
RKX( I) = COF( 1)* FEL( I)
RKY( I) = RKX( I)
C
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CLAR8 ( X, V, N )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION X( N)
C
DO 100 I = 1, N
C
X ( I) = V
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CLAR2 ( IX2, IV2, N )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 IX2( N)
C
DO 100 I = 1, N
C
IX2( I) = IV2
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DCIJ ( IA, JA, NJ, NIX, NJX )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 IA( NIX)
CHARACTER*4 ANIX, ANJX
C
DIMENSION JA( 2,NJX)
C
C ----- NI, NJ : Current length of IA, JA
C
C NIX, NJX : Max. length of IA, JA
C
DATA ANIX /'NIX '/, ANJX /'NJX '/
C
C ----- To count every ROW's element ---
C
NI = 1
C ------------------------------------------------
DO 100 I = 1, NJ
C
IF ( NI.EQ.JA( 1,I)) GO TO 101
C
C ----- To set IA( NI ) ---
C
IA( NI) = I - 1
NI = JA( 1,I)
IF ( NI.GT.NIX) GO TO 601
C
101 JA( 1,I) = I
C
100 CONTINUE
C ------------------------------------------------
C
C ----- To set last ROW's counter ---
C
IA( NI) = NJ
C
C ----- To set ( I,J) in JA for I > J ---
C
C ------------------------------------------------
DO 200 I = 2, NI
C
I1 = I - 1
IS = IA( I1)
C
C --------------------------------------
DO 300 J = 1, I1
C
JS = 1
IF ( J.GT.1) JS = IA( J-1) + 1
JL = IA( J)
C
C ----------------------------
DO 400 K = JS, JL
C
IF ( JA( 2, K).NE.I ) GO TO 400
K1 = K
GO TO 410
C
400 CONTINUE
C ----------------------------
C
C ------ ( I, J) not exist ---
C
GO TO 300
C
410 CONTINUE
C
C ----- To shift JA ---
C
NJIS = NJ - IS
C
C ----------------------------
DO 500 L = 1, NJIS
C
JA( 1,NJ-L+2) = JA( 1,NJ-L+1)
JA( 2,NJ-L+2) = JA( 2,NJ-L+1)
C
500 CONTINUE
C ----------------------------
NJ = NJ + 1
IF ( NJ.GT.NJX ) GO TO 602
C
C ----- To increase IA ---
C
DO 550 L = I, NI
C
IA( L) = IA( L) + 1
C
550 CONTINUE
C
C ----- To insert ( I,J) element ---
C
IS = IS + 1
JA( 1,IS) = JA( 1,K1)
JA( 2,IS) = J
C
300 CONTINUE
C --------------------------------------
C
200 CONTINUE
C ------------------------------------------------
RETURN
C
601 WRITE(6,2000) ANIX, NIX, ANIX, NI
STOP
C
602 WRITE(6,2000) ANJX, NJX, ANJX, NJ
2000 FORMAT(' ** Error ( DCIJ ) ',A3,'(',I5,') too small',A2,'=',I5)
C
END
C **********************************************************************
C
SUBROUTINE DCNL ( AV, I, J, A, L, K, NT )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION A( NT), L( NT)
C
C ----- Attention --- NO. OF NDE.LE.32768 ---
C
IF ( I.GT.J) RETURN
M = I* 10300 + J
IE = K
N1 = 1
C
IF ( K.EQ.0) N2 = K + 1
IF ( K.EQ.0) GOTO 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 600
IF ( M - L( N)) 200, 300, 400
C
200 CONTINUE
N2 = N
N = N2 - ( N2 - N1)/2
C
IF ( N.NE.N2) GO TO 100
IF ( M - L( N1)) 110, 115, 120
110 N = N1
120 CONTINUE
C
C -------------------------------------------
C
CALL SUBS ( AV, M, N, A, L, K, NT )
C
C -------------------------------------------
RETURN
C
115 N = N1
GO TO 300
C
550 N = N2
300 CONTINUE
A( N) = A( N) + AV
RETURN
C
400 CONTINUE
N1 = N
N = N1 + ( N2 - N1 )/ 2
IF ( N.NE.N1) GO TO 100
IF ( M - L( N2)) 500, 550, 600
500 N = N2
GO TO 120
C
600 K = K + 1
IF ( K.GT.NT) WRITE(6,2000)
2000 FORMAT(' ** Error ** ''NT'' IN ''SOR'' IS TOO SMALL **')
A( K) = AV
L( K) = M
RETURN
C
100 CONTINUE
C ------------------------------------------------
C
WRITE(6,2100) K
2100 FORMAT(' * Error --- Do LOOP Strange, K =',I10)
C
STOP
END
C **********************************************************************
C
SUBROUTINE DMAT ( X, Y, RKX, RKY, D, NDE, NCW, NNZ, NDM )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE
C
DIMENSION X( NNP), Y ( NNP), RKX( NEL), RKY( NEL),
1 D( NNZ), NDE( NEL,4), NCW( NNZ), A1 ( 4,4),
2 N( 4)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, PRN, RE
COMMON /FLW/ NVC, IVC( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 NTC, ITC( 50), TPC( 50), COF( 100)
C
CALL CLAR8 ( D, 0.D0, NNZ )
C
NZK = NNZ
C
C ----------------------------------------------------------
DO 100 IEL = 1, NEL
C
N( 1) = NDE( IEL,1)
N( 2) = NDE( IEL,2)
N( 3) = NDE( IEL,3)
N( 4) = NDE( IEL,4)
C
A = DABS( Y( NDE( IEL,1)) - Y( NDE( IEL,3)))/ 2.D0
B = DABS( X( NDE( IEL,1)) - X( NDE( IEL,3)))/ 2.D0
C
IF ( NDM-2) 2000, 2100, 2200
C
2000 CONTINUE
W1 = A/ B
W2 = B/ A
GO TO 2500
C
C ----- For Vorticity Equation ---
C
2100 CONTINUE
C
W1 = RKX( IEL)* A/ B ! RKX <==> Pr
W2 = RKY( IEL)* B/ A ! RKY <==> Pr
GOTO 2500
C
C ----- For Energy Equation -----------
C
2200 CONTINUE
TKAK = 1.D0 ! TKAK <==> 1
C =====================================
C
W1 = TKAK* A/ B
W2 = TKAK* B/ A
C
C ----- W1 = A / B, W2 = B / A ---
C
2500 CONTINUE
C
A1( 1,1) = W1/ 3.D0 + W2/ 3.D0
A1( 1,2) = - W1/ 3.D0 + W2/ 6.D0
A1( 1,3) = - W1/ 6.D0 - W2/ 6.D0
A1( 1,4) = W1/ 6.D0 - W2/ 3.D0
C
A1( 2,1) = - W1/ 3.D0 + W2/ 6.D0
A1( 2,2) = W1/ 3.D0 + W2/ 3.D0
A1( 2,3) = W1/ 6.D0 - W2/ 3.D0
A1( 2,4) = - W1/ 6.D0 - W2/ 6.D0
C
A1( 3,1) = - W1/ 6.D0 - W2/ 6.D0
A1( 3,2) = W1/ 6.D0 - W2/ 3.D0
A1( 3,3) = W1/ 3.D0 + W2/ 3.D0
A1( 3,4) = - W1/ 3.D0 + W2/ 6.D0
C
A1( 4,1) = W1/ 6.D0 - W2/ 3.D0
A1( 4,2) = - W1/ 6.D0 - W2/ 6.D0
A1( 4,3) = - W1/ 3.D0 + W2/ 6.D0
A1( 4,4) = W1/ 3.D0 + W2/ 3.D0
C
C ------------------------------------------------
DO 150 I = 1, 4
DO 150 J = 1, 4
C ------------------------------------------------
C
II = N ( I)
JJ = N ( J)
AV = A1( I,J)
C
CALL DCNL ( AV, II, JJ, D, NCW, NZK, NNZ )
C
150 CONTINUE
C ------------------------------------------------
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DSLV ( FVC, C, TSM, NCC, KCC, NBC, NNZ, NZ2,
1 OMG, ERR )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 KCC, NBC
C
DIMENSION FVC( NNP), C ( NNP), TSM( NNZ), NCC( 2,NZ2),
1 KCC( NNP), NBC( NNP)
C
COMMON /NNL/ NNP, NEL, LUV
C
C ----- S.O.R. Method ---
C
C --------------------------------------
CMAX = 0.D0
DO 100 I = 1, NNP
C
AC = DABS( C( I))
IF ( AC.GE.CMAX ) CMAX = AC
C
100 CONTINUE
C --------------------------------------
EPS = CMAX* 0.01D0
LMR = 3* NNP
KTL = 2* NNP
NTL = 0
C
900 NTL = NTL + 1
EMX = 0.D0
C
C ------------------------------------------------
DO 200 I = 1, NNP
C
IF ( NBC( I).EQ.9999) GO TO 200
IF ( I.EQ.1) GO TO 210
MS = KCC( I-1) + 1
GOTO 220
C
210 MS = 1
220 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) + OMG* DC
C
TEMP = DABS( C( I))
IF ( TEMP.LE.EPS) GO TO 200
TEMP = DABS( DC/ C( I))
IF ( TEMP.LE.EMX) GO TO 200
EMX = TEMP
C
200 CONTINUE
C ------------------------------------------------
C
IF ( EMX.LE.ERR) RETURN
IF ( NTL.GE.KTL) RETURN
IF ( NTL.GE.LMR) GO TO 900
GOTO 900
C
1000 WRITE(6,2000) I, MS, ML
2000 FORMAT ('1',' ** Error in "DSLV" ** ',3I6)
RETURN
END
C **********************************************************************
C
SUBROUTINE FVCTN ( D, CO, F, NC2, KCC, NNZ, NZ2 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 KCC
C
DIMENSION D( NNZ), CO( NNP), F( NNP), NC2( 2,NZ2), KCC( NNP)
C
COMMON /NNL/ NNP, NEL, LUV
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
JS = 1
JL = KCC( I)
IF ( I.GT.1) JS = KCC( I-1) + 1
FW = - F( I)
C
C --------------------------------------
DO 150 J = JS, JL
C
KPTR = NC2( 1,J)
KCOL = NC2( 2,J)
FW = FW + D( KPTR)* CO( KCOL)
C
150 CONTINUE
C --------------------------------------
C
F( I) = FW
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FVOR ( VOR, NDE, FVC, X, Y )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE
C
DIMENSION VOR( NNP), FVC( NNP), X( NNP), Y( NNP), NDE( NEL,4)
C
COMMON /NNL/ NNP, NEL, LUV
C
CALL CLAR8 ( FVC, 0.D0, NNP )
C
DO 100 NE = 1, NEL
C
A = DABS( Y( NDE( NE,1)) - Y( NDE( NE,3)))/ 2.D0
B = DABS( X( NDE( NE,1)) - X( NDE( NE,3)))/ 2.D0
C
S = A * B / 9.D0
C
N1 = NDE( NE,1)
N2 = NDE( NE,2)
N3 = NDE( NE,3)
N4 = NDE( NE,4)
C
FVC( N1) = FVC( N1)
1 + ( 4.D0* VOR( N1) + 2.D0* VOR( N2) + VOR( N3)
2 + 2.D0* VOR( N4))* S
C
FVC( N2) = FVC( N2)
1 + ( 2.D0* VOR( N1) + 4.D0* VOR( N2) + 2.D0* VOR( N3)
2 + VOR( N4))* S
C
FVC( N3) = FVC( N3)
1 + ( VOR( N1) + 2.D0* VOR( N2) + 4.D0* VOR( N3)
2 + 2.D0* VOR( N4))* S
C
FVC( N4) = FVC( N4)
1 + ( 2.D0* VOR( N1) + VOR( N2) + 2.D0* VOR( N3)
2 + 4.D0*VOR( N4))* S
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INPT ( X, Y, U, V, DSH, SH, NDE, ITB,
1 IEB, NTB, MTR, LBC, NNZ, NZ2, ICN, NE2,
2 HXX, HYY, HKE, FND )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, ITB, IEB, MTR, LBC, NTB
C
CHARACTER*3 HEAD
C
DIMENSION X ( NNP), Y ( NNP), U ( NEL), V ( NEL),
1 DSH( LUV), SH ( LUV), NDE( NEL,4), ITB( NNP),
2 IEB( NEL), NTB( LUV), MTR( NEL), LBC( LUV),
3 HXX( NE2), HYY( NE2), FND( NNP)
C
DIMENSION IWK( 10), WKX( 10), IWKX( 10,10), WK( 10,10)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /TIM/ DT, LPM, LOP, TSP
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, PRN, RE
COMMON /FLW/ NVC, IVC( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 NTC, ITC( 50), TPC( 50), COF(100)
C
C ----- Attention : In case of 121 - Nodes & 100-Elements ---
C
C Y
C |
C 121 111
C -------
C | |
C | | 10 X 10 Meshes
C ------- ===> X
C 11 1
C
C X( 1)= 1, Y( 1)= 0 ; X( 11)= 0, Y( 11)= 0
C X(111)= 1, Y(111)= 1 ; X(121)= 0, Y(121)= 1
C
NEX = 10
NEY = 10
NE2 = NEX + 2
C
NPX = 0
NLX = 0
C
NVC = 0
NPC = 0
NTC = 0
NLUV = 0
C
CALL CLAR8 ( COF, 1.D0, 100 ) ! == > Pr
CALL CLAR8 ( FND, 1.D0, NNP )
C
CALL CLAR8 ( U, 0.D0, NEL )
CALL CLAR8 ( V, 0.D0, NEL )
C
OMG( 1) = 1.6D0
OMG( 2) = 1.1D0
OMG( 3) = 1.1D0
C
ERR( 1) = 0.0001D0
ERR( 2) = 0.0001D0
ERR( 3) = 0.0001D0
C
DLT ( 1) = 0.001D0
DLT ( 2) = 0.01D0
DLT ( 3) = 0.01D0
C
C ----- Time scheme parameter ---
C
TSP = 1.D0
C
C ----- TIM ---
C
READ(7,1000) DT, LPM, LOP1, TSPA, DLT1
1000 FORMAT( 10X,F12.0,5X,I5,5X,I5,5X,F10.0,5X,F10.0)
C
IF ( DABS( DLT1).GE.0.0001D0) DLT( 1) = DLT1
C
C ----- MTR ---
C
READ(7,1100) J, COF( 1), RA, ICN, PRN
1100 FORMAT( 10X,I5,5X,F10.0,5X,F10.0,4X,I1,10X,F10.0)
C
C ----- Attention Heat & Fluid Only ---
C
WRITE(6,2000) NNP, NEL, LUV
2000 FORMAT(' ( NNP =',I5,' NEL =',I5,' LUV =',I4,' )')
WRITE(6,*) ' '
WRITE(6,*) ' ** ICN = 1 : Heat & Fluid, ICN = 0 : Fluid Only'
C
*** READ(6,*) ICN
ICN = 1
** ICN = 0
WRITE(6,2050) ICN
2050 FORMAT(' ** ICN =',I3)
C
IF ( ICN.EQ.1) THEN
C
C ----- Prandtle number ---
C
PRN = 0.71D0
C ------------------
COF( 1) = PRN
C ------------------
C
WRITE(6,*) ' '
WRITE(6,*) ' ** Ra = ?'
*** READ (5,*) RA
RA = 1000.D0
C
WRITE(6,2100) PRN, RA
2100 FORMAT(' * Pr =', F6.2,' * Ra =', F8.1)
C
END IF
C
C --------------------------------------
IF ( ICN.EQ.0 ) THEN
WRITE(6,*) ' '
WRITE(6,*) ' ** Re = ?'
** READ (5,*) RE
RE = 100.D0
COF( 1) = 1.D0/ RE
C ------------------------
WRITE(6,2200) RE
2200 FORMAT(/,' *** Re =', F7.1)
C
END IF
C --------------------------------------
C
WRITE(6,*) ' '
WRITE(6,*) ' ** DT = ? * LOP = ?'
*** READ( 5,*) DT, LPM
DT = 0.01D0
LPM = 100
C
** DT = 0.1D0
** LPM = 200
C
WRITE(6,2300) DT, LPM, DLT1
2300 FORMAT(' * DT =',F6.3,' LOP =',I4,' EPS =', F8.5,/ )
C
C ----- JTN ---
C
DO 100 II = 1, NNP
C
READ(7,1200) IWK( 1), WK( 1,1), WK( 2,1)
1200 FORMAT( 10X,I5,5X,2F10.0)
C
NPX = NPX + 1
ITB( NPX) = IWK( 1)
X ( NPX) = WK ( 1,1)
Y ( NPX) = WK ( 2,1)
C
100 CONTINUE
C ----- ELM --------------------------------------
C
DO 200 II = 1, NEL
C
READ(7,1300) ( IWKX( 1,J ), J = 1, 6 )
1300 FORMAT( 10X,6I5)
IWKX( 1,6) = MAX0( IWKX( 1,6),1)
NLX = NLX + 1
IEB( NLX) = IWKX( 1,1)
NDE( NLX,1) = IWKX( 1,2)
NDE( NLX,2) = IWKX( 1,3)
NDE( NLX,3) = IWKX( 1,4)
NDE( NLX,4) = IWKX( 1,5)
MTR( NLX) = IWKX( 1,6)
C
200 CONTINUE
C ------------------------------------------------
C
C ----- PSI ---
C
121 READ(7,1400) HEAD
1400 FORMAT(A3)
C
IF ( HEAD.NE.'PSI') GO TO 140
BACKSPACE 7
C
READ(7,1500) ( IWK( I), WKX( I), I=1, 3 )
1500 FORMAT( 10X,3( I5,5X,F10.0))
C
C ------------------------------------------------
DO 300 I = 1, 3
C
IF ( IWK( I).EQ.0 .AND. DABS(WKX( I)).LE.0.0001D0) GO TO 300
NPC = NPC + 1
IPC( NPC) = IWK( I)
PSC( NPC) = WKX( I)
C
300 CONTINUE
C ------------------------------------------------
GOTO 121
C
C ----- TEP ---
C
140 BACKSPACE 7
141 READ(7,1400) HEAD
C
IF ( HEAD.NE.'TEP' ) GO TO 150
BACKSPACE 7
C
READ(7,1500) ( IWK( I), WKX( I), I = 1 ,3 )
C
C --------------------------------------
DO 400 I = 1, 3
C
IF ( IWK( I).EQ.0 .AND. DABS( WKX( I)).LE.0.0001D0) GO TO 400
NTC = NTC + 1
ITC( NTC) = IWK( I)
TPC( NTC) = WKX( I)
C
400 CONTINUE
C --------------------------------------
C
GO TO 141
C
C ----- Wall ---
C
150 BACKSPACE 7
C
C ***************************************
C
151 READ(7,1400) HEAD
C
IF ( HEAD.NE.'WAL' ) GO TO 160
BACKSPACE 7
C
READ(7,1600) ( IWK( I),I=1,10)
1600 FORMAT( 10X,10I5)
C
C --------------------------------------
DO 500 I = 1, 10
C
IF ( IWK( I).EQ.0) GO TO 500
NLUV = NLUV + 1
LBC( NLUV)= IWK( I)
C
500 CONTINUE
C --------------------------------------
GOTO 151
C
160 CONTINUE
C
IF ( ICN.EQ.0) THEN
C
C --------------------------------------
DO 600 I = 1, NTC
C
TPC( I) = 0.D0
600 CONTINUE
C --------------------------------------
GOTO 900
C
END IF
C
900 CONTINUE
C
C ----- Attention ------------------------------------------
C
C ----- For Rectangular meshes ---
C
NZ2 = 4* 4 + ( NEX-1)* ( NEY-1)* 9 + ( NEX-1)* 6* 2
1 + ( NEY-1)* 6* 2
NNZ = ( NZ2 + NNP)/ 2
LUV = NLUV
WRITE(6,2500) NNZ, NZ2
2500 FORMAT(' ** NNZ =',I6,' * NZ2 =',I6)
WRITE(6,*) ' '
C
IF ( NNP.EQ.121) CALL BCS10 ( DSH, SH, NTB )
C -----------------------------------------------------
DO 700 I = 2, NE2 - 1
C
HXX( I) = X( I-1) - X( I)
HYY( I) = Y(( I-1)* ( NE2-1)+1) - Y(( I-2)* ( NE2-1)+1)
C
700 CONTINUE
C -----------------------------------------------------
C
HXX( 1) = HXX( 2)
HYY( 1) = HYY( 2)
HXX( NE2) = HXX( NE2-1)
HYY( NE2) = HYY( NE2-1)
C
* WRITE(6,*) ' * NS2=', NE2
* WRITE(6,*) ' * X=', ( X(I), I=1, NE2)
* WRITE(6,*) ' * HXX=', ( HXX(I), I=1, NE2)
* WRITE(6,*) ' * HYY=', ( HYY(I), I=1, NE2)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INSB ( X, Y, PSI, PSO, VOR, VRO, TEP,
1 TPO, U, V, UN, VN, MTR, RKX,
2 RKY, QE, NC1, NDE, IBP, IBV, IBT,
3 LBC, NCN, NEB, NNZ, ICN, NE2, HXX,
4 HYY, UO, VO, FND, FEL, ARW )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, MTR, IBP, IBV, IBT, LBC, NCN, NEB
C
DIMENSION X ( NNP), Y ( NNP), PSI( NNP), PSO( NNP),
1 VOR( NNP), VRO( NNP), TEP( NNP), TPO( NNP),
2 U ( NEL), V ( NEL), UN ( NNP), VN ( NNP),
3 RKX( NEL), RKY( NEL), QE ( NEL), NC1( NNZ),
4 NEB( NNP,10), HXX( NE2), HYY( NE2), FND( NNP),
5 FEL( NEL)
C
DIMENSION NDE( NEL,4), IBP( NNP), IBV( NNP), UO ( NEL),
1 VO ( NEL), IBT( NNP), LBC( LUV), NCN( NNP),
2 MTR( NEL), ARW( NNP)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /FLW/ NVC, IVC( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 NTC, ITC( 50), TPC( 50), COF(100)
C
CALL CLAR8 ( ARW, 0.D0, NNP )
CALL CLAR8 ( QE, 0.D0, NEL )
CALL CLAR8 ( UO, 0.D0, NEL )
CALL CLAR8 ( VO, 0.D0, NEL )
C
CALL CLAR2 ( IBP, 1, NNP )
CALL CLAR2 ( IBV, 1, NNP )
CALL CLAR2 ( IBT, 1, NNP )
C
CALL CLAR8 ( PSI, 0.D0, NNP )
CALL CLAR8 ( PSO, 0.D0, NNP )
CALL CLAR8 ( VOR, 0.D0, NNP )
CALL CLAR8 ( VRO, 0.D0, NNP )
C
CALL CLAR8 ( TEP, 0.D0, NNP )
CALL CLAR8 ( TPO, 0.D0, NNP )
C
IF ( NPC.LE.0) GO TO 110
C
C --------------------------------------
DO 100 I = 1, NPC
C
IBP( IPC( I)) = 9999
PSI( IPC( I)) = PSC( I)
C
100 CONTINUE
C --------------------------------------
110 IF ( NVC.LE.0) GO TO 121
C
DO 200 I = 1, NVC
C
IBV( IVC( I)) = 9999
VOR( IVC( I)) = VRC( I)
C
200 CONTINUE
C --------------------------------------
121 IF ( NTC.LE.0) GO TO 130
C
DO 300 I = 1, NTC
C
IBT( ITC( I)) = 9999
TEP( ITC( I)) = TPC( I)
C
300 CONTINUE
C --------------------------------------
130 IF ( LUV.LE.0) GO TO 150
C
DO 400 I = 1, LUV
C
IBV( LBC( I)) = 9999
400 CONTINUE
C --------------------------------------
C
150 CONTINUE
C
C ----- No. of meshes at each NDE ---
C
DO 500 I = 1, NNP
C
NCN( I) = 0
500 CONTINUE
C --------------------------------------
DO 600 I = 1, NEL
DO 600 J = 1, 4
C --------------------------------------
C
NE = NDE( I,J)
NCN ( NE) = NCN( NE) + 1
IF ( NCN ( NE).GT.10) WRITE(6,2000) NCN ( NE)
2000 FORMAT(' *** NO. OF ELEMENTS AT NDE=',I5,' GT 10 *** STOP' )
IF ( NCN ( NE).GT.10) STOP
NEB( NE, NCN(NE)) = I
C
600 CONTINUE
C ---------------------------------------------------------------------
C
CALL VELC ( PSI, U, V, NDE, UN, VN, X, Y, IBV, NCN, ICN, ARW )
C
C ---------------------------------------------------------------------
DO 700 I = 1, NEL
C
RKX( I) = COF( MTR( I)) ! COF : Pr
RKY( I) = COF( MTR( I)) ! COF : Pr
C
700 CONTINUE
C ------------------------------------------------
C
CALL MOVR( PSI, 1.D0, PSO, NNP)
C
CALL MOVR( VOR, 1.D0, VRO, NNP)
C
CALL MOVR( TEP, 1.D0, TPO, NNP)
C
CALL CKAK ( RKX, RKY, U, V, MTR, NE2, HXX, HYY, UN, VN,
1 NDE, FND, FEL )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE LPSI ( U, V, UN, VN, PSI, PSO, IBP,
1 VOR, DP, RKX, RKY, NCN, FVC, NC2,
2 KCC, NDE, MTR, X, Y, IBV, NNZ,
3 NZ2, NE2, HXX, HYY, FND, FEL, ARW )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 IBP, IBV, KCC, NDE, MTR, NCN
C
DIMENSION U ( NEL), V ( NEL), UN ( NNP), VN ( NNP),
1 PSI( NNP), PSO( NNP), IBP( NNP), VOR( NNP),
2 DP ( NNZ), RKX( NNP), RKY( NNP), NCN( NNP),
3 FVC( NNP), NC2( 2,NZ2), KCC( NNP), NDE( NEL,4),
4 MTR( NEL), X ( NNP), Y ( NNP), IBV( NNP),
5 HXX( NE2), HYY( NE2), FND( NNP), FEL( NEL),
6 ARW( NNP)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, PRN, RE
COMMON /FLW/ NVC, IVC( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 NTC, ITC( 50), TPC( 50), COF(100)
C
C ----- To construct F & solve equation ---
C
CALL FVOR ( VOR, NDE, FVC, X, Y )
C
CALL MOVR ( PSO, 1.D0, PSI, NNP)
C
CALL DSLV ( FVC, PSI, DP, NC2, KCC, IBP, NNZ, NZ2, OMG(1),
1 ERR(1) )
C
CALL VELC ( PSI, U, V, NDE, UN, VN, X, Y, IBV, NCN, ICN, ARW )
C
C ----- To compute RKX & RKY ---
C
CALL CKAK ( RKX, RKY, U, V, MTR, NE2, HXX, HYY, UN, VN,
1 NDE, FND, FEL )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE LPTP ( UN, VN, TEP, TPO, IBT, D, P2, QWK, FVC,
1 QE, NC2, KCC, NDE, X, Y, WNU, WNV, NNZ,
2 NZ2 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 KCC, IBT, NDE
C
DIMENSION UN( NNP), VN( NNP), TEP( NNP), TPO( NNP),
1 IBT( NNP), D ( NNZ), P2 ( NNZ), QWK( NEL,4),
2 FVC( NNP), QE( NEL), NC2( 2,NZ2), KCC( NNP),
3 NDE( NEL,4), X ( NNP), Y ( NNP), DS ( 121)
C
COMMON /NNL/ NNP,NEL,LUV
COMMON /TIM/ DT, LPM, LOP, TSP
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, PRN, RE
C
C ----- To construct F ---
C
DO 100 I = 1, NEL
C
QWK( I,1) = 0.D0
QWK( I,2) = 0.D0
QWK( I,3) = 0.D0
QWK( I,4) = 0.D0
C
100 CONTINUE
C
CALL MQV12 ( TEP, FVC, UN, VN, NDE, X, Y, QWK, DMX, DS,
1 WNU, WNV )
C
CALL FVCTN ( D, TPO, FVC, NC2, KCC, NNZ, NZ2 )
C
C ----- To solve time step equation ---
C
IF ( LOP.GT.1) CALL MOVR ( TPO, 1.D0 , TEP, NNP )
C
CALL DSLV ( FVC, TEP, P2, NC2, KCC, IBT, NNZ, NZ2,
1 OMG(3), ERR(3) )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE LPVR ( UN, VN, VOR, VRO, IBV, D, P2, FVC,
1 NC2, KCC, NDE, TEP, QWK, X, Y, NNZ,
2 NZ2, ICN, DMX, DN, WNU, WNV )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 KCC, IBV, NDE
C
DIMENSION UN ( NNP), VN ( NNP), VOR( NNP), VRO( NNP),
1 IBV( NNP), D ( NNZ), P2 ( NNZ), FVC( NNP),
2 NC2( 2,NZ2), KCC( NNP), NDE( NEL,4 ), TEP( NNP),
3 QWK( NEL,4), X ( NNP), Y ( NNP), DS ( 121),
4 WNU( NEL,16), WNV( NEL,16)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /TIM/ DT, LPM, LOP, TSP
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, PRN, RE
COMMON /FLW/ NVC, IVC( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 NTC, ITC( 50), TPC( 50), COF( 100)
C
C ---- To construct F ---
C
IF ( ICN.EQ.0) GO TO 90
C
DO 100 I = 1, NEL
C
TE1 = TEP( NDE( I,1))
TE2 = TEP( NDE( I,2))
TE3 = TEP( NDE( I,3))
TE4 = TEP( NDE( I,4))
C
C ----- Attention --- ( Coefficient of "Pr x Ra" ) ---
C
QWK( I,1) = RA* COF( 1)* COF( 1)* ( - 4.D0* TE1/3.D0
C
1 + 4.D0* TE2/3.D0 + 2.D0* TE3/3.D0 - 2.D0* TE4/3.D0)
2 / ( 2.D0* DABS( X( NDE( I,3)) - X( NDE( I,1))))
C
QWK( I,2) = RA* COF( 1)* COF( 1)* ( - 4.D0* TE1/3.D0
C
1 + 4.D0* TE2/3.D0 + 2.D0* TE3/3.D0 - 2.D0* TE4/3.D0)
2 / ( 2.D0* DABS( X( NDE( I,3)) - X( NDE( I,1))))
C
QWK( I,3 ) = RA* COF( 1)* COF( 1)* ( - 2.D0* TE1/3.D0
C
1 + 2.D0* TE2/3.D0 + 4.D0* TE3/3.D0 - 4.D0* TE4/3.D0)
2 / ( 2.D0* DABS( X( NDE( I,3)) - X( NDE( I,1))))
C
QWK( I,4 ) = RA* COF( 1)* COF( 1)* ( - 2.D0* TE1/3.D0
C
1 + 2.D0* TE2/3.D0 + 4.D0* TE3/3.D0 - 4.D0* TE4/3.D0)
2 / ( 2.D0* DABS( X( NDE( I,3 )) - X( NDE( I,1))))
C
100 CONTINUE
C
90 CONTINUE
C
CALL MQV12 ( VOR, FVC, UN, VN, NDE, X, Y, QWK, DMX, DS,
1 WNU, WNV )
C
DN = 0.D0
C
DO 150 I = 1, NNP
C
DN = DN + DABS( DS( I))
150 CONTINUE
C
DN = DN/ NNP
C
CALL FVCTN ( D, VRO, FVC, NC2, KCC, NNZ, NZ2 )
C
C ----- To solve simultaneous equation ---
C
DO 200 I = 1, NNP
C
IF ( IBV( I).EQ.9999) GO TO 200
VOR( I) = VRO( I)
C
200 CONTINUE
C
CALL DSLV ( FVC, VOR, P2, NC2, KCC, IBV, NNZ, NZ2, OMG( 2),
1 ERR( 2))
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MODD ( D, P1, P2, X, Y, RKX, RKY, NCW, NDE,
1 NNZ, NDM )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE
C
DIMENSION D( NNZ), P1 ( NNZ), P2 ( NNZ), X ( NNP), Y( NNP),
1 RKX( NNP), RKY( NNP), NCW( NNZ), NDE( NEL,4)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /TIM/ DT, LPM, LOP, TSP
C
C ----- To modify D ---
C
CALL DMAT ( X, Y, RKX, RKY, D, NDE, NCW, NNZ, NDM )
C
W1 = 1.D0/ DT
C
DO 100 I = 1, NNZ
C
P2( I) = TSP* D( I) + P1( I)* W1
D ( I) = P2( I) - D( I)
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MOVI ( IX, IA, IY, N )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION IX( N), IY( N)
C
DO 100 I = 1, N
C
IY ( I) = IX( I)* IA
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MOVR ( X, A, Y, N )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION X( N), Y( N)
C
DO 100 I = 1, N
C
Y ( I) = X( I)* A
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MQV12 ( C, F, UN, VN, NDE, X, Y, QWK, DMX, DS,
1 WNU, WNV )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE( NEL,4)
C
DIMENSION C( NNP), F( NNP), UN( NNP), VN( NNP),
1 X( NNP), Y( NNP), QWK( NEL,4 ), DS( NNP),
2 QD( 4), WNU( NEL,16), WNV( NEL,16)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /TIM/ DT, LPM, LOP, TSP
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, PRN, RE
C
CALL CLAR8 ( F, 0.D0, NNP)
C
CALL CLAR8 ( DS, 0.D0, NNP)
C
C ----- PSI & OMG : LINEAR ---
C
C --------------------------------------------------------------------
DO 100 NE = 1, NEL
C
C1 = C( NDE( NE,1))
C2 = C( NDE( NE,2))
C3 = C( NDE( NE,3))
C4 = C( NDE( NE,4))
C
A = DABS( Y( NDE( NE,1)) - Y( NDE( NE,3)))/ 2.D0
C
WNU( NE, 1) = - A/ 6.D0*C1 + A/ 6.D0*C2 + A/18.D0*C3 - A/18.D0*C4
WNU( NE, 2) = - A/12.D0*C1 + A/12.D0*C2 + A/36.D0*C3 - A/36.D0*C4
WNU( NE, 3) = - A/36.D0*C1 + A/36.D0*C2 + A/36.D0*C3 - A/36.D0*C4
WNU( NE, 4) = - A/18.D0*C1 + A/18.D0*C2 + A/18.D0*C3 - A/18.D0*C4
C
WNU( NE, 5) = WNU( NE, 2)
WNU( NE, 6) = - A/ 6.D0*C1 + A/ 6.D0*C2 + A/18.D0*C3 - A/18.D0*C4
WNU( NE, 7) = - A/18.D0*C1 + A/18.D0*C2 + A/18.D0*C3 - A/18.D0*C4
WNU( NE, 8) = - A/36.D0*C1 + A/36.D0*C2 + A/36.D0*C3 - A/36.D0*C4
C
WNU( NE, 9) = WNU( NE ,3)
WNU( NE,10) = WNU( NE ,7)
WNU( NE,11) = - A/18.D0*C1 + A/18.D0*C2 + A/ 6.D0*C3 - A/ 6.D0*C4
WNU( NE,12) = - A/36.D0*C1 + A/36.D0*C2 + A/12.D0*C3 - A/12.D0*C4
C
WNU( NE,13) = WNU( NE, 4)
WNU( NE,14) = WNU( NE, 8)
WNU( NE,15) = WNU( NE,12)
C
WNU( NE,16) = - A/18.D0*C1 + A/18.D0*C2 + A/ 6.D0*C3 - A/ 6.D0*C4
C
100 CONTINUE
C --------------------------------------------------------------------
DO 200 NE = 1, NEL
C
C1 = C( NDE( NE,1))
C2 = C( NDE( NE,2))
C3 = C( NDE( NE,3))
C4 = C( NDE( NE,4))
C
B = DABS( X( NDE( NE,1)) - X( NDE( NE,3)))/ 2.D0
C
WNV( NE, 1) = - B/ 6.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/ 6.D0*C4
WNV( NE, 2) = - B/18.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/18.D0*C4
WNV( NE, 3) = - B/36.D0*C1 - B/36.D0*C2 + B/36.D0*C3 + B/36.D0*C4
WNV( NE, 4) = - B/12.D0*C1 - B/36.D0*C2 + B/36.D0*C3 + B/12.D0*C4
C
WNV( NE, 5) = WNV( NE,2)
WNV( NE, 6) = - B/18.D0*C1 - B/ 6.D0*C2 + B/ 6.D0*C3 + B/18.D0*C4
WNV( NE, 7) = - B/36.D0*C1 - B/12.D0*C2 + B/12.D0*C3 + B/36.D0*C4
WNV( NE, 8) = - B/36.D0*C1 - B/36.D0*C2 + B/36.D0*C3 + B/36.D0*C4
C
WNV( NE, 9) = WNV( NE,3)
WNV( NE,10) = WNV( NE,7)
C
WNV( NE,11) = - B/18.D0*C1 - B/ 6.D0*C2 + B/ 6.D0*C3 + B/18.D0*C4
WNV( NE,12) = - B/18.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/18.D0*C4
C
WNV( NE,13) = WNV( NE,4)
WNV( NE,14) = WNV( NE,8)
WNV( NE,15) = WNV( NE,12)
C
WNV( NE,16) = - B/ 6.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/ 6.D0*C4
C
200 CONTINUE
C --------------------------------------------------------------------
DO 300 IEL = 1, NEL
C
C ----------------------------------------------------------
DO 320 I = 1, 4
C
II = NDE( IEL,I)
FU = 0.D0
FV = 0.D0
C
C ------------------------------------------------
DO 350 J = 1, 4
C
JJ = NDE( IEL,J)
C
FU = FU + WNU( IEL,4* ( I-1)+J)* UN( JJ)
FV = FV + WNV( IEL,4* ( I-1)+J)* VN( JJ)
C
350 CONTINUE
C ------------------------------------------------
C
F( II) = F( II) + FU + FV
C
320 CONTINUE
C ----------------------------------------------------------
C
300 CONTINUE
C --------------------------------------------------------------------
DO 400 NE = 1, NEL
C
C1 = UN ( NDE( NE,1))
C2 = UN ( NDE( NE,2))
C3 = UN ( NDE( NE,3))
C4 = UN ( NDE( NE,4))
C
A = DABS( Y( NDE( NE,1)) - Y( NDE( NE,3)))/ 2.D0
C
WNU( NE, 1) = - A/ 6.D0*C1 + A/ 6.D0*C2 + A/18.D0*C3 - A/18.D0*C4
WNU( NE, 2) = - A/12.D0*C1 + A/12.D0*C2 + A/36.D0*C3 - A/36.D0*C4
WNU( NE, 3) = - A/36.D0*C1 + A/36.D0*C2 + A/36.D0*C3 - A/36.D0*C4
WNU( NE, 4) = - A/18.D0*C1 + A/18.D0*C2 + A/18.D0*C3 - A/18.D0*C4
C
WNU( NE, 5) = WNU( NE,2)
WNU( NE, 6) = - A/ 6.D0*C1 + A/ 6.D0*C2 + A/18.D0*C3 - A/18.D0*C4
WNU( NE, 7) = - A/18.D0*C1 + A/18.D0*C2 + A/18.D0*C3 - A/18.D0*C4
WNU( NE, 8) = - A/36.D0*C1 + A/36.D0*C2 + A/36.D0*C3 - A/36.D0*C4
C
WNU( NE, 9) = WNU( NE,3)
WNU( NE,10) = WNU( NE,7)
WNU( NE,11) = - A/18.D0*C1 + A/18.D0*C2 + A/ 6.D0*C3 - A/ 6.D0*C4
WNU( NE,12) = - A/36.D0*C1 + A/36.D0*C2 + A/12.D0*C3 - A/12.D0*C4
C
WNU( NE,13) = WNU( NE,4)
WNU( NE,14) = WNU( NE,8)
WNU( NE,15) = WNU( NE,12)
WNU( NE,16) = - A/18.D0*C1 + A/18.D0*C2 + A/ 6.D0*C3 - A/ 6.D0*C4
C
400 CONTINUE
C --------------------------------------------------------------------
DO 500 NE = 1, NEL
C
C1 = VN ( NDE( NE,1))
C2 = VN ( NDE( NE,2))
C3 = VN ( NDE( NE,3))
C4 = VN ( NDE( NE,4))
C
B = DABS( X( NDE( NE,1)) - X( NDE( NE,3)))/ 2.D0
C
WNV( NE, 1) = - B/ 6.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/ 6.D0*C4
WNV( NE, 2) = - B/18.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/18.D0*C4
WNV( NE, 3) = - B/36.D0*C1 - B/36.D0*C2 + B/36.D0*C3 + B/36.D0*C4
WNV( NE, 4) = - B/12.D0*C1 - B/36.D0*C2 + B/36.D0*C3 + B/12.D0*C4
C
WNV( NE, 5) = WNV( NE,2)
WNV( NE, 6) = - B/18.D0*C1 - B/ 6.D0*C2 + B/ 6.D0*C3 + B/18.D0*C4
WNV( NE, 7) = - B/36.D0*C1 - B/12.D0*C2 + B/12.D0*C3 + B/36.D0*C4
WNV( NE, 8) = - B/36.D0*C1 - B/36.D0*C2 + B/36.D0*C3 + B/36.D0*C4
C
WNV( NE, 9) = WNV( NE,3)
WNV( NE,10) = WNV( NE,7)
WNV( NE,11) = - B/18.D0*C1 - B/ 6.D0*C2 + B/ 6.D0*C3 + B/18.D0*C4
WNV( NE,12) = - B/18.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/18.D0*C4
C
WNV( NE,13) = WNV( NE,4)
WNV( NE,14) = WNV( NE,8)
WNV( NE,15) = WNV( NE,12)
WNV( NE,16) = - B/ 6.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/ 6.D0*C4
C
500 CONTINUE
C --------------------------------------------------------------------
DO 600 IEL = 1, NEL
C
C ----------------------------------------------------------
DO 620 I = 1, 4
C
II = NDE( IEL,I)
FU = 0.D0
FV = 0.D0
DUV = 0.D0
C
C ------------------------------------------------
DO 650 J = 1, 4
C
JJ = NDE( IEL, J)
FU = FU + WNU( IEL,4*( I-1)+J)* C( JJ)
FV = FV + WNV( IEL,4*( I-1)+J)* C( JJ)
DUV = DUV + WNU( IEL,4*( I-1)+J) + WNV( IEL,4*( I-1)+J)
C
650 CONTINUE
C ------------------------------------------------
C
F ( II) = F ( II) + FU + FV
DS( II) = DS( II) + DUV
C
620 CONTINUE
C ----------------------------------------------------------
C
600 CONTINUE
C --------------------------------------------------------------------
DO 700 I = 1, NEL
C
A = DABS( Y( NDE( I,1)) - Y( NDE( I,3)))
B = DABS( X( NDE( I,1)) - X( NDE( I,3)))
C
N1 = NDE( I,1)
N2 = NDE( I,2)
N3 = NDE( I,3)
N4 = NDE( I,4)
C
QD( 1) = A* B* QWK( I,1)/ 4.D0
QD( 2) = A* B* QWK( I,2)/ 4.D0
QD( 3) = A* B* QWK( I,3)/ 4.D0
QD( 4) = A* B* QWK( I,4)/ 4.D0
C
F( N1) = F( N1) - QD( 1)
F( N2) = F( N2) - QD( 2)
F( N3) = F( N3) - QD( 3)
F( N4) = F( N4) - QD( 4)
C
700 CONTINUE
C --------------------------------------------------------------------
DMX = 0.D0
DO 800 I = 1, NNP
C
DMX = DMAX1( DMX, DABS( DS( I)))
800 CONTINUE
C --------------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE NBCN ( SH, DSH, PSI, VOR, LBC, NTB, NNP, LUV )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 LBC, NTB
C
DIMENSION SH ( LUV), DSH( LUV), PSI( NNP), VOR( NNP),
1 LBC( LUV), NTB( LUV,10)
C
C ----- Vorticity at wall ---
C
IF ( LUV.LE.0) RETURN
C
DO 100 I = 1, LUV
C
ND = LBC( I)
N1 = NTB( I,1)
C
APSI = PSI( N1)
AVOR = VOR( N1)
C
VOR( ND) = DSH( I)* ( PSI( ND) - APSI) - AVOR* 0.5D0
C
C ----- Attention : In case of 121/ 100 --------------------------
C
IF ( ND.GE.112.AND.ND.LE.120) VOR( ND) = VOR( ND) - SH( I)
C ----------------------------------------------------------------
C
100 CONTINUE
C
C ----- Attention : Vorticity = 0 at 4 corners ------
C
NEX = 10
NEY = 10
C
VOR( 1) = 0.D0
VOR( NEX+1) = 0.D0
C
VOR( 1+ NEY* 11) = 0.D0
VOR(( NEX+1)*( NEY+1)) = 0.D0
C
C ---------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLTD ( PSI, VOR, TEP, UN, VN, X, Y, NDE, ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /TIM/ DT, LPM, LOP, TSP
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, PRN, RE
COMMON /VPT/ VARV( 2000), VART( 2000), VRUV( 2000), VARP( 2000)
COMMON /FLW/ NVC, IVC( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 NTC, ITC( 50), TPC( 50), COF(100)
C
DIMENSION PSI( NNP), VOR( NNP), TEP( NNP), UN ( NNP),
1 VN ( NNP), X ( NNP), Y ( NNP), NDE( NEL,4)
C
IF ( ICN.EQ.1) RE = 1.D0 ! Dummy
C
WRITE(10,2000) ICN, LOP, DT, RE, RA, PRN, NNP, NEL
2000 FORMAT( I2, I4, F6.3, 3F9.2, 2I4)
WRITE(10,2100) ( X( I), Y( I), I = 1, NNP)
2100 FORMAT(10F8.4)
WRITE(10,2200) (( NDE( I,J), I = 1, NEL), J = 1, 4 )
2200 FORMAT(20I4)
WRITE(10,2300) ( VRUV( I), VARP( I), I = 1, LOP)
2300 FORMAT(10F11.5)
C
IF ( ICN.EQ.1) WRITE(10,2300) ( VARV( I), VART( I), I = 1, LOP)
IF ( ICN.EQ.0) WRITE(10,2300) ( VARV( I), I = 1, LOP)
C
WRITE(10,2100) ( UN( I), VN( I), I = 1, NNP)
IF ( ICN.EQ.1)
1WRITE(10,2100) ( PSI( I), VOR( I), TEP( I), I = 1, NNP)
IF ( ICN.EQ.0)
1WRITE(10,2100) ( PSI( I), VOR( I), I = 1, NNP)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PMAT ( P, NCC, NCW, X, Y, NDE, KCC, NC2, NNZ,
1 NZ2, ARS )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, KCC
C
DIMENSION P ( NNZ), NCC( NZ2), NCW( NNZ), X( NNP), Y( NNP),
1 NDE( NEL,4), KCC( NNP), NC2( 2,NZ2)
C
DIMENSION A1( 4,4), N( 4), ARS( NEL)
C
COMMON /NNL/ NNP, NEL, LUV
C
WD9 = 1.D0/ 9.D0
WD18 = 1.D0/ 18.D0
WD36 = 1.D0/ 36.D0
C
A1( 1,1) = WD9
A1( 1,2) = WD18
A1( 1,3) = WD36
A1( 1,4) = WD18
C
A1( 2,1) = WD18
A1( 2,2) = WD9
A1( 2,3) = WD18
A1( 2,4) = WD36
C
A1( 3,1) = WD36
A1( 3,2) = WD18
A1( 3,3) = WD9
A1( 3,4) = WD18
C
A1( 4,1) = WD18
A1( 4,2) = WD36
A1( 4,3) = WD18
A1( 4,4) = WD9
C
CALL CLAR8( P, 0.D0, NNZ )
C
NZK = 0
C
C ------------------------------------------------
DO 100 IEL = 1, NEL
C
N( 1) = NDE( IEL,1)
N( 2) = NDE( IEL,2)
N( 3) = NDE( IEL,3)
N( 4) = NDE( IEL,4)
C
ARS( IEL) = DABS( Y( NDE( IEL,1)) - Y( NDE( IEL,3)))
1 * DABS( X( NDE( IEL,1)) - X( NDE( IEL,3)))
C
C -------------------------------------
DO 150 I = 1, 4
DO 150 J = 1, 4
C -------------------------------------
C
II = N( I)
JJ = N( J)
C
A = DABS( Y( NDE( IEL,1)) - Y( NDE( IEL,3)))
B = DABS( X( NDE( IEL,1)) - X( NDE( IEL,3)))
AV = A* B* A1( I,J)
C
CALL DCNL ( AV, II, JJ, P, NCC, NZK, NNZ)
C
150 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
DO 200 II = 1, NZ2
C
NC2( 1,II) = NCC( II)/ 10300
NC2( 2,II) = NCC( II) - 10300* NC2( 1,II)
C
200 CONTINUE
C ----------------------------------------------------------
C
CALL MOVI ( NCC, 1, NCW, NNZ )
C
C ----------------------------------------------------------
C
CALL DCIJ ( KCC, NC2, NZK, NNP, NZ2 )
C
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SUBS ( AV, M, N, A, L, K, NT )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION A( NT), L( NT)
C
K = K + 1
IF ( K.GT.NT ) WRITE(6,2000)
2000 FORMAT(' * Error K.GT.NT *')
JE = K - N
C
DO 100 JJ = 1, JE
C
J1 = K - JJ + 1
A( J1) = A( J1-1)
L( J1) = L( J1-1)
C
100 CONTINUE
C
A( N) = AV
L( N) = M
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /NNL/ NNP, NEL, LUV
C
WRITE(6,*) ' '
WRITE(6,*) 'C ************************************************* C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C H F 2 D - F E R C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( Version 4.7 ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C 2D Thermal Fluid Analysis by FEM C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( TRectangular Linear 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
READ(7,1000) NNP, NEL, LUV
1000 FORMAT(10X,I5,5X,I5,5X,I5)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VELC ( PSI, U, V, NDE, UN, VN, X, Y, IBV,
1 NCN, ICN, ARW )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, NCN, IBV
C
DIMENSION PSI( NNP), U ( NEL), V ( NEL), NDE( NEL,4),
1 UN ( NNP), VN ( NNP), X ( NNP), Y ( NNP),
2 NCN( NNP), IBV( NNP), ARW( NNP)
C
COMMON /NNL/ NNP, NEL, LUV
C
CALL CLAR8 ( UN, 0.D0, NNP)
C
CALL CLAR8 ( VN, 0.D0, NNP)
C
CALL CLAR8 ( ARW, 0.D0, NNP)
C
C ----- To set of S, T & calculaTE for UN, VN ---
C
C ----------------------------------------------------------
DO 100 NE = 1, NEL
C
A = DABS( Y( NDE( NE,1)) - Y( NDE( NE,3)))
B = DABS( X( NDE( NE,1)) - X( NDE( NE,3)))
AB = A* B
C
N1 = NDE( NE,1)
N2 = NDE( NE,2)
N3 = NDE( NE,3)
N4 = NDE( NE,4)
C
ARW( N1) = ARW( N1) + AB
ARW( N2) = ARW( N2) + AB
ARW( N3) = ARW( N3) + AB
ARW( N4) = ARW( N4) + AB
C
SN1 = 0.D0
SN2 = X( N2) - X( N1)
SN3 = X( N2) - X( N1)
SN4 = 0.D0
C
TN1 = 0.D0
TN2 = 0.D0
TN3 = Y( N3) - Y( N2)
TN4 = Y( N3) - Y( N2)
C
UN( N1) = UN( N1) + (( SN1-B)* PSI( N1) - SN1* PSI( N2)
1 + SN1* PSI( N3) + ( B-SN1)* PSI( N4))
C
UN( N2) = UN( N2) + (( SN2-B)* PSI( N1) - SN2* PSI( N2)
1 + SN2* PSI( N3) + ( B-SN2)* PSI( N4))
C
UN( N3) = UN( N3) + (( SN3-B)* PSI( N1) - SN3* PSI( N2)
1 + SN3* PSI( N3) + ( B-SN3 )* PSI( N4))
C
UN( N4) = UN( N4) + (( SN4-B)* PSI( N1) - SN4* PSI( N2)
1 + SN4* PSI( N3) + ( B-SN4 )* PSI( N4))
C
C ---------------------------------------------------------------------
C
VN( N1) = VN( N1) - (( TN1-A)* PSI( N1) + ( A-TN1)* PSI( N2)
1 + TN1* PSI( N3) - TN1* PSI( N4))
C
VN( N2) = VN( N2) - (( TN2-A)* PSI( N1) + ( A-TN2)* PSI( N2)
1 + TN2* PSI( N3) - TN2* PSI( N4))
C
VN( N3) = VN( N3) - (( TN3-A)* PSI( N1) + ( A-TN3)* PSI( N2)
1 + TN3* PSI( N3) - TN3* PSI( N4))
C
VN( N4) = VN( N4) - (( TN4-A)* PSI( N1) + ( A-TN4)* PSI( N2)
1 + TN4* PSI( N3) - TN4* PSI( N4))
C
100 CONTINUE
C ----------------------------------------------------------
DO 200 I = 1, NNP
C
UN( I) = UN( I)/ ARW( I)
VN( I) = VN( I)/ ARW( I)
C
200 CONTINUE
C ----------------------------------------------------------
DO 300 NE = 1, NEL
C
N1 = NDE( NE,1)
N2 = NDE( NE,2)
N3 = NDE( NE,3)
N4 = NDE( NE,4)
C
A = DABS( Y( NDE( NE,1)) - Y( NDE( NE,3)))
B = DABS( X( NDE( NE,1)) - X( NDE( NE,3)))
C
S = B* 0.5D0
T = A* 0.5D0
AB = A* B
C
U( NE) = (( S - B )* PSI( N1) - S* PSI( N2)
1 + S* PSI( N3) + ( B-S)* PSI( N4))/ AB
C
V( NE) = - (( T - A )* PSI( N1) + ( A-T)* PSI( N2)
1 + T* PSI( N3) - T* PSI( N4))/ AB
C
300 CONTINUE
C ----------------------------------------------------------
DO 400 NE = 1, NEL
DO 400 II = 1, 4
C ----------------------------------------------------------
C
JJ = NDE( NE,II)
C
IF ( IBV( JJ).EQ.9999) UN( JJ) = 0.D0
IF ( IBV( JJ).EQ.9999) VN( JJ) = 0.D0
C
400 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************