diff options
Diffstat (limited to 'Trivac/src/TRIPXX.f')
| -rwxr-xr-x | Trivac/src/TRIPXX.f | 381 |
1 files changed, 381 insertions, 0 deletions
diff --git a/Trivac/src/TRIPXX.f b/Trivac/src/TRIPXX.f new file mode 100755 index 0000000..d801d3a --- /dev/null +++ b/Trivac/src/TRIPXX.f @@ -0,0 +1,381 @@ +*DECK TRIPXX + SUBROUTINE TRIPXX (IR,MAXKN,NEL,LL4,VOL,MAT,XSGD,XX,YY,ZZ,DD,KN, + 1 QFR,MUX,IPX,CYLIND,LC,T,TS,Q,QS,A11X) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Assembly of system matrices for a primal finite element method in 3-D. +* +*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 +* IR first dimension for matrix SGD. +* MAXKN first dimension for matrix KN. +* NEL total number of finite elements. +* LL4 order of system matrices. +* MAT mixture index assigned to each element. +* VOL volume of each element. +* XX X-directed mesh spacings. +* YY Y-directed mesh spacings. +* ZZ Z-directed mesh spacings. +* DD values used with a cylindrical geometry. +* KN element-ordered unknown list. +* QFR element-ordered boundary conditions. +* XSGD nuclear properties, derivatives or first variations of +* nuclear properties per material mixture: +* XSGD(L,1): X-oriented diffusion coefficients; +* XSGD(L,2): Y-oriented diffusion coefficients; +* XSGD(L,3): Z-oriented diffusion coefficients; +* XSGD(L,4): removal macroscopic cross section. +* MUX X-oriented compressed storage mode indices. +* MUY Y-oriented compressed storage mode indices. +* MUZ Z-oriented compressed storage mode indices. +* IPX X-oriented permutation matrices. +* IPY Y-oriented permutation matrices. +* IPZ Z-oriented permutation matrices. +* CYLIND cylinderization flag (=.true. for cylindrical geometry). +* LC order of the unit matrices. +* T cartesian linear product vector. +* TS cylindrical linear product vector. +* Q cartesian stiffness matrix. +* QS cylindrical stiffness matrix. +* +*Parameters: output +* A11X X-oriented matrix corresponding to the divergence (i.e +* leakage) and removal terms (should be initialized by the +* calling program). +* A11Y Y-oriented matrix corresponding to the divergence (i.e +* leakage) and removal terms (should be initialized by the +* calling program). +* A11Z Z-oriented matrix corresponding to the divergence (i.e +* leakage) and removal terms (should be initialized by the +* calling program). +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER IR,MAXKN,NEL,LL4,MAT(NEL),KN(MAXKN),MUX(LL4),IPX(LL4),LC + REAL VOL(NEL),XSGD(IR,4),XX(NEL),YY(NEL),ZZ(NEL),DD(NEL), + 1 QFR(6*NEL),T(LC),TS(LC),Q(LC,LC),QS(LC,LC),A11X(*) + LOGICAL CYLIND +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION VAR1,VOL1,VOL2,VOL3,QQX,QQY,QQZ + COMMON /ELEM2/LL,LCC,IJ1(125),IJ2(125),IJ3(125),ISR(6,25), + 1 Q3DP1(125,125),Q3DP2(125,125),Q3DP3(125,125),R3DP(125), + 2 Q3DC1(125,125),Q3DC2(125,125),Q3DC3(125,125),R3DC(125), + 3 R2DP(25),R2DC(25) +*---- +* X-DIRECTED COUPLINGS. +* +* ASSEMBLY OF MATRIX A11X. +*---- + CALL TRIPMA(LC,T,TS,Q,QS) + NUM1=0 + NUM2=0 + DO 90 K=1,NEL + L=MAT(K) + IF(L.EQ.0) GO TO 90 + VOL0=VOL(K) + IF(VOL0.EQ.0.0) GO TO 80 + DX=XX(K) + DY=YY(K) + DZ=ZZ(K) + VOL1=VOL0/(DX*DX) + VOL2=VOL0/(DY*DY) + VOL3=VOL0/(DZ*DZ) + DO 50 I=1,LL + IND1=KN(NUM1+I) + IF(IND1.EQ.0) GO TO 50 + INX1=IPX(IND1) + KEY0=MUX(INX1) + IF(CYLIND) THEN + RR=(R3DP(I)+R3DC(I)*DX/DD(K))*VOL0 + ELSE + RR=R3DP(I)*VOL0 + ENDIF + A11X(KEY0)=A11X(KEY0)+RR*XSGD(L,4) + KEY0=KEY0-INX1 + DO 40 J=1,LL + IND2=KN(NUM1+J) + IF(IND2.EQ.0) GO TO 40 + INX2=IPX(IND2) + IF(INX2.EQ.INX1) THEN + IF(CYLIND) THEN + QQX=(Q3DP1(I,J)+Q3DC1(I,J)*DX/DD(K))*VOL1 + QQY=(Q3DP2(I,J)+Q3DC2(I,J)*DX/DD(K))*VOL2 + QQZ=(Q3DP3(I,J)+Q3DC3(I,J)*DX/DD(K))*VOL3 + ELSE + QQX=Q3DP1(I,J)*VOL1 + QQY=Q3DP2(I,J)*VOL2 + QQZ=Q3DP3(I,J)*VOL3 + ENDIF + KEY=KEY0+INX2 + VAR1=QQX*XSGD(L,1)+QQY*XSGD(L,2)+QQZ*XSGD(L,3) + A11X(KEY)=REAL(A11X(KEY)+VAR1) + ELSE IF((INX2.LT.INX1).AND.(IJ2(I).EQ.IJ2(J)).AND. + 1 (IJ3(I).EQ.IJ3(J))) THEN + IF(CYLIND) THEN + QQX=(Q3DP1(I,J)+Q3DC1(I,J)*DX/DD(K))*VOL1 + ELSE + QQX=Q3DP1(I,J)*VOL1 + ENDIF + KEY=KEY0+INX2 + A11X(KEY)=REAL(A11X(KEY)+QQX*XSGD(L,1)) + ENDIF + 40 CONTINUE + 50 CONTINUE + DO 70 IC=1,6 + QFR1=QFR(NUM2+IC) + IF(QFR1.EQ.0.0) GO TO 70 + DO 60 I1=1,LCC + I=ISR(IC,I1) + IND1=KN(NUM1+I) + IF(IND1.EQ.0) GO TO 60 + INX1=IPX(IND1) + KEY=MUX(INX1) + IF(CYLIND) THEN + IF(IC.EQ.1) THEN + CRZ=-0.5*R2DP(I1) + ELSE IF(IC.EQ.2) THEN + CRZ=0.5*R2DP(I1) + ELSE + CRZ=R2DC(I1) + ENDIF + RR=(R2DP(I1)+CRZ*DX/DD(K)) + ELSE + RR=R2DP(I1) + ENDIF + A11X(KEY)=A11X(KEY)+RR*QFR1 + 60 CONTINUE + 70 CONTINUE + 80 NUM1=NUM1+LL + NUM2=NUM2+6 + 90 CONTINUE + RETURN + END +* + SUBROUTINE TRIPXY (IR,MAXKN,NEL,LL4,VOL,MAT,XSGD,XX,YY,ZZ,DD,KN, + 1 QFR,MUY,IPY,CYLIND,LC,T,TS,Q,QS,A11Y) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER IR,MAXKN,NEL,LL4,MAT(NEL),KN(MAXKN),MUY(LL4),IPY(LL4),LC + REAL VOL(NEL),XSGD(IR,4),XX(NEL),YY(NEL),ZZ(NEL),DD(NEL), + 1 QFR(6*NEL),T(LC),TS(LC),Q(LC,LC),QS(LC,LC),A11Y(*) + LOGICAL CYLIND +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION VAR1,VOL1,VOL2,VOL3,QQX,QQY,QQZ + COMMON /ELEM2/LL,LCC,IJ1(125),IJ2(125),IJ3(125),ISR(6,25), + 1 Q3DP1(125,125),Q3DP2(125,125),Q3DP3(125,125),R3DP(125), + 2 Q3DC1(125,125),Q3DC2(125,125),Q3DC3(125,125),R3DC(125), + 3 R2DP(25),R2DC(25) +*---- +* Y-DIRECTED COUPLINGS. +* +* ASSEMBLY OF MATRIX A11Y. +*---- + CALL TRIPMA(LC,T,TS,Q,QS) + NUM1=0 + NUM2=0 + DO 180 K=1,NEL + L=MAT(K) + IF(L.EQ.0) GO TO 180 + VOL0=VOL(K) + IF(VOL0.EQ.0.0) GO TO 170 + DX=XX(K) + DY=YY(K) + DZ=ZZ(K) + VOL1=VOL0/(DX*DX) + VOL2=VOL0/(DY*DY) + VOL3=VOL0/(DZ*DZ) + DO 140 I=1,LL + IND1=KN(NUM1+I) + IF(IND1.EQ.0) GO TO 140 + INY1=IPY(IND1) + KEY0=MUY(INY1) + IF(CYLIND) THEN + RR=(R3DP(I)+R3DC(I)*DX/DD(K))*VOL0 + ELSE + RR=R3DP(I)*VOL0 + ENDIF + A11Y(KEY0)=A11Y(KEY0)+RR*XSGD(L,4) + KEY0=KEY0-INY1 + DO 130 J=1,LL + IND2=KN(NUM1+J) + IF(IND2.EQ.0) GO TO 130 + INY2=IPY(IND2) + IF(INY2.EQ.INY1) THEN + IF(CYLIND) THEN + QQX=(Q3DP1(I,J)+Q3DC1(I,J)*DX/DD(K))*VOL1 + QQY=(Q3DP2(I,J)+Q3DC2(I,J)*DX/DD(K))*VOL2 + QQZ=(Q3DP3(I,J)+Q3DC3(I,J)*DX/DD(K))*VOL3 + ELSE + QQX=Q3DP1(I,J)*VOL1 + QQY=Q3DP2(I,J)*VOL2 + QQZ=Q3DP3(I,J)*VOL3 + ENDIF + KEY=KEY0+INY2 + VAR1=QQX*XSGD(L,1)+QQY*XSGD(L,2)+QQZ*XSGD(L,3) + A11Y(KEY)=REAL(A11Y(KEY)+VAR1) + ELSE IF((INY2.LT.INY1).AND.(IJ1(I).EQ.IJ1(J)).AND. + 1 (IJ3(I).EQ.IJ3(J))) THEN + IF(CYLIND) THEN + QQY=(Q3DP2(I,J)+Q3DC2(I,J)*DX/DD(K))*VOL2 + ELSE + QQY=Q3DP2(I,J)*VOL2 + ENDIF + KEY=KEY0+INY2 + A11Y(KEY)=REAL(A11Y(KEY)+QQY*XSGD(L,2)) + ENDIF + 130 CONTINUE + 140 CONTINUE + DO 160 IC=1,6 + QFR1=QFR(NUM2+IC) + IF(QFR1.EQ.0.0) GO TO 160 + DO 150 I1=1,LCC + I=ISR(IC,I1) + IND1=KN(NUM1+I) + IF(IND1.EQ.0) GO TO 150 + INY1=IPY(IND1) + KEY=MUY(INY1) + IF(CYLIND) THEN + IF(IC.EQ.1) THEN + CRZ=-0.5*R2DP(I1) + ELSE IF(IC.EQ.2) THEN + CRZ=0.5*R2DP(I1) + ELSE + CRZ=R2DC(I1) + ENDIF + RR=(R2DP(I1)+DX*CRZ/DD(K)) + ELSE + RR=R2DP(I1) + ENDIF + A11Y(KEY)=A11Y(KEY)+RR*QFR1 + 150 CONTINUE + 160 CONTINUE + 170 NUM1=NUM1+LL + NUM2=NUM2+6 + 180 CONTINUE + RETURN + END +* + SUBROUTINE TRIPXZ (IR,MAXKN,NEL,LL4,VOL,MAT,XSGD,XX,YY,ZZ,DD,KN, + 1 QFR,MUZ,IPZ,CYLIND,LC,T,TS,Q,QS,A11Z) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER IR,MAXKN,NEL,LL4,MAT(NEL),KN(MAXKN),MUZ(LL4),IPZ(LL4),LC + REAL VOL(NEL),XSGD(IR,4),XX(NEL),YY(NEL),ZZ(NEL),DD(NEL), + 1 QFR(6*NEL),T(LC),TS(LC),Q(LC,LC),QS(LC,LC),A11Z(*) + LOGICAL CYLIND +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION VAR1,VOL1,VOL2,VOL3,QQX,QQY,QQZ + COMMON /ELEM2/LL,LCC,IJ1(125),IJ2(125),IJ3(125),ISR(6,25), + 1 Q3DP1(125,125),Q3DP2(125,125),Q3DP3(125,125),R3DP(125), + 2 Q3DC1(125,125),Q3DC2(125,125),Q3DC3(125,125),R3DC(125), + 3 R2DP(25),R2DC(25) +*---- +* Z-DIRECTED COUPLINGS. +* +* ASSEMBLY OF MATRIX A11Z. +*---- + CALL TRIPMA(LC,T,TS,Q,QS) + NUM1=0 + NUM2=0 + DO 270 K=1,NEL + L=MAT(K) + IF(L.EQ.0) GO TO 270 + VOL0=VOL(K) + IF(VOL0.EQ.0.0) GO TO 260 + DX=XX(K) + DY=YY(K) + DZ=ZZ(K) + VOL1=VOL0/(DX*DX) + VOL2=VOL0/(DY*DY) + VOL3=VOL0/(DZ*DZ) + DO 230 I=1,LL + IND1=KN(NUM1+I) + IF(IND1.EQ.0) GO TO 230 + INZ1=IPZ(IND1) + KEY0=MUZ(INZ1) + IF(CYLIND) THEN + RR=(R3DP(I)+R3DC(I)*DX/DD(K))*VOL0 + ELSE + RR=R3DP(I)*VOL0 + ENDIF + A11Z(KEY0)=A11Z(KEY0)+RR*XSGD(L,4) + KEY0=KEY0-INZ1 + DO 220 J=1,LL + IND2=KN(NUM1+J) + IF(IND2.EQ.0) GO TO 220 + INZ2=IPZ(IND2) + IF(INZ2.EQ.INZ1) THEN + IF(CYLIND) THEN + QQX=(Q3DP1(I,J)+Q3DC1(I,J)*DX/DD(K))*VOL1 + QQY=(Q3DP2(I,J)+Q3DC2(I,J)*DX/DD(K))*VOL2 + QQZ=(Q3DP3(I,J)+Q3DC3(I,J)*DX/DD(K))*VOL3 + ELSE + QQX=Q3DP1(I,J)*VOL1 + QQY=Q3DP2(I,J)*VOL2 + QQZ=Q3DP3(I,J)*VOL3 + ENDIF + KEY=KEY0+INZ2 + VAR1=QQX*XSGD(L,1)+QQY*XSGD(L,2)+QQZ*XSGD(L,3) + A11Z(KEY)=REAL(A11Z(KEY)+VAR1) + ELSE IF((INZ2.LT.INZ1).AND.(IJ1(I).EQ.IJ1(J)).AND. + 1 (IJ2(I).EQ.IJ2(J))) THEN + IF(CYLIND) THEN + QQZ=(Q3DP3(I,J)+Q3DC3(I,J)*DX/DD(K))*VOL3 + ELSE + QQZ=Q3DP3(I,J)*VOL3 + ENDIF + KEY=KEY0+INZ2 + A11Z(KEY)=REAL(A11Z(KEY)+QQZ*XSGD(L,3)) + ENDIF + 220 CONTINUE + 230 CONTINUE + DO 250 IC=1,6 + QFR1=QFR(NUM2+IC) + IF(QFR1.EQ.0.0) GO TO 250 + DO 240 I1=1,LCC + I=ISR(IC,I1) + IND1=KN(NUM1+I) + IF(IND1.EQ.0) GO TO 240 + INZ1=IPZ(IND1) + KEY=MUZ(INZ1) + IF(CYLIND) THEN + IF(IC.EQ.1) THEN + CRZ=-0.5*R2DP(I1) + ELSE IF(IC.EQ.2) THEN + CRZ=0.5*R2DP(I1) + ELSE + CRZ=R2DC(I1) + ENDIF + RR=(R2DP(I1)+DX*CRZ/DD(K)) + ELSE + RR=R2DP(I1) + ENDIF + A11Z(KEY)=A11Z(KEY)+RR*QFR1 + 240 CONTINUE + 250 CONTINUE + 260 NUM1=NUM1+LL + NUM2=NUM2+6 + 270 CONTINUE + RETURN + END |