Initial commit from Polytechnique Montreal

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