summaryrefslogtreecommitdiff
path: root/Trivac/src/TRIHEX.f
diff options
context:
space:
mode:
Diffstat (limited to 'Trivac/src/TRIHEX.f')
-rwxr-xr-xTrivac/src/TRIHEX.f382
1 files changed, 382 insertions, 0 deletions
diff --git a/Trivac/src/TRIHEX.f b/Trivac/src/TRIHEX.f
new file mode 100755
index 0000000..c976e54
--- /dev/null
+++ b/Trivac/src/TRIHEX.f
@@ -0,0 +1,382 @@
+*DECK TRIHEX
+ SUBROUTINE TRIHEX (IOPT,LX,LZ,LL4,MAT,KN,NCODE,IPTRK)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Numbering corresponding to a mesh corner finite difference or
+* Lagrangian finite element discretization of hexagonal 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. Benaboud
+*
+*Parameters: input
+* IOPT type of hexagonal Lagrangian finite element:
+* = 1 for hexagonal element with 6 points;
+* = 2 for hexagonal element with 7 points and triangular element
+* (IOPT.le.2) for Bivac;
+* = 3 for hexagonal element with 6 points;
+* = 4 for hexagonal element with 7 points and triangular element
+* (IOPT.gt.2) for Trivac.
+* LX number of elements in the XY plane.
+* LZ number of elements along Z axis.
+* MAT mixture index assigned to each element.
+*
+*Parameters: output
+* LL4 total number of unknown (variational coefficients) per
+* energy group (order of system matrices).
+* KN element-ordered unknown list. Dimensionned to LC*LX*LZ
+* where LC= 7 for triangle/Bivac, 6 for hexagon/Bivac,
+* 14 for triangle/Trivac and 12 for hexagon/Trivac.
+* IPTRK L_TRACK pointer to the tracking information.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ INTEGER IOPT,LX,LZ,LL4,MAT(LX*LZ),KN(*),NCODE(6)
+ TYPE(C_PTR) IPTRK
+*----
+* LOCAL VARIABLES
+*----
+ LOGICAL LC1,LMD,LMG,LED,LEG,LIV,CADI,LNC1,LNC2
+ INTEGER IRO(2,7),ITAB(14),ICO(6,6)
+ INTEGER, DIMENSION(:), ALLOCATABLE :: IISS,IPER,NWXY,IENV,IDX,IDY
+ INTEGER, DIMENSION(:,:), ALLOCATABLE :: ICG
+ INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: NIK
+ DATA ITAB /0,2*1,0,2*-1,0,1,0,-1,0,1,0,-1/
+ DATA IRO /2,3,3,7,7,6,4,4,1,2,5,1,6,5/
+ DATA ICO /2,1,5,6,7,3,1,5,6,7,3,2,3,2,1,5,6,7,2,1,5,6,7,3,7,3,2,
+ > 1,5,6,3,2,1,5,6,7/
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(IISS(LX),IPER(LX),NWXY(LX),IENV(LX))
+ ALLOCATE(ICG(14,LX),NIK(3,7,LX))
+*
+ IZI = 0
+ IZE = 0
+ IZS = 0
+ KEL = 0
+ IPAR = IOPT
+ NC = INT((SQRT(REAL((4*LX-1)/3))+1.)/2.)
+ IF(IPAR.LT.1.OR.IPAR.GT.4) CALL XABORT('TRIHEX : INVALID DATA.')
+ M1 = 2 + 3*(NC-1)*(NC-2)
+ IF(NC.EQ.1) M1=1
+ DO 10 KX = 1,LX
+ IISS(KX) = 0
+ IENV(KX) = KX
+ IF(MAT(KX).LE.0) GO TO 10
+ KEL = KEL + 1
+ IISS(KEL) = KX
+ 10 CONTINUE
+ DO 20 IY = 1,14
+ ICG(IY,1) = ITAB(IY)
+ 20 CONTINUE
+ DO 45 J = 1,LX
+ DO 40 KF = 1, 6
+ N = NEIGHB(J,KF,9,LX,POIDS)
+ IF(N.GT.LX) GOTO 40
+ DO 30 I = 1,14
+ IF(I.LE.7) THEN
+ IF((KF.EQ.1).OR.(KF.EQ.5)) ICG(I,N)=ICG(I,J)+1
+ IF((KF.EQ.2).OR.(KF.EQ.4)) ICG(I,N)=ICG(I,J)-1
+ IF(KF.EQ.3) ICG(I,N)=ICG(I,J)-2
+ IF(KF.EQ.6) ICG(I,N)=ICG(I,J)+2
+ ELSE
+ IF((KF.EQ.1).OR.(KF.EQ.3)) ICG(I,N)=ICG(I,J)+1
+ IF((KF.EQ.4).OR.(KF.EQ.6)) ICG(I,N)=ICG(I,J)-1
+ IF(KF.EQ.5) ICG(I,N)=ICG(I,J)-2
+ IF(KF.EQ.2) ICG(I,N)=ICG(I,J)+2
+ ENDIF
+ 30 CONTINUE
+ 40 CONTINUE
+ 45 CONTINUE
+ ICHOIX = 1
+ CADI = .FALSE.
+ IF(IPAR.GT.2) THEN
+ CADI = .TRUE.
+ ICHOIX = 3
+ IPAR = IPAR - 2
+ ENDIF
+ NELEMW = -1
+ NELEMX = -1
+ NELEMY = -1
+ DO 190 IWXY = 1,ICHOIX
+ DO 51 IY = 1,7
+ DO 50 IX = 1,LX
+ NIK(IWXY,IY,IX) = 0
+ 50 CONTINUE
+ 51 CONTINUE
+ IF(NCODE(1).EQ.7) THEN
+ DO 60 J = 1,LX
+ IF(MAT(J).LE.0) GO TO 60
+ DO 55 KF = 1,6
+ N = NEIGHB(J,KF,9,LX,POIDS)
+ IF((N.GT.LX).OR.(MAT(N).LE.0)) THEN
+ NIK(IWXY,ICO(KF,2*IWXY-1),J) = -1
+ NIK(IWXY,ICO(KF,2*IWXY),J) = -1
+ ENDIF
+ 55 CONTINUE
+ 60 CONTINUE
+ ENDIF
+ IF(IWXY.EQ.2.OR.IWXY.EQ.3) M1 = M1 + NC - 1
+ CALL BIVPER(M1,IWXY,LX,LX,IPER,IENV)
+ DO 65 I = 1,LX
+ NWXY(IPER(I)) = I
+ 65 CONTINUE
+ NELEM = 1
+ DO 70 I = 1,NC
+ IV = NWXY(I)
+ IF(MAT(IV).LE.0) GO TO 70
+ IF(NIK(IWXY,1,IV).EQ.-1) GO TO 70
+ NIK(IWXY,1,IV) = NELEM
+ NELEM = NELEM + 1
+ 70 CONTINUE
+ DO 140 I = 1,LX
+ LIV = (I.GT.1)
+ IV1 = NWXY(I)
+ IF(MAT(IV1).LE.0) THEN
+ IFACE1 = IWXY + 5
+ IFACE2 = IWXY + 2
+ IF(IFACE1.GT.6) IFACE1 = IFACE1 - 6
+ IF(IFACE2.GT.6) IFACE2 = IFACE2 - 6
+ IND = NEIGHB(IV1,IFACE1,9,LX,P)
+ ING = NEIGHB(IV1,IFACE2,9,LX,P)
+ LED = (IND.GT.LX)
+ LEG = (ING.GT.LX)
+ LMD=.FALSE.
+ LMG=.FALSE.
+ IF(.NOT.LED) LMD = (MAT(IND).LE.0)
+ IF(.NOT.LEG) LMG = (MAT(ING).LE.0)
+ IF(LED.OR.LEG.OR.LMD.OR.LMG) THEN
+ IVK1 = -1
+ IVK2 = -1
+ IF((LMD.OR.LED).AND.(LMG.OR.LEG)) THEN
+ IVK1 = 1
+ IVK2 = 0
+ ELSE IF(LMD.OR.LED) THEN
+ IVK1 = 1
+ IVK2 = 3
+ ELSE IF(LMG.OR.LEG) THEN
+ IVK1 = 0
+ IVK2 = 0
+ ENDIF
+ IFACE3 = IWXY+IVK1+3
+ IFACE4 = IWXY+IVK2+2
+ IFACE5 = IWXY-IVK1+1
+ IF(IFACE3.GT.6) IFACE3 = IFACE3 - 6
+ IF(IFACE4.GT.6) IFACE4 = IFACE4 - 6
+ IF(IFACE5.GT.6) IFACE5 = IFACE5 - 6
+ IV3 = NEIGHB(IV1,IFACE3,9,LX,P)
+ IF((IV3.GT.LX).OR.(MAT(IV3).LE.0)) GOTO 90
+ 80 IF(NIK(IWXY,1,IV3).EQ.0) THEN
+ NIK(IWXY,1,IV3) = NELEM
+ NELEM = NELEM + 1
+ ENDIF
+ IV3 = NEIGHB(IV3,IFACE4,9,LX,P)
+ IF((IV3.LE.LX).AND.(MAT(IV3).GT.0)) THEN
+ IV5 = NEIGHB(IV3,IFACE4-1,9,LX,P)
+ IF((IV5.GT.LX).OR.(MAT(IV5).GT.0)) GO TO 140
+ GO TO 80
+ ENDIF
+ 90 IV4 = NEIGHB(IV1,IFACE5,9,LX,P)
+ IF((IV4.GT.LX).OR.(MAT(IV4).LE.0)) GOTO 140
+ 100 IF(NIK(IWXY,7,IV4).EQ.0) THEN
+ NIK(IWXY,7,IV4) = NELEM
+ NELEM = NELEM + 1
+ ENDIF
+ IV4 = NEIGHB(IV4,IFACE4,9,LX,P)
+ IF((IV4.LE.LX).AND.(MAT(IV4).GT.0)) GOTO 100
+ ENDIF
+ ELSE
+ IF(IWXY.EQ.1) THEN
+ DO 110 K = 0,4
+ IF((IPAR.EQ.1).AND.(K.EQ.2)) GO TO 110
+ IF(LIV) THEN
+ IV2 = NWXY(I-1)
+ IF((ICG(K+2,IV1).EQ.ICG(5,IV2))
+ > .AND.(MAT(IV2).GT.0)) THEN
+ NIK(1,2,IV1) = NIK(1,5,IV2)
+ NIK(1,3,IV1) = NIK(1,6,IV2)
+ GO TO 110
+ ENDIF
+ ENDIF
+ IF(NIK(1,K+2,IV1).EQ.-1) GO TO 110
+ NIK(1,K+2,IV1) = NELEM
+ NELEM = NELEM + 1
+ 110 CONTINUE
+ ELSE IF(IWXY.EQ.2) THEN
+ DO 120 K = 0,4
+ IF((IPAR.EQ.1).AND.(K.EQ.2)) GO TO 120
+ IF(LIV) THEN
+ IV2 = NWXY(I-1)
+ IF((ICG(K+1,IV1).EQ.ICG(6,IV2))
+ > .AND.(ICG(K+8,IV1).EQ.ICG(13,IV2))
+ > .AND.(MAT(IV2).GT.0)) THEN
+ NIK(2,2,IV1) = NIK(2,5,IV2)
+ GO TO 120
+ ELSE IF((ICG(K+1,IV1).EQ.ICG(7,IV2))
+ > .AND.(ICG(K+8,IV1).EQ.ICG(14,IV2))
+ > .AND.(MAT(IV2).GT.0)) THEN
+ NIK(2,3,IV1) = NIK(2,6,IV2)
+ GO TO 120
+ ENDIF
+ ENDIF
+ IF(NIK(2,K+2,IV1).EQ.-1) GO TO 120
+ NIK(2,K+2,IV1) = NELEM
+ NELEM = NELEM + 1
+ 120 CONTINUE
+ ELSE IF(IWXY.EQ.3) THEN
+ LK = 1
+ DO 130 K = 5,1,-1
+ LK = LK + 1
+ IF((IPAR.EQ.1).AND.(LK.EQ.4)) GO TO 130
+ IF(LIV) THEN
+ IV2 = NWXY(I-1)
+ IF((ICG(K,IV1).EQ.ICG(7,IV2))
+ > .AND.(ICG(K+7,IV1).EQ.ICG(14,IV2))
+ > .AND.(MAT(IV2).GT.0)) THEN
+ NIK(3,2,IV1) = NIK(3,5,IV2)
+ GO TO 130
+ ELSE IF((ICG(K-3,IV1).EQ.ICG(3,IV2))
+ > .AND.(ICG(K+4,IV1).EQ.ICG(10,IV2))
+ > .AND.(MAT(IV2).GT.0)) THEN
+ NIK(3,3,IV1) = NIK(3,6,IV2)
+ GO TO 130
+ ENDIF
+ ENDIF
+ IF(NIK(3,LK,IV1).EQ.-1) GO TO 130
+ NIK(3,LK,IV1) = NELEM
+ NELEM = NELEM + 1
+ 130 CONTINUE
+ ENDIF
+ ENDIF
+ 140 CONTINUE
+ DO 160 I = 1,LX
+ IV = NWXY(I)
+ IF(MAT(IV).LE.0) GO TO 160
+ DO 150 K = 1,2
+ IFACE = K+IWXY-1
+ IF(IFACE.GT.6) IFACE = IFACE - 6
+ INE = NEIGHB(IV,IFACE,9,LX,P)
+ IF((INE.GT.LX).OR.(MAT(INE).LE.0)) GO TO 150
+ IF(K.EQ.1) NIK(IWXY,1,IV) = NIK(IWXY,6,INE)
+ IF(K.EQ.2) NIK(IWXY,1,IV) = NIK(IWXY,3,INE)
+ 150 CONTINUE
+ 160 CONTINUE
+ DO 180 I = 1,LX
+ IV = NWXY(I)
+ IF(MAT(IV).LE.0) GO TO 180
+ LC1 = .FALSE.
+ LIV = .TRUE.
+ DO 170 K = 4,5
+ IFACE = K+IWXY-1
+ IF(IFACE.GT.6) IFACE = IFACE - 6
+ INE = NEIGHB(IV,IFACE,9,LX,P)
+ IF((INE.GT.LX).OR.(MAT(INE).LE.0)) THEN
+ IF(K.EQ.4) LC1 = .TRUE.
+ IF(NIK(IWXY,7,IV).EQ.-1) GO TO 170
+ IF(LC1.AND.K.EQ.5.AND.(NIK(IWXY,7,IV).EQ.0)) THEN
+ NIK(IWXY,7,IV) = NELEM
+ NELEM = NELEM + 1
+ ENDIF
+ ELSE
+ IF(K.EQ.4) THEN
+ LIV = .FALSE.
+ NIK(IWXY,7,IV) = NIK(IWXY,2,INE)
+ ENDIF
+ IF(K.EQ.5.AND.LIV) NIK(IWXY,7,IV)=NIK(IWXY,5,INE)
+ ENDIF
+ 170 CONTINUE
+ 180 CONTINUE
+ IF(IWXY.EQ.1) NELEMW = NELEM - 1
+ IF(IWXY.EQ.2) NELEMX = NELEM - 1
+ IF(IWXY.EQ.3) NELEMY = NELEM - 1
+ 190 CONTINUE
+ MEL = 0
+ NPT = NELEM -1
+ IF(ICHOIX.EQ.3) THEN
+ IF((NELEMW.NE.NELEMX).OR.(NELEMW.NE.NELEMY)) THEN
+ CALL XABORT('TRIHEX: ECHEC DE LA NUMEROTATION.')
+ ENDIF
+ ALLOCATE(IDX(NPT),IDY(NPT))
+ DO 200 I = 1,LX
+ IF(MAT(I).LE.0) GO TO 200
+ DO 195 J = 1,7
+ IF((IPAR.EQ.1).AND.(J.EQ.4)) GO TO 195
+ IF(NIK(1,J,I).EQ.-1) GO TO 195
+ IDX(NIK(1,J,I)) = NIK(2,IRO(1,J),I)
+ IDY(NIK(1,J,I)) = NIK(3,IRO(2,J),I)
+ 195 CONTINUE
+ 200 CONTINUE
+ CALL LCMPUT(IPTRK,'ILX',NPT,1,IDX)
+ CALL LCMPUT(IPTRK,'ILY',NPT,1,IDY)
+ DEALLOCATE(IDY,IDX)
+ ENDIF
+ NINC = 6
+ IF(CADI) NINC = 12
+ IF(IPAR.EQ.2) NINC = 7
+ IF((IPAR.EQ.2).AND.CADI) NINC = 14
+ LNC1 = .FALSE.
+ LNC2 = .FALSE.
+ KNZ = 0
+ IF(LZ.GT.1) KNZ = NPT
+ IF(CADI.AND.(LZ.GT.1)) CALL LCMPUT(IPTRK,'NCODE',6,1,NCODE)
+ IF((LZ.GT.1).AND.((NCODE(5).EQ.7).OR.(NCODE(6).EQ.7))) THEN
+ LNC1 = (NCODE(5).EQ.7)
+ LNC2 = (NCODE(6).EQ.7)
+ IF(LNC1) IZI = 1
+ IF(LNC2) IZS = 1
+ IF(LNC1.AND.LNC2) THEN
+ IZE = 1
+ IZS = 1
+ IZI = 0
+ ENDIF
+ ENDIF
+ DO 230 MZ = 1, LZ
+ KEL = 0
+ DO 220 MX = 1, LX
+ IF(MAT(MX).LE.0) GO TO 220
+ KEL = KEL + 1
+ ML = IISS(KEL)
+ IJX = 1
+ DO 210 JX = 1,7
+ KN(MEL+IJX) = 0
+ IF(CADI) KN(MEL+IJX+NINC/2) = 0
+ IF((IPAR.EQ.1).AND.(JX.EQ.4)) GO TO 210
+ IF(NIK(1,JX,ML).EQ.-1) GO TO 205
+ IF(MZ.EQ.1.AND.LNC1) THEN
+ KN(MEL+IJX+NINC/2) = NIK(1,JX,ML)
+ GO TO 205
+ ELSE IF(MZ.EQ.LZ.AND.LNC2) THEN
+ KN(MEL+IJX) = NIK(1,JX,ML)+(MZ-IZS-IZE)*KNZ
+ GO TO 205
+ ENDIF
+ KN(MEL+IJX) = NIK(1,JX,ML)+(MZ-1-IZE-IZI)*KNZ
+ IF(CADI) KN(MEL+IJX+NINC/2)=NIK(1,JX,ML)+(MZ-IZE-IZI)*KNZ
+ 205 IJX = IJX + 1
+ 210 CONTINUE
+ MEL = MEL + NINC
+ 220 CONTINUE
+ 230 CONTINUE
+ LL4 = NPT
+ IF(LZ.GT.1) THEN
+ LL4 = NPT*(LZ+1)
+ IF(LNC1.OR.LNC2) LL4 = NPT*LZ
+ IF(LNC1.AND.LNC2) LL4 = NPT*(LZ-1)
+ ENDIF
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(NIK,ICG,IENV,NWXY,IPER,IISS)
+ RETURN
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 13:22:45 +0000