diff options
Diffstat (limited to 'Dragon/src/MCTPTR.f')
| -rw-r--r-- | Dragon/src/MCTPTR.f | 178 |
1 files changed, 178 insertions, 0 deletions
diff --git a/Dragon/src/MCTPTR.f b/Dragon/src/MCTPTR.f new file mode 100644 index 0000000..daa2055 --- /dev/null +++ b/Dragon/src/MCTPTR.f @@ -0,0 +1,178 @@ +*DECK MCTPTR + SUBROUTINE MCTPTR(IPTRK,IPRINT,NDIM,MAXMSH,ITYPBC,NUCELL,MXGSUR, + 1 MXGREG,MAXPIN,IUNFLD,DGMESH,XYZL,NBIND,POS, + 2 LENGTH,VDIR,ODIR,IDS,IDSO,IREG,INDX,IDIRC, + 3 MESHC,NSURC,NREGC,NTPIN,CELLPO,PINCEN,INDEX, + 4 IDREG,DCMESH,ITPIN,DRAPIN) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Find region/surface index from an initial position and a path to +* travel. +* +*Copyright: +* Copyright (C) 2008 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): R. Le Tellier +* +*Parameters: input +* IPTRK pointer to the TRACKING data structure. +* IPRINT print level. +* NDIM problem dimensions. +* MAXMSH maximum number of elements in MESH array. +* ITYPBC type of boundary. +* NUCELL number of cell after unfolding in +* $X$, $Y$ and $Z$ directions. +* MXGSUR maximum number of surfaces for any geometry. +* MXGREG maximum number of region for any geometry. +* MAXPIN maximum number of pins in a cell. +* IUNFLD description of unfolded geometry. +* DGMESH meshing vector for global geometry. +* XYZL undefined. +* NBIND first dimension of INDX. +* VDIR travel direction (unit vector). +* +*Parameters: output +* IDS outer surface orientation index. ABS(IDSO). +* IDSO outer surface oreientation index with +-1:X+-; +-2:Y+-; +* +-3:Z+- faces. +* IREG region/surface index. +* ODIR search (octant) direction. +* +*Parameters: input/output +* POS initial/final position. +* LENGTH length to travel/remaining length on the path. +* INDX location index of the initial/final position in the +* geometry structure. +* +*Parameters: scratch +* IDIRC undefined. +* MESHC undefined. +* NSURC undefined. +* NREGC undefined. +* NTPIN undefined. +* CELLPO undefined. +* PINCEN undefined. +* INDEX undefined. +* IDREG undefined. +* DCMESH undefined. +* ITPIN undefined. +* DRAPIN undefined. +* +*----------------------------------------------------------------------- +* + USE GANLIB + IMPLICIT NONE +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPTRK + INTEGER IPRINT,NDIM,MAXMSH,ITYPBC,NUCELL(3),MXGSUR,MXGREG,MAXPIN, + 1 IUNFLD(2,NUCELL(1),NUCELL(2),NUCELL(3)),NBIND,ODIR(3),IREG, + 2 INDX(NBIND,0:2),IDIRC(2),MESHC(4,2),NSURC(2),NREGC(2),NTPIN, + 3 INDEX(5,-MXGSUR:MXGREG,2),IDREG(MXGREG,2),ITPIN(3,MAXPIN) + DOUBLE PRECISION DGMESH(-1:MAXMSH,4),XYZL(2,NDIM),POS(3),LENGTH, + 1 VDIR(3),CELLPO(3,2),PINCEN(3),DCMESH(-1:MAXMSH,4,2), + 2 DRAPIN(-1:4,MAXPIN) +*---- +* LOCAL VARIABLES +*---- + INTEGER IDIR,IBCO(3),IDS,IDSO,JJ + DOUBLE PRECISION VCOR(3),LEN,LENM + DOUBLE PRECISION EPS + PARAMETER(EPS=1.D-8) + INTEGER INDOS(2,3) + DATA INDOS / 2,3, + 1 3,1, + 2 1,2 / +*---- +* VERIFY IF A BOUNDARY IS REACHED +*---- + IDSO=0 + LENM=0.0D0 + IF(ITYPBC.EQ.0) THEN +* CARTESIAN BOUNDARY +* find corresponding geometry corner corresponding to travel direction +* check if a boundary is reached and if so, which one it is + IDS=0 + LENM=LENGTH + DO IDIR=1,NDIM + IF(VDIR(IDIR).GT.EPS) THEN + ODIR(IDIR)=1 + IBCO(IDIR)=2 + VCOR(IDIR)=XYZL(2,IDIR)-POS(IDIR) + LEN=VCOR(IDIR)/VDIR(IDIR) + IF(LEN.LT.LENM) THEN + LENM=LEN + IDS=IDIR + IF(IDIR.EQ.1) THEN + IREG=-2 ; + ELSEIF(IDIR.EQ.2) THEN + IREG=-4 ; + ELSEIF(IDIR.EQ.3) THEN + IREG=-6 ; + ENDIF + ENDIF + ELSEIF(VDIR(IDIR).LT.-EPS) THEN + ODIR(IDIR)=-1 + IBCO(IDIR)=1 + VCOR(IDIR)=XYZL(1,IDIR)-POS(IDIR) + LEN=VCOR(IDIR)/VDIR(IDIR) + IF(LEN.LT.LENM) THEN + LENM=LEN + IDS=IDIR + IF(IDIR.EQ.1) THEN + IREG=-1 ; + ELSEIF(IDIR.EQ.2) THEN + IREG=-3 ; + ELSEIF(IDIR.EQ.3) THEN + IREG=-5 ; + ENDIF + ENDIF + ELSE + ODIR(IDIR)=0 + ENDIF + ENDDO + ELSE +* CYLINDRICAL BOUNDARY + CALL XABORT('MCTPTR: CYLINDRICAL/HEXAGONAL BOUNDARY NOT IMPLEM' + 1 //'ENTED YET.') + ENDIF + IF(IDS.EQ.0) THEN +*---- +* NO: LOCATE POINT WITHIN THE GEOMETRY +*---- + DO IDIR=1,NDIM + POS(IDIR)=POS(IDIR)+LENGTH*VDIR(IDIR) + ENDDO + DO IDIR=NDIM+1,3 + POS(IDIR)=1.D0 + ENDDO +!!!!!! for the time being it is the same routine as for a starting point +!!!!!! an optimized version should use the info of the previous point +!!!!!! to start its search + CALL MCTCTR(IPTRK,IPRINT,NDIM,MAXMSH,NUCELL,MXGSUR,MXGREG, + 1 MAXPIN,IUNFLD,DGMESH,NBIND,ODIR,POS,IREG,INDX,IDIRC,MESHC, + 2 NSURC,NREGC,NTPIN,CELLPO,PINCEN,INDEX,IDREG,DCMESH,ITPIN, + 3 DRAPIN) + ELSE +*---- +* YES: LOCATE POINT ON THE BOUNDARY +*---- + POS(IDS)=XYZL(IBCO(IDS),IDS) + DO JJ=1,2 + IDIR=INDOS(JJ,IDS) + POS(IDIR)=POS(IDIR)+LENM*VDIR(IDIR) + ENDDO + LENGTH=LENGTH-LENM + IDSO=IDS + IF(IBCO(IDS).EQ.1) IDSO=-IDSO + ENDIF +* + RETURN + END |