Initial commit from Polytechnique Montreal

1 files changed, 210 insertions, 0 deletions

diff --git a/Trivac/src/BIVA02.f b/Trivac/src/BIVA02.f
new file mode 100755
index 0000000..da999ea
--- /dev/null
+++ b/Trivac/src/BIVA02.f
@@ -0,0 +1,210 @@
+*DECK BIVA02
+      SUBROUTINE BIVA02(ITY,SGD,CYLIND,IELEM,ICOL,NREG,LL4,NBMIX,IIMAX,
+     1 XX,YY,DD,MAT,KN,QFR,VOL,MU,LC,R,V,SYS)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Assembly of a within-group (leakage and removal) or out-of-group
+* system matrix in mixed-dual finite element diffusion approximation
+* (Cartesian 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.
+* SGD     nuclear properties. SGD(:,1) and SGD(:,2) are diffusion
+*         coefficients. SGD(:,3) are removal macroscopic cross sections.
+* CYLIND  cylinderization flag (=.true. for cylindrical geometry).
+* IELEM   degree of the Lagrangian finite elements: =1 (linear);
+*         =2 (parabolic); =3 (cubic); =4 (quartic).
+* ICOL    type of quadrature: =1 (analytical integration);
+*         =2 (Gauss-Lobatto); =3 (Gauss-Legendre).
+* NREG    number of elements in BIVAC.
+* LL4     number of unknowns per group in BIVAC.
+* NBMIX   number of macro-mixtures.
+* IIMAX   allocated dimension of array SYS.
+* XX      X-directed mesh spacings.
+* YY      Y-directed mesh spacings.
+* DD      value used used with a cylindrical geometry.
+* MAT     mixture index per region.
+* KN      element-ordered unknown list.
+* QFR     element-ordered boundary conditions.
+* VOL     volume of regions.
+* MU      indices used with compressed diagonal storage mode matrix SYS.
+* LC      number of polynomials in a complete 1-D basis.
+* R       cartesian mass matrix.
+* V       nodal coupling matrix.
+*
+*Parameters: output
+* SYS     system matrix.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER ITY,IELEM,ICOL,NREG,LL4,NBMIX,IIMAX,MAT(NREG),KN(5*NREG),
+     1 MU(LL4),LC
+      REAL SGD(NBMIX,3),XX(NREG),YY(NREG),DD(NREG),QFR(4*NREG),
+     1 VOL(NREG),R(LC,LC),V(LC,LC-1),SYS(IIMAX)
+      LOGICAL CYLIND
+*----
+*  LOCAL VARIABLES
+*----
+      REAL QQ(5,5)
+*
+      IF((CYLIND).AND.((IELEM.GT.1).OR.(ICOL.NE.2)))
+     1 CALL XABORT('BIVA02: TYPE OF DISCRETIZATION NOT IMPLEMENTED.')
+*----
+*  ASSEMBLY OF A SYSTEM MATRIX.
+*----
+      IF(ITY.EQ.0) THEN
+*        LEAKAGE-REMOVAL SYSTEM MATRIX ASSEMBLY.
+         DO 12 I0=1,IELEM
+         DO 11 J0=1,IELEM
+         QQ(I0,J0)=0.0
+         DO 10 K0=2,IELEM
+         QQ(I0,J0)=QQ(I0,J0)+V(K0,I0)*V(K0,J0)/R(K0,K0)
+   10    CONTINUE
+   11    CONTINUE
+   12    CONTINUE
+         NUM1=0
+         NUM2=0
+         DO 80 K=1,NREG
+         L=MAT(K)
+         IF(L.EQ.0) GO TO 80
+         VOL0=VOL(K)
+         IF(VOL0.EQ.0.0) GO TO 70
+         DX=XX(K)
+         DY=YY(K)
+         IF(CYLIND) THEN
+            DIN=1.0-0.5*DX/DD(K)
+            DOT=1.0+0.5*DX/DD(K)
+         ELSE
+            DIN=1.0
+            DOT=1.0
+         ENDIF
+*
+         DO 60 I0=1,IELEM
+         INX1=ABS(KN(NUM1+2))+I0-1
+         INX2=ABS(KN(NUM1+3))+I0-1
+         INY1=ABS(KN(NUM1+4))+I0-1
+         INY2=ABS(KN(NUM1+5))+I0-1
+         DO 50 J0=1,IELEM
+         JND1=KN(NUM1+1)+(I0-1)*IELEM+J0-1
+         KEY=MU(JND1)
+         SYS(KEY)=SYS(KEY)+VOL0*SGD(L,3)
+         DO 20 K0=1,J0
+         IF(QQ(J0,K0).EQ.0.0) GO TO 20
+         KND1=KN(NUM1+1)+(I0-1)*IELEM+K0-1
+         KEY=MU(JND1)-JND1+KND1
+         SYS(KEY)=SYS(KEY)+VOL0*QQ(J0,K0)*SGD(L,1)/(DX*DX)
+   20    CONTINUE
+         IF(KN(NUM1+2).NE.0) THEN
+            IF(JND1.GT.INX1) KEY=MU(JND1)-JND1+INX1
+            IF(JND1.LT.INX1) KEY=MU(INX1)-INX1+JND1
+            SG=REAL(SIGN(1,KN(NUM1+2)))
+            SYS(KEY)=SYS(KEY)+SG*(VOL0/DX)*DIN*V(1,J0)
+         ENDIF
+         IF(KN(NUM1+3).NE.0) THEN
+            IF(INX2.GT.JND1) KEY=MU(INX2)-INX2+JND1
+            IF(INX2.LT.JND1) KEY=MU(JND1)-JND1+INX2
+            SG=REAL(SIGN(1,KN(NUM1+3)))
+            SYS(KEY)=SYS(KEY)+SG*(VOL0/DX)*DOT*V(IELEM+1,J0)
+         ENDIF
+         JND1=KN(NUM1+1)+(J0-1)*IELEM+I0-1
+         DO 30 K0=1,J0
+         IF(QQ(J0,K0).EQ.0.0) GO TO 30
+         KND1=KN(NUM1+1)+(K0-1)*IELEM+I0-1
+         KEY=MU(JND1)-JND1+KND1
+         SYS(KEY)=SYS(KEY)+VOL0*QQ(J0,K0)*SGD(L,2)/(DY*DY)
+   30    CONTINUE
+         IF(KN(NUM1+4).NE.0) THEN
+            IF(JND1.GT.INY1) KEY=MU(JND1)-JND1+INY1
+            IF(JND1.LT.INY1) KEY=MU(INY1)-INY1+JND1
+            SG=REAL(SIGN(1,KN(NUM1+4)))
+            SYS(KEY)=SYS(KEY)+SG*(VOL0/DY)*V(1,J0)
+         ENDIF
+         IF(KN(NUM1+5).NE.0) THEN
+            IF(INY2.GT.JND1) KEY=MU(INY2)-INY2+JND1
+            IF(INY2.LT.JND1) KEY=MU(JND1)-JND1+INY2
+            SG=REAL(SIGN(1,KN(NUM1+5)))
+            SYS(KEY)=SYS(KEY)+SG*(VOL0/DY)*V(IELEM+1,J0)
+         ENDIF
+   50    CONTINUE
+         IF(KN(NUM1+2).NE.0) THEN
+            KEY=MU(INX1)
+            SYS(KEY)=SYS(KEY)-DIN*(VOL0*R(1,1)/SGD(L,1)+QFR(NUM2+1))
+         ENDIF
+         IF(KN(NUM1+3).NE.0) THEN
+            KEY=MU(INX2)
+            SYS(KEY)=SYS(KEY)-DOT*(VOL0*R(IELEM+1,IELEM+1)/SGD(L,1)
+     1      +QFR(NUM2+2))
+         ENDIF
+         IF(KN(NUM1+4).NE.0) THEN
+            KEY=MU(INY1)
+            SYS(KEY)=SYS(KEY)-VOL0*R(1,1)/SGD(L,2)-QFR(NUM2+3)
+         ENDIF
+         IF(KN(NUM1+5).NE.0) THEN
+            KEY=MU(INY2)
+            SYS(KEY)=SYS(KEY)-VOL0*R(IELEM+1,IELEM+1)/SGD(L,2)
+     1      -QFR(NUM2+4)
+         ENDIF
+         IF(ICOL.NE.2) THEN
+            IF((KN(NUM1+2).NE.0).AND.(KN(NUM1+3).NE.0)) THEN
+               IF(INX2.GT.INX1) KEY=MU(INX2)-INX2+INX1
+               IF(INX2.LE.INX1) KEY=MU(INX1)-INX1+INX2
+               SG=REAL(SIGN(1,KN(NUM1+2))*SIGN(1,KN(NUM1+3)))
+               IF(INX1.EQ.INX2) SG=2.0*SG
+               SYS(KEY)=SYS(KEY)-SG*VOL0*R(IELEM+1,1)/SGD(L,1)
+            ENDIF
+            IF((KN(NUM1+4).NE.0).AND.(KN(NUM1+5).NE.0)) THEN
+               IF(INY2.GT.INY1) KEY=MU(INY2)-INY2+INY1
+               IF(INY2.LE.INY1) KEY=MU(INY1)-INY1+INY2
+               SG=REAL(SIGN(1,KN(NUM1+4))*SIGN(1,KN(NUM1+5)))
+               IF(INY1.EQ.INY2) SG=2.0*SG
+               SYS(KEY)=SYS(KEY)-SG*VOL0*R(IELEM+1,1)/SGD(L,2)
+            ENDIF
+         ENDIF
+   60   CONTINUE
+   70   NUM1=NUM1+5
+        NUM2=NUM2+4
+   80   CONTINUE
+      ELSE
+*        CROSS SECTION SYSTEM MATRIX ASSEMBLY. COMPONENTS WITH 1E-10
+*        FACTORS ARE INTRODUCED TO MAKE THE MATRIX INVERTIBLE.
+         NUM1=0
+         DO 110 K=1,NREG
+         L=MAT(K)
+         IF(L.EQ.0) GO TO 110
+         VOL0=VOL(K)
+         IF(VOL0.EQ.0.0) GO TO 100
+         DO 95 I0=1,IELEM
+         INX1=ABS(KN(NUM1+2))+I0-1
+         INX2=ABS(KN(NUM1+3))+I0-1
+         INY1=ABS(KN(NUM1+4))+I0-1
+         INY2=ABS(KN(NUM1+5))+I0-1
+         IF(KN(NUM1+2).NE.0) SYS(MU(INX1))=SYS(MU(INX1))+1.0E-30
+         IF(KN(NUM1+3).NE.0) SYS(MU(INX2))=SYS(MU(INX2))+1.0E-30
+         IF(KN(NUM1+4).NE.0) SYS(MU(INY1))=SYS(MU(INY1))+1.0E-30
+         IF(KN(NUM1+5).NE.0) SYS(MU(INY2))=SYS(MU(INY2))+1.0E-30
+         DO 90 J0=1,IELEM
+         JND1=KN(NUM1+1)+(I0-1)*IELEM+J0-1
+         KEY=MU(JND1)
+         SYS(KEY)=SYS(KEY)+VOL0*SGD(L,1)
+   90    CONTINUE
+   95    CONTINUE
+  100    NUM1=NUM1+5
+  110    CONTINUE
+      ENDIF
+      RETURN
+      END