Initial commit from Polytechnique Montreal

1 files changed, 405 insertions, 0 deletions

diff --git a/Trivac/src/BIVDKN.f b/Trivac/src/BIVDKN.f
new file mode 100755
index 0000000..ac927d4
--- /dev/null
+++ b/Trivac/src/BIVDKN.f
@@ -0,0 +1,405 @@
+*DECK BIVDKN
+      SUBROUTINE BIVDKN (MAXEV,IMPX,LX,LY,CYLIND,IELEM,ICOL,L4,NCODE,
+     1 ICODE,ZCODE,MAT,VOL,XXX,YYY,XX,YY,DD,KN,QFR,IQFR,BFR,IDL,MU)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Numbering corresponding to a mixed-dual formulation of the finite-
+* element discretization in a 2-D geometry. This version does not
+* support diagonal symmetries.
+*
+*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
+* MAXEV   allocated storage for vector MU.
+* IMPX    print parameter.
+* LX      number of elements along the X axis.
+* LY      number of elements along the Y axis.
+* 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).
+* NCODE   type of boundary condition applied on each side
+*         (i=1: X-  i=2: X+  i=3: Y-  i=4: Y+):
+*         NCODE(I)=1: VOID;   NCODE(I)=2: REFL;   NCODE(I)=4: TRAN;
+*         NCODE(I)=5: SYME;   NCODE(I)=7: ZERO.
+* ICODE   physical albedo index on each side of the domain.
+* ZCODE   ZCODE(I) is the albedo corresponding to boundary condition
+*         'VOID' on each side (ZCODE(I)=0.0 by default).
+* MAT     mixture index assigned to each element.
+* XXX     Cartesian coordinates along the X axis.
+* YYY     Cartesian coordinates along the Y axis.
+*
+*Parameters: output
+* L4      total number of unknown (variational coefficients) per
+*         energy group (order of system matrices).
+* VOL     volume of each element.
+* XX      X-directed mesh spacings.
+* YY      Y-directed mesh spacings.
+* DD      value used with a cylindrical geometry.
+* KN      element-ordered unknown list.
+* QFR     element-ordered boundary conditions.
+* IQFR    element-ordered physical albedo indices.
+* BFR     element-ordered surface fractions.
+* IDL     position of integrated fluxes into unknown vector.
+* MU      compressed storage mode indices.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER MAXEV,IMPX,LX,LY,IELEM,ICOL,L4,NCODE(4),ICODE(4),
+     1 MAT(LX*LY),KN(5*LX*LY),IQFR(4*LX*LY),IDL(LX*LY),MU(MAXEV)
+      REAL ZCODE(4),VOL(LX*LY),XXX(LX+1),YYY(LY+1),XX(LX*LY),YY(LX*LY),
+     1 DD(LX*LY),QFR(4*LX*LY),BFR(4*LX*LY)
+      LOGICAL CYLIND
+*----
+*  LOCAL VARIABLES
+*----
+      LOGICAL COND,LOG1,LOG2,LOG3,LOG4
+      CHARACTER TEXT8*8
+      REAL ZALB(4)
+      INTEGER, DIMENSION(:), ALLOCATABLE :: IP
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(IP(MAXEV))
+*----
+*  IDENTIFICATION OF THE GEOMETRY. MAIN LOOP OVER THE ELEMENTS.
+*----
+      DO 10 I=1,4
+      IF(ZCODE(I).NE.1.0) THEN
+         ZALB(I)=2.0*(1.0+ZCODE(I))/(1.0-ZCODE(I))
+      ELSE
+         ZALB(I)=1.0E20
+      ENDIF
+   10 CONTINUE
+      IF(IMPX.GT.0) WRITE(6,700) LX,LY
+      KN(:5*LX*LY)=0
+      SURFTOT=0.0
+      NUM1=0
+      NUM2=0
+      KEL=0
+      DO 151 K1=1,LY
+      DO 150 K2=1,LX
+      KEL=KEL+1
+      XX(KEL)=0.0
+      YY(KEL)=0.0
+      VOL(KEL)=0.0
+      IF(MAT(KEL).EQ.0) GO TO 150
+      XX(KEL)=XXX(K2+1)-XXX(K2)
+      YY(KEL)=YYY(K1+1)-YYY(K1)
+      IF(CYLIND) DD(KEL)=0.5*(XXX(K2)+XXX(K2+1))
+      IND1=(K1-1)*(3*LX+1)
+      KN(NUM1+1)=IND1+LX+2*K2
+      KN(NUM1+2)=IND1+LX+2*K2-1
+      KN(NUM1+3)=IND1+LX+2*K2+1
+      KN(NUM1+4)=IND1+K2
+      KN(NUM1+5)=IND1+3*LX+K2+1
+      QFR(NUM2+1:NUM2+4)=0.0
+      IQFR(NUM2+1:NUM2+4)=0
+      BFR(NUM2+1:NUM2+4)=0.0
+      FRX=1.0
+      FRY=1.0
+*----
+*  VOID, REFL OR ZERO BOUNDARY CONTITION.
+*----
+      IF(K2.EQ.1) THEN
+         LOG1=.TRUE.
+      ELSE
+         LOG1=(MAT(KEL-1).EQ.0)
+      ENDIF
+      IF(LOG1) THEN
+         COND=(NCODE(1).EQ.2).OR.((NCODE(1).EQ.1).AND.(ZCODE(1).EQ.1.0))
+         IF(COND) THEN
+            KN(NUM1+2)=0
+         ELSE IF(NCODE(1).EQ.1) THEN
+            IF(ICODE(1).EQ.0) THEN
+              QFR(NUM2+1)=ZALB(1)
+            ELSE
+              QFR(NUM2+1)=1.0
+              IQFR(NUM2+1)=ICODE(1)
+            ENDIF
+         ENDIF
+      ENDIF
+*
+      IF(K2.EQ.LX) THEN
+         LOG2=.TRUE.
+      ELSE
+         LOG2=(MAT(KEL+1).EQ.0)
+      ENDIF
+      IF(LOG2) THEN
+         COND=(NCODE(2).EQ.2).OR.((NCODE(2).EQ.1).AND.(ZCODE(2).EQ.1.0))
+         IF(COND) THEN
+            KN(NUM1+3)=0
+         ELSE IF(NCODE(2).EQ.1) THEN
+            IF(ICODE(2).EQ.0) THEN
+              QFR(NUM2+2)=ZALB(2)
+            ELSE
+              QFR(NUM2+2)=1.0
+              IQFR(NUM2+2)=ICODE(2)
+            ENDIF
+         ENDIF
+      ENDIF
+*
+      IF(K1.EQ.1) THEN
+         LOG3=.TRUE.
+      ELSE
+         LOG3=(MAT(KEL-LX).EQ.0)
+      ENDIF
+      IF(LOG3) THEN
+         COND=(NCODE(3).EQ.2).OR.((NCODE(3).EQ.1).AND.(ZCODE(3).EQ.1.0))
+         IF(COND) THEN
+            KN(NUM1+4)=0
+         ELSE IF(NCODE(3).EQ.1) THEN
+            IF(ICODE(3).EQ.0) THEN
+              QFR(NUM2+3)=ZALB(3)
+            ELSE
+              QFR(NUM2+3)=1.0
+              IQFR(NUM2+3)=ICODE(3)
+            ENDIF
+         ENDIF
+      ENDIF
+*
+      IF(K1.EQ.LY) THEN
+         LOG4=.TRUE.
+      ELSE
+         LOG4=(MAT(KEL+LX).EQ.0)
+      ENDIF
+      IF(LOG4) THEN
+         COND=(NCODE(4).EQ.2).OR.((NCODE(4).EQ.1).AND.(ZCODE(4).EQ.1.0))
+         IF(COND) THEN
+            KN(NUM1+5)=0
+         ELSE IF(NCODE(4).EQ.1) THEN
+            IF(ICODE(4).EQ.0) THEN
+              QFR(NUM2+4)=ZALB(4)
+            ELSE
+              QFR(NUM2+4)=1.0
+              IQFR(NUM2+4)=ICODE(4)
+            ENDIF
+         ENDIF
+      ENDIF
+*----
+*  TRAN BOUNDARY CONDITION.
+*----
+      IF((K2.EQ.LX).AND.(NCODE(2).EQ.4)) THEN
+         KN(NUM1+3)=KN(NUM1+3)-2*LX
+      ENDIF
+      IF((K1.EQ.LY).AND.(NCODE(4).EQ.4)) THEN
+         KN(NUM1+5)=K2
+      ENDIF
+*----
+*  SYME BOUNDARY CONDITION.
+*----
+      IF((NCODE(1).EQ.5).AND.(K2.EQ.1)) THEN
+         QFR(NUM2+1)=QFR(NUM2+2)
+         IQFR(NUM2+1)=IQFR(NUM2+2)
+         FRX=0.5
+         KN(NUM1+2)=-KN(NUM1+3)
+      ELSE IF((NCODE(2).EQ.5).AND.(K2.EQ.LX)) THEN
+         QFR(NUM2+2)=QFR(NUM2+1)
+         IQFR(NUM2+2)=IQFR(NUM2+1)
+         FRX=0.5
+         KN(NUM1+3)=-KN(NUM1+2)
+      ENDIF
+      IF((NCODE(3).EQ.5).AND.(K1.EQ.1)) THEN
+         QFR(NUM2+3)=QFR(NUM2+4)
+         FRY=0.5
+         KN(NUM1+4)=-KN(NUM1+5)
+      ELSE IF((NCODE(4).EQ.5).AND.(K1.EQ.LY)) THEN
+         QFR(NUM2+4)=QFR(NUM2+3)
+         IQFR(NUM2+4)=IQFR(NUM2+3)
+         FRY=0.5
+         KN(NUM1+5)=-KN(NUM1+4)
+      ENDIF
+*
+      VOL0=XX(KEL)*YY(KEL)*FRX*FRY
+      IF(CYLIND) THEN
+         VOL0=6.2831853072*DD(KEL)*VOL0
+      ENDIF
+      VOL(KEL)=VOL0
+      QFR(NUM2+1)=QFR(NUM2+1)*VOL0/XX(KEL)
+      QFR(NUM2+2)=QFR(NUM2+2)*VOL0/XX(KEL)
+      QFR(NUM2+3)=QFR(NUM2+3)*VOL0/YY(KEL)
+      QFR(NUM2+4)=QFR(NUM2+4)*VOL0/YY(KEL)
+*
+      IF(((NCODE(1).EQ.1).OR.(NCODE(1).EQ.7)).AND.LOG1)
+     1 BFR(NUM2+1)=VOL0/XX(KEL)
+      IF(((NCODE(2).EQ.1).OR.(NCODE(2).EQ.7)).AND.LOG2)
+     1 BFR(NUM2+2)=VOL0/XX(KEL)
+      IF(((NCODE(3).EQ.1).OR.(NCODE(3).EQ.7)).AND.LOG3)
+     1 BFR(NUM2+3)=VOL0/YY(KEL)
+      IF(((NCODE(4).EQ.1).OR.(NCODE(4).EQ.7)).AND.LOG4)
+     1 BFR(NUM2+4)=VOL0/YY(KEL)
+      SURFTOT=SURFTOT+BFR(NUM2+1)+BFR(NUM2+2)+BFR(NUM2+3)+BFR(NUM2+4)
+      NUM1=NUM1+5
+      NUM2=NUM2+4
+  150 CONTINUE
+  151 CONTINUE
+* END OF THE MAIN LOOP OVER ELEMENTS.
+*
+* COMPUTE THE SURFACE FRACTIONS.
+      IF(SURFTOT.GT.0.0) THEN
+         DO 155 I=1,4*LX*LY
+         BFR(I)=BFR(I)/SURFTOT
+  155    CONTINUE
+      ENDIF
+*----
+*  REMOVING THE UNUSED UNKNOWNS INDICES FROM KN.
+*----
+      LL4=LY*(3*LX+1)+LX
+      WRITE (TEXT8,'(I8)') LL4
+      IF(LL4.GT.MAXEV) CALL XABORT('BIVDKN: MAXEV SHOULD BE INCREASED '
+     1 //'TO'//TEXT8//'.')
+      DO 160 IND=1,LL4
+      IP(IND)=0
+  160 CONTINUE
+      DO 170 NUM1=1,5*LX*LY
+      IF(KN(NUM1).NE.0) IP(ABS(KN(NUM1)))=1
+  170 CONTINUE
+      L4=0
+      DO 180 IND=1,LL4
+      IF(IP(IND).EQ.1) THEN
+         L4=L4+1
+         IP(IND)=L4
+      ENDIF
+  180 CONTINUE
+      DO 190 NUM1=1,5*LX*LY
+      IF(KN(NUM1).NE.0) KN(NUM1)=SIGN(IP(ABS(KN(NUM1))),KN(NUM1))
+  190 CONTINUE
+*----
+*  PROCESS CASES WITH IELEM.GT.1.
+*----
+      IF(IELEM.GT.1) THEN
+         LL4=0
+         DO 220 IND=1,L4
+         IP(IND)=LL4+1
+         NUM1=0
+         DO 210 KEL=1,LX*LY
+         IF(MAT(KEL).EQ.0) GO TO 210
+         IF(ABS(KN(NUM1+1)).EQ.IND) THEN
+            LL4=LL4+IELEM**2
+            GO TO 220
+         ELSE
+            DO 200 I=2,5
+            IF(ABS(KN(NUM1+I)).EQ.IND) THEN
+               LL4=LL4+IELEM
+               GO TO 220
+            ENDIF
+  200       CONTINUE
+         ENDIF
+         NUM1=NUM1+5
+  210    CONTINUE
+         CALL XABORT('BIVDKN: FAILURE OF THE RENUMBERING ALGORITHM.')
+  220    CONTINUE
+         L4=LL4
+         DO 230 NUM1=1,5*LX*LY
+         IF(KN(NUM1).NE.0) KN(NUM1)=SIGN(IP(ABS(KN(NUM1))),KN(NUM1))
+  230    CONTINUE
+      ENDIF
+      NUM1=0
+      DO 235 KEL=1,LX*LY
+      IDL(KEL)=0
+      IF(MAT(KEL).EQ.0) GO TO 235
+      IDL(KEL)=KN(NUM1+1)
+      NUM1=NUM1+5
+  235 CONTINUE
+      WRITE (TEXT8,'(I8)') L4
+      IF(L4.GT.MAXEV) CALL XABORT('BIVDKN: MAXEV SHOULD BE INCREASED TO'
+     1 //TEXT8//'.')
+      IF(IMPX.GT.2) WRITE(6,710) (VOL(I),I=1,LX*LY)
+*----
+*  COMPUTE THE SYSTEM MATRIX BANDWIDTH.
+*----
+      DO 240 I=1,L4
+      MU(I)=1
+  240 CONTINUE
+      NUM1=0
+      DO 270 KEL=1,LX*LY
+      IF(MAT(KEL).EQ.0) GO TO 270
+      DO 260 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 250 J0=1,IELEM
+      JND1=KN(NUM1+1)+(I0-1)*IELEM+J0-1
+      IF(IELEM.GE.4) MU(JND1)=MAX(MU(JND1),J0)
+      IF(KN(NUM1+2).NE.0) THEN
+         MU(JND1)=MAX(MU(JND1),JND1-INX1+1)
+         MU(INX1)=MAX(MU(INX1),INX1-JND1+1)
+      ENDIF
+      IF(KN(NUM1+3).NE.0) THEN
+         MU(INX2)=MAX(MU(INX2),INX2-JND1+1)
+         MU(JND1)=MAX(MU(JND1),JND1-INX2+1)
+      ENDIF
+      JND1=KN(NUM1+1)+(J0-1)*IELEM+I0-1
+      IF(IELEM.GE.4) MU(JND1)=MAX(MU(JND1),(J0-1)*IELEM+1)
+      IF(KN(NUM1+4).NE.0) THEN
+         MU(JND1)=MAX(MU(JND1),JND1-INY1+1)
+         MU(INY1)=MAX(MU(INY1),INY1-JND1+1)
+      ENDIF
+      IF(KN(NUM1+5).NE.0) THEN
+         MU(INY2)=MAX(MU(INY2),INY2-JND1+1)
+         MU(JND1)=MAX(MU(JND1),JND1-INY2+1)
+      ENDIF
+  250 CONTINUE
+      IF(ICOL.NE.2) THEN
+         IF((KN(NUM1+2).NE.0).AND.(KN(NUM1+3).NE.0)) THEN
+            MU(INX2)=MAX(MU(INX2),INX2-INX1+1)
+            MU(INX1)=MAX(MU(INX1),INX1-INX2+1)
+         ENDIF
+         IF((KN(NUM1+4).NE.0).AND.(KN(NUM1+5).NE.0)) THEN
+            MU(INY2)=MAX(MU(INY2),INY2-INY1+1)
+            MU(INY1)=MAX(MU(INY1),INY1-INY2+1)
+         ENDIF
+      ENDIF
+  260 CONTINUE
+      NUM1=NUM1+5
+  270 CONTINUE
+      IIMAX=0
+      DO 280 I=1,L4
+      IIMAX=IIMAX+MU(I)
+      MU(I)=IIMAX
+  280 CONTINUE
+*
+      IF(IMPX.GT.2) THEN
+         WRITE (6,720) IIMAX
+         NUM1=0
+         NUM2=0
+         WRITE (6,750)
+         DO 500 K=1,LX*LY
+         IF(MAT(K).EQ.0) GO TO 500
+         WRITE (6,755) K,(KN(NUM1+I),I=1,5),(QFR(NUM2+I),I=1,4),
+     1   (BFR(NUM2+I),I=1,4)
+         NUM1=NUM1+5
+         NUM2=NUM2+4
+  500    CONTINUE
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(IP)
+      RETURN
+*
+  700 FORMAT(/42H BIVDKN: MIXED-DUAL FINITE ELEMENT METHOD.//7H NUMBER,
+     1 27H OF ELEMENTS ALONG X AXIS =,I3/26H NUMBER OF ELEMENTS ALONG ,
+     2 8HY AXIS =,I3)
+  710 FORMAT(/20H VOLUMES PER ELEMENT/(1X,1P,10E13.4))
+  720 FORMAT(/52H NUMBER OF TERMS IN THE COMPRESSED SYSTEM MATRICES =,
+     1 I7)
+  750 FORMAT(/22H NUMBERING OF UNKNOWNS/1X,21(1H-)//
+     1 8H ELEMENT,2X,7HNUMBERS,30X,23HVOID BOUNDARY CONDITION,23X,
+     2 17HSURFACE FRACTIONS)
+  755 FORMAT (1X,I4,2X,5I7,2X,1P,4E11.2,3X,4E10.2)
+      END