diff options
Diffstat (limited to 'Trivac/src/BIVA03.f')
| -rwxr-xr-x | Trivac/src/BIVA03.f | 176 |
1 files changed, 176 insertions, 0 deletions
diff --git a/Trivac/src/BIVA03.f b/Trivac/src/BIVA03.f new file mode 100755 index 0000000..f9c3f71 --- /dev/null +++ b/Trivac/src/BIVA03.f @@ -0,0 +1,176 @@ +*DECK BIVA03 + SUBROUTINE BIVA03(ITY,MAXKN,MAXQF,SGD,NREG,LL4,ISPLH,NELEM,NBMIX, + 1 IIMAX,SIDE,MAT,KN,QFR,VOL,MU,R,RH,QH,RT,QT,SYS) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Assembly of a within-group (leakage and removal) or out-of-group +* system matrix in mesh-corner finite-difference diffusion +* approximation (hexagonal geometry). +* +*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 +* ITY type of assembly: =0: leakage-removal matrix assembly; +* =1: cross section matrix assembly. +* MAXKN dimension of array KN. +* MAXQF dimension of array QFR. +* SGD nuclear properties. SGD(:,1) and SGD(:,2) are diffusion +* coefficients. SGD(:,3) are removal macroscopic cross sections. +* NREG number of hexagons in BIVAC. +* LL4 order of the matrix SYS. +* ISPLH hexagonal geometry flag: +* =1: hexagonal elements; >1: triangular elements. +* NELEM number of finite elements (hexagons or triangles) excluding +* the virtual elements. +* NBMIX number of macro-mixtures. +* IIMAX allocated dimension of array SYS. +* SIDE side of the hexagons. +* MAT mixture index per hexagon. +* KN element-ordered unknown list. +* QFR element-ordered information. +* VOL volume of the hexagons. +* MU indices used with the compressed diagonal storage mode matrix +* SYS. +* R unit matrix. +* RH unit matrix. +* QH unit matrix. +* RT unit matrix. +* QT unit matrix. +* +*Parameters: output +* SYS system matrix. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER ITY,MAXKN,MAXQF,NREG,LL4,ISPLH,NELEM,NBMIX,IIMAX, + 1 MAT(NREG),KN(MAXKN),MU(LL4) + REAL SGD(NBMIX,3),SIDE,QFR(MAXQF),VOL(NREG),R(2,2),RH(6,6), + 1 QH(6,6),RT(3,3),QT(3,3),SYS(IIMAX) +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION RR,RRH,QQH + INTEGER ISR(6,2),ISRH(6,2),ISRT(3,2) + REAL RH2(6,6),QH2(6,6) + DATA ISRH/2,1,4,5,6,3,1,4,5,6,3,2/ + DATA ISRT/1,2,3,2,3,1/ +*---- +* RECOVER THE HEXAGONAL MASS (RH2) AND STIFFNESS (QH2) MATRICES. +*---- + IF(ISPLH.EQ.1) THEN +* HEXAGONAL BASIS. + LH=6 + DO 15 I=1,6 + DO 10 J=1,2 + ISR(I,J)=ISRH(I,J) + 10 CONTINUE + 15 CONTINUE + DO 25 I=1,6 + DO 20 J=1,6 + RH2(I,J)=RH(I,J) + QH2(I,J)=QH(I,J) + 20 CONTINUE + 25 CONTINUE + CONST=1.5*SQRT(3.0) + CONSB=2.0*SQRT(3.0)/3.0 + AA=SIDE + ELSE +* TRIANGULAR BASIS. + LH=3 + DO 35 I=1,3 + DO 30 J=1,2 + ISR(I,J)=ISRT(I,J) + 30 CONTINUE + 35 CONTINUE + DO 45 I=1,3 + DO 40 J=1,3 + RH2(I,J)=RT(I,J) + QH2(I,J)=QT(I,J) + 40 CONTINUE + 45 CONTINUE + CONST=0.25*SQRT(3.0) + CONSB=2.0*SQRT(3.0) + AA=SIDE/REAL(ISPLH-1) + ENDIF +*---- +* ASSEMBLY OF A SYSTEM MATRIX. +*---- + IF(ITY.EQ.0) THEN +* LEAKAGE-REMOVAL SYSTEM MATRIX ASSEMBLY. + NUM1=0 + DO 105 K=1,NELEM + KHEX=KN(NUM1+LH+1) + IF(VOL(KHEX).EQ.0.0) GO TO 100 + L=MAT(KHEX) + VOL0=QFR(NUM1+LH+1) + DO 60 I=1,LH + IND1=KN(NUM1+I) + IF(IND1.EQ.0) GO TO 60 + KEY1=MU(IND1)-IND1 + DO 50 J=1,LH + IND2=KN(NUM1+J) + IF((IND2.EQ.0).OR.(IND2.GT.IND1)) GO TO 50 + QQH=QH2(I,J)/(CONST*AA*AA) + RRH=RH2(I,J)/CONST + IF((QQH.EQ.0.0).AND.(RRH.EQ.0.0)) GO TO 50 + KEY=KEY1+IND2 + SYS(KEY)=SYS(KEY)+REAL(QQH*SGD(L,1)+RRH*SGD(L,3))*VOL0 + 50 CONTINUE + 60 CONTINUE + DO 90 IC=1,LH + QFR1=QFR(NUM1+IC) + IF(QFR1.EQ.0.0) GO TO 90 + DO 80 I1=1,2 + IND1=KN(NUM1+ISR(IC,I1)) + IF(IND1.EQ.0) GO TO 80 + KEY1=MU(IND1)-IND1 + DO 70 J1=1,2 + IND2=KN(NUM1+ISR(IC,J1)) + IF((IND2.EQ.0).OR.(IND2.GT.IND1)) GO TO 70 + RR=R(I1,J1) + IF(RR.EQ.0.0) GO TO 70 + KEY=KEY1+IND2 + SYS(KEY)=SYS(KEY)+REAL(RR)*QFR1 + 70 CONTINUE + 80 CONTINUE + 90 CONTINUE + 100 NUM1=NUM1+LH+1 + 105 CONTINUE + ELSE +* CROSS SECTION SYSTEM MATRIX ASSEMBLY + NUM1=0 + DO 135 K=1,NELEM + KHEX=KN(NUM1+LH+1) + IF(VOL(KHEX).EQ.0.0) GO TO 130 + L=MAT(KHEX) + VOL0=QFR(NUM1+LH+1) + DO 120 I=1,LH + IND1=KN(NUM1+I) + IF(IND1.EQ.0) GO TO 120 + KEY1=MU(IND1)-IND1 + DO 110 J=1,LH + IND2=KN(NUM1+J) + IF((IND2.EQ.0).OR.(IND2.GT.IND1)) GO TO 110 + RRH=RH2(I,J)/CONST + IF(RRH.EQ.0.0) GO TO 110 + KEY=KEY1+IND2 + SYS(KEY)=SYS(KEY)+REAL(RRH)*SGD(L,1)*VOL0 + 110 CONTINUE + 120 CONTINUE + 130 NUM1=NUM1+LH+1 + 135 CONTINUE + ENDIF + RETURN + END |