Initial commit from Polytechnique Montreal

1 files changed, 254 insertions, 0 deletions

diff --git a/Trivac/src/TRICH1.f b/Trivac/src/TRICH1.f
new file mode 100755
index 0000000..6aa237d
--- /dev/null
+++ b/Trivac/src/TRICH1.f
@@ -0,0 +1,254 @@
+*DECK TRICH1
+      SUBROUTINE TRICH1(IELEM,IDIM,LX,LY,LZ,L4,MAT,KN,MUX,MUY,MUZ,IPY,
+     1 IPZ,IMPX)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* compute the compressed diagonal storage indices (MUX, MUY and MUZ)
+* and the permutation vectors (IPY and IPZ) for an ADI splitting of
+* the nodal collocation leakage 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
+* IELEM   degree of the polynomial expansion: =1 (linear);
+*         =2 (parabolic); =3 (cubic); =4 (quartic).
+* IDIM    number of dimensions.
+* LX      number of mesh along the X axis.
+* LY      number of mesh along the Y axis.
+* LZ      number of mesh along the Z axis.
+* L4      order of system matrices
+* MAT     mixture index assigned to each element.
+* KN      element-ordered unknown list.
+* IMPX    print parameter (equal to zero for no print).
+*
+*Parameters: output
+* MUX     X-oriented compressed storage mode indices.
+* MUY     Y-oriented compressed storage mode indices.
+* MUZ     Z-oriented compressed storage mode indices.
+* IPY     Y-oriented permutation matrices.
+* IPZ     Z-oriented permutation matrices.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER IELEM,IDIM,LX,LY,LZ,L4,MAT(LX*LY*LZ),KN(7*LX*LY*LZ),
+     1 MUX(L4),MUY(L4),MUZ(L4),IPY(L4),IPZ(L4),IMPX
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER, DIMENSION(:), ALLOCATABLE :: IWRK
+*
+      IORD(J,K,L,LL,IEL,IW)=(IEL*L+K)*LL*IEL+(1+IEL*(IW-1))+J
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(IWRK(LX*LY*LZ))
+*
+      IWRK(:LX*LY*LZ)=0
+      LL4=0
+      KEL=0
+      DO 22 K0=1,LZ
+      DO 21 K1=1,LY
+      DO 20 K2=1,LX
+      KEL=KEL+1
+      IF(MAT(KEL).EQ.0) GO TO 20
+      LL4=LL4+1
+      IWRK((K0-1)*LX*LY+(K1-1)*LX+K2)=LL4
+   20 CONTINUE
+   21 CONTINUE
+   22 CONTINUE
+*----
+*  COMPUTE THE PERMUTATION VECTORS IPY AND IPZ
+*----
+      IF(IDIM.GE.2) THEN
+         INX1=0
+         DO 52 K0=1,LZ
+         DO 51 K2=1,LX
+         DO 50 K1=1,LY
+         INX2=IWRK((K0-1)*LX*LY+(K1-1)*LX+K2)
+         IF(INX2.LE.0) GO TO 50
+         INX1=INX1+1
+         IF(IDIM.EQ.2) THEN
+            DO 31 K=0,IELEM-1
+            DO 30 J=0,IELEM-1
+            I=IORD(J,K,0,LL4,IELEM,INX1)
+            IPY(IORD(K,J,0,LL4,IELEM,INX2))=I
+   30       CONTINUE
+   31       CONTINUE
+         ELSE IF(IDIM.EQ.3) THEN
+            DO 42 L=0,IELEM-1
+            DO 41 K=0,IELEM-1
+            DO 40 J=0,IELEM-1
+            I=IORD(J,K,L,LL4,IELEM,INX1)
+            IPY(IORD(K,J,L,LL4,IELEM,INX2))=I
+   40       CONTINUE
+   41       CONTINUE
+   42       CONTINUE
+         ENDIF
+   50    CONTINUE
+   51    CONTINUE
+   52    CONTINUE
+         IF(INX1.NE.LL4) CALL XABORT('TRICH1: FAILURE OF THE RENUMBERI'
+     1   //'NG ALGORITHM(1)')
+         IF(IDIM.EQ.3) THEN
+            INX1=0
+            DO 72 K1=1,LY
+            DO 71 K2=1,LX
+            DO 70 K0=1,LZ
+            INX2=IWRK((K0-1)*LX*LY+(K1-1)*LX+K2)
+            IF(INX2.LE.0) GO TO 70
+            INX1=INX1+1
+            DO 62 L=0,IELEM-1
+            DO 61 K=0,IELEM-1
+            DO 60 J=0,IELEM-1
+            I=IORD(J,K,L,LL4,IELEM,INX1)
+            IPZ(IORD(K,L,J,LL4,IELEM,INX2))=I
+   60       CONTINUE
+   61       CONTINUE
+   62       CONTINUE
+   70       CONTINUE
+   71       CONTINUE
+   72       CONTINUE
+            IF(INX1.NE.LL4) CALL XABORT('TRICH1: FAILURE OF THE RENUMB'
+     1      //'ERING ALGORITHM(2)')
+         ENDIF
+      ENDIF
+*
+      L2M=0
+      DO 80 KEL=1,LX*LY*LZ
+      IF(MAT(KEL).EQ.0) GO TO 80
+      L2M=L2M+1
+      IWRK(KEL)=L2M
+   80 CONTINUE
+      DO 90 I=1,L4
+      MUY(I)=0
+      MUZ(I)=0
+   90 CONTINUE
+      LL5=L4/IELEM**(IDIM-1)
+*----
+*  COMPUTE VECTOR MUX
+*----
+      NUM1=0
+      DO 130 KEL=1,LL4
+      KK1=KN(NUM1+1)
+      KK2=KN(NUM1+2)
+      DO 100 J=0,IELEM-1
+      INX1=IORD(J,0,0,LL4,IELEM,KEL)
+      MUX(INX1)=J+1
+*     X- SIDE:
+      IF(KK1.GT.0) THEN
+         INX2=IORD(0,0,0,LL4,IELEM,IWRK(KK1))
+         MUX(INX1)=MAX0(MUX(INX1),INX1-INX2+1)
+      ENDIF
+*     X+ SIDE:
+      IF(KK2.GT.0) THEN
+         INX2=IORD(0,0,0,LL4,IELEM,IWRK(KK2))
+         MUX(INX1)=MAX0(MUX(INX1),INX1-INX2+1)
+      ENDIF
+  100 CONTINUE
+      NUM1=NUM1+6
+  130 CONTINUE
+*----
+*  COMPUTE VECTOR MUY
+*----
+      IF(IDIM.GE.2) THEN
+         NUM1=0
+         DO 160 KEL=1,LL4
+         KK3=KN(NUM1+3)
+         KK4=KN(NUM1+4)
+         DO 140 K=0,IELEM-1
+         INY1=IPY(IORD(0,K,0,LL4,IELEM,KEL))
+         MUY(INY1)=K+1
+*        Y- SIDE:
+         IF(KK3.GT.0) THEN
+            INY2=IPY(IORD(0,0,0,LL4,IELEM,IWRK(KK3)))
+            MUY(INY1)=MAX0(MUY(INY1),INY1-INY2+1)
+         ENDIF
+*        Y+ SIDE:
+         IF(KK4.GT.0) THEN
+            INY2=IPY(IORD(0,0,0,LL4,IELEM,IWRK(KK4)))
+            MUY(INY1)=MAX0(MUY(INY1),INY1-INY2+1)
+         ENDIF
+  140    CONTINUE
+         NUM1=NUM1+6
+  160    CONTINUE
+*----
+*  COMPUTE VECTOR MUZ
+*----
+         IF(IDIM.EQ.3) THEN
+            NUM1=0
+            DO 180 KEL=1,LL4
+            KK5=KN(NUM1+5)
+            KK6=KN(NUM1+6)
+            DO 170 L=0,IELEM-1
+            INZ1=IPZ(IORD(0,0,L,LL4,IELEM,KEL))
+            MUZ(INZ1)=L+1
+*           Z- SIDE:
+            IF(KK5.GT.0) THEN
+               INZ2=IPZ(IORD(0,0,0,LL4,IELEM,IWRK(KK5)))
+               MUZ(INZ1)=MAX0(MUZ(INZ1),INZ1-INZ2+1)
+            ENDIF
+*           Z+ SIDE:
+            IF(KK6.GT.0) THEN
+               INZ2=IPZ(IORD(0,0,0,LL4,IELEM,IWRK(KK6)))
+               MUZ(INZ1)=MAX0(MUZ(INZ1),INZ1-INZ2+1)
+            ENDIF
+  170       CONTINUE
+            NUM1=NUM1+6
+  180       CONTINUE
+            DO 195 J=1,IELEM-1
+            DO 190 I=1,LL5
+            MUX(I+J*LL5)=MUX(I)
+            MUY(I+J*LL5)=MUY(I)
+            MUZ(I+J*LL5)=MUZ(I)
+  190       CONTINUE
+  195       CONTINUE
+            LL5=IELEM*LL5
+         ENDIF
+         DO 205 J=1,IELEM-1
+         DO 200 I=1,LL5
+         MUX(I+J*LL5)=MUX(I)
+         MUY(I+J*LL5)=MUY(I)
+         MUZ(I+J*LL5)=MUZ(I)
+  200    CONTINUE
+  205    CONTINUE
+      ENDIF
+*
+      MUXMAX=0
+      MUYMAX=0
+      MUZMAX=0
+      IIMAXX=0
+      IIMAXY=0
+      IIMAXZ=0
+      DO 210 I=1,L4
+      MUXMAX=MAX(MUXMAX,MUX(I))
+      MUYMAX=MAX(MUYMAX,MUY(I))
+      MUZMAX=MAX(MUZMAX,MUZ(I))
+      IIMAXX=IIMAXX+MUX(I)
+      MUX(I)=IIMAXX
+      IIMAXY=IIMAXY+MUY(I)
+      MUY(I)=IIMAXY
+      IIMAXZ=IIMAXZ+MUZ(I)
+      MUZ(I)=IIMAXZ
+  210 CONTINUE
+      IF(IMPX.GT.0) WRITE (6,230) MUXMAX,MUYMAX,MUZMAX
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(IWRK)
+      RETURN
+*
+  230 FORMAT(/41H TRICH1: MAXIMUM BANDWIDTH ALONG X AXIS =,I5/
+     1 27X,14HALONG Y AXIS =,I5/27X,14HALONG Z AXIS =,I5)
+      END