Initial commit from Polytechnique Montreal

1 files changed, 457 insertions, 0 deletions

diff --git a/Dragon/src/XELBIN.f b/Dragon/src/XELBIN.f
new file mode 100644
index 0000000..f84be2f
--- /dev/null
+++ b/Dragon/src/XELBIN.f
@@ -0,0 +1,457 @@
+*DECK XELBIN
+      SUBROUTINE XELBIN( IPGEOM,   NDIM, NGEOME, L1CELL, NTYPES,  NGIDL,
+     >                    NTIDL, NBLOCK, MAXGRI,  NUNKO,   IPRT,  CELLG,
+     >                    NSURO,  NVOLO, IDLGEO, MATGEO, KEYGEO, IDLTYP,
+     >                   IDLBLK, KEYTYP, MATTYP, KEYINT)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Identify every zone of every type  to its material and 
+* interface all internal surfaces for cells present in the supercell. 
+*
+*Copyright:
+* Copyright (C) 1987 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): R. Roy
+*
+*Parameters: input
+* IPGEOM  pointer to the geometry (l_geom).             
+* NDIM    number of dimensions (2 or 3).                    
+* NGEOME  number of geometries.                             
+* L1CELL  .true. if only one cell.                     
+* NTYPES  number of types.                                  
+* NGIDL   lenght of geometric numbering.               
+* NTIDL   lenght of type numbering.                    
+* NBLOCK  number of blocks.                                 
+* MAXGRI  number of cells along each axis.                  
+* NUNKO   old number of unknowns.                           
+* IPRT    intermediate printing level for output.      
+* CELLG   to keep geometry  names.                     
+* NSURO   number of surfaces of each geometry.              
+* NVOLO   number of zones of each geometry.                 
+* IDLGEO  position of each geometry in the             
+*         geometry numbering scheme.        
+* MATGEO  material numbers corresponding to geometries.     
+* KEYGEO  geometric key for each type.                 
+* IDLTYP  position of each type in numbering scheme.   
+* IDLBLK  position of each block in numbering scheme.  
+* KEYTYP  type key for each block.                     
+*
+*Parameters: output
+* MATTYP  material numbers for zones of every type.         
+* KEYINT  interface key (giving the connected surface).
+*
+*-----------------------------------------------------------------------
+*
+      USE                GANLIB
+      IMPLICIT           NONE
+*
+      TYPE(C_PTR)        IPGEOM 
+      INTEGER            NDIM, NGEOME, NTYPES, NGIDL, NTIDL, NBLOCK,
+     >                   NUNKO, IPRT
+      INTEGER             NSURO(NGEOME),  NVOLO(NGEOME), IDLGEO(NGEOME),
+     >                   MATGEO( NGIDL), KEYGEO(NTYPES), IDLTYP(NTYPES),
+     >                   MATTYP( NTIDL), KEYTYP(NBLOCK), IDLBLK(NBLOCK),
+     >                   KEYINT(NUNKO ), MAXGRI(NDIM)  , CELLG(3*NTYPES)
+*
+      INTEGER            ILO(3,2), NO(2), KTYP(2),
+     >                     KMAT(2),    KSUR(2), KABSO(2), KSID(2),
+     >                   ICOORD(3),   NCODE(6)
+      CHARACTER          GEOCEL*12, TEDATA*12, TEMESH(4)*7
+      LOGICAL            SWKILL, L1CELL, LL1, LL2
+      INTEGER            NSTATE, IOUT, MAXSPL
+      PARAMETER        ( NSTATE=40, IOUT=6, MAXSPL=100 )
+      INTEGER            ISTATE(NSTATE),ISPLT(MAXSPL)
+      INTEGER            NUMGEO, NUMTYP, NUMBLK, I, K
+      INTEGER            NBMD
+      INTEGER            IMYT, IUNK, ITYP, IMYG, IGEO, NSUX, NVOX, ICYL,
+     >                   ICX, ICY, ICZ, LR, LX, LY, LZ, KOLD, ITYPG,
+     >                   ISUR, IX, IY, IZ, IOFF, KNEW, ILEN, ITYLCM,
+     >                   ISX, ISY, ISZ, ISR, KIOFX, KIOFY, KIOFZ,
+     >                   J0, J1, J2, JC, JR, IP0, IP1, IP2, N, NP1, NP2,
+     >                   K0, K1, K2, K3, KR, IBLK, ISUX
+      EQUIVALENCE      ( ICOORD(1),LX ),(ICOORD(2),LY),(ICOORD(3),LZ )
+      DATA        TEMESH / 'X', 'Y', 'Z', 'R'/
+*
+      NUMGEO(I,K)= I + IDLGEO(K)
+      NUMTYP(I,K)= I + IDLTYP(K)
+      NUMBLK(I,K)= I + IDLBLK(K)
+*
+      SWKILL= .FALSE.
+      LL1= .FALSE.
+      LL2= .FALSE.
+      DO 10 IMYT= 1,  NTIDL
+          MATTYP(IMYT)=  0
+   10 CONTINUE
+      DO 20 IUNK= 1, NUNKO
+          KEYINT(IUNK)=  0
+   20 CONTINUE
+      DO 40 ITYP= 1, NTYPES
+         IGEO  = KEYGEO( ITYP   )
+         NVOX  = NVOLO( IGEO )
+         NSUX  = NSURO( IGEO ) 
+         IF( .NOT.L1CELL )THEN
+            WRITE(GEOCEL( 1: 4), '(A4)') CELLG(3*ITYP-2)
+            WRITE(GEOCEL( 5: 8), '(A4)') CELLG(3*ITYP-1)
+            WRITE(GEOCEL( 9:12), '(A4)') CELLG(3*ITYP  )
+            CALL LCMSIX(IPGEOM, GEOCEL, 1)
+         ELSE
+            CALL LCMGET(IPGEOM,'NCODE', NCODE)
+            LL1=((NCODE(2).EQ.3).AND.(NCODE(3).EQ.3))
+            LL2=((NCODE(1).EQ.3).AND.(NCODE(4).EQ.3))
+         ENDIF
+         ISTATE(:NSTATE)=0
+         CALL LCMGET(IPGEOM,'STATE-VECTOR',ISTATE)
+         ITYPG= ISTATE(1)
+         IF( ITYPG.EQ.20) THEN
+*           FOR *CARCEL* GEOMETRIES
+            ICYL= 1
+            ICX=  1
+            ICY=  2
+            ICZ=  3
+         ELSEIF(ITYPG.EQ.3.OR.ITYPG.EQ.6 )THEN
+*           FOR *CARCEL*, *TUBE* OR *TUBEZ* GEOMETRIES
+            ICYL= 1
+            ICX=  1
+            ICY=  2
+            ICZ=  3
+            IF( LL1.OR.LL2 )THEN
+               CALL XABORT( 'XELBIN: DIAGONAL SYMETRIES NOT POSSIBLE')
+            ENDIF
+         ELSEIF( ITYPG.GT.20 )THEN
+*           FOR *CARCELX*, *CARCELY* OR *CARCELZ*
+            ICYL= 1
+            ICZ= ITYPG-20
+            ICX= MOD(ICZ  , 3) + 1
+            ICY= MOD(ICZ+1, 3) + 1
+         ELSE
+*           FOR *CAR2D* OR *CAR3D*
+            ICYL= 0
+            ICX=  1
+            ICY=  2
+            ICZ=  3
+         ENDIF
+         LR=    ISTATE(2)
+         LX=    MAX(1,ISTATE(3))
+         LY=    MAX(1,ISTATE(4))
+         LZ=    MAX(1,ISTATE(5))
+         KOLD=  ISTATE(6)
+         DO 30 ISUR= NSUX, -1
+            MATTYP(NUMTYP(ISUR,ITYP))= MATGEO(NUMGEO(ISUR,IGEO))
+   30    CONTINUE
+*
+*        GET MIXTURE NUMBERS
+         CALL LCMLEN(IPGEOM, 'MIX', ILEN, ITYLCM)
+         IF( ILEN.NE.KOLD )THEN
+            WRITE(IOUT,*) 'LENGHT(MIX)=  ',ILEN
+            WRITE(IOUT,*) '# OF VOLUMES= ',KOLD
+            CALL LCMLIB(IPGEOM)
+            CALL XABORT( 'XELBIN: INVALID NUMBER OF MIXTURES')
+         ENDIF
+         CALL LCMGET(IPGEOM,'MIX',MATTYP(NUMTYP(1,ITYP)))
+*
+*        IN THE CASE OF DIAGONAL SYMMETRY IN 'ONE-CELL'
+*        CAR2D AND CAR3D  GEOMETRY UNFOLD MIXTURES 
+*
+      IF(ITYPG .LT. 20) THEN
+         K3=ISTATE(6)
+         NBMD=(LZ*LY*(LX+1))/2
+         IF(K3 .EQ. NBMD) THEN
+*----
+* MIXTURE ENTERED IN DIAGONAL FORM
+*----
+         IF( LL1 )THEN
+            DO 72 IZ=LZ,1,-1
+            IOFF=(IZ-1)*LX*LY
+            DO 71 IY=LY,1,-1
+            DO 60 IX=LX,IY+1,-1
+            MATTYP(NUMTYP(IOFF+(IY-1)*LX+IX,ITYP))=
+     >                    MATTYP(NUMTYP(IOFF+(IX-1)*LY+IY,ITYP))
+   60       CONTINUE
+            DO 70 IX=IY,1,-1
+            MATTYP(NUMTYP(IOFF+(IY-1)*LX+IX,ITYP))=
+     >                    MATTYP(NUMTYP(K3,ITYP))
+            K3=K3-1
+   70       CONTINUE
+   71       CONTINUE
+   72       CONTINUE
+            KOLD= LX*LY*LZ
+         ELSEIF( LL2 )THEN
+            DO 82 IZ=LZ,1,-1
+            IOFF=(IZ-1)*LX*LY
+            DO 81 IY=LY,1,-1
+            DO 80 IX=LX,IY,-1
+            MATTYP(NUMTYP(IOFF+(IY-1)*LX+IX,ITYP))=
+     >                    MATTYP(NUMTYP(K3,ITYP))
+            K3=K3-1
+   80       CONTINUE
+   81       CONTINUE
+   82       CONTINUE
+            DO 92 IZ=1,LZ
+            IOFF=(IZ-1)*LX*LY
+            DO 91 IY=1,LY
+            DO 90 IX=1,IY-1
+            MATTYP(NUMTYP(IOFF+(IY-1)*LX+IX,ITYP))=
+     >               MATTYP(NUMTYP(IOFF+(IX-1)*LY+IY,ITYP))
+   90       CONTINUE
+   91       CONTINUE
+   92       CONTINUE
+            KOLD= LX*LY*LZ
+         ENDIF
+         ENDIF
+      ENDIF
+*
+*        FOR THE PARTICULAR CASE OF *TUBE* OR *TUBEZ* GEOMETRIES
+         IF( ITYPG.EQ.3.OR.ITYPG.EQ.6 )THEN
+            DO 39 IZ= 1, LZ
+               MATTYP(NUMTYP(KOLD+IZ,ITYP))= -2
+   39       CONTINUE
+            KOLD= KOLD+LZ
+         ENDIF
+*
+*        FILL UP MATTYP ACCORDING TO SPLITTING VALUES.
+         KNEW= NVOX
+         ISR= 0
+         ISX= 0
+         ISY= 0
+         ISZ= 0
+         DO 308 K0= ICOORD(ICZ),1,-1
+            KIOFZ= KOLD
+            TEDATA= 'SPLIT'//TEMESH(ICZ)
+            CALL LCMLEN(IPGEOM,TEDATA,ILEN,ITYLCM)
+            IF( ILEN.GT.MAXSPL )THEN
+               CALL XABORT('XELBIN: SPLIT OVERFLOW ('//TEDATA//')')
+            ELSEIF( ILEN.EQ.0 )THEN
+               ISZ= 1
+            ELSE
+               CALL LCMGET(IPGEOM,TEDATA,ISPLT)
+               ISZ= ISPLT(K0)
+            ENDIF
+         DO 307 J0=ISZ,1,-1
+            KOLD= KIOFZ
+         DO 306 K1= ICOORD(ICY),1,-1
+            KIOFY= KOLD
+            TEDATA= 'SPLIT'//TEMESH(ICY)
+            CALL LCMLEN(IPGEOM,TEDATA,ILEN,ITYLCM)
+            IF( ILEN.GT.MAXSPL )THEN
+               CALL XABORT('XELBIN: SPLIT OVERFLOW ('//TEDATA//')')
+            ELSEIF( ILEN.EQ.0 )THEN
+               ISY= 1
+            ELSE
+               CALL LCMGET(IPGEOM,TEDATA,ISPLT)
+               ISY= ISPLT(K1)
+            ENDIF
+         DO 305 J1=ISY,1,-1
+            KOLD= KIOFY
+         DO 304 K2= ICOORD(ICX),1,-1
+            KIOFX= KOLD
+            TEDATA= 'SPLIT'//TEMESH(ICX)
+            CALL LCMLEN(IPGEOM,TEDATA,ILEN,ITYLCM)
+            IF( ILEN.GT.MAXSPL )THEN
+               CALL XABORT('XELBIN: SPLIT OVERFLOW ('//TEDATA//')')
+            ELSEIF( ILEN.EQ.0 )THEN
+               ISX= 1
+            ELSE
+               CALL LCMGET(IPGEOM,TEDATA,ISPLT)
+               ISX= ISPLT(K2)
+            ENDIF
+         DO 303 J2=ISX,1,-1
+            KOLD= KIOFX
+*           FOR RECTANGULAR OUTER REGIONS.
+            IMYT= MATTYP(NUMTYP(KOLD,ITYP))
+            MATTYP(NUMTYP(KNEW,ITYP))= IMYT
+            KNEW= KNEW-1
+            KOLD= KOLD-1
+            IF( ICYL.EQ.1 )THEN
+*              FOR CYLINDRICAL INNER REGIONS.
+               DO 302 KR= LR,1,-1
+                  TEDATA= 'SPLIT'//TEMESH(4)
+                  CALL LCMLEN(IPGEOM,TEDATA,ILEN,ITYLCM)
+                  IF( ILEN.GT.MAXSPL )THEN
+                    CALL XABORT('XELBIN: SPLIT OVERFLOW ('//TEDATA//')')
+                  ELSEIF( ILEN.EQ.0 )THEN
+                     ISR= 1
+                  ELSE
+                     CALL LCMGET(IPGEOM,TEDATA,ISPLT)
+                     ISR= ABS(ISPLT(KR))
+                  ENDIF
+                  IMYT= MATTYP(NUMTYP(KOLD,ITYP))
+               DO 301 JR=ISR,1,-1
+                  MATTYP(NUMTYP(KNEW,ITYP))= IMYT
+                  KNEW= KNEW-1
+  301          CONTINUE
+                  KOLD= KOLD-1
+  302          CONTINUE
+            ENDIF
+  303    CONTINUE
+  304    CONTINUE
+  305    CONTINUE
+  306    CONTINUE
+  307    CONTINUE
+  308    CONTINUE
+         IF( KNEW.NE.0 )THEN
+            WRITE(IOUT,*) 'XELBIN: KNEW.NE.0 = PROBLEM WITH SPLITTING'
+            SWKILL= .TRUE.
+         ENDIF
+         IF( KOLD.NE.0 )THEN
+            WRITE(IOUT,*) 'XELBIN: KOLD.NE.0 = PROBLEM WITH SPLITTING'
+            SWKILL= .TRUE.
+         ENDIF
+*
+         IF( .NOT.L1CELL ) CALL LCMSIX(IPGEOM, ' ',    2)
+   40 CONTINUE
+*
+*     RECOMPOSE INTERNAL SURFACES COUPLING (INTERFACES)
+*        THIS ASSUMES THAT AN ORDERING OF SURFACES IS DONE
+*        BECAUSE:  SIDE-BY-SIDE INTERFACES
+*                  ARE SUPPOSED IN INCREASING POSITION.
+      DO 220 N= 1, NDIM
+*
+*        DEFINITION OF THE SIDE NUMBER TO COUPLE.
+         KSID(1)=  -2*N
+         KSID(2)= (-2*N) + 1
+         NP1   = MOD(N  ,NDIM) + 1
+         IF( NDIM.EQ.3 )THEN
+            NP2   = MOD(N+1,NDIM) + 1
+            DO 112 IP1= 1, MAXGRI(NP1)
+               ILO(NP1,1)= IP1
+               ILO(NP1,2)= IP1
+            DO 111 IP2= 1, MAXGRI(NP2)
+               ILO(NP2,1)= IP2
+               ILO(NP2,2)= IP2
+            DO 110 IP0= 1, MAXGRI(N)-1
+               ILO(N  ,1)= IP0
+               ILO(N  ,2)= IP0 + 1
+               DO 100  JC= 1, 2
+                  NO(JC)= MAXGRI(1)*(MAXGRI(2)*ILO(3,JC)+ILO(2,JC)-
+     >                               MAXGRI(2))+ILO(1,JC)-MAXGRI(1)
+                  KTYP(JC)= KEYTYP( NO(JC) )
+                  IF( KTYP(JC).EQ.0 ) GO TO 110
+                  IGEO  = KEYGEO( KTYP(JC) )
+*                 SEARCH FROM THE END
+                  KSUR(JC)= NSURO(IGEO)
+                  KMAT(JC)= MATTYP( NUMTYP(KSUR(JC),KTYP(JC)) )
+  100          CONTINUE
+*
+*              ORDERING INTERFACING OF THE TWO BLOCKS.
+  101          CONTINUE
+                  IF( KMAT(1).EQ.KSID(1).AND.KMAT(2).EQ.KSID(2) )THEN
+                     IF( KSUR(1).EQ.0 .OR. KSUR(2).EQ.0 ) GO TO 109
+                     KABSO(1)= NUMBLK( KSUR(1),NO(1) )
+                     KABSO(2)= NUMBLK( KSUR(2),NO(2) )
+                     KEYINT( KABSO(1) )= KABSO(2)
+                     KEYINT( KABSO(2) )= KABSO(1)
+                     KSUR(1)= KSUR(1)+1
+                     KSUR(2)= KSUR(2)+1
+                  ELSE
+                     IF( KMAT(1).NE.KSID(1) ) KSUR(1)= KSUR(1)+1
+                     IF( KMAT(2).NE.KSID(2) ) KSUR(2)= KSUR(2)+1
+                  ENDIF
+                  IF( KSUR(1).NE.0 )THEN
+                     KMAT(1)= MATTYP( NUMTYP(KSUR(1),KTYP(1)) )
+                  ELSE
+                     KMAT(1)= KSID(1)
+                  ENDIF
+                  IF( KSUR(2).NE.0 )THEN
+                     KMAT(2)= MATTYP( NUMTYP(KSUR(2),KTYP(2)) )
+                  ELSE
+                     KMAT(2)= KSID(2)
+                  ENDIF
+               GO TO 101
+  109          IF( KSUR(1).NE.0 .OR. KSUR(2).NE.0 )THEN
+                  WRITE(IOUT,'(1H ,I8,4H OF ,I8,5H <=> ,I8,4H OF ,I8)')
+     >                          KSUR(1),  NO(1),     KSUR(2),  NO(2)
+                  SWKILL=.TRUE.
+               ENDIF
+  110       CONTINUE
+  111       CONTINUE
+  112       CONTINUE
+         ELSEIF( NDIM.EQ.2 )THEN
+            DO 215 IP1= 1, MAXGRI(NP1)
+               ILO(NP1,1)= IP1
+               ILO(NP1,2)= IP1
+            DO 210 IP0= 1, MAXGRI(N)-1
+               ILO(N  ,1)= IP0
+               ILO(N  ,2)= IP0 + 1
+               DO 200  JC= 1, 2
+                  NO(JC)= MAXGRI(1) * (ILO(2,JC) - 1) + ILO(1,JC)
+                  KTYP(JC)= KEYTYP( NO(JC) )
+                  IF( KTYP(JC).EQ.0 ) GO TO 210
+                  IGEO  = KEYGEO( KTYP(JC) )
+*                 SEARCH FROM THE END
+                  KSUR(JC)= NSURO(IGEO)
+                  KMAT(JC)= MATTYP( NUMTYP(KSUR(JC),KTYP(JC)) ) 
+  200          CONTINUE
+*
+*              ORDERING INTERFACING OF THE TWO BLOCKS.
+  201          CONTINUE
+                  IF( KMAT(1).EQ.KSID(1).AND.KMAT(2).EQ.KSID(2) )THEN
+                     IF( KSUR(1).EQ.0 .OR. KSUR(2).EQ.0 ) GO TO 209
+                     KABSO(1)= NUMBLK( KSUR(1),NO(1) )
+                     KABSO(2)= NUMBLK( KSUR(2),NO(2) )
+                     KEYINT( KABSO(1) )= KABSO(2)
+                     KEYINT( KABSO(2) )= KABSO(1)
+                     KSUR(1)= KSUR(1)+1
+                     KSUR(2)= KSUR(2)+1
+                  ELSE
+                     IF( KMAT(1).NE.KSID(1) ) KSUR(1)= KSUR(1)+1
+                     IF( KMAT(2).NE.KSID(2) ) KSUR(2)= KSUR(2)+1
+                  ENDIF
+                  IF( KSUR(1).NE.0 )THEN
+                     KMAT(1)= MATTYP( NUMTYP(KSUR(1),KTYP(1)) )
+                  ELSE
+                     KMAT(1)= KSID(1)
+                  ENDIF
+                  IF( KSUR(2).NE.0 )THEN
+                     KMAT(2)= MATTYP( NUMTYP(KSUR(2),KTYP(2)) )
+                  ELSE
+                     KMAT(2)= KSID(2)
+                  ENDIF
+               GO TO 201
+  209          IF( KSUR(1).NE.0 .OR. KSUR(2).NE.0 )THEN
+                  WRITE(IOUT,'(1H ,I8,4H OF ,I8,5H <=> ,I8,4H OF ,I8)')
+     >                          KSUR(1),  NO(1),     KSUR(2),  NO(2)
+                  SWKILL=.TRUE.
+               ENDIF
+  210       CONTINUE
+  215       CONTINUE
+         ELSE
+            CALL LCMLIB(IPGEOM)
+            CALL XABORT( 'XELBIN: *** FALSE NDIM VALUE')
+         ENDIF
+  220 CONTINUE
+*
+      IF( IPRT.GE.100 .OR. SWKILL )THEN
+         IUNK= 0
+         WRITE(IOUT,'(/40H       KEYINT       COUPLE    MATERIAL  )')
+         DO 250 IBLK= 1, NBLOCK
+            ITYP= KEYTYP(IBLK)
+            IGEO= KEYGEO(ITYP)
+            NVOX= NVOLO(IGEO)
+            NSUX= NSURO(IGEO)
+            DO 240 ISUX= NSUX, NVOX
+               IUNK= IUNK+1
+               IMYT= MATTYP( NUMTYP(ISUX,ITYP) )
+               IF( ISUX.LT.0 )THEN
+                  WRITE(IOUT,
+     >            '(5H SUR(,I8,5H) => ,I8,4H OF ,    I8)')
+     >                      IUNK,      KEYINT(IUNK), IMYT
+               ELSEIF( ISUX.GT.0 )THEN
+                  IMYG= MATGEO( NUMGEO(ISUX,IGEO) )
+                  WRITE(IOUT,
+     >            '(5H VOL(,I8,5H) => ,I8,4H OF ,    I8,1H(,I8,1H))')
+     >                      IUNK,      KEYINT(IUNK), IMYT,  IMYG
+               ENDIF
+  240       CONTINUE
+            WRITE(IOUT,'(/1X)')
+  250    CONTINUE
+      ENDIF
+      IF( SWKILL ) CALL XABORT( 'XELBIN: IMPOSSIBLE TO INTERFACE')
+*
+      RETURN
+      END