diff options
Diffstat (limited to 'Dragon/src/TLMDIR.f')
| -rw-r--r-- | Dragon/src/TLMDIR.f | 407 |
1 files changed, 407 insertions, 0 deletions
diff --git a/Dragon/src/TLMDIR.f b/Dragon/src/TLMDIR.f new file mode 100644 index 0000000..3903a2c --- /dev/null +++ b/Dragon/src/TLMDIR.f @@ -0,0 +1,407 @@ +*DECK TLMDIR + SUBROUTINE TLMDIR(IPMAT ,IFTRK ,IPRINT,ISPEC,NSKTRK,NBTR , + > NDIM ,NSOUT, NREG ,MXSUB,MXSEG ,NANGL , + > NBDR ,NPLOTS,IPLOT ,IPLP ,DANGLT,DVNOR , + > MATALB,LMIX ) +* +*----------------------------------------------------------------------- +* +*Purpose: +* To generate the Matlab instruction for drawing the +* lines for the directions selected. +* +*Copyright: +* Copyright (C) 2006 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): +* C. Plamondon, G. Marleau +* +*Parameters: input +* IPMAT pointer to Matlab-m file. +* IFTRK pointer to the TRACKING file. +* IPRINT print level. +* ISPEC type of boundary (=0 TISO; =1 TSPC). +* NSKTRK number of records to skip on track file before tracking +* lines can be extracted. +* NBTR number of tracks. +* NDIM number of dimensions for problem. +* NSOUT number of surfaces for problem. +* NREG number of regions for problem. +* MXSUB maximum number of subtracks in a line. +* MXSEG maximum number of segments in a line. +* NANGL number of direction for tracking. +* NBDR number of direction for volume normalization. +* NPLOTS number of plots. +* IPLOT plot number being processed. +* IPLP integer plot parameters. +* DANGLT track directions. +* DVNOR track normalization factor for regional volumes. +* MATALB surface direction and region material identification array. +* LMIX flag set to .true. to draw mixture lines. +* +*---------- +* + IMPLICIT NONE +*---- +* Subroutine arguments +*---- + INTEGER IPMAT,IFTRK + INTEGER IPRINT,ISPEC,NSKTRK,NBTR,NDIM,NSOUT,NREG,MXSUB, + > MXSEG,NANGL,NBDR,NPLOTS,IPLOT + INTEGER IPLP(6,NPLOTS),MATALB(-NSOUT:NREG) + DOUBLE PRECISION DANGLT(NDIM,NANGL),DVNOR(NREG,NBDR) + LOGICAL LMIX +*---- +* Local parameters +*---- + INTEGER IOUT + LOGICAL LNEW + CHARACTER NAMSBR*6 + PARAMETER (IOUT=6,NAMSBR='TLMDIR') +*---- +* Local variables for tracking file +*---- + INTEGER ILINE,IDUM,NSUB,NBSEG,NTLINE,ISEG,IPLANE,IPTA2, + > IPTA3 + DOUBLE PRECISION WEIGHT +*---- +* Other local variables +*---- + INTEGER IREG,ILREG,ILSUR,IDIR,ISV,IDL,IPL,IU,IV,IPM, + > ITRACE,II,ISUB + DOUBLE PRECISION DXYZ(3,2),FLEN + CHARACTER TITLE*36 +*---- +* Allocatable arrays +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: NUMERO,KANGL + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: LENGTH + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: TORIG +*---- +* Scratch storage allocation +* NUMERO region/surface identification number for segment. +* LENGTH segment length. +*---- + ALLOCATE(NUMERO(MXSEG),LENGTH(MXSEG),KANGL(MXSUB), + > TORIG(NDIM,MXSUB)) +*---- +* Processing starts: +* print routine opening header if required +* and initialize various parameters. +*---- + IF(IPRINT .GE. 1) WRITE(IOUT,6000) NAMSBR +*---- +* Print IPMAT header +*---- + IDL=IPLP(2,IPLOT) + IPL=IPLP(3,IPLOT) + IU=IPLP(4,IPLOT) + IV=IPLP(5,IPLOT) + WRITE(TITLE,'(A4,I4,A7,I1,A4,I6,A4,I6)') + >'Dir=',IDL,' Plane=',IPL,' IU=',IU,' IV=',IV + IF(LMIX) THEN + WRITE(IPMAT,7000) NAMSBR,TITLE,MAXVAL(MATALB(1:NREG))+1 + ELSE + WRITE(IPMAT,7000) NAMSBR,TITLE,NREG + ENDIF +*---- +* Print matlab instructions for line segment +*---- + DO IREG=1,NREG + IF(IPRINT .GE. 10) THEN + WRITE(IOUT,6002) IREG + ENDIF + DO ILINE=1,NSKTRK + READ(IFTRK) IDUM + ENDDO + ILREG=0 +*---- +* Scan over lines +*---- + DO ILINE=1,NBTR + READ(IFTRK) NSUB,NBSEG,WEIGHT, + > (KANGL(II),II=1,NSUB), + > (NUMERO(ISEG),ISEG=1,NBSEG), + > (LENGTH(ISEG),ISEG=1,NBSEG), + > NTLINE,IPLANE,IPTA2,IPTA3, + > ((TORIG(IDIR,ISUB),IDIR=1,NDIM),ISUB=1,NSUB) + IF(ISPEC.EQ.0) THEN + IF(NSUB.NE.1) CALL XABORT(NAMSBR//': NSUB.NE.1') + DO IDIR=1,NDIM + DXYZ(IDIR,1)=TORIG(IDIR,1) + ENDDO + ITRACE=0 + IF(IDL .EQ. 0 .OR. IDL .EQ. KANGL(1) ) THEN + IF(NDIM .EQ. 2) THEN + ITRACE=1 + ELSE + IF(IPL .EQ. 0 .OR. IPL .EQ. IPLANE) THEN + IF(IU .GT. 0) THEN + IF(IU .EQ. IPTA2) ITRACE=1 + ELSE IF(IV .GT. 0) THEN + IF(IV .EQ. IPTA3) ITRACE=1 + ELSE + ITRACE=1 + ENDIF + ENDIF + ENDIF + ENDIF + DO ISEG=1,NBSEG + ISV=NUMERO(ISEG) + IF(ISV .GT. 0) THEN + IF(NBDR.EQ.1) THEN + FLEN=LENGTH(ISEG)/DVNOR(ISV,1) + ELSE + FLEN=LENGTH(ISEG)/DVNOR(ISV,KANGL(1)+1) + ENDIF + DO IDIR=1,NDIM + DXYZ(IDIR,2)=DXYZ(IDIR,1)+DANGLT(IDIR,KANGL(1))*FLEN + ENDDO + IF(IREG .EQ. ISV .AND. ITRACE .EQ. 1) THEN + ILREG=ILREG+1 + IF(ILREG .EQ. 1) WRITE(IPMAT,7002) + WRITE(IPMAT,7004) + > ((DXYZ(IDIR,IPM),IPM=1,2),IDIR=1,NDIM), + > FLOAT(ILINE),FLOAT(ISEG),LENGTH(ISEG) + ENDIF + DO IDIR=1,NDIM + DXYZ(IDIR,1)=DXYZ(IDIR,2) + ENDDO + ENDIF + ENDDO + ELSE IF(ISPEC.EQ.1) THEN + ISUB=0 + LNEW=.TRUE. + DO ISEG=1,NBSEG + ISV=NUMERO(ISEG) + IF(ISV.GT.0) THEN + IF(LNEW) THEN + ISUB=ISUB+1 + IF(ISUB.GT.NSUB) CALL XABORT(NAMSBR// + > ': NSUB overflow') + DO IDIR=1,NDIM + DXYZ(IDIR,1)=TORIG(IDIR,ISUB) + ENDDO + LNEW=.FALSE. + ENDIF + IF(NBDR.EQ.1) THEN + FLEN=LENGTH(ISEG)/DVNOR(ISV,1) + ELSE + FLEN=LENGTH(ISEG)/DVNOR(ISV,KANGL(ISUB)+1) + ENDIF + DO IDIR=1,NDIM + DXYZ(IDIR,2)=DXYZ(IDIR,1)+DANGLT(IDIR,KANGL(ISUB))*FLEN + ENDDO + IF(IREG .EQ. ISV) THEN + ILREG=ILREG+1 + IF(ILREG .EQ. 1) WRITE(IPMAT,7002) + IF(IDL .EQ. 0 .OR. IDL .EQ. KANGL(ISUB)) THEN + WRITE(IPMAT,7004) + > ((DXYZ(IDIR,IPM),IPM=1,2),IDIR=1,NDIM), + > FLOAT(ILINE),FLOAT(ISEG),LENGTH(ISEG) + ENDIF + ENDIF + DO IDIR=1,NDIM + DXYZ(IDIR,1)=DXYZ(IDIR,2) + ENDDO + ELSE + LNEW=.TRUE. + ENDIF + ENDDO + IF(ISUB.NE.NSUB) CALL XABORT(NAMSBR//': Algorithm failure') + ENDIF + ENDDO +*---- +* Write Matlab commands to trace lines +*---- + IF(ILREG .GE. 1) THEN + WRITE(IPMAT,7003) + IF(LMIX) THEN + IF(NDIM .EQ. 2) THEN + WRITE(IPMAT,7012) MATALB(IREG)+1 + ELSE + WRITE(IPMAT,7013) MATALB(IREG)+1 + ENDIF + ELSE + IF(NDIM .EQ. 2) THEN + WRITE(IPMAT,7012) IREG + ELSE + WRITE(IPMAT,7013) IREG + ENDIF + ENDIF +*---- +* Change colour for next region +*---- + WRITE(IPMAT,7090) + IF(IREG .NE. NREG) THEN + IF(IPLP(1,IPLOT) .GT. 0) WRITE(IPMAT,7091) + ENDIF + ENDIF + REWIND IFTRK + ENDDO + IF(IPLP(6,IPLOT) .EQ. 0) WRITE(IPMAT,7093) + WRITE(IPMAT,7091) +*---- +* Print matlab instructions for surface points +*---- + IF(IPLP(6,IPLOT) .EQ. 1) THEN + DO ILINE=1,NSKTRK + READ(IFTRK) IDUM + ENDDO + ILSUR=0 +*---- +* Scan over lines +*---- + DO ILINE=1,NBTR + READ(IFTRK) NSUB,NBSEG,WEIGHT, + > (KANGL(II),II=1,NSUB), + > (NUMERO(ISEG),ISEG=1,NBSEG), + > (LENGTH(ISEG),ISEG=1,NBSEG), + > NTLINE,IPLANE,IPTA2,IPTA3, + > ((TORIG(IDIR,ISUB),IDIR=1,NDIM),ISUB=1,NSUB) + IF(ISPEC.EQ.0) THEN + IF(NSUB.NE.1) CALL XABORT(NAMSBR//': NSUB.NE.1') + DO IDIR=1,NDIM + DXYZ(IDIR,1)=TORIG(IDIR,1) + ENDDO + ITRACE=0 + IF(IDL .EQ. 0 .OR. IDL .EQ. KANGL(1) ) THEN + IF(NDIM .EQ. 2) THEN + ITRACE=1 + ELSE + IF(IPL .EQ. 0 .OR. IPL .EQ. IPLANE) THEN + IF(IU .GT. 0) THEN + IF(IU .EQ. IPTA2) ITRACE=1 + ELSE IF(IV .GT. 0) THEN + IF(IV .EQ. IPTA3) ITRACE=1 + ELSE + ITRACE=1 + ENDIF + ENDIF + ENDIF + ENDIF + DO ISEG=1,NBSEG + ISV=NUMERO(ISEG) + IF(ISV .LT. 0) THEN + IF(ITRACE .EQ. 1) THEN + ILSUR=ILSUR+1 + IF(ILSUR .EQ. 1) WRITE(IPMAT,7005) + WRITE(IPMAT,7004) + > (DXYZ(IDIR,1),IDIR=1,NDIM), + > FLOAT(ILINE),FLOAT(ISEG),LENGTH(ISEG) + ENDIF + ELSE IF(ISV .GT. 0) THEN + IF(NBDR.EQ.1) THEN + FLEN=LENGTH(ISEG)/DVNOR(ISV,1) + ELSE + FLEN=LENGTH(ISEG)/DVNOR(ISV,KANGL(1)+1) + ENDIF + DO IDIR=1,NDIM + DXYZ(IDIR,1)=DXYZ(IDIR,1)+DANGLT(IDIR,KANGL(1))*FLEN + ENDDO + ENDIF + ENDDO + ELSE IF(ISPEC.EQ.1) THEN + ISUB=0 + LNEW=.TRUE. + DO ISEG=1,NBSEG + ISV=NUMERO(ISEG) + IF(ISV.GT.0) THEN + IF(LNEW) THEN + ISUB=ISUB+1 + IF(ISUB.GT.NSUB) CALL XABORT(NAMSBR// + > ': NSUB overflow') + DO IDIR=1,NDIM + DXYZ(IDIR,1)=TORIG(IDIR,ISUB) + ENDDO + LNEW=.FALSE. + ENDIF + IF(NBDR.EQ.1) THEN + FLEN=LENGTH(ISEG)/DVNOR(ISV,1) + ELSE + FLEN=LENGTH(ISEG)/DVNOR(ISV,KANGL(ISUB)+1) + ENDIF + DO IDIR=1,NDIM + DXYZ(IDIR,1)=DXYZ(IDIR,1)+DANGLT(IDIR,KANGL(ISUB))*FLEN + ENDDO + ELSE + LNEW=.TRUE. + ENDIF + ENDDO + IF(ISUB.NE.NSUB) CALL XABORT(NAMSBR//': Algorithm failure') + ENDIF + ENDDO +*---- +* Write Matlab commands to trace lines +*---- + WRITE(IPMAT,7003) + IF(NDIM .EQ. 2) THEN + WRITE(IPMAT,7010) + ELSE + WRITE(IPMAT,7011) + ENDIF + WRITE(IPMAT,7092) + WRITE(IPMAT,7093) + WRITE(IPMAT,7091) + REWIND IFTRK + ENDIF +*---- +* Processing finished, return +*---- + IF(IPRINT .GE. 1) THEN + WRITE(IOUT,6001) NAMSBR + ENDIF +*---- +* Scratch storage deallocation +*---- + DEALLOCATE(TORIG,KANGL,LENGTH,NUMERO) + RETURN +*---- +* Output formats +*---- + 6000 FORMAT('(* Output from --',A6,'-- follows ') + 6001 FORMAT(' Output from --',A6,'-- completed *)') + 6002 FORMAT(' Processing lines for region = ',I8) +*---- +* Matlab .m file format +*---- + 7000 FORMAT('%'/'% Output from ',A6/'%'/ + >'hold on;'/7Htitle(',A36,3H');/9Hxcol=jet(,i5,2H);) + 7002 FORMAT('TLMIntegrationLines=[') + 7003 FORMAT(12X,'];') + 7004 FORMAT(9F18.10) + 7005 FORMAT('TLMSurfacePoints=[') + 7010 FORMAT('plot(TLMSurfacePoints(:,1),', + > 'TLMSurfacePoints(:,2),',6H'k.');) + 7011 FORMAT('plot3(TLMSurfacePoints(:,1),', + > 'TLMSurfacePoints(:,2),', + > 'TLMSurfacePoints(:,3),',6H'k.');) + 7012 FORMAT('[m,n]=size(TLMIntegrationLines);'/ + > 'idreg=',I5,';'/ + > 'for i=1:m'/ + > ' TLMcolorset=line([TLMIntegrationLines(i,1),', + > 'TLMIntegrationLines(i,2)],', + > '[TLMIntegrationLines(i,3),', + > 'TLMIntegrationLines(i,4)]);'/ + > ' set(TLMcolorset,',8H'Color',,'xcol(idreg,:));'/ + > 'end;') + 7013 FORMAT('[m,n]=size(TLMIntegrationLines);'/ + > 'idreg=',I5,';'/ + > 'for i=1:m'/ + > ' TLMcolorset=line([TLMIntegrationLines(i,1),', + > 'TLMIntegrationLines(i,2)],', + > '[TLMIntegrationLines(i,3),', + > 'TLMIntegrationLines(i,4)],', + > '[TLMIntegrationLines(i,5),', + > 'TLMIntegrationLines(i,6)]);'/ + > ' set(TLMcolorset,',8H'Color',,'xcol(idreg,:));'/ + > 'end;') + 7090 FORMAT('clear TLMIntegrationLines TLMcolorset ;') + 7091 FORMAT('pause ;') + 7092 FORMAT('clear TLMSurfacePoints ;') + 7093 FORMAT('hold off;') + END |