summaryrefslogtreecommitdiff
path: root/Dragon/src/SNTT3D.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/SNTT3D.f')
-rw-r--r--Dragon/src/SNTT3D.f609
1 files changed, 609 insertions, 0 deletions
diff --git a/Dragon/src/SNTT3D.f b/Dragon/src/SNTT3D.f
new file mode 100644
index 0000000..3162d44
--- /dev/null
+++ b/Dragon/src/SNTT3D.f
@@ -0,0 +1,609 @@
+*DECK SNTT3D
+ SUBROUTINE SNTT3D(IGE,IMPX,LX,LY,LZ,SIDE,IELEM,NLF,NPQ,NSCT,IQUAD,
+ 1 NCODE,ZCODE,MAT,XXX,YYY,ZZZ,VOL,IDL,DU,DE,DZ,W,MRMX,MRMY,MRMZ,DC,
+ 2 DB,DA,PL,LL4,NUN,EELEM,WX,WE,CST,IBFP,ISCHM,ESCHM,IGLK,MN,DN,IL,
+ 3 IM,ISCAT)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Numbering corresponding to a 3-D Cartesian with discrete ordinates
+* approximation of the flux.
+*
+*Copyright:
+* Copyright (C) 2007 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): N. Martin and C. Bienvenue
+*
+*Parameters: input
+* IGE geometry type with (=0 for Cartesian, =2 for hexagonal).
+* IMPX print parameter.
+* LX number of elements along the X axis;
+* OR, number of elements in X-Y plane in hexagonal geometry,
+* including lozenges and further submeshing of lozenges
+* e.g. domain of 7 hex. by 5 levels with a submeshing of 2
+* will have LX=7*3*2*2=84
+* LY number of elements along the Y axis.
+* LZ number of elements along the Z axis.
+* SIDE side of hexagon.
+* IELEM measure of order of the spatial approximation polynomial:
+* =1 constant - only for HODD, classical diamond scheme
+* (default for HODD);
+* =2 linear - default for DG;
+* =3 parabolic;
+* =4 cubic - only for DG.
+* NLF SN order for the flux (even number).
+* NPQ number of SN directions in eight octants.
+* NSCT maximum number of spherical harmonics moments of the flux.
+* IQUAD type of SN quadrature (1 Level symmetric, type IQUAD;
+* 4 Legendre-Chebyshev; 5 symmetric Legendre-Chebyshev;
+* 6 quadruple range).
+* NCODE type of boundary condition applied on each side
+* (i=1 X-; i=2 X+; i=3 Y-; i=4 Y+; i=5 Z-; i=6 Z+):
+* =1 VOID; =2 REFL; =4 TRAN.
+* 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.
+* ZZZ Cartesian coordinates along the Z axis.
+* EELEM measure of order of the energy approximation polynomial:
+* =1 constant - default for HODD;
+* =2 linear - default for DG;
+* >3 higher orders.
+* IBFP type of energy proparation relation:
+* =0 no Fokker-Planck term;
+* =1 Galerkin type;
+* =2 heuristic Przybylski and Ligou type.
+* ISCHM method of spatial discretisation:
+* =1 High-Order Diamond Differencing (HODD) - default;
+* =2 Discontinuous Galerkin finite element method (DG);
+* =3 Adaptive weighted method (AWD).
+* ESCHM method of energy discretisation:
+* =1 High-Order Diamond Differencing (HODD) - default;
+* =2 Discontinuous Galerkin finite element method (DG);
+* =3 Adaptive weighted method (AWD).
+* IGLK angular interpolation type:
+* =0 classical SN method.
+* =1 Galerkin quadrature method (M = inv(D))
+* =2 Galerkin quadrature method (D = inv(M))
+* ISCAT maximum number of spherical harmonics moments of the flux.
+*
+*Parameters: output
+* VOL volume of each element.
+* IDL isotropic flux indices.
+* DU first direction cosines ($\\mu$).
+* DE second direction cosines ($\\eta$).
+* DZ third direction cosines ($\\xi$).
+* W weights.
+* MRMX quadrature index.
+* MRMY quadrature index.
+* MRMZ quadrature index.
+* DC diamond-scheme parameter.
+* DB diamond-scheme parameter.
+* DA diamond-scheme parameter.
+* PL discrete values of the spherical harmonics corresponding
+* to the 3D SN quadrature.
+* LL4 number of unknowns being solved for, over the domain. This
+* includes the various moments of the isotropic (and if present,
+* anisotropic) flux.
+* NUN total number of unknowns stored in the FLUX vector per group.
+* This includes LL4 (see above) as well as any surface boundary
+* fluxes, if present.
+* WX spatial closure relation weighting factors.
+* WE energy closure relation weighting factors.
+* CST constants for the polynomial approximations.
+* MN moment-to-discrete matrix.
+* DN discrete-to-moment matrix.
+* IL indexes (l) of each spherical harmonics in the
+* interpolation basis.
+* IM indexes (m) of each spherical harmonics in the
+* interpolation basis.
+*
+*-----------------------------------------------------------------------
+*
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ INTEGER IMPX,LX,LY,LZ,IELEM,NLF,NPQ,NSCT,IQUAD,NCODE(6),
+ 1 MAT(LX,LY,LZ),IDL(LX*LY*LZ),MRMX(NPQ),MRMY(NPQ),MRMZ(NPQ),EELEM,
+ 2 IBFP,ISCHM,ESCHM,IL(NSCT),IM(NSCT),ISCAT,IGLK
+ REAL ZCODE(6),VOL(LX,LY,LZ),XXX(LX+1),YYY(LY+1),ZZZ(LZ+1),
+ 1 DU(NPQ),DE(NPQ),DZ(NPQ),W(NPQ),DC(LX,LY,NPQ),DB(LX,LZ,NPQ),
+ 2 DA(LY,LZ,NPQ),PL(NSCT,NPQ),WX(IELEM+1),
+ 3 WE(EELEM+1),CST(MAX(IELEM,EELEM)),MN(NPQ,NSCT),DN(NSCT,NPQ)
+*----
+* LOCAL VARIABLES
+*----
+ CHARACTER HSMG*131
+ LOGICAL L1,L2,L3,L4,L5,L6
+ PARAMETER(RLOG=1.0E-8,PI=3.141592654)
+ REAL PX,PE
+ DOUBLE PRECISION MND(NPQ,NPQ),NORM,IPROD
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: JOP
+ REAL, ALLOCATABLE, DIMENSION(:) :: XX,YY,ZZ,UU,WW,TPQ,UPQ,VPQ,WPQ
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: V,V2
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: U
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: RLM
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(XX(LX),YY(LY),ZZ(LZ))
+*----
+* UNFOLD HEIGHT-OCTANT QUADRATURES.
+*----
+ IF(MOD(NLF,2).EQ.1) CALL XABORT('SNTT3D: EVEN NLF EXPECTED.')
+ IF(IQUAD.EQ.10) THEN
+ NPQ0=NLF**2/4
+ ELSE
+ NPQ0=NLF*(NLF/2+1)/4
+ ENDIF
+ ALLOCATE(JOP(NLF/2),UU(NLF/2),WW(NLF/2),TPQ(NPQ0),UPQ(NPQ0),
+ 1 VPQ(NPQ0),WPQ(NPQ0))
+ IF(IQUAD.EQ.1) THEN
+!Level-symmetric quadrature of type 1
+ CALL SNQU01(NLF,JOP,UU,WW,TPQ,UPQ,VPQ,WPQ)
+ ELSE IF(IQUAD.EQ.2) THEN
+!Level-symmetric quadrature of type 2
+ CALL SNQU02(NLF,JOP,UU,WW,TPQ,UPQ,VPQ,WPQ)
+ ELSE IF(IQUAD.EQ.3) THEN
+!Snow Level-symmetric type quadrature
+ CALL SNQU03(NLF,JOP,UU,WW,TPQ,UPQ,VPQ,WPQ)
+ ELSE IF(IQUAD.EQ.4) THEN
+!Legendre-Chebyshev quadrature
+ CALL SNQU04(NLF,JOP,UU,WW,TPQ,UPQ,VPQ,WPQ)
+ ELSE IF(IQUAD.EQ.5) THEN
+!Symmetric Legendre-Chebyshev quadrature
+ CALL SNQU05(NLF,TPQ,UPQ,VPQ,WPQ)
+ ELSE IF(IQUAD.EQ.6) THEN
+!Quadruple Range quadrature
+ CALL SNQU06(NLF,TPQ,UPQ,VPQ,WPQ)
+ ELSE IF(IQUAD.EQ.10) THEN
+ CALL SNQU10(NLF,JOP,UU,WW,TPQ,UPQ,VPQ,WPQ)
+ ELSE
+ CALL XABORT('SNTT3D: UNKNOWN QUADRATURE TYPE.')
+ ENDIF
+ N=0
+ IOF=0
+ DO 320 I=1,NLF/2
+ JOF=IOF+NLF-2*I+2
+ DO 330 J=0,NLF/2-I
+ IOF=IOF+1
+ KOF=IOF+(NLF+2)*NLF/4
+ LOF=IOF+(NLF+2)*NLF/2
+ MRMX(IOF)=JOF
+ MRMY(IOF)=KOF
+ MRMZ(IOF)=LOF
+ DU(IOF)=-UPQ(N+J+1)
+ DE(IOF)=-VPQ(N+J+1)
+ DZ(IOF)=-TPQ(N+J+1)
+ W(IOF)=WPQ(N+J+1)
+ JOF=JOF-1
+ 330 CONTINUE
+ DO 340 J=NLF/2-I,0,-1
+ IOF=IOF+1
+ KOF=IOF+(NLF+2)*NLF/4
+ LOF=IOF+(NLF+2)*NLF/2
+ MRMX(IOF)=JOF
+ MRMY(IOF)=KOF
+ MRMZ(IOF)=LOF
+ DU(IOF)=UPQ(N+J+1)
+ DE(IOF)=-VPQ(N+J+1)
+ DZ(IOF)=-TPQ(N+J+1)
+ W(IOF)=WPQ(N+J+1)
+ JOF=JOF-1
+ 340 CONTINUE
+ N=N+NLF/2-I+1
+ 320 CONTINUE
+ N=0
+ DO 350 I=1,NLF/2
+ JOF=IOF+NLF-2*I+2
+ DO 360 J=0,NLF/2-I
+ IOF=IOF+1
+ KOF=IOF-(NLF+2)*NLF/4
+ LOF=IOF+(NLF+2)*NLF/2
+ MRMX(IOF)=JOF
+ MRMY(IOF)=KOF
+ MRMZ(IOF)=LOF
+ DU(IOF)=-UPQ(N+J+1)
+ DE(IOF)=VPQ(N+J+1)
+ DZ(IOF)=-TPQ(N+J+1)
+ W(IOF)=WPQ(N+J+1)
+ JOF=JOF-1
+ 360 CONTINUE
+ DO 370 J=NLF/2-I,0,-1
+ IOF=IOF+1
+ KOF=IOF-(NLF+2)*NLF/4
+ LOF=IOF+(NLF+2)*NLF/2
+ MRMX(IOF)=JOF
+ MRMY(IOF)=KOF
+ MRMZ(IOF)=LOF
+ DU(IOF)=UPQ(N+J+1)
+ DE(IOF)=VPQ(N+J+1)
+ DZ(IOF)=-TPQ(N+J+1)
+ W(IOF)=WPQ(N+J+1)
+ JOF=JOF-1
+ 370 CONTINUE
+ N=N+NLF/2-I+1
+ 350 CONTINUE
+ N=0
+ DO 380 I=1,NLF/2
+ JOF=IOF+NLF-2*I+2
+ DO 390 J=0,NLF/2-I
+ IOF=IOF+1
+ KOF=IOF+(NLF+2)*NLF/4
+ LOF=IOF-(NLF+2)*NLF/2
+ MRMX(IOF)=JOF
+ MRMY(IOF)=KOF
+ MRMZ(IOF)=LOF
+ DU(IOF)=-UPQ(N+J+1)
+ DE(IOF)=-VPQ(N+J+1)
+ DZ(IOF)=TPQ(N+J+1)
+ W(IOF)=WPQ(N+J+1)
+ JOF=JOF-1
+ 390 CONTINUE
+ DO 400 J=NLF/2-I,0,-1
+ IOF=IOF+1
+ KOF=IOF+(NLF+2)*NLF/4
+ LOF=IOF-(NLF+2)*NLF/2
+ MRMX(IOF)=JOF
+ MRMY(IOF)=KOF
+ MRMZ(IOF)=LOF
+ DU(IOF)=UPQ(N+J+1)
+ DE(IOF)=-VPQ(N+J+1)
+ DZ(IOF)=TPQ(N+J+1)
+ W(IOF)=WPQ(N+J+1)
+ JOF=JOF-1
+ 400 CONTINUE
+ N=N+NLF/2-I+1
+ 380 CONTINUE
+ N=0
+ DO 410 I=1,NLF/2
+ JOF=IOF+NLF-2*I+2
+ DO 420 J=0,NLF/2-I
+ IOF=IOF+1
+ KOF=IOF-(NLF+2)*NLF/4
+ LOF=IOF-(NLF+2)*NLF/2
+ MRMX(IOF)=JOF
+ MRMY(IOF)=KOF
+ MRMZ(IOF)=LOF
+ DU(IOF)=-UPQ(N+J+1)
+ DE(IOF)=VPQ(N+J+1)
+ DZ(IOF)=TPQ(N+J+1)
+ W(IOF)=WPQ(N+J+1)
+ JOF=JOF-1
+ 420 CONTINUE
+ DO 430 J=NLF/2-I,0,-1
+ IOF=IOF+1
+ KOF=IOF-(NLF+2)*NLF/4
+ LOF=IOF-(NLF+2)*NLF/2
+ MRMX(IOF)=JOF
+ MRMY(IOF)=KOF
+ MRMZ(IOF)=LOF
+ DU(IOF)=UPQ(N+J+1)
+ DE(IOF)=VPQ(N+J+1)
+ DZ(IOF)=TPQ(N+J+1)
+ W(IOF)=WPQ(N+J+1)
+ JOF=JOF-1
+ 430 CONTINUE
+ N=N+NLF/2-I+1
+ 410 CONTINUE
+ DEALLOCATE(WPQ,VPQ,UPQ,TPQ,WW,UU,JOP)
+ IF(IMPX.GE.4) THEN
+ WRITE(6,'(/41H SNTT3D:HEIGHT-OCTANT ANGULAR QUADRATURES:/26X,
+ 1 2HMU,9X,3HETA,10X,2HXI,6X,6HWEIGHT)')
+ SUM=0.0
+ DO 70 N=1,NPQ
+ SUM=SUM+W(N)
+ WRITE(6,'(1X,4I5,1P,4E12.4)') N,MRMX(N),MRMY(N),MRMZ(N),DU(N),
+ 1 DE(N),DZ(N),W(N)
+ 70 CONTINUE
+ WRITE(6,'(54X,10(1H-)/52X,1P,E12.4)') SUM
+ ENDIF
+*----
+* IDENTIFICATION OF THE GEOMETRY.
+*----
+ IF(IGE.EQ.0) THEN
+* ----------
+* 3D CARTESIAN
+* ----------
+ DO 83 N=1,NPQ
+ VU=DU(N)
+ VE=DE(N)
+ VZ=DZ(N)
+ DO 82 I=1,LX
+ DO 81 J=1,LY
+ DO 80 K=1,LZ
+ XX(I)=XXX(I+1)-XXX(I)
+ YY(J)=YYY(J+1)-YYY(J)
+ ZZ(K)=ZZZ(K+1)-ZZZ(K)
+ DA(J,K,N)=VU*YY(J)*ZZ(K)
+ DB(I,K,N)=VE*XX(I)*ZZ(K)
+ DC(I,J,N)=VZ*XX(I)*YY(J)
+ VOL(I,J,K)=XX(I)*YY(J)*ZZ(K)
+ 80 CONTINUE
+ 81 CONTINUE
+ 82 CONTINUE
+ 83 CONTINUE
+ ELSEIF(IGE.EQ.2) THEN
+* ----------
+* 3D HEXAGONAL
+* ----------
+ DET = SQRT(3.0)*(SIDE**2)/2.0
+ DO 93 N=1,NPQ
+ VU=DU(N)
+ VE=DE(N)
+ VZ=DZ(N)
+ DO 92 K=1,LZ
+ DO 91 J=1,LY
+ DO 90 I=1,LX
+ ZZ(K)=ZZZ(K+1)-ZZZ(K)
+ DA(J,K,N)=VU*ZZ(K)
+ DB(I,K,N)=VE*ZZ(K)
+ DC(I,J,N)=VZ*DET
+ VOL(I,J,K)=DET*ZZ(K)
+ 90 CONTINUE
+ 91 CONTINUE
+ 92 CONTINUE
+ 93 CONTINUE
+ ENDIF
+*----
+* GENERATE SPHERICAL HARMONICS FOR SCATTERING SOURCE.
+*----
+ IOF=0
+ DO 211 L=0,ISCAT-1
+ DO 210 M=-L,L
+ IOF=IOF+1
+ IF(IOF.GT.ISCAT) GO TO 211
+ DO 200 N=1,NPQ
+ PL(IOF,N)=PNSH(L,M,DU(N),DE(N),DZ(N))
+ 200 CONTINUE
+ 210 CONTINUE
+ 211 CONTINUE
+*----
+* GENERATE MAPPING MATRIX FOR GALERKIN QUADRATURE METHOD
+*----
+ MN(:NPQ,:NSCT)=0.0
+ DN(:NSCT,:NPQ)=0.0
+ IL(:NSCT)=0
+ IM(:NSCT)=0
+ IF(IGLK.NE.0) THEN
+ ALLOCATE(U(NPQ,NPQ),RLM(NPQ,ISCAT,2*ISCAT-1),V(NPQ),V2(NPQ))
+ RLM(:NPQ,:ISCAT-1,:2*ISCAT-1)=0.0
+ DO L=0,ISCAT-1
+ DO M=-L,L
+ DO N=1,NPQ
+ RLM(N,L+1,M+L+1)=PNSH(L,M,DU(N),DE(N),DZ(N))
+ ENDDO
+ ENDDO
+ ENDDO
+ ! GRAM-SCHMIDT PROCEDURE TO FIND INDEPENDANT SET
+ ! OF SPHERICAL HARMONICS WITH ANY QUADRATURE
+ U(:NPQ,:NPQ)=0.0D0
+ NORM=0.0D0
+ DO N=1,NPQ
+ NORM=NORM+RLM(N,1,1)**2
+ ENDDO
+ NORM=SQRT(NORM)
+ DO N=1,NPQ
+ IF(IGLK.EQ.1) THEN
+ MND(1,N)=W(N)*RLM(N,1,1)
+ ELSEIF(IGLK.EQ.2) THEN
+ MND(N,1)=(2.0*L+1.0)/(4.0*PI)*RLM(N,1,1)
+ ELSE
+ CALL XABORT('UNKNOWN GALERKIN QUADRATURE METHOD.')
+ ENDIF
+ U(N,1)=RLM(N,1,1)/NORM
+ ENDDO
+ IND=1
+ ! ITERATE OVER THE SPHERICAL HARMONICS
+ DO 212 L=1,ISCAT-1
+ DO 213 M=-L,L
+ V2(:NPQ)=0.0D0
+ DO N=1,IND
+ IPROD=0.0D0
+ DO N2=1,NPQ
+ IPROD=IPROD+U(N2,N)*RLM(N2,L+1,M+L+1)
+ ENDDO
+ DO N2=1,NPQ
+ V2(N2)=V2(N2)+IPROD*U(N2,N)
+ ENDDO
+ ENDDO
+ V(:NPQ)=0.0D0
+ DO N=1,NPQ
+ V(N)=RLM(N,L+1,M+L+1)-V2(N)
+ ENDDO
+ NORM=0.0D0
+ DO N=1,NPQ
+ NORM=NORM+V(N)**2
+ ENDDO
+ NORM=SQRT(NORM)
+ ! KEEP THE SPHERICAL HARMONICS IF IT IS INDEPENDANT
+ IF(NORM.GE.1.0E-5) THEN
+ IND=IND+1
+ DO N=1,NPQ
+ U(N,IND)=V(N)/NORM
+ IF(IGLK.EQ.1) THEN
+ MND(IND,N)=W(N)*RLM(N,L+1,M+L+1)
+ ELSEIF(IGLK.EQ.2) THEN
+ MND(N,IND)=(2.0*L+1.0)/(4.0*PI)*RLM(N,L+1,M+L+1)
+ ELSE
+ CALL XABORT('UNKNOWN GALERKIN QUADRATURE METHOD.')
+ ENDIF
+ ENDDO
+ IL(IND)=L
+ IM(IND)=M
+ ENDIF
+ IF(IND.EQ.NPQ) GOTO 217
+ 213 ENDDO
+ 212 ENDDO
+ CALL XABORT('SNTT3D: THE'//
+ 1 ' GRAM-SCHMIDTH PROCEDURE TO FIND A SUITABLE INTERPOLATION'//
+ 2 ' BASIS REQUIRE HIGHER LEGENDRE ORDER.')
+ ! FIND INVERSE MATRIX
+ 217 IF(IGLK.EQ.1) THEN
+ DN=REAL(MND)
+ CALL ALINVD(NPQ,MND,NPQ,IER)
+ IF(IER.NE.0) CALL XABORT('SNTT3D: SINGULAR MATRIX.')
+ MN=REAL(MND)
+ ELSEIF(IGLK.EQ.2) THEN
+ MN=REAL(MND)
+ CALL ALINVD(NPQ,MND,NPQ,IER)
+ IF(IER.NE.0) CALL XABORT('SNTT3D: SINGULAR MATRIX.')
+ DN=REAL(MND)
+ ELSE
+ CALL XABORT('UNKNOWN GALERKIN QUADRATURE METHOD.')
+ ENDIF
+ DEALLOCATE(U,RLM,V,V2)
+ ELSE
+ IND=1
+ DO L=0,ISCAT-1
+ DO 218 M=-L,L
+ IL(IND)=L
+ IM(IND)=M
+ DO N=1,NPQ
+ DN(IND,N)=W(N)*PNSH(L,M,DU(N),DE(N),DZ(N))
+ MN(N,IND)=(2.0*L+1.0)/(4.0*PI)
+ 1 *PNSH(L,M,DU(N),DE(N),DZ(N))
+ ENDDO
+ IND=IND+1
+ 218 ENDDO
+ ENDDO
+ ENDIF
+*----
+* GENERATE THE WEIGHTING PARAMETERS OF THE CLOSURE RELATION.
+*----
+ PX=1
+ PE=1
+ IF(ISCHM.EQ.1.OR.ISCHM.EQ.3) THEN
+ PX=1
+ ELSEIF(ISCHM.EQ.2) THEN
+ PX=0
+ ELSE
+ CALL XABORT('SNTT3D: UNKNOWN TYPE OF SPATIAL CLOSURE RELATION.')
+ ENDIF
+ IF(MOD(IELEM,2).EQ.1) THEN
+ WX(1)=-PX
+ WX(2:IELEM+1:2)=1+PX
+ IF(IELEM.GE.2) WX(3:IELEM+1:2)=1-PX
+ ELSE
+ WX(1)=PX
+ WX(2:IELEM+1:2)=1-PX
+ IF(IELEM.GE.2) WX(3:IELEM+1:2)=1+PX
+ ENDIF
+ IF(IBFP.NE.0) THEN
+ IF(ESCHM.EQ.1.OR.ESCHM.EQ.3) THEN
+ PE=1
+ ELSEIF(ESCHM.EQ.2) THEN
+ PE=0
+ ELSE
+ CALL XABORT('SNTT3D: UNKNOWN TYPE OF ENERGY CLOSURE RELATION.')
+ ENDIF
+ IF(MOD(EELEM,2).EQ.1) THEN
+ WE(1)=-PE
+ WE(2:EELEM+1:2)=1+PE
+ IF(EELEM.GE.2) WE(3:EELEM+1:2)=1-PE
+ ELSE
+ WE(1)=PE
+ WE(2:EELEM+1:2)=1-PE
+ IF(EELEM.GE.2) WE(3:EELEM+1:2)=1+PE
+ ENDIF
+ ENDIF
+ ! NORMALIZED LEGENDRE POLYNOMIAL CONSTANTS
+ DO IEL=1,MAX(IELEM,EELEM)
+ CST(IEL)=SQRT(2.0*IEL-1.0)
+ ENDDO
+*----
+* COMPUTE ISOTROPIC FLUX INDICES.
+*----
+ NM=IELEM**3*EELEM
+ NMX=IELEM**2*EELEM
+ NMY=NMX
+ NMZ=NMX
+ NME=IELEM**3
+ LL4=LX*LY*LZ*NSCT*NM
+ IF(IGE.EQ.0)THEN
+ NUN=LL4+(NMX*LY*LZ+NMY*LX*LZ+NMZ*LX*LY)*NPQ
+ ELSEIF(IGE.EQ.2)THEN
+ IF((NCODE(5)==1).and.(NCODE(5)==NCODE(6)))THEN
+ NUN=LL4
+ ELSE
+ NUN=LL4+(LX*LY*LZ)*NMZ*NPQ
+ ENDIF
+ ELSE
+ CALL XABORT('SNTT3D: CHECK SPATIAL SCHEME DISCRETISATION '//
+ 1 'PARAMETER.')
+ ENDIF
+ DO I=1,LX*LY*LZ
+ IDL(I)=(I-1)*NSCT*NM+1
+ ENDDO
+*----
+* SET BOUNDARY CONDITIONS.
+*----
+ DO 240 I=1,6
+ IF(NCODE(I).NE.1) ZCODE(I)=1.0
+ IF(NCODE(I).EQ.5) CALL XABORT('SNTT3D: SYME BC NOT ALLOWED.')
+ IF(NCODE(I).EQ.7) CALL XABORT('SNTT3D: ZERO FLUX BC NOT ALLOWED.')
+ 240 CONTINUE
+*----
+* CHECK FOR INVALID VIRTUAL ELEMENTS.
+*----
+ DO 292 I=2,LX-1
+ DO 291 J=2,LY-1
+ DO 290 K=2,LZ-1
+ IF(MAT(I,J,K).EQ.0) THEN
+ L1=(NCODE(1).NE.1)
+ DO 251 J1=1,J-1
+ DO 250 K1=1,K-1
+ L1=L1.OR.(MAT(I,J1,K1).NE.0)
+ 250 CONTINUE
+ 251 CONTINUE
+ L2=(NCODE(2).NE.1)
+ DO 261 J1=J+1,LY
+ DO 260 K1=K+1,LZ
+ L2=L2.OR.(MAT(I,J1,K1).NE.0)
+ 260 CONTINUE
+ 261 CONTINUE
+ L3=(NCODE(3).NE.1)
+ DO 271 I1=1,I-1
+ DO 270 K1=1,K-1
+ L3=L3.OR.(MAT(I1,J,K1).NE.0)
+ 270 CONTINUE
+ 271 CONTINUE
+ L4=(NCODE(4).NE.1)
+ DO 281 I1=I+1,LX
+ DO 280 K1=K+1,LZ
+ L4=L4.OR.(MAT(I1,J,K1).NE.0)
+ 280 CONTINUE
+ 281 CONTINUE
+ L5=(NCODE(5).NE.1)
+ DO 301 I1=1,I-1
+ DO 300 J1=1,J-1
+ L5=L5.OR.(MAT(I1,J1,K).NE.0)
+ 300 CONTINUE
+ 301 CONTINUE
+ L6=(NCODE(6).NE.0)
+ DO 311 I1=I+1,LX
+ DO 310 J1=I+1,LY
+ L6=L6.OR.(MAT(I1,J1,K).NE.0)
+ 310 CONTINUE
+ 311 CONTINUE
+ IF(L1.AND.L2.AND.L3.AND.L4.AND.L5.AND.L6) THEN
+ WRITE(HSMG,'(17HSNTT3D: ELEMENT (,I3,1H,,I3,11H) CANNOT BE,
+ 1 9H VIRTUAL.)') I,J,K
+ CALL XABORT(HSMG)
+ ENDIF
+ ENDIF
+ 290 CONTINUE
+ 291 CONTINUE
+ 292 CONTINUE
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(ZZ,YY,XX)
+ RETURN
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 15:15:25 +0000