1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
|
*DECK TRICHD
SUBROUTINE TRICHD(IMPX,LX,LY,LZ,CYLIND,IELEM,L4,LL4F,LL4X,
1 LL4Y,LL4Z,MAT,VOL,XX,YY,ZZ,DD,KN,V,MUX,MUY,MUZ,IPBBX,IPBBY,IPBBZ,
2 BBX,BBY,BBZ)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Thomas-Raviart (dual) finite element unknown numbering for ADI
* solution in a 3D domain. Compute the storage info for ADI matrices
* in compressed diagonal storage mode. Compute the ADI permutation
* vectors. Compute the group-independent XB, YB and ZB matrices.
*
*Copyright:
* Copyright (C) 2002 Ecole Polytechnique de Montreal
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version
*
*Author(s): A. Hebert
*
*Parameters: input
* IMPX print parameter.
* LX number of elements along the X axis.
* LY number of elements along the Y axis.
* LZ number of elements along the Z axis.
* CYLIND cylindrical geometry flag (set with CYLIND=.true.).
* IELEM degree of the Lagrangian finite elements: =1 (linear);
* =2 (parabolic); =3 (cubic).
* L4 total number of unknown (variational coefficients) per
* energy group (order of system matrices).
* LL4F exact number of flux unknowns.
* LL4X exact number of X-directed current unknowns.
* LL4Y exact number of Y-directed current unknowns.
* LL4Z exact number of Z-directed current unknowns.
* MAT mixture index assigned to each element.
* VOL volume of each element.
* XX X-directed mesh spacings.
* YY Y-directed mesh spacings.
* ZZ Z-directed mesh spacings.
* DD used with cylindrical geometry.
* KN element-ordered unknown list.
* V finite element unit matrix.
*
*Parameters: output
* MUX X-directed compressed diagonal mode indices.
* MUY Y-directed compressed diagonal mode indices.
* MUZ Z-directed compressed diagonal mode indices.
* IPBBX X-directed perdue storage indices.
* IPBBY Y-directed perdue storage indices.
* IPBBZ Z-directed perdue storage indices.
* BBX X-directed flux-current matrices.
* BBY Y-directed flux-current matrices.
* BBZ Z-directed flux-current matrices.
*
*-----------------------------------------------------------------------
*
*----
* SUBROUTINE ARGUMENTS
*----
LOGICAL CYLIND
INTEGER IMPX,LX,LY,LZ,IELEM,L4,LL4F,LL4X,LL4Y,LL4Z,
1 MAT(LX*LY*LZ),KN(LX*LY*LZ*(1+6*IELEM**2)),MUX(L4),MUY(L4),
2 MUZ(L4),IPBBX(2*IELEM,LL4X),IPBBY(2*IELEM,LL4Y),
3 IPBBZ(2*IELEM,LL4Z)
REAL VOL(LX*LY*LZ),XX(LX*LY*LZ),YY(LX*LY*LZ),ZZ(LX*LY*LZ),
1 DD(LX*LY*LZ),V(IELEM+1,IELEM),BBX(2*IELEM,LL4X),
2 BBY(2*IELEM,LL4Y),BBZ(2*IELEM,LL4Z)
*
IF(IELEM.GT.4) CALL XABORT('TRICHD: 1 .LE. IELEM .LE. 3.')
IF(L4.NE.LL4F+LL4X+LL4Y+LL4Z) CALL XABORT('TRICHD: INVALID L4.')
*----
* COMPUTE THE X-ORIENTED SYSTEM BANDWIDTH VECTOR
*----
MUX(:L4)=1
IPBBX(:2*IELEM,:LL4X)=0
NUM1=0
DO 20 KEL=1,LX*LY*LZ
IF(MAT(KEL).EQ.0) GO TO 20
DO 12 K3=0,IELEM-1
DO 11 K2=0,IELEM-1
KN1=KN(NUM1+2+K3*IELEM+K2)
KN2=KN(NUM1+2+IELEM**2+K3*IELEM+K2)
INX1=ABS(KN1)-LL4F
INX2=ABS(KN2)-LL4F
IF((KN1.NE.0).AND.(KN2.NE.0)) THEN
MUX(INX2)=MAX(MUX(INX2),INX2-INX1+1)
MUX(INX1)=MAX(MUX(INX1),INX1-INX2+1)
ENDIF
DO 10 K1=0,IELEM-1
JND1=KN(NUM1+1)+(K3*IELEM+K2)*IELEM+K1
IF(KN1.NE.0) CALL TRINDX(JND1,IPBBX(1,INX1),2*IELEM)
IF(KN2.NE.0) CALL TRINDX(JND1,IPBBX(1,INX2),2*IELEM)
10 CONTINUE
11 CONTINUE
12 CONTINUE
NUM1=NUM1+1+6*IELEM**2
20 CONTINUE
*----
* COMPUTE THE Y-ORIENTED SYSTEM BANDWIDTH VECTOR
*----
MUY(:L4)=1
IPBBY(:2*IELEM,:LL4Y)=0
NUM1=0
DO 50 KEL=1,LX*LY*LZ
IF(MAT(KEL).EQ.0) GO TO 50
DO 42 K3=0,IELEM-1
DO 41 K1=0,IELEM-1
KN1=KN(NUM1+2+2*IELEM**2+K3*IELEM+K1)
KN2=KN(NUM1+2+3*IELEM**2+K3*IELEM+K1)
INY1=ABS(KN1)-LL4F-LL4X
INY2=ABS(KN2)-LL4F-LL4X
IF((KN1.NE.0).AND.(KN2.NE.0)) THEN
MUY(INY2)=MAX(MUY(INY2),INY2-INY1+1)
MUY(INY1)=MAX(MUY(INY1),INY1-INY2+1)
ENDIF
DO 40 K2=0,IELEM-1
JND1=KN(NUM1+1)+(K3*IELEM+K2)*IELEM+K1
IF(KN1.NE.0) CALL TRINDX(JND1,IPBBY(1,INY1),2*IELEM)
IF(KN2.NE.0) CALL TRINDX(JND1,IPBBY(1,INY2),2*IELEM)
40 CONTINUE
41 CONTINUE
42 CONTINUE
NUM1=NUM1+1+6*IELEM**2
50 CONTINUE
*----
* COMPUTE THE Z-ORIENTED SYSTEM BANDWIDTH VECTOR
*----
MUZ(:L4)=1
IPBBZ(:2*IELEM,:LL4Z)=0
NUM1=0
DO 70 KEL=1,LX*LY*LZ
IF(MAT(KEL).EQ.0) GO TO 70
DO 62 K2=0,IELEM-1
DO 61 K1=0,IELEM-1
KN1=KN(NUM1+2+4*IELEM**2+K2*IELEM+K1)
KN2=KN(NUM1+2+5*IELEM**2+K2*IELEM+K1)
INZ1=ABS(KN1)-LL4F-LL4X-LL4Y
INZ2=ABS(KN2)-LL4F-LL4X-LL4Y
IF((KN1.NE.0).AND.(KN2.NE.0)) THEN
MUZ(INZ2)=MAX(MUZ(INZ2),INZ2-INZ1+1)
MUZ(INZ1)=MAX(MUZ(INZ1),INZ1-INZ2+1)
ENDIF
DO 60 K3=0,IELEM-1
JND1=KN(NUM1+1)+(K3*IELEM+K2)*IELEM+K1
IF(KN1.NE.0) CALL TRINDX(JND1,IPBBZ(1,INZ1),2*IELEM)
IF(KN2.NE.0) CALL TRINDX(JND1,IPBBZ(1,INZ2),2*IELEM)
60 CONTINUE
61 CONTINUE
62 CONTINUE
NUM1=NUM1+1+6*IELEM**2
70 CONTINUE
*
MUXMAX=0
IIMAXX=0
DO 80 I=1,LL4X
MUXMAX=MAX(MUXMAX,MUX(I))
IIMAXX=IIMAXX+MUX(I)
MUX(I)=IIMAXX
80 CONTINUE
*
MUYMAX=0
IIMAXY=0
DO 90 I=1,LL4Y
MUYMAX=MAX(MUYMAX,MUY(I))
IIMAXY=IIMAXY+MUY(I)
MUY(I)=IIMAXY
90 CONTINUE
*
MUZMAX=0
IIMAXZ=0
DO 100 I=1,LL4Z
MUZMAX=MAX(MUZMAX,MUZ(I))
IIMAXZ=IIMAXZ+MUZ(I)
MUZ(I)=IIMAXZ
100 CONTINUE
IF(IMPX.GT.0) THEN
WRITE (6,600) MUXMAX,MUYMAX,MUZMAX
WRITE (6,610) IIMAXX,IIMAXY,IIMAXZ
ENDIF
*----
* COMPUTE THE FLUX-CURRENT COUPLING MATRICES XB, YB AND ZB.
*----
BBX(:2*IELEM,:LL4X)=0.0
BBY(:2*IELEM,:LL4Y)=0.0
BBZ(:2*IELEM,:LL4Z)=0.0
NUM1=0
DO 270 IE=1,LX*LY*LZ
L=MAT(IE)
IF(L.EQ.0) GO TO 270
VOL0=VOL(IE)
IF(VOL0.EQ.0.0) GO TO 260
DX=XX(IE)
DY=YY(IE)
DZ=ZZ(IE)
IF(CYLIND) THEN
DIN=1.0-0.5*DX/DD(IE)
DOT=1.0+0.5*DX/DD(IE)
ELSE
DIN=1.0
DOT=1.0
ENDIF
*
DO 152 K3=0,IELEM-1
DO 151 K2=0,IELEM-1
INX1=ABS(KN(NUM1+2+K3*IELEM+K2))-LL4F
INX2=ABS(KN(NUM1+2+IELEM**2+K3*IELEM+K2))-LL4F
DO 150 K1=0,IELEM-1
JND1=KN(NUM1+1)+(K3*IELEM+K2)*IELEM+K1
IF(KN(NUM1+2+K3*IELEM+K2).NE.0) THEN
KK=0
DO 110 I=1,2*IELEM
IF(IPBBX(I,INX1).EQ.JND1) THEN
KK=I
GO TO 120
ENDIF
110 CONTINUE
CALL XABORT('TRICHD: BUG1.')
120 SG=REAL(SIGN(1,KN(NUM1+2+K3*IELEM+K2)))
BBX(KK,INX1)=BBX(KK,INX1)+SG*(VOL0/DX)*DIN*V(1,K1+1)
ENDIF
IF(KN(NUM1+2+IELEM**2+K3*IELEM+K2).NE.0) THEN
KK=0
DO 130 I=1,2*IELEM
IF(IPBBX(I,INX2).EQ.JND1) THEN
KK=I
GO TO 140
ENDIF
130 CONTINUE
CALL XABORT('TRICHD: BUG2.')
140 SG=REAL(SIGN(1,KN(NUM1+2+IELEM**2+K3*IELEM+K2)))
BBX(KK,INX2)=BBX(KK,INX2)+SG*(VOL0/DX)*DOT*V(IELEM+1,K1+1)
ENDIF
150 CONTINUE
151 CONTINUE
152 CONTINUE
*
DO 202 K3=0,IELEM-1
DO 201 K1=0,IELEM-1
INY1=ABS(KN(NUM1+2+2*IELEM**2+K3*IELEM+K1))-LL4F-LL4X
INY2=ABS(KN(NUM1+2+3*IELEM**2+K3*IELEM+K1))-LL4F-LL4X
DO 200 K2=0,IELEM-1
JND1=KN(NUM1+1)+(K3*IELEM+K2)*IELEM+K1
IF(KN(NUM1+2+2*IELEM**2+K3*IELEM+K1).NE.0) THEN
KK=0
DO 160 I=1,2*IELEM
IF(IPBBY(I,INY1).EQ.JND1) THEN
KK=I
GO TO 170
ENDIF
160 CONTINUE
CALL XABORT('TRICHD: BUG3.')
170 SG=REAL(SIGN(1,KN(NUM1+2+2*IELEM**2+K3*IELEM+K1)))
BBY(KK,INY1)=BBY(KK,INY1)+SG*(VOL0/DY)*V(1,K2+1)
ENDIF
IF(KN(NUM1+2+3*IELEM**2+K3*IELEM+K1).NE.0) THEN
KK=0
DO 180 I=1,2*IELEM
IF(IPBBY(I,INY2).EQ.JND1) THEN
KK=I
GO TO 190
ENDIF
180 CONTINUE
CALL XABORT('TRICHD: BUG4.')
190 SG=REAL(SIGN(1,KN(NUM1+2+3*IELEM**2+K3*IELEM+K1)))
BBY(KK,INY2)=BBY(KK,INY2)+SG*(VOL0/DY)*V(IELEM+1,K2+1)
ENDIF
200 CONTINUE
201 CONTINUE
202 CONTINUE
*
DO 252 K2=0,IELEM-1
DO 251 K1=0,IELEM-1
INZ1=ABS(KN(NUM1+2+4*IELEM**2+K2*IELEM+K1))-LL4F-LL4X-LL4Y
INZ2=ABS(KN(NUM1+2+5*IELEM**2+K2*IELEM+K1))-LL4F-LL4X-LL4Y
DO 250 K3=0,IELEM-1
JND1=KN(NUM1+1)+(K3*IELEM+K2)*IELEM+K1
IF(KN(NUM1+2+4*IELEM**2+K2*IELEM+K1).NE.0) THEN
KK=0
DO 210 I=1,2*IELEM
IF(IPBBZ(I,INZ1).EQ.JND1) THEN
KK=I
GO TO 220
ENDIF
210 CONTINUE
CALL XABORT('TRICHD: BUG5.')
220 SG=REAL(SIGN(1,KN(NUM1+2+4*IELEM**2+K2*IELEM+K1)))
BBZ(KK,INZ1)=BBZ(KK,INZ1)+SG*(VOL0/DZ)*V(1,K3+1)
ENDIF
IF(KN(NUM1+2+5*IELEM**2+K2*IELEM+K1).NE.0) THEN
KK=0
DO 230 I=1,2*IELEM
IF(IPBBZ(I,INZ2).EQ.JND1) THEN
KK=I
GO TO 240
ENDIF
230 CONTINUE
CALL XABORT('TRICHD: BUG6.')
240 SG=REAL(SIGN(1,KN(NUM1+2+5*IELEM**2+K2*IELEM+K1)))
BBZ(KK,INZ2)=BBZ(KK,INZ2)+SG*(VOL0/DZ)*V(IELEM+1,K3+1)
ENDIF
250 CONTINUE
251 CONTINUE
252 CONTINUE
260 NUM1=NUM1+1+6*IELEM**2
270 CONTINUE
RETURN
*
600 FORMAT(/52H TRICHD: MAXIMUM BANDWIDTH FOR X-ORIENTED MATRICES =,
1 I4/27X,25HFOR Y-ORIENTED MATRICES =,I4/27X,16HFOR Z-ORIENTED M,
2 9HATRICES =,I4)
610 FORMAT(/40H TRICHD: LENGTH OF X-ORIENTED MATRICES =,I10/16X,
1 24HOF Y-ORIENTED MATRICES =,I10/16X,24HOF Z-ORIENTED MATRICES =,
2 I10)
END
|
