diff options
Diffstat (limited to 'Trivac/src/TRIDIG.f')
| -rwxr-xr-x | Trivac/src/TRIDIG.f | 171 |
1 files changed, 171 insertions, 0 deletions
diff --git a/Trivac/src/TRIDIG.f b/Trivac/src/TRIDIG.f new file mode 100755 index 0000000..479353e --- /dev/null +++ b/Trivac/src/TRIDIG.f @@ -0,0 +1,171 @@ +*DECK TRIDIG + SUBROUTINE TRIDIG(HNAME,IPTRK,IPSYS,IMPX,MAXMIX,NEL,IPR,MAT,VOL, + 1 SGD) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Driver for the assembly of a cross section diagonal matrix. +* +*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 +* HNAME name of the diagonal matrix. +* IPTRK L_TRACK pointer to the TRIVAC tracking information. +* IPSYS L_SYSTEM pointer to system matrices. +* IMPX print parameter. Equal to zero for no print. +* MAXMIX dimension of matrix SGD. +* NEL total number of finite elements. +* IPR type of assembly: +* =3: the new contribution is added to existing matrix. +* MAT index-number of the mixture type assigned to each volume. +* VOL volumes. +* SGD cross section per material mixture. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + CHARACTER HNAME*(*) + TYPE(C_PTR) IPTRK,IPSYS + INTEGER IMPX,MAXMIX,NEL,IPR,MAT(NEL) + REAL VOL(NEL),SGD(MAXMIX) +*---- +* LOCAL VARIABLES +*---- + PARAMETER (NSTATE=40) + CHARACTER TEXT12*12 + LOGICAL CYLIND,CHEX + INTEGER ISTATE(NSTATE) + INTEGER, DIMENSION(:), ALLOCATABLE :: KN,IPERT,IPW + REAL, DIMENSION(:), ALLOCATABLE :: T,TS,FRZ + REAL, DIMENSION(:,:), ALLOCATABLE :: R,RH,RT + REAL, DIMENSION(:), ALLOCATABLE :: XORZ,DD + REAL, DIMENSION(:), POINTER :: VEC + TYPE(C_PTR) VEC_PTR +*---- +* RECOVER TRIVAC SPECIFIC TRACKING INFORMATION +*---- + CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE) + ITYPE=ISTATE(6) + CYLIND=(ITYPE.EQ.3).OR.(ITYPE.EQ.6) + CHEX=(ITYPE.EQ.8).OR.(ITYPE.EQ.9) + IELEM=ABS(ISTATE(9)) + LL4=ISTATE(11) + ICHX=ISTATE(12) + ISPLH=ISTATE(13) + LX=ISTATE(14) + LY=ISTATE(15) + LZ=ISTATE(16) + LL4F=ISTATE(25) + CALL LCMLEN(IPTRK,'KN',MAXKN,ITYLCM) + ALLOCATE(KN(MAXKN),XORZ(LX*LY*LZ)) + CALL LCMGET(IPTRK,'KN',KN) + IF(CHEX) THEN + CALL LCMGET(IPTRK,'ZZ',XORZ) + CALL LCMGET(IPTRK,'SIDE',SIDE) + ELSE + CALL LCMGET(IPTRK,'XX',XORZ) + ALLOCATE(DD(LX*LY*LZ)) + CALL LCMGET(IPTRK,'DD',DD) + ENDIF + TEXT12=HNAME + IF(IMPX.GT.0) WRITE(6,'(/37H TRIDIG: ASSEMBLY OF DIAGONAL MATRIX , + 1 1H'',A12,2H''.)') TEXT12 +*---- +* INITIALIZATION OF A DIAGONAL SYSTEM MATRIX +*---- + IF(ICHX.EQ.2) THEN + IF(IPR.EQ.3) THEN + CALL LCMGPD(IPSYS,TEXT12,VEC_PTR) + CALL C_F_POINTER(VEC_PTR,VEC,(/ LL4F /)) + ELSE + VEC_PTR=LCMARA(LL4F) + CALL C_F_POINTER(VEC_PTR,VEC,(/ LL4F /)) + VEC(:LL4F)=0.0 + ENDIF + ELSE + IF(IPR.EQ.3) THEN + CALL LCMGPD(IPSYS,TEXT12,VEC_PTR) + CALL C_F_POINTER(VEC_PTR,VEC,(/ LL4 /)) + ELSE + VEC_PTR=LCMARA(LL4) + CALL C_F_POINTER(VEC_PTR,VEC,(/ LL4 /)) + VEC(:LL4)=0.0 + ENDIF + ENDIF +*---- +* COMPUTE THE DIAGONAL SYSTEM MATRIX +*---- + IF((ICHX.EQ.1).AND.(.NOT.CHEX)) THEN +* VARIATIONAL COLLOCATION METHOD. + CALL LCMSIX(IPTRK,'BIVCOL',1) + CALL LCMLEN(IPTRK,'T',LC,ITYLCM) + ALLOCATE(T(LC),TS(LC)) + CALL LCMGET(IPTRK,'T',T) + CALL LCMGET(IPTRK,'TS',TS) + CALL LCMSIX(IPTRK,' ',2) + CALL TRIASP(IELEM,MAXMIX,NEL,LL4,CYLIND,SGD,XORZ,DD,VOL,MAT, + 1 KN,LC,T,TS,VEC) + DEALLOCATE(T,TS) + ELSE IF((ICHX.EQ.1).AND.CHEX) THEN +* MESH CORNER FINITE DIFFERENCES IN HEXAGONAL GEOMETRY. + CALL LCMSIX(IPTRK,'BIVCOL',1) + ALLOCATE(R(2,2),RH(6,6),RT(3,3)) + CALL LCMGET(IPTRK,'R',R) + CALL LCMGET(IPTRK,'RH',RH) + CALL LCMGET(IPTRK,'RT',RT) + CALL LCMSIX(IPTRK,' ',2) + CALL TRIAHP(MAXKN,ISPLH,MAXMIX,NEL,LL4,SGD,SIDE,XORZ,VOL,MAT, + 1 KN,R,RH,RT,VEC) + DEALLOCATE(RT,RH,R) + ELSE IF((ICHX.EQ.2).AND.CHEX) THEN +* DUAL (THOMAS-RAVIART-SCHNEIDER) FINITE ELEMENT METHOD. + NBLOS=LX*LZ/3 + ALLOCATE(IPERT(NBLOS),FRZ(NBLOS)) + CALL LCMGET(IPTRK,'IPERT',IPERT) + CALL LCMGET(IPTRK,'FRZ',FRZ) + CALL TRIASH(IELEM,MAXMIX,LL4,NBLOS,MAT,SIDE,XORZ,FRZ,SGD,KN, + 1 IPERT,VEC) + DEALLOCATE(FRZ,IPERT) + ELSE IF(.NOT.CHEX) THEN +* DUAL FINITE ELEMENT METHOD. + IDIM=1 + IF((ITYPE.EQ.5).OR.(ITYPE.EQ.6).OR.(ITYPE.EQ.8)) IDIM=2 + IF((ITYPE.EQ.7).OR.(ITYPE.EQ.9)) IDIM=3 + CALL TRIASD(MAXKN,IELEM,ICHX,IDIM,MAXMIX,NEL,LL4,SGD,VOL,MAT, + 1 KN,VEC) + ELSE IF(CHEX.AND.(ISPLH.EQ.1)) THEN +* MESH CENTERED FINITE DIFFERENCES IN HEXAGONAL GEOMETRY. + CALL TRIAHD(MAXMIX,NEL,LL4,SGD,VOL,MAT,VEC) + ELSE IF(CHEX.AND.(ISPLH.GT.1)) THEN +* MESH CENTERED FINITE DIFFERENCES IN TRIANGULAR GEOMETRY. + ALLOCATE(IPW(LL4)) + CALL LCMGET(IPTRK,'IPW',IPW) + CALL TRIMTD(ISPLH,MAXMIX,NEL,LL4,VOL,MAT,SGD,KN,IPW,VEC) + DEALLOCATE(IPW) + ENDIF +*---- +* STORAGE OF THE DIAGONAL SYSTEM MATRIX +*---- + IF(ICHX.EQ.2) THEN + CALL LCMPPD(IPSYS,TEXT12,LL4F,2,VEC_PTR) + ELSE + CALL LCMPPD(IPSYS,TEXT12,LL4,2,VEC_PTR) + ENDIF +*---- +* RELEASE TRIVAC SPECIFIC TRACKING INFORMATION +*---- + IF(.NOT.CHEX) DEALLOCATE(DD) + DEALLOCATE(XORZ,KN) + RETURN + END |