From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Trivac/src/PN3HWW.f | 560 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 560 insertions(+) create mode 100755 Trivac/src/PN3HWW.f (limited to 'Trivac/src/PN3HWW.f') diff --git a/Trivac/src/PN3HWW.f b/Trivac/src/PN3HWW.f new file mode 100755 index 0000000..bc49e10 --- /dev/null +++ b/Trivac/src/PN3HWW.f @@ -0,0 +1,560 @@ +*DECK PN3HWW + SUBROUTINE PN3HWW(NBMIX,NBLOS,IELEM,ICOL,NLF,NVD,NAN,LL4F,LL4W, + 1 MAT,SIGT,SIGTI,SIDE,ZZ,FRZ,QFR,IPERT,KN,MUW,IPBBW,LC,R,V,BBW, + 2 TTF,AW,C11W) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Assembly of system matrices for a Thomas-Raviart-Schneider (dual) +* finite element method in hexagonal 3-D simplified PN approximation. +* Note: system matrices should be initialized by the calling program. +* +*Copyright: +* Copyright (C) 2009 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 +* NBMIX number of mixtures. +* NBLOS number of lozenges per direction, taking into account +* mesh-splitting. +* IELEM degree of the Lagrangian finite elements: =1 (linear); +* =2 (parabolic); =3 (cubic). +* ICOL type of quadrature: =1 (analytical integration); +* =2 (Gauss-Lobatto); =3 (Gauss-Legendre). +* NLF number of Legendre orders for the flux (even number). +* NVD type of void boundary condition if NLF>0 and ICOL=3. +* NAN number of Legendre orders for the cross sections. +* LL4F number of flux components. +* LL4W number of W-directed currents. +* LL4X number of X-directed currents. +* LL4Y number of Y-directed currents. +* LL4Z number of Z-directed currents. +* MAT mixture index assigned to each lozenge. +* SIGT total minus self-scattering macroscopic cross sections. +* SIGT(:,NAN) generally contains the total cross section only. +* SIGTI inverse macroscopic cross sections ordered by mixture. +* SIGTI(:,NAN) generally contains the inverse total cross +* section only. +* SIDE side of an hexagon. +* ZZ Z-directed mesh spacings. +* FRZ volume fractions for the axial SYME boundary condition. +* QFR element-ordered boundary conditions. +* IPERT mixture permutation index. +* KN ADI permutation indices for the volumes and currents. +* MUW W-directed compressed storage mode indices. +* MUX X-directed compressed storage mode indices. +* MUY Y-directed compressed storage mode indices. +* MUZ Z-directed compressed storage mode indices. +* IPBBW W-directed perdue storage indices. +* IPBBX X-directed perdue storage indices. +* IPBBY Y-directed perdue storage indices. +* IPBBZ Z-directed perdue storage indices. +* LC order of the unit matrices. +* R unit matrix. +* V unit matrix. +* BBW W-directed flux-current matrices. +* BBX X-directed flux-current matrices. +* BBY Y-directed flux-current matrices. +* BBZ Z-directed flux-current matrices. +* +*Parameters: output +* TTF flux-flux matrices. +* AW W-directed main current-current matrices. Dimensionned to +* MUW(LL4W)*NLF/2. +* AX X-directed main current-current matrices. Dimensionned to +* MUX(LL4X)*NLF/2. +* AY Y-directed main current-current matrices. Dimensionned to +* MUY(LL4Y)*NLF/2. +* AZ Z-directed main current-current matrices. Dimensionned to +* MUZ(LL4Z)*NLF/2. +* C11W W-directed main current-current matrices to be factorized. +* Dimensionned to MUW(LL4W)*NLF/2. +* C11X X-directed main current-current matrices to be factorized. +* Dimensionned to MUX(LL4X)*NLF/2. +* C11Y Y-directed main current-current matrices to be factorized. +* Dimensionned to MUY(LL4Y)*NLF/2. +* C11Z Z-directed main current-current matrices to be factorized. +* Dimensionned to MUZ(LL4Z)*NLF/2. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NBMIX,NBLOS,IELEM,ICOL,NLF,NVD,NAN,LL4F,LL4W, + 1 MAT(3,NBLOS),MUW(LL4W),IPBBW(2*IELEM,LL4W),LC,IPERT(NBLOS), + 2 KN(NBLOS,3+6*(IELEM+2)*IELEM**2) + REAL SIGT(NBMIX,NAN),SIGTI(NBMIX,NAN),SIDE,ZZ(3,NBLOS),FRZ(NBLOS), + 1 QFR(NBLOS,8),R(LC,LC),V(LC,LC-1),BBW(2*IELEM,LL4W), + 2 TTF(LL4F*NLF/2),AW(*),C11W(*) +*---- +* LOCAL VARIABLES +*---- + REAL QQ(5,5) + DOUBLE PRECISION FFF,TTTT,VOL0,GARS,GARSI,FACT,VAR1 +*---- +* W-ORIENTED COUPLINGS +*---- + ZMARS=0.0 + IF(ICOL.EQ.3) THEN + IF(NVD.EQ.0) THEN + NZMAR=NLF+1 + ELSE IF(NVD.EQ.1) THEN + NZMAR=NLF + ELSE IF(NVD.EQ.2) THEN + NZMAR=65 + ENDIF + ELSE + NZMAR=65 + ENDIF + DO 25 I0=1,IELEM + DO 20 J0=1,IELEM + FFF=0.0D0 + DO 10 K0=2,IELEM + FFF=FFF+V(K0,I0)*V(K0,J0)/R(K0,K0) + 10 CONTINUE + IF(ABS(FFF).LE.1.0E-6) FFF=0.0D0 + QQ(I0,J0)=REAL(FFF) + 20 CONTINUE + 25 CONTINUE + MUMAX=MUW(LL4W) + DO 120 IL=0,NLF-1 + IF(MOD(IL,2).EQ.1) ZMARS=PNMAR2(NZMAR,IL,IL) + FACT=REAL(2*IL+1) +*---- +* ASSEMBLY OF THE W-ORIENTED COEFFICIENT MATRICES AT ORDER IL. +*---- + NELEH=(IELEM+1)*IELEM**2 + TTTT=0.5D0*SQRT(3.D00)*SIDE*SIDE + NUM=0 + DO 70 KEL=1,NBLOS + IF(IPERT(KEL).EQ.0) GO TO 70 + IBM=MAT(1,IPERT(KEL)) + NUM=NUM+1 + IF(IBM.EQ.0) GO TO 70 + DZ=ZZ(1,IPERT(KEL)) + VOL0=TTTT*DZ*FRZ(KEL) + GARS=SIGT(IBM,MIN(IL+1,NAN)) + IF(MOD(IL,2).EQ.0) THEN +* EVEN PARITY EQUATION. + VAR1=FACT*VOL0*GARS + DO 32 K3=0,IELEM-1 + DO 31 K2=0,IELEM-1 + DO 30 K1=0,IELEM-1 + IOF=(IL/2)*LL4F + JND1=IOF+(((NUM-1)*IELEM+K3)*IELEM+K2)*IELEM+K1+1 + JND2=IOF+(((KN(NUM,1)-1)*IELEM+K3)*IELEM+K2)*IELEM+K1+1 + JND3=IOF+(((KN(NUM,2)-1)*IELEM+K3)*IELEM+K2)*IELEM+K1+1 + TTF(JND1)=TTF(JND1)+REAL(VAR1) + TTF(JND2)=TTF(JND2)+REAL(VAR1) + TTF(JND3)=TTF(JND3)+REAL(VAR1) + 30 CONTINUE + 31 CONTINUE + 32 CONTINUE + ELSE +* PARTIAL INVERSION OF THE ODD PARITY EQUATION. MODIFICATION OF +* THE EVEN PARITY EQUATION. + GARSI=SIGTI(IBM,MIN(IL+1,NAN)) + IF(IELEM.GT.1) THEN + KOFF=((IL-1)/2)*LL4F + DO 42 K3=0,IELEM-1 + DO 41 K2=0,IELEM-1 + DO 40 K1=0,IELEM-1 + JND1=KOFF+(((NUM-1)*IELEM+K3)*IELEM+K2)*IELEM+K1+1 + JND2=KOFF+(((KN(NUM,1)-1)*IELEM+K3)*IELEM+K2)*IELEM+K1+1 + JND3=KOFF+(((KN(NUM,2)-1)*IELEM+K3)*IELEM+K2)*IELEM+K1+1 + VAR1=(REAL(IL)**2)*VOL0*QQ(K3+1,K3+1)*GARSI/(FACT*DZ*DZ) + TTF(JND1)=TTF(JND1)+REAL(VAR1) + TTF(JND2)=TTF(JND2)+REAL(VAR1) + TTF(JND3)=TTF(JND3)+REAL(VAR1) + IF(IL.LE.NLF-3) THEN + JND1=JND1+LL4F + JND2=JND2+LL4F + JND3=JND3+LL4F + VAR1=(REAL(IL+1)**2)*VOL0*QQ(K3+1,K3+1)*GARSI/(FACT*DZ*DZ) + TTF(JND1)=TTF(JND1)+REAL(VAR1) + TTF(JND2)=TTF(JND2)+REAL(VAR1) + TTF(JND3)=TTF(JND3)+REAL(VAR1) + ENDIF + 40 CONTINUE + 41 CONTINUE + 42 CONTINUE + ENDIF +* +* ODD PARITY EQUATION. + DO 63 K5=0,1 ! TWO LOZENGES PER HEXAGON + DO 62 K4=0,IELEM-1 + DO 61 K3=0,IELEM-1 + DO 60 K2=1,IELEM+1 + KNW1=KN(NUM,3+K5*NELEH+(K4*IELEM+K3)*(IELEM+1)+K2) + INW1=ABS(KNW1) + DO 50 K1=1,IELEM+1 + KNW2=KN(NUM,3+K5*NELEH+(K4*IELEM+K3)*(IELEM+1)+K1) + INW2=ABS(KNW2) + IF((KNW2.NE.0).AND.(KNW1.NE.0).AND.(INW1.GE.INW2)) THEN + KEY=((IL-1)/2)*MUMAX+MUW(INW1)-INW1+INW2 + SG=REAL(SIGN(1,KNW1)*SIGN(1,KNW2)) + VAR1=(4./3.)*SG*FACT*VOL0*GARS*R(K2,K1) + AW(KEY)=AW(KEY)-REAL(VAR1) + ENDIF + 50 CONTINUE + IF(KNW1.NE.0) THEN + KEY=((IL-1)/2)*MUMAX+MUW(INW1) + IF((K2.EQ.1).AND.(K5.EQ.0)) THEN + VAR1=0.5*FACT*QFR(NUM,1)*ZMARS + AW(KEY)=AW(KEY)-REAL(VAR1) + ELSE IF((K2.EQ.IELEM+1).AND.(K5.EQ.1)) THEN + VAR1=0.5*FACT*QFR(NUM,2)*ZMARS + AW(KEY)=AW(KEY)-REAL(VAR1) + ENDIF + ENDIF + 60 CONTINUE + 61 CONTINUE + 62 CONTINUE + 63 CONTINUE + ENDIF + 70 CONTINUE +* + IF(MOD(IL,2).EQ.1) THEN + DO 80 I0=1,MUMAX + C11W(((IL-1)/2)*MUMAX+I0)=-AW(((IL-1)/2)*MUMAX+I0) + 80 CONTINUE + MUIM1=0 + DO 110 I=1,LL4W + MUI=MUW(I) + DO 100 J=I-(MUI-MUIM1)+1,I + KEY=((IL-1)/2)*MUMAX+(MUI-I+J) + DO 95 I0=1,2*IELEM + II=IPBBW(I0,I) + IF(II.EQ.0) GO TO 100 + DO 90 J0=1,2*IELEM + JJ=IPBBW(J0,J) + IF(II.EQ.JJ) C11W(KEY)=C11W(KEY)+REAL(IL**2)*BBW(I0,I)* + 1 BBW(J0,J)/TTF(((IL-1)/2)*LL4F+II) + 90 CONTINUE + 95 CONTINUE + 100 CONTINUE + MUIM1=MUI + 110 CONTINUE + ENDIF + 120 CONTINUE + RETURN + END +* + SUBROUTINE PN3HWX(NBMIX,NBLOS,IELEM,ICOL,NLF,NVD,NAN,LL4F,LL4W, + 1 LL4X,MAT,SIGT,SIDE,ZZ,FRZ,QFR,IPERT,KN,MUX,IPBBX,LC,R,BBX,TTF, + 2 AX,C11X) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NBMIX,NBLOS,IELEM,ICOL,NLF,NVD,NAN,LL4F,LL4W,LL4X, + 1 MAT(3,NBLOS),MUX(LL4X),IPBBX(2*IELEM,LL4X),LC,IPERT(NBLOS), + 2 KN(NBLOS,3+6*(IELEM+2)*IELEM**2) + REAL SIGT(NBMIX,NAN),SIDE,ZZ(3,NBLOS),FRZ(NBLOS),QFR(NBLOS,8), + 1 R(LC,LC),BBX(2*IELEM,LL4X),TTF(LL4F*NLF/2),AX(*),C11X(*) +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION TTTT,VOL0,GARS,FACT,VAR1 +*---- +* X-ORIENTED COUPLINGS +*---- + ZMARS=0.0 + IF(ICOL.EQ.3) THEN + IF(NVD.EQ.0) THEN + NZMAR=NLF+1 + ELSE IF(NVD.EQ.1) THEN + NZMAR=NLF + ELSE IF(NVD.EQ.2) THEN + NZMAR=65 + ENDIF + ELSE + NZMAR=65 + ENDIF + MUMAX=MUX(LL4X) + DO 200 IL=1,NLF-1,2 + ZMARS=PNMAR2(NZMAR,IL,IL) + FACT=REAL(2*IL+1) +*---- +* ASSEMBLY OF THE X-ORIENTED COEFFICIENT MATRICES AT ODD ORDER IL. +*---- + NELEH=(IELEM+1)*IELEM**2 + TTTT=0.5D0*SQRT(3.D00)*SIDE*SIDE + NUM=0 + DO 150 KEL=1,NBLOS + IF(IPERT(KEL).EQ.0) GO TO 150 + IBM=MAT(1,IPERT(KEL)) + NUM=NUM+1 + IF(IBM.EQ.0) GO TO 150 + VOL0=TTTT*ZZ(1,IPERT(KEL))*FRZ(KEL) + GARS=SIGT(IBM,MIN(IL+1,NAN)) +* + DO 143 K5=0,1 ! TWO LOZENGES PER HEXAGON + DO 142 K4=0,IELEM-1 + DO 141 K3=0,IELEM-1 + DO 140 K2=1,IELEM+1 + KNX1=KN(NUM,3+(K5+2)*NELEH+(K4*IELEM+K3)*(IELEM+1)+K2) + INX1=ABS(KNX1)-LL4W + DO 130 K1=1,IELEM+1 + KNX2=KN(NUM,3+(K5+2)*NELEH+(K4*IELEM+K3)*(IELEM+1)+K1) + INX2=ABS(KNX2)-LL4W + IF((KNX2.NE.0).AND.(KNX1.NE.0).AND.(INX1.GE.INX2)) THEN + KEY=((IL-1)/2)*MUMAX+MUX(INX1)-INX1+INX2 + SG=REAL(SIGN(1,KNX1)*SIGN(1,KNX2)) + VAR1=(4./3.)*SG*FACT*VOL0*GARS*R(K2,K1) + AX(KEY)=AX(KEY)-REAL(VAR1) + ENDIF + 130 CONTINUE + IF(KNX1.NE.0) THEN + KEY=((IL-1)/2)*MUMAX+MUX(INX1) + IF((K2.EQ.1).AND.(K5.EQ.0)) THEN + VAR1=0.5*FACT*QFR(NUM,3)*ZMARS + AX(KEY)=AX(KEY)-REAL(VAR1) + ELSE IF((K2.EQ.IELEM+1).AND.(K5.EQ.1)) THEN + VAR1=0.5*FACT*QFR(NUM,4)*ZMARS + AX(KEY)=AX(KEY)-REAL(VAR1) + ENDIF + ENDIF + 140 CONTINUE + 141 CONTINUE + 142 CONTINUE + 143 CONTINUE + 150 CONTINUE +* + DO 160 I0=1,MUMAX + C11X(((IL-1)/2)*MUMAX+I0)=-AX(((IL-1)/2)*MUMAX+I0) + 160 CONTINUE + MUIM1=0 + DO 190 I=1,LL4X + MUI=MUX(I) + DO 180 J=I-(MUI-MUIM1)+1,I + KEY=((IL-1)/2)*MUMAX+(MUI-I+J) + DO 175 I0=1,2*IELEM + II=IPBBX(I0,I) + IF(II.EQ.0) GO TO 180 + DO 170 J0=1,2*IELEM + JJ=IPBBX(J0,J) + IF(II.EQ.JJ) THEN + VAR1=REAL(IL**2)*BBX(I0,I)*BBX(J0,J)/TTF(((IL-1)/2)*LL4F+II) + C11X(KEY)=C11X(KEY)+REAL(VAR1) + ENDIF + 170 CONTINUE + 175 CONTINUE + 180 CONTINUE + MUIM1=MUI + 190 CONTINUE + 200 CONTINUE + RETURN + END +* + SUBROUTINE PN3HWY(NBMIX,NBLOS,IELEM,ICOL,NLF,NVD,NAN,LL4F,LL4W, + 1 LL4X,LL4Y,MAT,SIGT,SIDE,ZZ,FRZ,QFR,IPERT,KN,MUY,IPBBY,LC,R,BBY, + 2 TTF,AY,C11Y) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NBMIX,NBLOS,IELEM,ICOL,NLF,NVD,NAN,LL4F,LL4W,LL4X,LL4Y, + 1 MAT(3,NBLOS),MUY(LL4Y),IPBBY(2*IELEM,LL4Y),LC,IPERT(NBLOS), + 2 KN(NBLOS,3+6*(IELEM+2)*IELEM**2) + REAL SIGT(NBMIX,NAN),SIDE,ZZ(3,NBLOS),FRZ(NBLOS),QFR(NBLOS,8), + 1 R(LC,LC),BBY(2*IELEM,LL4Y),TTF(LL4F*NLF/2),AY(*),C11Y(*) +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION TTTT,VOL0,GARS,FACT,VAR1 +*---- +* Y-ORIENTED COUPLINGS +*---- + ZMARS=0.0 + IF(ICOL.EQ.3) THEN + IF(NVD.EQ.0) THEN + NZMAR=NLF+1 + ELSE IF(NVD.EQ.1) THEN + NZMAR=NLF + ELSE IF(NVD.EQ.2) THEN + NZMAR=65 + ENDIF + ELSE + NZMAR=65 + ENDIF + MUMAX=MUY(LL4Y) + DO 280 IL=1,NLF-1,2 + ZMARS=PNMAR2(NZMAR,IL,IL) + FACT=REAL(2*IL+1) +*---- +* ASSEMBLY OF THE Y-ORIENTED COEFFICIENT MATRICES AT ODD ORDER IL. +*---- + NELEH=(IELEM+1)*IELEM**2 + TTTT=0.5D0*SQRT(3.D00)*SIDE*SIDE + NUM=0 + DO 230 KEL=1,NBLOS + IF(IPERT(KEL).EQ.0) GO TO 230 + IBM=MAT(1,IPERT(KEL)) + NUM=NUM+1 + IF(IBM.EQ.0) GO TO 230 + VOL0=TTTT*ZZ(1,IPERT(KEL))*FRZ(KEL) + GARS=SIGT(IBM,MIN(IL+1,NAN)) +* + DO 223 K5=0,1 ! TWO LOZENGES PER HEXAGON + DO 222 K4=0,IELEM-1 + DO 221 K3=0,IELEM-1 + DO 220 K2=1,IELEM+1 + KNY1=KN(NUM,3+(K5+4)*NELEH+(K4*IELEM+K3)*(IELEM+1)+K2) + INY1=ABS(KNY1)-LL4W-LL4X + DO 210 K1=1,IELEM+1 + KNY2=KN(NUM,3+(K5+4)*NELEH+(K4*IELEM+K3)*(IELEM+1)+K1) + INY2=ABS(KNY2)-LL4W-LL4X + IF((KNY2.NE.0).AND.(KNY1.NE.0).AND.(INY1.GE.INY2)) THEN + KEY=((IL-1)/2)*MUMAX+MUY(INY1)-INY1+INY2 + SG=REAL(SIGN(1,KNY1)*SIGN(1,KNY2)) + VAR1=(4./3.)*SG*FACT*VOL0*GARS*R(K2,K1) + AY(KEY)=AY(KEY)-REAL(VAR1) + ENDIF + 210 CONTINUE + IF(KNY1.NE.0) THEN + KEY=((IL-1)/2)*MUMAX+MUY(INY1) + IF((K2.EQ.1).AND.(K5.EQ.0)) THEN + VAR1=0.5*FACT*QFR(NUM,5)*ZMARS + AY(KEY)=AY(KEY)-REAL(VAR1) + ELSE IF((K2.EQ.IELEM+1).AND.(K5.EQ.1)) THEN + VAR1=0.5*FACT*QFR(NUM,6)*ZMARS + AY(KEY)=AY(KEY)-REAL(VAR1) + ENDIF + ENDIF + 220 CONTINUE + 221 CONTINUE + 222 CONTINUE + 223 CONTINUE + 230 CONTINUE +* + DO 240 I0=1,MUMAX + C11Y(((IL-1)/2)*MUMAX+I0)=-AY(((IL-1)/2)*MUMAX+I0) + 240 CONTINUE + MUIM1=0 + DO 270 I=1,LL4Y + MUI=MUY(I) + DO 260 J=I-(MUI-MUIM1)+1,I + KEY=((IL-1)/2)*MUMAX+(MUI-I+J) + DO 255 I0=1,2*IELEM + II=IPBBY(I0,I) + IF(II.EQ.0) GO TO 260 + DO 250 J0=1,2*IELEM + JJ=IPBBY(J0,J) + IF(II.EQ.JJ) C11Y(KEY)=C11Y(KEY)+REAL(IL**2)*BBY(I0,I)* + 1 BBY(J0,J)/TTF(((IL-1)/2)*LL4F+II) + 250 CONTINUE + 255 CONTINUE + 260 CONTINUE + MUIM1=MUI + 270 CONTINUE + 280 CONTINUE + RETURN + END +* + SUBROUTINE PN3HWZ(NBMIX,NBLOS,IELEM,ICOL,NLF,NVD,NAN,LL4F,LL4W, + 1 LL4X,LL4Y,LL4Z,MAT,SIGT,SIDE,ZZ,FRZ,QFR,IPERT,KN,MUZ,IPBBZ,LC, + 2 R,BBZ,TTF,AZ,C11Z) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NBMIX,NBLOS,IELEM,ICOL,NLF,NVD,NAN,LL4F,LL4W,LL4X, + 1 LL4Y,LL4Z,MAT(3,NBLOS),MUZ(LL4Z),IPBBZ(2*IELEM,LL4Z),LC, + 2 IPERT(NBLOS),KN(NBLOS,3+6*(IELEM+2)*IELEM**2) + REAL SIGT(NBMIX,NAN),SIDE,ZZ(3,NBLOS),FRZ(NBLOS),QFR(NBLOS,8), + 1 R(LC,LC),BBZ(2*IELEM,LL4Z),TTF(LL4F*NLF/2),AZ(*),C11Z(*) +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION TTTT,VOL0,GARS,FACT,VAR1 +*---- +* Z-ORIENTED COUPLINGS +*---- + ZMARS=0.0 + IF(ICOL.EQ.3) THEN + IF(NVD.EQ.0) THEN + NZMAR=NLF+1 + ELSE IF(NVD.EQ.1) THEN + NZMAR=NLF + ELSE IF(NVD.EQ.2) THEN + NZMAR=65 + ENDIF + ELSE + NZMAR=65 + ENDIF + MUMAX=MUZ(LL4Z) + DO 360 IL=1,NLF-1,2 + ZMARS=PNMAR2(NZMAR,IL,IL) + FACT=REAL(2*IL+1) +*---- +* ASSEMBLY OF THE Z-ORIENTED COEFFICIENT MATRICES AT ODD ORDER IL. +*---- + NELEH=(IELEM+1)*IELEM**2 + TTTT=0.5D0*SQRT(3.D00)*SIDE*SIDE + NUM=0 + DO 310 KEL=1,NBLOS + IF(IPERT(KEL).EQ.0) GO TO 310 + IBM=MAT(1,IPERT(KEL)) + IF(IBM.EQ.0) GO TO 310 + NUM=NUM+1 + VOL0=TTTT*ZZ(1,IPERT(KEL))*FRZ(KEL) + GARS=SIGT(IBM,MIN(IL+1,NAN)) +* + DO 302 K5=0,2 ! THREE LOZENGES PER HEXAGON + DO 301 K2=0,IELEM-1 + DO 300 K1=0,IELEM-1 + KNZ1=KN(NUM,3+6*NELEH+2*K5*IELEM**2+K2*IELEM+K1+1) + INZ1=ABS(KNZ1)-LL4W-LL4X-LL4Y + KNZ2=KN(NUM,3+6*NELEH+(2*K5+1)*IELEM**2+K2*IELEM+K1+1) + INZ2=ABS(KNZ2)-LL4W-LL4X-LL4Y + IF(KNZ1.NE.0) THEN + KEY=((IL-1)/2)*MUMAX+MUZ(INZ1) + VAR1=FACT*VOL0*GARS*R(1,1)+0.5*FACT*QFR(NUM,7)*ZMARS + AZ(KEY)=AZ(KEY)-REAL(VAR1) + ENDIF + IF(KNZ2.NE.0) THEN + KEY=((IL-1)/2)*MUMAX+MUZ(INZ2) + VAR1=FACT*VOL0*GARS*R(IELEM+1,IELEM+1)+0.5*FACT*QFR(NUM,8)*ZMARS + AZ(KEY)=AZ(KEY)-REAL(VAR1) + ENDIF + IF((ICOL.NE.2).AND.(KNZ1.NE.0).AND.(KNZ2.NE.0)) THEN + IF(INZ2.GT.INZ1) KEY=((IL-1)/2)*MUMAX+MUZ(INZ2)-INZ2+INZ1 + IF(INZ2.LE.INZ1) KEY=((IL-1)/2)*MUMAX+MUZ(INZ1)-INZ1+INZ2 + SG=REAL(SIGN(1,KNZ1)*SIGN(1,KNZ2)) + IF(INZ1.EQ.INZ2) SG=2.0*SG + VAR1=SG*FACT*VOL0*GARS*R(IELEM+1,1) + AZ(KEY)=AZ(KEY)-REAL(VAR1) + ENDIF + 300 CONTINUE + 301 CONTINUE + 302 CONTINUE + 310 CONTINUE +* + DO 320 I0=1,MUMAX + C11Z(((IL-1)/2)*MUMAX+I0)=-AZ(((IL-1)/2)*MUMAX+I0) + 320 CONTINUE + MUIM1=0 + DO 350 I=1,LL4Z + MUI=MUZ(I) + DO 340 J=I-(MUI-MUIM1)+1,I + KEY=((IL-1)/2)*MUMAX+(MUI-I+J) + DO 335 I0=1,2*IELEM + II=IPBBZ(I0,I) + IF(II.EQ.0) GO TO 340 + DO 330 J0=1,2*IELEM + JJ=IPBBZ(J0,J) + IF(II.EQ.JJ) C11Z(KEY)=C11Z(KEY)+REAL(IL**2)*BBZ(I0,I)* + 1 BBZ(J0,J)/TTF(((IL-1)/2)*LL4F+II) + 330 CONTINUE + 335 CONTINUE + 340 CONTINUE + MUIM1=MUI + 350 CONTINUE + 360 CONTINUE + RETURN + END -- cgit v1.2.3