戻る

PROGRAM MAIN
C
C ******************************************************************** C
C C
C CVDF3D - FEM C
C C
C Three-Dimensional Convection-Diffusion Analysis C
C C
C by Finite Element Method C
C C
C ( Rectangular Parallelepiped Elements ) C
C C
C ( V.3.8 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ******************************************************************** C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 KCC, NCC, IBC
C
C ----- NNP ---
DIMENSION X ( 2541), Y ( 2541), Z ( 2541), U ( 2541),
1 V ( 2541), W ( 2541), C ( 2541), CC ( 2541),
2 ICC( 2541), F ( 2541), B ( 2541), ITB( 2541),
3 CAA( 2541), KCC( 2541), IBC( 2541), CA ( 2541),
4 CN ( 2541)
C
C ----- NEL ---
DIMENSION UE ( 2000), VE ( 2000), WE ( 2000), RKX( 2000),
1 RKY( 2000), RKZ( 2000), DLX( 2000), DLY( 2000),
2 DLZ( 2000), AV ( 2000), VS ( 2000), COF( 2000)
C
C ----- NEL,8 ---
C
DIMENSION NDE( 2000,8)
C
C ----- NNZ ---
C
DIMENSION P( 50000), D( 50000), A( 50000), NCW( 50000)
C
C ----- NZ2* 2 ---
C
DIMENSION NCC( 2,100000)
C
COMMON /CON2/ COEF, NPB, PHI( 20), SMH, DTD
C
OPEN ( 9, FILE='CVDF3D_FEM.DAT' ) ! DT = 0.05 LOP = 20
OPEN ( 10, FILE='CVDF3D_FEM_FIG.DATA' )
C
C --------------------------------------------------------------------
C
CALL TTLE
C
900 CONTINUE
C
CALL INPT ( X, Y, Z, NDE, ITB, AV, DLX, DLY, DLZ,
1 NCL, CC, ICC, DT, LPE, NNP, NEL, NNZ, NZ2 )
C
CALL INIT ( U, V, W, C, UE, VE, WE, VS, NCL,
1 CC, ICC, IBC, RKX, RKY, RKZ, NDE, COF, DT,
2 NNP, NEL )
C
CALL PMAT ( P, AV, NDE, NCW, KCC, NCC, NNP, NEL, NNZ,
1 NZ2 )
C
CALL DMAT ( D, AV, DLX, DLY, DLZ, COF, NDE, RKX, RKY,
1 RKZ, NCW, NEL, NNZ )
C
C ----------------------------------------------------------
DO 100 I = 1, NPB
C
PHA = PHI( I)
LOP = 0
SEM = 0.D0
C
CALL TOUT ( CAA, C, ITB, LOP, CA, CN, DT, LPE, SEM, ERR,
1 NNP, *999 )
C
999 LOP = LOP + 1
C
CALL FVCT ( F, C, UE, VE, WE, AV, DLX, DLY, DLZ, VS,
1 NDE, NNP, NEL )
C
CALL ABVC ( A, B, P, D, F, C, PHA, NCC, KCC, DT,
1 NNP, NNZ, NZ2 )
C
CALL SOLV ( NNP, NNZ, NZ2, B, C, A, NCC, KCC, IBC, LOP)

C
CALL TOUT ( C, CAA, ITB, LOP, CA, CN, DT, LPE, SEM, ERR,
1 NNP, *999 )
C
100 CONTINUE
C ----------------------------------------------------------
C
CALL FNSB ( SEM, ERR, *900 )
C
STOP
END
C **********************************************************************
C
SUBROUTINE ABVC ( A, B, P, D, F, C, PHI, NCC, KCC,
1 DT, NNP, NNZ, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NCC, KCC
C
DIMENSION B( NNP), F( NNP), C ( NNP), KCC( NNP), P( NNZ),
1 D( NNZ), A( NNZ), NCC( 2,NZ2)
C
C --------------------------------------
C
CALL CLER ( A, 0.D0, NNZ )
C
C --------------------------------------
C
CALL CLER ( B, 0.D0, NNP )
C
C --------------------------------------
C
W1 = 1.D0/ DT
W2 = 1.D0 - PHI
C
C ----- { B } = ( [ P ] / DT - (1 - PHI)* [ D ] ){C} ---
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
JS = 1
IF ( I.GT.1) JS = KCC( I-1) + 1
JL = KCC( I)
FW = 0.D0
C
C --------------------------------------
DO 150 J = JS, JL
C
KPR = NCC( 1,J)
KCL = NCC( 2,J)
BBB = P( KPR)* W1 - D( KPR)* W2
FW = FW + BBB* C( KCL)
C
150 CONTINUE
C --------------------------------------
B( I) = FW
C
100 CONTINUE
C ------------------------------------------------
C
C ----- [ A ] = [ P ]/DT + PHI* [ D ] ---
C
DO 200 I = 1, NNZ
C
A( I) = D( I)* PHI + P( I)* W1
200 CONTINUE
C
C ----- { B } = { B } - {F} -------------
C
DO 300 I = 1, NNP
C
B( I) = B( I) - F( I)
300 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE AGCL ( NNP, NEL, X, Y, Z, NDE, AV, DLX, DLY,
1 DLZ, BV, RV, DT )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CON2/ COEF, NPB, PHI( 20), SMH, DTD
COMMON /COND/ UC, VC, WC, UIT, VIT, WIT, CIT
C
DIMENSION X ( NNP), Y ( NNP), Z ( NNP), AV ( NEL),
1 DLX( NEL), DLY( NEL), DLZ( NEL), NDE( NEL,8)
C
C --------------------------------------
DO 100 I = 1, NEL
C
N1 = NDE( I,1)
N2 = NDE( I,2)
N4 = NDE( I,4)
N5 = NDE( I,5)
C
DLX( I) = DABS( X( N2) - X( N1)) ! = (2a)
DLY( I) = DABS( Y( N4) - Y( N1)) ! = (2b)
DLZ( I) = DABS( Z( N5) - Z( N1)) ! = (2c)
C
C ----- Attention ----------------------------------------------
C
AV( I) = DLX( I)* DLY( I)* DLZ( I) ! = (2a)* (2b)* (2c)
C
C --------------------------------------------------------------
IF ( AV( I).LE.0.D0) THEN
WRITE(6,*) ' ** DATA: Strange ** STOP ** '
STOP
END IF
C
100 CONTINUE
C --------------------------------------
C
DLM = 0.D0
C --------------------------------------
DO 200 I = 1, NEL
C
IF ( DLX( I).GT.DLM) DLM = DLX( I)
200 CONTINUE
C --------------------------------------
C
BV = UC* DT/ DLM
RV = COEF* DT/(DLM* DLM)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ANAS ( NNP, C, LOP, DT )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CON2/ COEF, NPB, PHI( 20), SMH, DTD
COMMON /COND/ UC, VC, WC, UIT, VIT, WIT, CIT
COMMON /XYZ / XP( 1500), YP( 1500), ZP( 1500)
C
DIMENSION C( NNP)
C
C ----- Analytical solution ---
C
RKX = COEF
RKY = COEF
RKZ = COEF
C
RKT = DTD
C ===============
C
DX = SMH
DY = SMH
DZ = SMH
C
T = FLOAT( LOP)* DT
RR = RKX + RKY + RKZ
C
C ----- Attention ---
C
IF ( RR.GT.0.00001) RKT = ( RKT/( RKX + RKY + RKZ))* 3.D0
C
C --------------------------------------------------------------------
DO 100 I = 1, NNP
C
CF1 = ( DX/2.D0 - ( XP(I)+0.5D0-UC*T))/( 2.D0*DSQRT( RKX*T + RKT))
C =========
CF1 = ERF( CF1)
C
CF2 = ( DX/2.D0+( XP(I)+0.5D0 - UC*T))/( 2.D0*DSQRT( RKX*T + RKT))
CF2 = ERF( CF2)
C
CF3 = ( DY/2.D0 - ( YP(I) - VC*T))/( 2.D0*DSQRT( RKY*T + RKT ))
CF3 = ERF( CF3)
C
CF4 = ( DY/2.D0 + ( YP(I) - VC*T))/( 2.D0*DSQRT( RKY*T + RKT ))
CF4 = ERF( CF4)
C
CF5 = ( DZ/2.D0 - ( ZP(I) - WC*T))/( 2.D0*DSQRT( RKZ*T + RKT ))
CF5 = ERF( CF5)
C
CF6 = ( DZ/ 2.D0 + ( ZP(I) - WC*T))/( 2.D0* DSQRT( RKZ*T + RKT))
CF6 = ERF( CF6 )
C
C( I)= 0.125D0* ( CF1+ CF2)* ( CF3+ CF4)* ( CF5+ CF6 )* 100.D0
C
100 CONTINUE
C --------------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ARST ( NNP, NEL, NDE, NNZ, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION NDE( NEL,8), INDE2( 35000), M( 28), KZ3( 2,28)
C
DATA KZ3 / 1,2, 1,3, 1,4, 1,5, 1,6, 1,7, 1,8, 2,3, 2,4, 2,5, 2,6,
1 2,7, 2,8, 3,4, 3,5, 3,6, 3,7, 3,8, 4,5, 4,6, 4,7, 4,8,
2 5,6, 5,7, 5,8, 6,7, 6,8, 7,8 /
C
L = 0
C ------------------------------------------------
DO 100 I = 1, NEL
C
C --------------------------------------
DO 120 J = 1, 28
C
N1 = NDE( I, KZ3( 1,J))
N2 = NDE( I, KZ3( 2,J))
M1 = MIN0( N1, N2)
M2 = MAX0( N1, N2)
M( J) = M1* 10000 + M2
C
120 CONTINUE
C --------------------------------------
DO 140 J = 1, 28
C
IL = 1
20 IF ( IL.GT.L) GO TO 10
IF ( M( J).EQ.INDE2( IL)) GO TO 140
IL = IL + 1
GO TO 20
10 L = L + 1
INDE2( L) = M( J)
C
140 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
NNZ = NNP + L
NZ2 = NNP + L* 2
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CFDF ( NEL, RKX, RKY, RKZ, UE, VE, WE, COF, DT )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /RTIM/ KEY
COMMON /CON2/ COEF, NPB, PHI( 20), SMH, DTD
C
DIMENSION UE ( NEL), VE ( NEL), WE ( NEL), COF( NEL),
1 RKX( NEL), RKY( NEL), RKZ( NEL)
C
DT6 = DT/ 6.D0
C --------------------------------------------------------------------
DO 100 I = 1, NEL
C
RKX( I) = COEF
RKY( I) = COEF
RKZ( I) = COEF
C
IF( KEY.EQ.1) GO TO 150
C
C ----- With Corrected Diffusivity -------------------------
C
RKX( I) = COEF + ( UE( I) + VE( I) + WE( I))
1 *( UE( I) + VE( I) + WE( I))* DT6
C
RKY( I) = COEF + ( UE( I) + VE( I) + WE( I))
1 *( UE( I) + VE( I) + WE( I))* DT6
C
RKZ( I) = COEF + ( UE( I) + VE( I) + WE( I))
1 *( UE( I) + VE( I) + WE( I))* DT6
C ----------------------------------------------------------
150 CONTINUE
C
COF( I) = RKX( I) + RKY( I) + RKZ( I)
C
100 CONTINUE
C --------------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CLER ( RIX, RIV, N)
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION RIX( N)
C
DO 100 I = 1, N
C
RIX( I) = RIV
100 CONTINUE
C
RETURN
END
C **********************************************************************

SUBROUTINE CLR2 ( IX2, IV2, N)
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION IX2( N)
C
DO 100 I = 1, N
C
IX2 ( I) = IV2
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CNAL ( PDV, I, J, PD, NCC, K, NNZ )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 M2( 2)
EQUIVALENCE ( M, M2( 1))
C
DIMENSION PD( NNZ), NCC( NNZ)
C
C ----- Attention --- No. of Nodes.LE.32768 ---
C
IF ( I.GT.J) GO TO 99
C
M2( 1) = J
M2( 2) = I
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 ( M - NCC( N)) 21, 22, 23
21 CONTINUE
N2 = N
N = N2 - ( N2 - N1)/ 2
C
IF ( N.NE.N2) GO TO 100
IF ( M - NCC( N1)) 31, 32, 33
31 N= N1
33 CONTINUE
C
C ------------------------------------------------
C
CALL SUBS ( PDV, M, N, PD, NCC, K, NNZ )
C
C ------------------------------------------------
GO TO 99
C
32 N= N1
GO TO 22
C
42 N= N2
22 CONTINUE
PD( N) = PD( N) + PDV
GO TO 99
C
23 CONTINUE
N1 = N
N = N1 + ( N2 - N1)/ 2
IF ( N.NE.N1) GO TO 100
IF ( M - NCC( N2)) 41, 42, 43
41 N= N2
GO TO 33
C
43 K = K + 1
IF ( K.GT.NNZ) WRITE(6,2000)
2000 FORMAT(' ** Error ** ''NNZ'' IN ''SOR'' IS TOO SMALL **')
PD ( K) = PDV
NCC( K) = M
GO TO 99
C
100 CONTINUE
C ----------------------------------------------------------
C
WRITE(6,2100) K
2100 FORMAT(' * Error -- DO LOP strange, K =',I10 )
STOP
C
99 RETURN
END
C **********************************************************************
C
SUBROUTINE CNIJ ( KCC, NCC, NZ, NNP, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 KCC, NCC
C
DIMENSION KCC( NNP), NCC( 2, NZ2)
C
C NP, NZ : Current length of KCC, NCC
C
C NNP, NZ2 : Max. length of KCC, NCC
C
CHARACTER IANP * 4, IANZ2 * 4
DATA IANP /'NNP '/, IANZ2 /'NZ2'/
C
NP = 1
C --------------------------------------
DO 100 I = 1, NZ
C
IF ( NP.EQ.NCC( 1,I)) GO TO 10
C
C ----- KCC ( NP) ---
C
KCC( NP) = I - 1
NP = NCC( 1,I)
IF ( NP.GT.NNP) GO TO 30
10 NCC( 1,I) = I
C
100 CONTINUE
C --------------------------------------
C
C ---- Last row's counter ---
KCC( NP) = NZ
C
C ---- ( I,J) in NCC for I > J ---
C
C ------------------------------------------------
DO 200 I = 2, NP
C
I1 = I - 1
IS = KCC( I1)
C
C --------------------------------------
DO 220 J = 1, I1
C
JS = 1
IF ( J.GT.1) JS = KCC( J-1) + 1
JL = KCC( J)
C
C ----------------------------
DO 240 K = JS, JL
C
IF ( NCC( 2, K ).NE.I) GO TO 240
K1 = K
GO TO 20
C
240 CONTINUE
C ----------------------------
C
C ---- ( I,J) not exist ---
GO TO 220
C
20 CONTINUE
C ---- To shift NCC ---
C
NZIS = NZ - IS
C
C ----------------------------
DO 260 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
260 CONTINUE
C ----------------------------
C
NZ = NZ + 1
IF ( NZ.GT.NZ2 ) GO TO 40
C ------ To increase KCC ---
C
C ----------------------------
DO 280 L = I, NP
C
KCC( L) = KCC( L) + 1
280 CONTINUE
C ----------------------------
C ------- To insert ( I,J) element ---
C
IS = IS + 1
NCC( 1,IS) = NCC( 1,K1)
NCC( 2,IS) = J
C
220 CONTINUE
C --------------------------------------
C
200 CONTINUE
C ------------------------------------------------
RETURN
C
30 WRITE(6,2000) IANP, NNP, IANP, NP
STOP
40 WRITE(6,2000) IANZ2, NZ2, IANZ2, NZ
STOP
2000 FORMAT(' ** Error( CNIJ) ',A3,'(',I5,') too small',A2,'=',I5 )
C
END
C **********************************************************************
C
SUBROUTINE CPER ( NNP, C, CAA, LOP, SEM, DT )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION C( NNP), CAA( NNP)
C
T = FLOAT( LOP)* DT
SUM = 0.D0
C
C ----- Relative Error ---
C
C ------------------------------------------------
CM = 0.D0
DO 100 I = 1, NNP
C
IF ( C( I).GT.CM ) CM = C( I)
100 CONTINUE
C ------------------------------------------------
CAM = 0.D0
DO 200 I = 1, NNP
C
IF ( CAA( I).GT.CAM ) CAM = CAA( I)
200 CONTINUE
C ------------------------------------------------
DO 300 I = 1, NNP
C
IF ( CAM.LT.1.D-8 ) GO TO 300
SUM = SUM + ( DABS( C( I) - CAA( I))/ DABS( CAM))
C
300 CONTINUE
C ------------------------------------------------
C
SUM = ( SUM* 100.D0/ DFLOAT( NNP))
SEM = SEM + SUM
C ======================
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DMAT ( D, AV, DLX, DLY, DLZ, COF, NDE,
1 RKX, RKY, RKZ, NCW, NEL, NNZ )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION DLX( NEL), DLY( NEL), DLZ( NEL), AV ( NEL),
1 RKX( NEL), RKY( NEL), RKZ( NEL), COF( NEL),
2 NDE( NEL,8), D ( NNZ), NCW( NNZ), DX ( 4,4),
3 DY ( 4,4), DZ ( 4,4), DDX( 8,8), DDY( 8,8),
4 DDZ( 8,8)
C
CALL CLER ( D, 0.D0, NNZ)
C
NZWK = NNZ
C
CXX = 0.D0
C
C ------------------------------------------------
DO 100 I = 1, NEL
C
IF ( COF( I).GE.CXX) CXX = COF( I)
100 CONTINUE
C ------------------------------------------------
C
CXX = CXX* 0.001D0
C
C ------------------------------------------------
DO 200 I = 1, 4
DO 200 J = 1, 4
C
IJ = I + J
C
DX( I,J) = - 4.D0
DY( I,J) = 2.D0
DZ( I,J) = 2.D0
C
IF ( IJ.EQ.5) DX( I,J) = 2.D0
IF ( IJ.EQ.5) DY( I,J) = - 4.D0
IF ( IJ.EQ.5) DZ( I,J) = 2.D0
C
IF ( MOD( IJ,2).EQ.0) DX( I,J) = - 2.D0
IF ( MOD( IJ,2).EQ.0) DY( I,J) = - 2.D0
IF ( MOD( IJ,2).EQ.0) DZ( I,J) = 1.
C
IF ( I.EQ.J) DX( I,J) = 4.D0
IF ( I.EQ.J) DY( I,J) = 4.D0
IF ( I.EQ.J) DZ( I,J) = 4.D0
C
200 CONTINUE
C ----------------------------------------------------------
DO 300 IEL = 1, NEL
C
AKX = RKX( IEL)* AV( IEL)/ ( 2.D0* DLX( IEL)* DLX( IEL))
C
C ----- Attention ----------------------------------------------
C
C ==> RKX* 8* a* b* c/( 2* ( 2*a)* ( 2*a) ) ==> RKX* b* c/ a
C ===============
C --------------------------------------------------------------
C
AKY = RKY( IEL)* AV( IEL)/ ( 2.D0* DLY( IEL)* DLY( IEL))
AKZ = RKZ( IEL)* AV( IEL)/ ( 2.D0* DLZ( IEL)* DLZ( IEL))
C
C ------------------------------------------------
DO 350 I = 1, 8
DO 350 J = 1, 8
C ------------------------------------------------
C
II = NDE( IEL,I)
JJ = NDE( IEL,J)
C
IF ( I.LE.4 .AND. J.LE.4) DDX( I,J) = DX( I,J)
IF ( I.LE.4 .AND. J.LE.4) DDY( I,J) = DY( I,J)
IF ( I.LE.4 .AND. J.LE.4) DDZ( I,J) = DZ( I,J)
C
IF ( I.LE.4 .AND. J.GT.4) DDX( I,J) = 0.5D0* DX( I,J-4)
IF ( I.LE.4 .AND. J.GT.4) DDY( I,J) = 0.5D0* DY( I,J-4)
IF ( I.LE.4 .AND. J.GT.4) DDZ( I,J) = - DZ( I,J-4)
C
IF ( I.GT.4 .AND. J.LE.4) DDX( I,J) = 0.5D0* DX( I-4,J)
IF ( I.GT.4 .AND. J.LE.4) DDY( I,J) = 0.5D0* DY( I-4,J)
IF ( I.GT.4 .AND. J.LE.4) DDZ( I,J) = - DZ( I-4,J)
C
IF ( I.GT.4 .AND. J.GT.4) DDX( I,J) = DX( I-4,J-4)
IF ( I.GT.4 .AND. J.GT.4) DDY( I,J) = DY( I-4,J-4)
IF ( I.GT.4 .AND. J.GT.4) DDZ( I,J) = DZ( I-4,J-4)
C
AVQQ = DDX( I,J)* AKX + DDY( I,J)* AKY + DDZ( I,J)* AKZ
C
C ----------------------------------------------------
C
CALL CNAL ( AVQQ, II,JJ, D, NCW, NZWK, NNZ )
C
C ----------------------------------------------------
C
350 CONTINUE
C ------------------------------------------------
C
300 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DMST ( NNP, NEL, NLAY, NFIX )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION CID( 3), DLAY( 3), NLAY( 3)
C
CHARACTER CID*4, DLAY*4, CIDX*4, DIR*4
C
DATA CID /'LAYE','FIXC','FIXU'/
DATA DLAY /'X ','Y ','Z '/
C
NFIX = 0
C
REWIND 9
10 READ(9,1000,END=90) CIDX
1000 FORMAT(A4)
C
C --------------------------------------
DO 100 I = 1, 2
C
IF ( CIDX.EQ.CID( I)) GO TO 20
100 CONTINUE
C --------------------------------------
GO TO 10
C
20 BACKSPACE 9
GO TO ( 30, 40), I
C
30 CONTINUE
READ(9,1100) DIR
1100 FORMAT(10X,A1)
C
C --------------------------------------
DO 200 I = 1, 3
C
IF ( DIR.NE.DLAY( I)) GO TO 200
NLAY( I) = NLAY( I) + 1
GO TO 35
C
200 CONTINUE
C --------------------------------------
35 GO TO 10
40 CONTINUE
READ(9,1200) DIR
1200 FORMAT(10X,A4)
C
NFIX = NFIX + 1
GO TO 10
C
90 CONTINUE
NNP = NLAY( 1)* NLAY( 2)* NLAY( 3)
NEL = ( NLAY( 1)-1)* ( NLAY( 2)-1)* ( NLAY( 3)-1)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DTGN ( LYO, CRD, NNP, NEL, X, Y, Z, NDE, ITB,
1 IEB, SCL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- 3-Dimensional Node Data, Element Data Generation ---
C
DIMENSION X ( NNP), Y ( NNP), Z ( NNP), ITB( NNP),
1 IEB( NEL), NDE( NEL,8), LYO( 50,3), CRD( 50,3),
2 SCL ( 3), IWK( 8)
C
COMMON /MSZE/ NLAY( 3)
COMMON / XYZ/ XP( 1500), YP( 1500), ZP( 1500)
C
C --- NDE point DATA generation ---
C
J = 0
C
NLY1 = NLAY( 1)
NLY2 = NLAY( 2)
NLY3 = NLAY( 3)
C
C ----------------------------------------------------------
DO 100 I3 = 1, NLY3
DO 100 I2 = 1, NLY2
DO 100 I1 = 1, NLY1
C ----------------------------------------------------------
C
J = J + 1
NO = LYO( I3,3)* 10000 + LYO( I2,2)* 100 + LYO( I1,1)
ITB( J) = NO
C
X ( J) = CRD( I1,1)* SCL( 1)
Y ( J) = CRD( I2,2)* SCL( 2)
Z ( J) = CRD( I3,3)* SCL( 3)
C
XP( J) = CRD( I1,1)* SCL( 1)
YP( J) = CRD( I2,2)* SCL( 2)
ZP( J) = CRD( I3,3)* SCL( 3)
C
100 CONTINUE
C ----------------------------------------------------------
C
NNP = J
C
C ----- Node Data Generation ---
C
J = 0
NX1 = NLAY( 1) - 1
NY1 = NLAY( 2) - 1
NZ1 = NLAY( 3) - 1
C
J3 = NLAY( 1)* NLAY( 2)
J2 = NLAY( 1)
J1 = 1
C
C ----------------------------------------------------------
DO 200 K3 = 1, NZ1
DO 200 K2 = 1, NY1
DO 200 K1 = 1, NX1
C ----------------------------------------------------------
C
J = J + 1
IEB( J)= J
C
IWK( 1) = J3* ( K3-1) + J2* ( K2-1) + J1* ( K1-1) + 1
IWK( 2) = IWK( 1) + 1
IWK( 4) = IWK( 1) + J2
IWK( 3) = IWK( 4) + 1
IWK( 5) = IWK( 1) + J3
IWK( 6) = IWK( 5) + 1
IWK( 8) = IWK( 5) + J2
IWK( 7) = IWK( 8) + 1
C
C --------------------------------------
DO 250 I = 1, 8
C
NDE( J,I) = IWK( I)
250 CONTINUE
C --------------------------------------
C
200 CONTINUE
C ----------------------------------------------------------
NEL = J
C
RETURN
END
C **********************************************************************
C
FUNCTION ERF ( A )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
L = 100
C
DERF = 0.D0
DH = DBLE( A /DFLOAT( L))
C
C ---------------------------------------------------------------
DO 100 I = 1, L
C
DU = DH* DBLE( I)
DV = DH* DBLE( I-1)
C
DERF = DERF + ( DEXP( - DU* DU) + DEXP( - DV* DV))/ 2.D0* DH
C
100 CONTINUE
C ---------------------------------------------------------------
C
ERF = REAL( DERF* 2.D0/ DSQRT( 3.14159 265358 D0))
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FNSB ( SEM, ERR, * )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
WRITE(6,2000) SEM, ERR
2000 FORMAT(/' *** Final Error =',F9.3,' (%)',' Error/Loop =',F9.3)
WRITE(6,2100)
2100 FORMAT(/' Computation to be continued ? ( 0: No 1: Yes )')
C
ICTK = 0
IF ( ICTK.EQ.1) RETURN 1
C
CLOSE ( 9)
CLOSE (10)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FVCT ( F, C, UE, VE, WE, AV, DLX, DLY, DLZ,
1 VS, NDE, NNP, NEL)
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION C ( NNP), F ( NNP), DLX( NEL), DLY( NEL),
1 DLZ( NEL), AV( NEL), UE ( NEL), VE ( NEL),
2 WE ( NEL), VS( NEL), NDE( NEL,8), F1 ( 4,4),
3 F2 ( 4,4), F3( 4,4), FX ( 8,8), FY ( 8,8),
4 FZ ( 8,8)
C
CALL CLER ( F, 0.D0, NNP )
C
C ------------------------------------------------
DO 100 I = 1, 4
DO 100 J = 1, 4
C ------------------------------------------------
C
IJ = I + J
C
F1( I,J) = 2.D0
IF ( J.EQ.1.OR.J.EQ.4) F1( I,J) = - 2.D0
IF ( I.LE.2.AND.J.EQ.1.OR.I.GE.3.AND.J.EQ.4) F1( I,J) = - 4.D0
IF ( I.LE.2.AND.J.EQ.2.OR.I.GE.3.AND.J.EQ.3) F1( I,J) = 4.D0
C
F2( I,J) = 2.D0
IF ( J.EQ.1 .OR. J.EQ.2) F2( I,J) = - 2.D0
IF ( I.EQ.J.AND.J.LE.2.OR.IJ.EQ.5.AND.J.LE.2) F2( I,J) = - 4.D0
IF ( I.EQ.J.AND.J.GE.3.OR.IJ.EQ.5.AND.J.GE.3) F2( I,J) = 4.D0
C
F3( I,J) = - 2.D0
IF ( MOD( IJ, 2).EQ.0) F3( I,J) = -1.D0
IF ( I.EQ.J) F3( I,J) = -4.D0
C
100 CONTINUE
C ------------------------------------------------
DO 200 I = 1, 8
DO 200 J = 1, 8
C ------------------------------------------------
C
IF ( I.LE.4 .AND. J.LE.4) FX( I,J) = F1( I,J)
IF ( I.LE.4 .AND. J.LE.4) FY( I,J) = F2( I,J)
IF ( I.LE.4 .AND. J.LE.4) FZ( I,J) = F3( I,J)
C
IF ( I.LE.4 .AND. J.GT.4) FX( I,J) = 0.5D0* F1( I,J-4)
IF ( I.LE.4 .AND. J.GT.4) FY( I,J) = 0.5D0* F2( I,J-4)
IF ( I.LE.4 .AND. J.GT.4) FZ( I,J) = - F3( I,J-4)
C
IF ( I.GT.4 .AND. J.LE.4) FX( I,J) = 0.5D0* F1( I-4, J)
IF ( I.GT.4 .AND. J.LE.4) FY( I,J) = 0.5D0* F2( I-4, J)
IF ( I.GT.4 .AND. J.LE.4) FZ( I,J) = F3( I-4, J)
C
IF ( I.GT.4 .AND. J.GT.4) FX( I,J) = F1( I-4, J-4)
IF ( I.GT.4 .AND. J.GT.4) FY( I,J) = F2( I-4, J-4)
IF ( I.GT.4 .AND. J.GT.4) FZ( I,J) = - F3( I-4, J-4)
C
200 CONTINUE
C ----------------------------------------------------------
DO 300 IEL = 1, NEL
C
IF ( VS( IEL).LT.1.E-8 ) GO TO 300
VQX = UE( IEL)* AV( IEL)/ ( 4.D0* DLX( IEL))
C
C ----- Attention ---------------------------------
C
C ==> UE* 8* a*b*c/( 4* (2*a) ) ==> UE* b* c
C ==========
C -------------------------------------------------
C
VQY = VE( IEL)* AV( IEL)/ ( 4.D0* DLY( IEL))
VQZ = WE( IEL)* AV( IEL)/ ( 4.D0* DLZ( IEL))
C
C ------------------------------------------------
DO 320 I = 1, 8
C
NI = NDE( IEL, I)
C
FU = 0.D0
FV = 0.D0
FW = 0.D0
C
C --------------------------------------
DO 340 J = 1, 8
C
NJ = NDE ( IEL,J)
FU = FU + FX( I,J)* C( NJ)
FV = FV + FY( I,J)* C( NJ)
FW = FW + FZ( I,J)* C( NJ)
C
340 CONTINUE
C --------------------------------------
C
F ( NI) = F ( NI) + ( VQX* FU + VQY* FV + VQZ* FW )
C
320 CONTINUE
C ------------------------------------------------
C
300 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INIT ( U, V, W, C, UE, VE, WE, VS, NCL,
1 CC, ICC, IBC, RKX, RKY, RKZ, NDE, COF, DT,
2 NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 IBC
C
COMMON /COND/ UC, VC, WC, UIT, VIT, WIT, CIT
C
DIMENSION U ( NNP), V ( NNP), W ( NNP), C ( NNP),
1 CC( NNP), ICC( NNP), IBC( NNP), UE ( NEL),
2 VE( NEL), WE ( NEL), VS ( NEL), NDE( NEL,8),
3 RKX( NEL), RKY( NEL), RKZ( NEL), COF( NEL)
C
CALL CLER ( U, UIT, NNP )
C
CALL CLER ( V, VIT, NNP )
C
CALL CLER ( W, WIT, NNP )
C
CALL CLER ( C, CIT, NNP )
C
CALL CLR2 ( IBC, 1, NNP )
C
C --------------------------------------
DO 100 I = 1, NNP
C
U( I) = UC
V( I) = VC
W( I) = WC
C
100 CONTINUE
C --------------------------------------
C
IF ( NCL.LE.0 ) GO TO 900
C
C --------------------------------------
DO 200 I = 1, NCL
C
C ( ICC( I)) = CC( I)
IBC( ICC( I)) = 9999
C
200 CONTINUE
C --------------------------------------
C
900 CONTINUE
C
C ------------------------------------------------------------
C
CALL VLEM ( NNP, NEL, NDE, U, V, W, UE, VE, WE, VS )
C
C ------------------------------------------------------------
C
CALL CFDF ( NEL, RKX, RKY, RKZ, UE, VE, WE, COF, DT )
C
C ------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INPT ( X, Y, Z, NDE, ITB, AV, DLX, DLY,
1 DLZ, NCL, CC, ICC, DT, LPE, NNP, NEL,
2 NNZ, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /RTIM/ KEY
COMMON /MSZE/ NLAY( 3)
COMMON /CON1/ OMG, EROR, NTR
COMMON /CON2/ COEF, NPB, PHI( 20), SMH, DTD
COMMON /COND/ UC, VC, WC, UIT, VIT, WIT, CIT
C
C ----- NNP ---
C
DIMENSION X ( 2541), Y ( 2541), Z( 2541), CC( 2541),
1 ICC( 2541), ITB( 2541)
C
C ----- NEL ---
C
DIMENSION AV ( 2000), IEB( 2000), DLX( 2000), DLY( 2000),
1 DLZ( 2000)
C
C ----- NEL*8 ---
C
DIMENSION NDE( 2000,8)
C
C ---------------------------------------------------------------
DIMENSION LYO( 50,3), CRD( 50,3), SCL( 3), NLYE( 3),
1 IVK( 2), VK ( 2)
C ---------------------------------------------------------------
C
WRITE(6,*) ' '
WRITE(6,*) ' ** DELT, LOP = ?'
*** READ (5,*) DT, LPE
C
DT = 0.05D0
LPE = 20
C ==================
C
C -------------------------------------------
C
CALL ZRVC ( 3, NLYE )
C
C -------------------------------------------
C
CALL DMST ( NNP, NEL, NLYE, NFIX )
C
C -------------------------------------------
REWIND 9
C
KEY = 0 ! With Correction for Difuusivity
C ====================================
OMG = 1.1D0
EROR = 0.0001D0
NTR = 100
C
NPB = 1
PHI( 1)= 1.D0
C
WRITE(6,2000) KEY, OMG, EROR, NTR, NPB, ( PHI( I),I = 1, NPB )
2000 FORMAT(/,' * KEY =',I2,' OMG =',F5.2,' EROR =',F7.4,' NTR =',
1 I4,' NPB =',I2,' PHI =',10F5.2)
C
C ----- Attention ---
C
COEF = 0.025D0
C ===================
C
SMH = 0.1D0
C ===================
C
DTD = 0.0005D0
C ===================
C
SCL( 1) = 1.D0
SCL( 2) = 1.D0
SCL( 3) = 1.D0
C
UIT = 0.D0
VIT = 0.D0
WIT = 0.D0
C
CIT = 0.D0
C
CALL ZRVC ( 150, LYO )
C
CALL ZRVC ( 3, NLAY )
C
REWIND 9
C ------------------------------------------------
DO 100 I = 1, 3
DO 100 L = 1, NLYE( I)
C ------------------------------------------------
C
READ(9,1000) NO, COR
1000 FORMAT(13X,I2,5X,F10.0)
C
C ----- To store sequentially ---
C
NLA = NLAY( I)
C
C --------------------------------------
DO 120 J = 1, NLA
C
IF ( NO.GT.LYO( J,I)) GO TO 120
C ----------------------------
DO 140 K = J, NLA
C
K1 = NLAY( I) + J - K
LYO( K1+1,I) = LYO( K1,I)
CRD( K1+1,I) = CRD( K1,I)
C
140 CONTINUE
C ----------------------------
LYO( J,I) = NO
CRD( J,I) = COR
GO TO 70
C
120 CONTINUE
C --------------------------------------
C
J = NLAY( I) + 1
LYO( J,I) = NO
CRD( J,I) = COR
70 NLAY( I) = NLAY( I) + 1
C
100 CONTINUE
C --------------------------------------------------------------------
C
CALL DTGN ( LYO, CRD, NNP, NEL, X, Y, Z, NDE, ITB, IEB, SCL )
C
C --------------------------------------------------------------------
REWIND 9
C
IF ( NFIX.EQ.0) GO TO 800
C
C ------------------------------------------------
DO 200 L = 1, NFIX
C
READ(9,1100) ( IVK( J), VK( J), J = 1, 2)
1100 FORMAT(19X,I6,5X,F10.0,4X,I6,5X,F10.0)
C
C ------------------------------------------------
DO 200 J = 1, 2
C
IF ( IVK( J).EQ.0) GO TO 200
K = NSER( IVK( J),ITB, NNP )
NCL = NCL + 1
ICC( NCL) = K
CC ( NCL) = VK( J)
C
200 CONTINUE
C ------------------------------------------------
C
800 CONTINUE
C ---------------------------------------------------------------------
C
UC = 0.5D0
C ===============
VC = 0.D0
WC = 0.D0
C
C ---------------------------------------------------------------------
C
CALL AGCL ( NNP, NEL, X, Y, Z, NDE, AV, DLX, DLY,
1 DLZ, BV, RV, DT )
C
C ---------------------------------------------------------------------
C
CALL ARST ( NNP, NEL, NDE, NNZ, NZ2 )
C
C ---------------------------------------------------------------------
WRITE(6,2100) NNP, NEL, NNZ, NZ2, NLAY
2100 FORMAT(' ** NNP =',I5,' NEL =',I5,' NNZ =',I7,' NZ2 =',I7,
1 ' NLAY=',3I3)
WRITE(6,2200) COEF, SMH, DTD, UC, VC, WC
2200 FORMAT(/,' ** COEF, SMH, DTD =',3F7.3,' * UC =',F6.3,' VC =',F6.3,
1 ' WC =',F6.3)
WRITE(6,2400) DT, LPE, BV, RV
2400 FORMAT(/,' * DT =',F7.3,' LPE =',I5,' B =',F7.3,2X,
1 ' R =', F8.4,/)
WRITE(6,*) '** Computation started **'
WRITE(6,*) ' '
C --------------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MOVEI ( IX, IA, IY, NQ )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X( NQ), Y( NQ), IX( NQ), IY( NQ)
C
C ----------------------------
DO 100 I = 1, NQ
C
IY( I) = IX( I)* IA
100 CONTINUE
C ----------------------------
RETURN
C
C --------------------------------------
ENTRY MOVER ( X, A, Y, N)
C
C ----------------------------
DO 200 I = 1, N
C
Y( I) = X( I)* A
200 CONTINUE
C ----------------------------
C
RETURN
END
C **********************************************************************
C
FUNCTION NSER ( I, IT, N)
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- Table by binary search:
C IT: Table N: Size of table I: Search of value ---
C ----------------------------
C
DIMENSION IT( 1)
C
I1 = 0
I2 = N + 1
10 K = ( I1 + I2) / 2
IF ( K.EQ.I1 .OR. K.EQ.I2) GO TO 50
IF ( I - IT( K)) 20, 40, 30
20 I2 = K
GO TO 10
C
30 I1 = K
GO TO 10
C
C ----- Found ---
C
40 NSER = K
RETURN
C
C ----- Not found ---
C
50 NSER = 0
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PMAT ( P, AV, NDE, NCW, KCC, NCC, NNP, NEL,
1 NNZ, NZ2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 KCC, NCC
C
DIMENSION KCC( NNP), AV ( NEL), NDE( NEL,8), P ( NNZ),
1 NCW( NNZ), NCC( 2,NZ2), P1 ( 4,4), PP( 8,8)
C
CALL ZARY ( 4, 4, P1 )
C
C --------------------------------------
DO 100 I = 1, 4
DO 100 J = 1, 4
C --------------------------------------
C
IJ = I + J
P1( I,J) = 4.D0
IF ( MOD( IJ,2).EQ.0) P1( I,J) = 2.D0
IF ( I.EQ.J) P1( I,J) = 8.D0
C
100 CONTINUE
C --------------------------------------
C
CALL ZARY ( 8, 8, PP )
C
C --------------------------------------
C
CALL CLER ( P, 0.D0, NNZ )
C
C --------------------------------------
NZWK = 0
C --------------------------------------
DO 200 IEL = 1, NEL
C
AVQ = AV( IEL)/ 12.D0
C
C ----- Attention ---
C
C ==> 8* a* b* c/12 ==> 2* a* b* c/3
C ===============
C
C ---------------------------
DO 250 I = 1, 8
C
DO 250 J = 1, 8
C ---------------------------
C
II = NDE( IEL,I)
JJ = NDE( IEL,J)
C
IF( I.LE.4.AND.J.LE.4) PP( I,J) = P1( I,J)
IF( I.LE.4.AND.J.GT.4) PP( I,J) = P1( I,J-4)* 0.5D0
C
IF( I.GT.4.AND.J.LE.4) PP( I,J) = P1( I-4,J)* 0.5D0
IF( I.GT.4.AND.J.GT.4) PP( I,J) = P1( I-4,J-4)
C
AVQQ = AVQ* PP( I,J)
C ----------------------------------------------------
C
CALL CNAL ( AVQQ, II,JJ, P, NCC, NZWK, NNZ )
C
C ----------------------------------------------------
C
250 CONTINUE
C ----------------------------
C
200 CONTINUE
C --------------------------------------
C
CALL MOVEI ( NCC, 1, NCW, NNZ )
C
C --------------------------------------
DO 300 II = 1, NZ2
C
JJ = NCC( 2,II)
NCC( 2,II) = NCC( 1,II)
NCC( 1,II) = JJ
C
300 CONTINUE
C ----------------------------------------------
C
CALL CNIJ ( KCC, NCC, NZWK, NNP, NZ2 )
C
C ----------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SOLV ( NNP, NNZ, NZ2, FVC, C, TSM, NCC, KCC, NBC,
1 LOP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NCC, KCC, NBC
C
COMMON /CON1/ OMG, EROR, NTR
C
DIMENSION FVC( NNP), C ( NNP), KCC( NNP), NBC( NNP),
1 TSM( NNZ), NCC( 2, NZ2)
C
C ----- S.O.R. method --- Attention --- TSM ( NNZ)---
C
CXX = 0.D0
C
C --------------------------------------
DO 100 I = 1, NNP
C
AC = DABS( C( I))
IF ( AC.GE.CXX) CXX = AC
C
100 CONTINUE
C --------------------------------------
C
EPS = CXX* 0.01D0
LMR = 3* NNP
KTL = 2* NNP
C
IF ( NTR.GE.1) KTL = NTR
C
NTL = 0
990 NTL = NTL + 1
ERX = 0.D0
C
C --------------------------------------
DO 200 I = 1, NNP
C
IF ( NBC( I).EQ.9999 ) GO TO 200
IF ( I.EQ.1) GO TO 210
MS = KCC( I-1) + 1
GO TO 220
C
210 MS = 1
220 ML = KCC( I)
IF ( MS.GT.ML) GO TO 1000
DK = FVC( I)
C
C ----------------------------
DO 250 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 250
C
DK = DK - TSM( M1)* C( M2)
C
250 CONTINUE
C ----------------------------
C
DC = - C( I) + DSM* DK
C( I) = C( I) + OMG* DC
C
TEMP = DABS ( C( I))
C
IF ( TEMP.LE.EPS) GO TO 200
TEMP = DABS( DC/C( I))
IF ( TEMP.LE.ERX) GO TO 200
ERX = TEMP
NERX = I
C
200 CONTINUE
C --------------------------------------
C
IF ( ERX.LE.EROR ) GO TO 700
IF ( NTL.GE.KTL ) RETURN
GO TO 990
C
700 CONTINUE
IF ( MOD( LOP,5).EQ.0 )
C ============================
C
1WRITE(6,2000) LOP, EPS, NTL
2000 FORMAT(' * Loop =',I4,' EPS =',F9.5,' No.of Iter.=',I4 )
RETURN
C
1000 NER = 201
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SUBS ( PDV, M, N, PD, NCC, K, NNZ )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION PD( NNZ), NCC( NNZ)
C
K = K + 1
C
IF ( K.GT.NNZ) WRITE(6,2000)
2000 FORMAT(' * Error K.GT.NNZ *')
C
JE = K - N
C ---------------------------------
DO 100 JJ = 1, JE
C
J1 = K - JJ + 1
PD ( J1) = PD ( J1-1)
NCC( J1) = NCC( J1-1)
C
100 CONTINUE
C ---------------------------------
C
PD ( N) = PDV
NCC( N) = M
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TOUT ( C, CAA, ITB, LOP, CA, CN, DT, LPE,
1 SEM, ERR, NNP, * )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /MSZE/ NLAY( 3)
COMMON /XYZ / XP( 1500), YP( 1500), ZP( 1500)
COMMON /CON2/ COEF, NPB, PHI( 20), SMH, DTD
C
DIMENSION C ( NNP), ITB( NNP), CAA( NNP), CA( NNP), CN( NNP),
1 ARA( 12), CRA( 12)
C
IF ( LOP.EQ.0) THEN
C --------------------------------------
C
CALL CLER ( C, 0.D0, NNP )
C
C --------------------------------------
C
CALL ANAS ( NNP, CAA, LOP, DT )
C
C --------------------------------------
END IF
C
C ----- Analytical Solution --------------
C
CALL ANAS ( NNP, CAA, LOP, DT )
C
C ----------------------------------------
C
TM = DT* DFLOAT( LOP)
CAM = 0.D0
C ------------------------------------------------
DO 100 I = 1, NNP
C
IF ( CAA( I).GT.CAM) CAM = CAA( I)
100 CONTINUE
C ------------------------------------------------
IF ( CAM.LT.1.D-8) CAM = 1.D0
C
L = 0
C ------------------------------------------------
DO 200 I = 1, NNP
C
II = ITB( I) - 10117
C ******* Attention ********
C
IF ( II.LE.0) THEN
L = L + 1
CA( L) = CAA( I)/ CAM
CN( L) = C ( I)/ CAM
END IF
C
200 CONTINUE
C ------------------------------------------------
DO 300 I = 1, 11
C
ARA( I+1) = CA( I)
CRA( I+1) = CN( I)
C
300 CONTINUE
C ------------------------------------------------
C
IF ( LOP.NE.LPE) GO TO 333
C
* WRITE(6,*) ' '
* WRITE(6,*) '*** Ratio for Analytical Solution **'
* WRITE(6,2300) ( ARA( I), I= 2, 12)
* 2300 FORMAT(2X,11F7.3)
* WRITE(6,*) ' '
* WRITE(6,*) '*** Ratio for Numerical Solution ** '
* WRITE(6,2300) ( CRA( I), I = 2, 12 )
C
333 CONTINUE
IF ( LOP.EQ.0 ) THEN
WRITE(10,2000) COEF, DT, LPE
2000 FORMAT(2F7.4,I3)
WRITE(10,2100) ( NLAY(I), I=1, 3 )
2100 FORMAT(3I3)
WRITE(10,2200) ( ARA( I), CRA( I), I = 2, 12 )
2200 FORMAT(10F8.4)
END IF
C
C -------------------------------------------------------
C
IF ( LOP.EQ.LPE/2) WRITE(10,2200) ( ARA(I), CRA(I), I=2,12 )
C
C -------------------------------------------------------
C
IF ( LOP.EQ.LPE) WRITE(10,2200) ( ARA(I), CRA(I), I=2,12 )
C
C -------------------------------------------------------
C
IF ( LOP.NE.0) THEN
C ------------------------------------------------
C
CALL CPER ( NNP, C, CAA, LOP, SEM, DT )
C
C ------------------------------------------------
ERR = SEM/ DFLOAT( LOP)
END IF
C
IF ( LOP.LT.LPE) RETURN 1
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,100)
WRITE(6,101) 'C *********************************************** C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C CVDF3D - FEM C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C Three-Dimensional Convection-Diffusion C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C Analysis by FEM C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C ( Modified Galerkin Method ) C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C ( Rectangular Parallelepiped Elements ) C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C ( V. 3.8 ) C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C Copyright 2011 : Yasuhiro MATSUDA C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C *********************************************** C'
100 FORMAT(/)
101 FORMAT(2X,A55)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VLEM ( NNP, NEL, NDE, U, V, W, UE, VE, WE, VS )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION U ( NNP), V ( NNP), W ( NNP), UE( NEL), VE( NEL),
1 WE( NEL), VS( NEL), NDE( NEL,8)
C
C ------------------------------------------------
DO 100 I = 1, NEL
C
UT = 0.D0
VT = 0.D0
WT = 0.D0
C --------------------------------------
DO 150 J = 1, 8
C
N = NDE( I,J)
UT = UT + U( N)
VT = VT + V( N)
WT = WT + W( N)
C
150 CONTINUE
C -------------------------------------
UE( I) = UT/ 8.D0
VE( I) = VT/ 8.D0
WE( I) = WT/ 8.D0
C
100 CONTINUE
C ------------------------------------------------
C
CALL CLER ( VS, 0.D0, NEL )
C
C ------------------------------------------------
C
VMAX = 0.D0
C
C ------------------------------------------------
DO 200 I = 1, NEL
C
VS( I) = DSQRT( UE(I)* UE(I) + VE(I)* VE(I) + WE(I)* WE(I))
IF ( VS( I).GE.VMAX) VMAX = VS( I)
C
200 CONTINUE
C ------------------------------------------------
VMT = VMAX* 0.01D0
C
C ------------------------------------------------
DO 300 I = 1, NEL
C
IF ( VS( I).LT.VMT) VS( I) = 0.D0
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ZARY ( M, N, A )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION A( M, N)
C
C -------------------------
DO 100 I = 1, M
DO 100 J = 1, N
C
A( I,J)= 0.D0
100 CONTINUE
C -------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ZRVC ( M, IB )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION IB( M)
C
C ------------------------
DO 100 I = 1, M
C
IB( I) = 0
100 CONTINUE
C ------------------------
C
RETURN
END
C *******************************************************************
C ******************** DATA : CVDF3D_FEM.DAT***************
LAYER X 1 -1.000
LAYER X 2 -0.875
LAYER X 3 -0.75
LAYER X 4 -0.625
LAYER X 5 -0.5
LAYER X 6 -0.375
LAYER X 7 -0.25
LAYER X 8 -0.125
LAYER X 9 0.
LAYER X 10 0.125
LAYER X 11 0.25
LAYER X 12 0.375
LAYER X 13 0.5
LAYER X 14 0.625
LAYER X 15 0.75
LAYER X 16 0.875
LAYER X 17 1.000
LAYER Y 1 0.000
LAYER Y 2 0.125
LAYER Y 3 0.25
LAYER Y 4 0.375
LAYER Y 5 0.5
LAYER Y 6 0.625
LAYER Y 7 0.75
LAYER Y 8 0.875
LAYER Y 9 1.000
LAYER Z 1 0.
LAYER Z 2 0.125
LAYER Z 3 0.25
LAYER Z 4 0.375
LAYER Z 5 0.5
LAYER Z 6 0.625
LAYER Z 7 0.75
LAYER Z 8 0.875
LAYER Z 9 1.000
END
C *************************************************************