Initial commit from Polytechnique Montreal

1 files changed, 249 insertions, 0 deletions

diff --git a/Trivac/src/OUTHOM.f b/Trivac/src/OUTHOM.f
new file mode 100755
index 0000000..45341bf
--- /dev/null
+++ b/Trivac/src/OUTHOM.f
@@ -0,0 +1,249 @@
+*DECK OUTHOM
+      SUBROUTINE OUTHOM(MAXNEL,IPGEOM,IMPX,NEL,IELEM,ICOL,HTRACK,MAT,
+     1 NZS,IHOM)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Read an modify the merge indices.
+*
+*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
+* MAXNEL  maximum number of elements.
+* IPGEOM  L_GEOM pointer to the geometry.
+* IMPX    print parameter.
+* NEL     total number of finite elements.
+* IELEM   degree of the Lagrangian finite elements:
+* ICOL    type of quadrature used to integrate the mass matrix
+* HTRACK  type of tracking (equal to 'BIVAC' or 'TRIVAC').
+* MAT     index-number of the mixture type assigned to each volume.
+*
+*Parameters: output
+* NZS     number of merged regions.
+* IHOM    merge indices.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPGEOM
+      CHARACTER HTRACK*12
+      INTEGER MAXNEL,IMPX,NEL,IELEM,ICOL,MAT(NEL),NZS,IHOM(NEL)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (NSTATE=40)
+      INTEGER ISTATE(NSTATE),NCODE(6),ICODE(6)
+      REAL ZCODE(6)
+      LOGICAL ILK,LDIAG,CHEX,LFOLD1,LFOLD2,LFOLD3
+      DOUBLE PRECISION DFLOTT
+      CHARACTER TEXT4*4,HSMG*131
+      INTEGER, DIMENSION(:), ALLOCATABLE :: MAT2,DPP,MX,XXX,YYY,ZZZ
+      INTEGER, DIMENSION(:), ALLOCATABLE :: ISPLX,ISPLY,ISPLZ
+      EQUIVALENCE (ITYPE,ISTATE(1)),(LR1,ISTATE(2)),(LX1,ISTATE(3)),
+     1 (LY1,ISTATE(4)),(LZ1,ISTATE(5))
+*----
+*  DETERMINE THE MESH SPLITTING INFO FROM THE GEOMETRY.
+*----
+      ALLOCATE(ISPLX(MAXNEL),ISPLY(MAXNEL),ISPLZ(MAXNEL))
+      ALLOCATE(XXX(MAXNEL+1),YYY(MAXNEL+1),ZZZ(MAXNEL+1))
+*
+      ALLOCATE(MAT2(MAXNEL))
+      CALL READ3D(MAXNEL,MAXNEL,MAXNEL,MAXNEL,IPGEOM,IHEX,IR,ILK,SIDE,
+     1 XXX,YYY,ZZZ,IMPX,LX,LY,LZ,MAT2,IPAS,NCODE,ICODE,ZCODE,ISPLX,
+     2 ISPLY,ISPLZ,ISPLH,ISPLL)
+      DEALLOCATE(MAT2)
+*
+      CALL LCMGET(IPGEOM,'STATE-VECTOR',ISTATE)
+      LYOLD=1
+      LZOLD=1
+      NELOLD=0
+      IF(ITYPE.EQ.2) THEN
+*        1-D CARTESIAN GEOMETRY.
+         LXOLD=LX1
+         NELOLD=LXOLD
+      ELSE IF((ITYPE.EQ.3).OR.(ITYPE.EQ.4)) THEN
+*        1-D CYLINDRICAL/SPHERICAL GEOMETRY.
+         LXOLD=LR1
+         NELOLD=LXOLD
+      ELSE IF(ITYPE.EQ.5) THEN
+*        2-D CARTESIAN GEOMETRY.
+         LXOLD=LX1
+         LYOLD=LY1
+         NELOLD=LXOLD*LYOLD
+         LDIAG=.FALSE.
+         DO 30 IC=1,4
+         LDIAG=LDIAG.OR.(NCODE(IC).EQ.3)
+   30    CONTINUE
+         IF(LDIAG) NELOLD=(LXOLD+1)*LXOLD/2
+      ELSE IF(ITYPE.EQ.6) THEN
+*        2-D CYLINDRICAL GEOMETRY.
+         LXOLD=LR1
+         LZOLD=LZ1
+         NELOLD=LXOLD*LZOLD
+      ELSE IF(ITYPE.EQ.7) THEN
+*        3-D CARTESIAN GEOMETRY.
+         LXOLD=LX1
+         LYOLD=LY1
+         LZOLD=LZ1
+         NELOLD=LXOLD*LYOLD*LZOLD
+         LDIAG=.FALSE.
+         DO 40 IC=1,4
+         LDIAG=LDIAG.OR.(NCODE(IC).EQ.3)
+   40    CONTINUE
+         IF(LDIAG) NELOLD=(LXOLD+1)*LXOLD*LZOLD/2
+      ELSE IF(ITYPE.EQ.8) THEN
+*        2-D HEXAGONAL GEOMETRY.
+         LXOLD=LX1
+         NELOLD=LXOLD
+      ELSE IF(ITYPE.EQ.9) THEN
+*        3-D HEXAGONAL GEOMETRY.
+         LXOLD=LX1
+         LZOLD=LZ1
+         NELOLD=LXOLD*LZOLD
+      ENDIF
+*----
+*  READ THE MERGE INDICES.
+*----
+      CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DFLOTT)
+      IF(INDIC.EQ.1) THEN
+        DO 160 K=1,NELOLD
+       IHOM(K)=0
+  160   CONTINUE
+        IHOM(1)=NITMA
+        NZS=NITMA
+        DO 170 K=2,NELOLD
+        CALL REDGET(INDIC,IHOM(K),FLOTT,TEXT4,DFLOTT)
+        IF(INDIC.NE.1) CALL XABORT('OUTHOM: INTEGER EXPECTED.')
+        NZS=MAX(NZS,IHOM(K))
+  170   CONTINUE
+        IF((ITYPE.EQ.8).OR.(ITYPE.EQ.9)) CALL LCMGET(IPGEOM,'IHEX',IHEX)
+      ELSE IF((INDIC.EQ.3).AND.(TEXT4.EQ.'NONE')) THEN
+        NZS=NEL
+        DO 180 K=1,NEL
+        IHOM(K)=K
+  180   CONTINUE
+        GO TO 270
+      ELSE IF((INDIC.EQ.3).AND.(TEXT4.EQ.'IN')) THEN
+        DO 190 K=1,NELOLD
+        IHOM(K)=K
+  190   CONTINUE
+        NZS=NELOLD
+        IF((ITYPE.EQ.8).OR.(ITYPE.EQ.9)) CALL LCMGET(IPGEOM,'IHEX',IHEX)
+      ELSE IF((INDIC.EQ.3).AND.(TEXT4.EQ.'MIX')) THEN
+        CALL LCMLEN(IPGEOM,'MIX',ILONG,ITYLCM)
+        IF(ILONG.NE.NELOLD) THEN
+          WRITE(HSMG,'(42HOUTHOM: INCONSISTENT INTG MIX OPTION (EXPE,
+     1    24HCTED NUMBER OF MIXTURES=,I5,24H; VALUE FOUND IN L_GEOM ,
+     2    7HOBJECT=,I5,2H).)') NELOLD,ILONG
+          CALL XABORT(HSMG)
+        ENDIF
+        CALL LCMGET(IPGEOM,'MIX',IHOM)
+        NZS=0
+        DO 200 K=1,NELOLD
+        NZS=MAX(NZS,IHOM(K))
+  200   CONTINUE
+        IF((ITYPE.EQ.8).OR.(ITYPE.EQ.9)) CALL LCMGET(IPGEOM,'IHEX',IHEX)
+      ELSE
+        CALL XABORT('OUTHOM: INVALID KEY WORD.')
+      ENDIF
+*----
+*  UNFOLD HEXAGONAL GEOMETRY IN BIVAC AND TRIVAC CASES.
+*----
+      CHEX=(ITYPE.EQ.8).OR.(ITYPE.EQ.9)
+      LFOLD1=CHEX.AND.(IHEX.NE.9).AND.(HTRACK.EQ.'TRIVAC')
+      LFOLD2=CHEX.AND.(IHEX.NE.9).AND.(HTRACK.EQ.'BIVAC').AND.
+     1       (IELEM.GT.0).AND.(ICOL.LE.3)
+      LFOLD3=CHEX.AND.(IHEX.NE.9).AND.((HTRACK.EQ.'MCCG').OR.
+     1       (HTRACK.EQ.'EXCELL'))
+      IF(LFOLD1.OR.LFOLD2.OR.LFOLD3) THEN
+         IF(NELOLD.NE.LXOLD*LZOLD) CALL XABORT('OUTHOM: HEXAGONAL SPLI'
+     1   //'T ERROR.')
+         ALLOCATE(DPP(MAXNEL),MX(NELOLD))
+         DO 205 I=1,NELOLD
+         MX(I)=IHOM(I)
+  205    CONTINUE
+         LXOLD=LX1
+         CALL BIVALL(MAXNEL,IHEX,LXOLD,LX,DPP)
+         DO 215 KZ=1,LZOLD
+         DO 210 KX=1,LX
+         IHOM(KX+(KZ-1)*LX)=0
+         KEL=DPP(KX)+(KZ-1)*LXOLD
+         IF(KEL.GT.LXOLD*LZOLD) CALL XABORT('OUTHOM: MX OVERFLOW.')
+         IHOM(KX+(KZ-1)*LX)=MX(KEL)
+  210    CONTINUE
+  215    CONTINUE
+         DEALLOCATE(MX,DPP)
+         LXOLD=LX
+         IHEX=9
+      ENDIF
+*----
+*  MESH-SPLITTING FOR THE IHOM VECTOR.
+*----
+      IF(NZS.GT.NELOLD) CALL XABORT('OUTHOM: FAILURE 1.')
+      IF(ISTATE(11).NE.0) THEN
+         CALL SPLIT0(MAXNEL,ITYPE,NCODE,LXOLD,LYOLD,LZOLD,ISPLX,ISPLY,
+     1   ISPLZ,0,ISPLL,NEL2,LX,LY,LZ,SIDE,XXX,YYY,ZZZ,IHOM,.FALSE.,IMPX)
+      ENDIF
+*----
+*  FORCE DIAGONAL SYMMETRY AND UNFOLD THE IHOM VECTOR.
+*----
+      IF((NCODE(2).EQ.3).AND.(NCODE(3).EQ.3)) THEN
+         IF(NEL.EQ.LX*LY*LZ) THEN
+            K=(LX*(LX+1)/2)*LZ
+            DO 232 IZ=LZ,1,-1
+            IOFF=(IZ-1)*LX*LY
+            DO 231 IY=LY,1,-1
+            DO 220 IX=LX,IY+1,-1
+            IHOM(IOFF+(IY-1)*LX+IX)=IHOM(IOFF+(IX-1)*LY+IY)
+  220       CONTINUE
+            DO 230 IX=IY,1,-1
+            IHOM(IOFF+(IY-1)*LX+IX)=IHOM(K)
+            K=K-1
+  230       CONTINUE
+  231       CONTINUE
+  232       CONTINUE
+            IF(K.NE.0) CALL XABORT('OUTHOM: FAILURE 2.')
+         ENDIF
+      ELSE IF((NCODE(1).EQ.3).AND.(NCODE(4).EQ.3)) THEN
+         IF(NEL.EQ.LX*LY*LZ) THEN
+            K=(LX*(LX+1)/2)*LZ
+            DO 242 IZ=LZ,1,-1
+            IOFF=(IZ-1)*LX*LY
+            DO 241 IY=LY,1,-1
+            DO 240 IX=LX,IY,-1
+            IHOM(IOFF+(IY-1)*LX+IX)=IHOM(K)
+            K=K-1
+  240       CONTINUE
+  241       CONTINUE
+  242       CONTINUE
+            DO 252 IZ=1,LZ
+            IOFF=(IZ-1)*LX*LY
+            DO 251 IY=1,LY
+            DO 250 IX=1,IY-1
+            IHOM(IOFF+(IY-1)*LX+IX)=IHOM(IOFF+(IX-1)*LY+IY)
+  250       CONTINUE
+  251       CONTINUE
+  252       CONTINUE
+            IF(K.NE.0) CALL XABORT('OUTHOM: FAILURE 3.')
+         ENDIF
+      ENDIF
+      DEALLOCATE(ZZZ,YYY,XXX,ISPLZ,ISPLY,ISPLX)
+      DO 260 K=1,NEL
+      IF(MAT(K).EQ.0) IHOM(K)=0
+  260 CONTINUE
+  270 IF(IMPX.GT.0) THEN
+        WRITE(6,'(/15H MERGING INDEX:/(1X,14I5))') (IHOM(K),K=1,NEL)
+      ENDIF
+      RETURN
+      END