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
|
*DECK BIVA05
SUBROUTINE BIVA05(ITY,SGD,IELEM,NBLOS,LL4,NBMIX,IIMAX,SIDE,MAT,
1 IPERT,KN,QFR,MU,LC,R,V,H,SYS)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Assembly of a within-group (leakage and removal) or out-of-group
* system matrix in a Thomas-Raviart-Schneider (dual) finite element
* diffusion approximation (hexagonal geometry).
*
*Copyright:
* Copyright (C) 2006 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
* ITY type of assembly: =0: leakage-removal matrix assembly;
* =1: cross section matrix assembly.
* SGD nuclear properties. SGD(:,1) and SGD(:,2) are diffusion
* coefficients. SGD(:,3) are removal macroscopic cross sections.
* IELEM degree of the Lagrangian finite elements: =1 (linear);
* =2 (parabolic); =3 (cubic); =4 (quartic).
* NBLOS number of lozenges per direction, taking into account
* mesh-splitting.
* LL4 number of unknowns per group in BIVAC.
* NBMIX number of macro-mixtures.
* IIMAX allocated dimension of array SYS.
* SIDE side of the hexagons.
* MAT mixture index per lozenge.
* IPERT mixture permutation index.
* KN element-ordered unknown list.
* QFR element-ordered boundary conditions.
* MU indices used with compressed diagonal storage mode matrix SYS.
* LC order of the unit matrices.
* R Cartesian mass matrix.
* V nodal coupling matrix.
* H Piolat (hexagonal) coupling matrix.
*
*Parameters: output
* SYS system matrix.
*
*-----------------------------------------------------------------------
*
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER ITY,IELEM,NBLOS,LL4,NBMIX,IIMAX,MAT(3,NBLOS),IPERT(NBLOS),
1 KN(NBLOS,4+6*IELEM*(IELEM+1)),MU(LL4),LC
REAL SGD(NBMIX,3),SIDE,QFR(NBLOS,6),R(LC,LC),V(LC,LC-1),
1 H(LC,LC-1),SYS(IIMAX)
*----
* LOCAL VARIABLES
*----
PARAMETER(MAXIEL=3)
DOUBLE PRECISION CTRAN(MAXIEL*(MAXIEL+1),MAXIEL*(MAXIEL+1))
*----
* ASSEMBLY OF A SYSTEM MATRIX.
*----
TTTT=0.5*SQRT(3.0)*SIDE*SIDE
IF(IELEM.GT.MAXIEL) CALL XABORT('BIVA05: MAXIEL OVERFLOW.')
IF(ITY.EQ.0) THEN
* COMPUTE THE TRANVERSE COUPLING PIOLAT UNIT MATRIX
CTRAN(:MAXIEL*(MAXIEL+1),:MAXIEL*(MAXIEL+1))=0.0D0
CNORM=SIDE*SIDE/SQRT(3.0)
I=0
DO 22 JS=1,IELEM
DO 21 JT=1,IELEM+1
J=0
I=I+1
SSS=1.0
DO 20 IT=1,IELEM
DO 10 IS=1,IELEM+1
J=J+1
CTRAN(I,J)=SSS*CNORM*H(IS,JS)*H(JT,IT)
10 CONTINUE
SSS=-SSS
20 CONTINUE
21 CONTINUE
22 CONTINUE
*
* LEAKAGE-REMOVAL SYSTEM MATRIX ASSEMBLY
NELEM=IELEM*(IELEM+1)
COEF=2.0*SIDE*SIDE/SQRT(3.0)
NUM=0
DO 70 KEL=1,NBLOS
IF(IPERT(KEL).EQ.0) GO TO 70
IBM=MAT(1,IPERT(KEL))
IF(IBM.EQ.0) GO TO 70
NUM=NUM+1
DINV=1.0/SGD(IBM,1)
SIG=SGD(IBM,3)
DO 43 K4=0,1
DO 42 K3=0,IELEM-1
DO 41 K2=1,IELEM+1
KNW1=KN(NUM,4+K4*NELEM+K3*(IELEM+1)+K2)
KNX1=KN(NUM,4+(K4+2)*NELEM+K3*(IELEM+1)+K2)
KNY1=KN(NUM,4+(K4+4)*NELEM+K3*(IELEM+1)+K2)
INW1=ABS(KNW1)
INX1=ABS(KNX1)
INY1=ABS(KNY1)
DO 30 K1=1,IELEM+1
KNW2=KN(NUM,4+K4*NELEM+K3*(IELEM+1)+K1)
KNX2=KN(NUM,4+(K4+2)*NELEM+K3*(IELEM+1)+K1)
KNY2=KN(NUM,4+(K4+4)*NELEM+K3*(IELEM+1)+K1)
INW2=ABS(KNW2)
INX2=ABS(KNX2)
INY2=ABS(KNY2)
IF((KNW2.NE.0).AND.(KNW1.NE.0).AND.(INW1.GE.INW2)) THEN
L=MU(INW1)-INW1+INW2
SG=REAL(SIGN(1,KNW1)*SIGN(1,KNW2))
SYS(L)=SYS(L)-SG*COEF*DINV*R(K2,K1)
IF(INW1.EQ.INW2) THEN
IF((K1.EQ.1).AND.(K4.EQ.0)) SYS(L)=SYS(L)-QFR(NUM,1)
IF((K1.EQ.IELEM+1).AND.(K4.EQ.1)) SYS(L)=SYS(L)-QFR(NUM,2)
ENDIF
ENDIF
IF((KNX2.NE.0).AND.(KNX1.NE.0).AND.(INX1.GE.INX2)) THEN
L=MU(INX1)-INX1+INX2
SG=REAL(SIGN(1,KNX1)*SIGN(1,KNX2))
SYS(L)=SYS(L)-SG*COEF*DINV*R(K2,K1)
IF(INX1.EQ.INX2) THEN
IF((K1.EQ.1).AND.(K4.EQ.0)) SYS(L)=SYS(L)-QFR(NUM,3)
IF((K1.EQ.IELEM+1).AND.(K4.EQ.1)) SYS(L)=SYS(L)-QFR(NUM,4)
ENDIF
ENDIF
IF((KNY2.NE.0).AND.(KNY1.NE.0).AND.(INY1.GE.INY2)) THEN
L=MU(INY1)-INY1+INY2
SG=REAL(SIGN(1,KNY1)*SIGN(1,KNY2))
SYS(L)=SYS(L)-SG*COEF*DINV*R(K2,K1)
IF(INY1.EQ.INY2) THEN
IF((K1.EQ.1).AND.(K4.EQ.0)) SYS(L)=SYS(L)-QFR(NUM,5)
IF((K1.EQ.IELEM+1).AND.(K4.EQ.1)) SYS(L)=SYS(L)-QFR(NUM,6)
ENDIF
ENDIF
30 CONTINUE
DO 40 K1=0,IELEM-1
IF(V(K2,K1+1).EQ.0.0) GO TO 40
IF(K4.EQ.0) THEN
SSS=(-1.0)**K1
JND1=KN(NUM,1)+K3*IELEM+K1
JND2=KN(NUM,2)+K3*IELEM+K1
JND3=KN(NUM,3)+K3*IELEM+K1
ELSE
SSS=1.0
JND1=KN(NUM,2)+K1*IELEM+K3
JND2=KN(NUM,3)+K1*IELEM+K3
JND3=KN(NUM,4)+K1*IELEM+K3
ENDIF
IF(KNW1.NE.0) THEN
L=MU(JND1)-JND1+INW1
IF(JND1.LT.INW1) L=MU(INW1)-INW1+JND1
SG=REAL(SIGN(1,KNW1))
SYS(L)=SYS(L)+SG*SSS*SIDE*V(K2,K1+1)
ENDIF
IF(KNX1.NE.0) THEN
L=MU(JND2)-JND2+INX1
IF(JND2.LT.INX1) L=MU(INX1)-INX1+JND2
SG=REAL(SIGN(1,KNX1))
SYS(L)=SYS(L)+SG*SSS*SIDE*V(K2,K1+1)
ENDIF
IF(KNY1.NE.0) THEN
L=MU(JND3)-JND3+INY1
IF(JND3.LT.INY1) L=MU(INY1)-INY1+JND3
SG=REAL(SIGN(1,KNY1))
SYS(L)=SYS(L)+SG*SSS*SIDE*V(K2,K1+1)
ENDIF
40 CONTINUE
41 CONTINUE
42 CONTINUE
43 CONTINUE
ITRS=0
DO I=1,NBLOS
IF(KN(I,1).EQ.KN(NUM,4)) THEN
ITRS=I
GO TO 45
ENDIF
ENDDO
CALL XABORT('BIVA05: ITRS FAILURE.')
45 DO 55 I=1,NELEM
KNW1=KN(ITRS,4+I)
KNX1=KN(NUM,4+2*NELEM+I)
KNY1=KN(NUM,4+4*NELEM+I)
INW1=ABS(KNW1)
INX1=ABS(KNX1)
INY1=ABS(KNY1)
DO 50 J=1,NELEM
KNW2=KN(NUM,4+NELEM+J)
KNX2=KN(NUM,4+3*NELEM+J)
KNY2=KN(NUM,4+5*NELEM+J)
INW2=ABS(KNW2)
INX2=ABS(KNX2)
INY2=ABS(KNY2)
IF((KNY2.NE.0).AND.(KNW1.NE.0).AND.(INW1.LT.INY2)) THEN
L=MU(INY2)-INY2+INW1
SG=REAL(SIGN(1,KNW1)*SIGN(1,KNY2))
SYS(L)=SYS(L)-SG*DINV*REAL(CTRAN(I,J)) ! y w
ELSE IF((KNY2.NE.0).AND.(KNW1.NE.0).AND.(INW1.GT.INY2)) THEN
L=MU(INW1)-INW1+INY2
SG=REAL(SIGN(1,KNW1)*SIGN(1,KNY2))
SYS(L)=SYS(L)-SG*DINV*REAL(CTRAN(I,J)) ! w y
ENDIF
IF((KNW2.NE.0).AND.(KNX1.NE.0).AND.(INW2.LT.INX1)) THEN
L=MU(INX1)-INX1+INW2
SG=REAL(SIGN(1,KNX1)*SIGN(1,KNW2))
SYS(L)=SYS(L)-SG*DINV*REAL(CTRAN(I,J)) ! x w
ELSE IF((KNW2.NE.0).AND.(KNX1.NE.0).AND.(INW2.GT.INX1)) THEN
L=MU(INW2)-INW2+INX1
SG=REAL(SIGN(1,KNX1)*SIGN(1,KNW2))
SYS(L)=SYS(L)-SG*DINV*REAL(CTRAN(I,J)) ! w x
ENDIF
IF((KNX2.NE.0).AND.(KNY1.NE.0).AND.(INX2.LT.INY1)) THEN
L=MU(INY1)-INY1+INX2
SG=REAL(SIGN(1,KNY1)*SIGN(1,KNX2))
SYS(L)=SYS(L)-SG*DINV*REAL(CTRAN(I,J)) ! y x
ELSE IF((KNX2.NE.0).AND.(KNY1.NE.0).AND.(INX2.GT.INY1)) THEN
L=MU(INX2)-INX2+INY1
SG=REAL(SIGN(1,KNY1)*SIGN(1,KNX2))
SYS(L)=SYS(L)-SG*DINV*REAL(CTRAN(I,J)) ! x y
ENDIF
50 CONTINUE
55 CONTINUE
DO 65 K2=0,IELEM-1
DO 60 K1=0,IELEM-1
JND1=KN(NUM,1)+K2*IELEM+K1
JND2=KN(NUM,2)+K2*IELEM+K1
JND3=KN(NUM,3)+K2*IELEM+K1
L=MU(JND1)
SYS(L)=SYS(L)+TTTT*SIG
L=MU(JND2)
SYS(L)=SYS(L)+TTTT*SIG
L=MU(JND3)
SYS(L)=SYS(L)+TTTT*SIG
60 CONTINUE
65 CONTINUE
70 CONTINUE
ELSE
* CROSS SECTION SYSTEM MATRIX ASSEMBLY
NUM=0
DO 90 KEL=1,NBLOS
IF(IPERT(KEL).EQ.0) GO TO 90
IBM=MAT(1,IPERT(KEL))
IF(IBM.EQ.0) GO TO 90
NUM=NUM+1
SIG=SGD(IBM,1)
DO 85 K2=0,IELEM-1
DO 80 K1=0,IELEM-1
JND1=KN(NUM,1)+K2*IELEM+K1
JND2=KN(NUM,2)+K2*IELEM+K1
JND3=KN(NUM,3)+K2*IELEM+K1
L=MU(JND1)
SYS(L)=SYS(L)+TTTT*SIG
L=MU(JND2)
SYS(L)=SYS(L)+TTTT*SIG
L=MU(JND3)
SYS(L)=SYS(L)+TTTT*SIG
80 CONTINUE
85 CONTINUE
90 CONTINUE
ENDIF
RETURN
END
|
