戻る

PROGRAM MAIN
C
C ******************************************************************** C
C C
C DIFF3D - FEM61 C
C C
C Three-Dimensional Diffusion Analysis by FEM C
C C
C ( Linear-Hexahedra Element ) C
C C
C ( V.3.0 ) 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 C ( 2541), CC ( 2541), ICC( 2541), F ( 2541),
1 B ( 2541), ITB( 2541), CAA( 2541), KCC( 2541),
2 IBC( 2541), CA ( 2541), CN ( 2541)
C
DIMENSION X( 2541), Y( 2541), Z( 2541)
C
C ----- NEL ---
DIMENSION RKX( 2000), RKY( 2000), RKZ( 2000), DLX( 2000),
1 DLY( 2000), DLZ( 2000), AV ( 2000), VS ( 2000),
2 COF( 2000)
C
C ----- NEL,8 ---
DIMENSION NDE( 2000,8)
C
C ----- NNZ ---
DIMENSION P( 50000), D( 50000), A( 50000), NCW( 500000)
C
C ----- NNZ2*2 ---
DIMENSION NCC( 2,100000)
C
COMMON /CN2/ COEF, NMP, PHI( 20), DND, DTD
C
OPEN ( 9, FILE='DIFF3D.DAT' )
OPEN (10, FILE='DIFF3D_FIG_OUT.DATA' )
C
CALL TTLE
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 ( C, VS, NCL, CC, ICC, IBC, RKX, RKY, RKZ,
1 NDE, COF, DT, 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, NMP
C
PHA = PHI( I)
C
LOP = 0
SEM = 0.D0
C
CALL ROUT ( CAA, C, ITB, LOP, CA, CN, DT, LPE, SEM, ERR,
1 NNP, *999 )
C
999 LOP = LOP + 1
C
CALL FVCT ( F, C, AV, DLX, DLY, DLZ, VS, 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 ROUT ( C, CAA, ITB, LOP, CA, CN, DT, LPE, SEM, ERR,
1 NNP, *999 )
C
100 CONTINUE
C ----------------------------------------------------------
C
CLOSE ( 9)
CLOSE (10)
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 CLAR ( A, 0.D0, NNZ )
C
C --------------------------------------
C
CALL CLAR ( B, 0.D0, NNP )
C
C --------------------------------------
C
W1 = 1.D0/ DT
W2 = 1.D0 - PHI
C ----------------------------------------------
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
C --------------------------------------
DO 200 I = 1, NNZ
C
A( I) = D( I)* PHI + P( I)* W1
200 CONTINUE
C --------------------------------------
C
C {B} = {B} - {F}
C
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, B, R, DT )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CND/ CIT
COMMON /CN2/ COEF, NMP, PHI( 20), DND, DTD
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)
C
X1 = X( N1)
Y1 = Y( N1)
Z1 = Z( N1)
C
N2 = NDE( I,2)
X2 = X( N2)
C
N4 = NDE( I,4)
Y4 = Y( N4)
C
N5 = NDE( I,5)
Z5 = Z( N5)
C
DLX( I) = DABS( X2 - X1)
DLY( I) = DABS( Y4 - Y1)
DLZ( I) = DABS( Z5 - Z1)
C
AV( I) = DLX( I)* DLY( I)* DLZ( I)
C
IF ( AV( I).GT.0.D0) GO TO 50
IF ( AV( I).LE.0.D0)
1 WRITE(6,*) ' ** DATA: Strange ** STOP ** '
IF ( AV( I).LE.0.D0) STOP
C
50 CONTINUE
C
100 CONTINUE
C ------------------------------------------------
DLM = 0.D0
DO 200 I = 1, NEL
C
IF ( DLX( I).GT.DLM ) DLM = DLX( I)
200 CONTINUE
C ------------------------------------------------
C
UC = 0.D0
*** B = UC* DT / DLM
R = 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 /CND/ CIT
COMMON /CN2/ COEF, NMP, PHI( 20), DND, DTD
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
DND = 0.1D0
DX = DND
DY = DND
DZ = DND
C
T = FLOAT( LOP)* DT
C
RR = DABS( RKX) + DABS( RKY) + DABS( RKZ)
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))
1 / ( 2.D0* DSQRT( RKX* T + RKT))
CF1 = ERF( CF1)
C
CF2 = ( DX/ 2.D0 + ( XP( I) + 0.5D0))
1 / ( 2.D0* DSQRT( RKX* T + RKT ))
CF2 = ERF( CF2)
C
CF3 = ( DY/ 2.D0 - ( YP( I)))/ ( 2.D0* DSQRT( RKY* T + RKT ))
CF3 = ERF( CF3)
C
CF4 = ( DY/ 2.D0 + ( YP( I)))/ ( 2.D0* DSQRT( RKY* T + RKT ))
CF4 = ERF( CF4)
C
CF5 = ( DZ/ 2.D0 - ( ZP( I)))/ ( 2.D0* DSQRT( RKZ* T + RKT ))
CF5 = ERF( CF5)
C
CF6 = ( DZ/ 2.D0 + ( ZP( I)))/ ( 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
FUNCTION ERF ( A )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
L = 100
C
DRF = 0.D0
DH = DBLE( A/ DFLOAT( L))
C
C ----------------------------------------------------------
DO 100 I = 1, L
C
DU = DH* DBLE( I)
DV = DH* DBLE( I-1)
DRF= DRF + ( DEXP( - DU* DU) + DEXP( - DV* DV))/ 2.D0* DH
C
100 CONTINUE
C ----------------------------------------------------------
C
ERF = DRF* 2.D0/ DSQRT( 3.1415 92653 58 D0)
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
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
C
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 CLAR ( 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 CONAL ( 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. ofO NDE.LE.32768 ---
C
IF ( I.GT.J) GO TO 99
C
M2( 1) = J
M2( 2) = I
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
IF ( N.NE.N2) GO TO 100
IF ( M - NCC( N1)) 31, 32, 33
31 N = N1
C
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 **')
C
PD ( K) = PDV
NCC( K) = M
GO TO 99
C
100 CONTINUE
C ------------------------------------------------
C
WRITE(6,2100) K
2100 FORMAT(' * Error -- DO Loop Strange, K =',I10)
C
STOP
C
99 RETURN
C
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 NNP, NZ2 : Max. length of KCC, NCC
C
CHARACTER IANNP * 4, IANZ2 * 4
C
DATA IANNP /'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)
C
IF ( NP.GT.NNP) GO TO 30
10 NCC( 1,I) = I
C
100 CONTINUE
C ------------------------------------------------
C
C ---- Last row's counter ---
C
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 ---
C
GO TO 220
C
20 CONTINUE
C
C ---- To shift NCC ---
C
NZIS = NZ - IS
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
C ------ To increase KCC ---
C
C ----------------------------
DO 280 L = I, NP
C
KCC( L) = KCC( L) + 1
280 CONTINUE
C ----------------------------
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) IANNP, NNP, IANNP, NP
2000 FORMAT(' ** ErroR( CNIJ) ',A3,'(',I5,') too small',A2,'=',I5)
STOP
C
40 WRITE(6,2000) IANZ2, NZ2, IANZ2, NZ
STOP
C
END
C **********************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,100)
WRITE(6,101) 'C ********************************************** C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C DIFF3D - FEM_61 C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C Three-Dimensional Diffusion Analysis by FEM C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C ( Linear Hexahedra Element ) C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C ( V.3.0) 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'
WRITE(6,*) ' '
C
100 FORMAT(/)
101 FORMAT( 12X,A55)
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
C -------------------
NFIX = 0
C ----- Attention ---
C
REWIND 9
C
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
C
40 CONTINUE
READ(9,1200) DIR
1200 FORMAT(10X,A4)
C
NFIX = NFIX + 1
GO TO 10
C
90 NNP = NLAY( 1)* NLAY( 2)* NLAY( 3)
NEL = ( NLAY( 1)-1)* ( NLAY( 2)-1)* ( NLAY( 3)-1)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DMAT ( D, AV, DLX, DLY, DLZ, COF, NDE, RKX,
1 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 CLAR ( D, 0.D0, NNZ)
C
NZWK = NNZ
CMX = 0.D0
C
C -------------------------------------------
DO 100 I = 1, NEL
C
IF ( COF( I).GE.CMX) CMX = COF( I)
100 CONTINUE
CMX = CMX * 0.001D0
C -------------------------------------------
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))
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
CALL CONAL ( AVQQ, II,JJ, D, NCW, NZWK, NNZ )
C
350 CONTINUE
C --------------------------------------
C
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SOLV ( NNP, NNZ, NZ2, FVEC, C, TSM, NCC, KCC,
1 NBC, LOP)
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NCC, KCC, NBC
C
COMMON /CN1/ OMG, EROR, NITR
C
DIMENSION FVEC( NNP), C ( NNP), KCC( NNP), NBC( NNP),
1 TSM ( NNZ), NCC( 2,NZ2)
C
C ----- S.O.R. method --- Attention --- TSM ( NNZ) ---
C
CMX = 0.D0
C
DO 100 I = 1, NNP
C
AC = DABS( C( I))
IF ( AC.GE.CMX) CMX = AC
C
100 CONTINUE
C
EPS = CMX* 0.01D0
LMPR = 3* NNP
KITL = 2* NNP
C
IF ( NITR.GE.1) KITL = NITR
NITL = 0
990 NITL = NITL + 1
C
ERMX = 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 = FVEC( 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
DK = DK - TSM( M1)* C( M2)
C
250 CONTINUE
C --------------------------------------
DC = - C( I) + DSM * DK
C( I) = C( I) + OMG * DC
TEMP = DABS ( C( I))
IF ( TEMP.LE.EPS) GO TO 200
TEMP = DABS ( DC/ C( I))
IF ( TEMP.LE.ERMX) GO TO 200
ERMX = TEMP
NERMX = I
C
200 CONTINUE
C ------------------------------------------------
C
IF ( ERMX.LE.EROR) GO TO 700
IF ( NITL.GE.KITL) GO TO 701
IF ( NITL.GE.LMPR) GO TO 800
GO TO 990
C
800 CONTINUE
GO TO 990
C
700 IF ( MOD( LOP,10).EQ.0)
1WRITE(6,2220) LOP, EPS, NITL
2220 FORMAT(' * LOP =',I5,' EPS =',F10.5,' No.of Iter.=',I5 )
RETURN
C
701 RETURN
C
1000 NER = 201
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DTGN ( LYO, CRD, NNP, NEL, X, Y, Z, NDE, ITB, IETB,
1 SCL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C 3- Dimensional NDE Data, Element Data Generation
C
DIMENSION X ( NNP), Y ( NNP), Z ( NNP), ITB( NNP),
1 IETB( NEL), NDE( NEL,8), LYO( 50,3), CRD( 50,3),
2 SCL ( 3), IWK( 8)
C
COMMON /MSZ/ 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
C
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
IETB( J)= J
IWK ( 1) = J3* ( K3-1) + J2* ( K2-1) + J1* ( K1-1) + 1
C
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
SUBROUTINE FVCT ( F, C, AV, DLX, DLY, DLZ, VS, NDE, NNP,
1 NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION C ( NNP), F ( NNP), DLX( NEL), DLY( NEL),
1 DLZ( NEL), AV( NEL), VS ( NEL), NDE( NEL,8),
2 F1 ( 4,4), F2( 4,4), F3 ( 4,4), FX ( 8,8),
3 FY ( 8,8), FZ( 8,8)
C
CALL CLAR ( F, 0.D0, NNP)
C
C ----------------------------------------------------------
DO 100 I = 1, 4
DO 100 J = 1, 4
C ----------------------------------------------------------
C
IJ = I + J
F1( I,J) = 2.D0
C
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)
1 F1( I,J) = - 4.D0
IF ( I.LE.2.AND.J.EQ.2 .OR. I.GE.3.AND.J.EQ.3)
1 F1( I,J) = 4.D0
F2( I,J) = 2.D0
C
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)
1 F2( I,J) = - 4.D0
IF ( I.EQ.J.AND.J.GE.3 .OR. IJ.EQ.5.AND.J.GE.3)
1 F2( I,J) = 4.D0
F3( I,J) = - 2.D0
C
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 ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INIT ( C, VS, NCL, CC, ICC, IBC, RKX, RKY,
1 RKZ, NDE, COF, DT, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 IBC
C
COMMON /CND/ CIT
C
DIMENSION C ( NNP), CC( NNP), ICC( NNP), IBC( NNP),
1 VS ( NEL), NDE( NEL,8), RKX( NEL), RKY( NEL),
2 RKZ( NEL), COF( NEL)
C
CALL CLAR ( C, CIT, NNP)
C
CALL CLR2 ( IBC, 1, NNP)
C
IF ( NCL.LE.0) GO TO 900
C
C ---------------------------------
DO 100 I = 1, NCL
C
C ( ICC( I)) = CC( I)
IBC( ICC( I)) = 9999
C
100 CONTINUE
C ---------------------------------
C
900 CONTINUE
C
CALL CLAR ( VS, 0.D0, NEL )
C
CALL CFDF ( NEL, RKX, RKY, RKZ, COF, DT )
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 /TIM/ KEY
COMMON /CND/ CIT
COMMON /MSZ/ NLAY( 3)
COMMON /CN1/ OMG, EROR, NITR
COMMON /CN2/ COEF, NMP, PHI( 20), DND, DTD
C
C ----- NNP ---
DIMENSION X ( 2541), Y ( 2541), Z( 2541), CC( 2541),
1 ICC( 2541), ITB( 2541)
C
C ----- NEL ---
DIMENSION AV ( 2000), IETB( 2000), DLX( 2000), DLY( 2000),
1 DLZ( 2000)
C
C ----- NEL*8 ---
DIMENSION NDE( 2000,8)
C
DIMENSION LYO( 50,3), CRD( 50,3), SCL( 3), NLYE( 3),
1 IVK( 2), VK( 2)
C
WRITE(6, *) ' ** DELT, LOP = ?'
*** READ(5,*) DT, LPE
C
DT = 0.005D0
LPE = 40
C
CALL ZRVC ( 3, NLYE )
C
CALL DMST ( NNP, NEL, NLYE, NFIX )
C
REWIND 9
C
KEY = 0 ! Correction for Difuusivity
C
OMG = 1.1D0
EROR = 0.0001D0
NITR = 100
C
NMP = 1
PHI(1) = 1.D0
C
WRITE(6,2000) KEY, OMG, EROR, NITR, NMP,( PHI( I),I = 1, NMP)
2000 FORMAT(/,' * KEY =',I2,' OMG =',F6.2,' EROR =',F8.5,' NITR =',I4,
1 ' NMP =',I2,' PHI=',10F6.2 )
C
COEF = 0.2D0
C
DND = 0.1D0
DTD = 0.0005D0
C
SCL( 1) = 1.D0
SCL( 2) = 1.D0
SCL( 3) = 1.D0
C
CIT = 0.D0
C
WRITE(6,2100) COEF, DND, DTD
2100 FORMAT(/,' ** COEF, DND, DTD =',3F7.3)
C
WRITE(6,2200) SCL( 1), SCL( 2), SCL( 3), CIT
2200 FORMAT(' * SCL=',3F7.2,' CIT =',4F7.2 )
C
CALL ZRVC ( 150, LYO )
C
CALL ZRVC ( 3, NLAY )
C
C ----------------------------------------------------------
DO 100 I = 1, 3
DO 100 L = 1, NLYE( I)
C ----------------------------------------------------------
C
READ(9,1700) NO, COR
1700 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
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 ------------------------------------------------
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, IETB, SCL )
C
C --------------------------------------------------------------------
*** UC = 0.5D0
UC = 0.D0
VC = 0.D0
WC = 0.D0
C
IF ( NFIX.EQ.0) GO TO 800
C
C ----------------------------------------------------------
DO 200 L = 1, NFIX
C
READ(9,1900) ( IVK( J), VK( J), J = 1, 2)
1900 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 = NSEAR ( IVK( J), ITB, NNP)
NCL = NCL + 1
ICC( NCL) = K
CC ( NCL) = VK( J)
C
200 CONTINUE
C ----------------------------------------------------------
800 CONTINUE
C
C --------------------------------------------------------------------
C
CALL AGCL ( NNP, NEL, X, Y, Z, NDE, AV, DLX, DLY, DLZ,
1 B, R, DT )
C
C --------------------------------------------------------------------
C
CALL ARST ( NNP, NEL, NDE, NNZ, NZ2 )
C
C ------------------------------------------------
C
WRITE(6,2300) NNP, NEL, NNZ, NZ2, NLAY, DT, LPE, R
2300 FORMAT(/,' *** NNP =',I5,' * NEL =',I5,' * NNZ =',I6,' * NZ2 =',I6
1 ,' * NLY=',3I4,//,' * DT =',F7.3,' * LPE =',I5,
2 5X,' * R =', F8.4 )
C
WRITE(6,*) ' '
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CFDF ( NEL, RKX, RKY, RKZ, COF, DT )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /TIM/ KEY
COMMON /CN2/ COEF, NMP, PHI( 20), DND, DTD
C
DIMENSION COF( NEL), 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
COF( I) = RKX( I) + RKY( I) + RKZ( I)
C
100 CONTINUE
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
DO 100 I = 1, NQ
C
IY( I) = IX( I)*IA
100 CONTINUE
RETURN
C
ENTRY MOVER ( X, A, Y, N )
C
DO 200 I = 1, N
C
Y( I) = X( I)*A
200 CONTINUE
RETURN
END
C **********************************************************************
C
FUNCTION NSEAR ( I, IT, N )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- Table by binary search ---
C
C IT : Table N : Size of table I : Search of value
C
DIMENSION IT( 1)
C
I1 = 0
I2 = N + 1
10 K = ( I1 + I2)/ 2
C
IF ( K.EQ.I1 .OR. K.EQ.I2) GO TO 50
IF ( I - IT( K)) 20, 40, 30
C
20 I2 = K
GO TO 10
C
30 I1 = K
GO TO 10
C
C ----- Found ---
C
40 NSEAR = K
RETURN
C
C ----- Not found ---
C
50 NSEAR = 0
RETURN
END
C **********************************************************************
C
SUBROUTINE NMER ( 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
C
SUMS = 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 -------------------------------------------
CAAM = 0.D0
DO 200 I = 1, NNP
C
IF ( CAA( I).GT.CAAM) CAAM = CAA( I)
200 CONTINUE
C -------------------------------------------
DO 300 I = 1, NNP
C
IF ( CAAM.LT.1.D-8) GO TO 300
SUMS = SUMS + ( DABS( C( I) - CAA( I))/ DABS( CAAM))
C
300 CONTINUE
C -------------------------------------------
SUMS = ( SUMS* 100.D0/ DFLOAT( NNP))
SEM = SEM + SUMS
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
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
CALL CLAR ( P, 0.D0, NNZ )
C
NZWK = 0
C
C ----------------------------------------------------------
DO 200 IEL = 1, NEL
C
AVQ = AV( IEL)/ 12.D0
C ------------------------------------------------
DO 250 I = 1, 8
DO 250 J = 1, 8
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) = 0.5D0 * P1( I,J-4)
C
IF( I.GT.4.AND.J.LE.4) PP( I,J) = 0.5D0 * P1( I-4, J)
IF( I.GT.4.AND.J.GT.4) PP( I,J) = P1( I-4, J-4)
C
AVQQ = AVQ * PP( I,J)
C
CALL CONAL ( AVQQ, II,JJ, P, NCC, NZWK, NNZ )
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 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
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
PD ( N ) = PDV
NCC( N ) = M
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ROUT ( C, CAA, ITB, LOP, CA, CN, DT, LPE, SEM,
1 ERR, NNP, * )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /MSZ/ NLAY( 3)
COMMON /XYZ/ XP( 1500), YP( 1500), ZP( 1500)
COMMON /CN2/ COEF, NMP, PHI( 20), DND, DTD
C
DIMENSION C ( NNP), ITB( NNP), CAA( NNP), CA( NNP), CN( NNP),
1 ARA( 12), CRA( 12)
C
C --------------------------------------
IF ( LOP.EQ.0) THEN
C
CALL CLAR( C, 0.D0, NNP)
C
CALL ANAS ( NNP, CAA, LOP, DT )
C
END IF
C --------------------------------------
C
C ----- Analytical Solution --------------
C
CALL ANAS ( NNP, CAA, LOP, DT )
C
C ----------------------------------------
C
C ------------------------------------------------
CAAM = 0.D0
DO 100 I = 1, NNP
C
IF ( CAA( I).GT.CAAM) CAAM = CAA( I)
100 CONTINUE
C ------------------------------------------------
IF ( CAAM.LT.1.D-8) CAAM = 1.D0
IF ( LOP.EQ.0.OR.LOP.EQ.LPE) GO TO 120
GO TO 140
C
120 CONTINUE
C
140 CONTINUE
C --------------------------------------
L = 0
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)/ CAAM
CN( L) = C( I) / CAAM
END IF
C
200 CONTINUE
C --------------------------------------
DO 250 I= 1, NNP
C
IF ( DABS( CN( I)).GE.1.5D0) GO TO 900
250 CONTINUE
C --------------------------------------
GO TO 270
C
900 CONTINUE
WRITE(6,*) ' *** Diverged ! *** CN( I) > = 1.5 ***'
STOP
C
270 CONTINUE
C --------------------------------------
DO 300 I = 1, 11
C
ARA( I+1) = CA( I)
CRA( I+1) = CN( I)
C
300 CONTINUE
C --------------------------------------
IF ( LOP.NE.LPE) GO TO 333
C
WRITE(6,*) ' '
WRITE(6,*) ' *** Ratio for Analytical Solution **'
WRITE(6,2000) ( ARA( I), I= 2, 12)
2000 FORMAT(2X,11F7.3)
WRITE(6,*) ' '
WRITE(6,*) ' *** Ratio for Numerical Solution ** '
WRITE(6,2000) ( CRA( I), I = 2, 12 )
C
333 CONTINUE
C -------------------------------------------------------
IF ( LOP.EQ.0) THEN
WRITE(10,2100) COEF, DT, LPE
2100 FORMAT(2F7.4,I3)
WRITE(10,2200) ( NLAY(I), I=1, 3 )
2200 FORMAT(3I3)
WRITE(10,2300) ( ARA( I), CRA( I), I = 2, 12 )
END IF
C -------------------------------------------------------
C
IF ( LOP.EQ. LPE/2) WRITE(10,2300) ( ARA( I), CRA( I), I = 2, 12)
IF ( LOP.EQ. LPE) WRITE(10,2300) ( ARA( I), CRA( I), I = 2, 12)
2300 FORMAT(10F8.4)
C
C -------------------------------------------------------
IF ( LOP.NE.0) THEN
C
CALL NMER ( NNP, C, CAA, LOP, SEM, DT )
C
ERR = SEM/ DFLOAT( LOP)
C
END IF
C ------------------------------------
IF ( LOP.LT.LPE) RETURN 1
C
WRITE(6,2400) SEM, ERR
2400 FORMAT(/' *** Final Error =',F9.3,' (%)',' Error / LOP =',F9.3)

RETURN
END
C **********************************************************************
C
SUBROUTINE ZARY ( M, N, A )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION A( M, N )
C
DO 100 I = 1, M
DO 100 J = 1, N
C
A( I,J) = 0.D0
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ZRVC ( M, IB )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION IB( M)
C
DO 100 I = 1, M
C
IB( I) = 0
100 CONTINUE
C
RETURN
END
C *******************************************************************