Initial commit from Polytechnique Montreal

1 files changed, 310 insertions, 0 deletions

diff --git a/Dragon/src/NXTVHC.f b/Dragon/src/NXTVHC.f
new file mode 100644
index 0000000..60e18c7
--- /dev/null
+++ b/Dragon/src/NXTVHC.f
@@ -0,0 +1,310 @@
+*DECK NXTVHC
+      SUBROUTINE NXTVHC(IPRINT,NDIM  ,MXMESH,MAXSUR,MAXREG,
+     >                  MESH  ,DMESH ,NBSUR ,NBREG ,INDXSR,SURVOL,
+     >                  POSTRI)
+*
+*----------
+*
+*Purpose:
+* Compute the volume and area associated with each region
+* or surface for a annular/hexagon with triangular mesh 
+* in 2-D or 3-D geometry.
+*
+*Copyright:
+* Copyright (C) 2005 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):
+* G. Marleau.
+*
+*Parameters: input
+* IPRINT  print level.
+* NDIM    dimension of problem.
+* MXMESH  maximum number of spatial subdivision in
+*         $X$, $Y$ and $Z$.
+* MAXSUR  maximum number of surfaces in the geometry.
+* MAXREG  maximum number of regions in the geometry.
+* MESH    effective number of spatial subdivision in
+*         each direction ($X$, $Y$ and $Z$).
+* DMESH   spatial description of the parallepiped.
+* POSTRI  triangle position:
+*         POSTRI(1,*,*,*) is X position;
+*         POSTRI(2,*,*,*) is Y position;
+*         POSTRI(*,1,*,*) is location of first corner;
+*         POSTRI(*,2,*,*) is location of second corner;
+*         POSTRI(*,3,*,*) is location of third corner;
+*         POSTRI(*,*,i,j) is location of triangle i in cector j.
+*
+*Parameters: output
+* NBSUR   number of surfaces in the geometry.
+* NBREG   number of regions in the geometry.
+* INDXSR  local indexing of surfaces/regions.
+* SURVOL  area/volume of regions.
+*
+*Comments:
+*  1- Contents of the DMESH array:
+*     hexagonal mesh is DMESH(i,1) for i=0,MESH(1);
+*     mesh in $Z$  is z(k)=DMESH(k,3) for k=0,MESH(3);
+*     annular regions in the $X-Y$ plane
+*     centre of cylinder in (x,y)=(DMESH(-1,1),DMESH(-1,2))
+*     radius of shells      r(l)=DMESH(l,4), l=1,MESH(4)
+*  2- Contents of the INDXSR array:
+*     For i>0
+*       INDXSR(1,i)= iu is the $U$ location of region i
+*       INDXSR(2,i)= iv is the $V$ location of region i
+*       INDXSR(3,i)= iz is the $Z$ location of region i
+*       INDXSR(4,i)= ir is the $R$ location of region i
+*       INDXSR(5,i)= iw is the $W$ location of region i
+*     For i<0
+*       INDXSR(1,i)= iu is the $U$ location of surface i
+*       INDXSR(2,i)= iv is the $V$ location of surface i
+*       INDXSR(3,i)= iz is the $Z$ location of surface i
+*       INDXSR(4,i)= ir is the $R$ location of surface i
+*       INDXSR(5,i)= iw is the $W$ location of surface i
+*       with INDXSR(n,i)=-1 for surface associated with
+*                           location 0 in direction n.
+*       with INDXSR(n,i)=-2 for surface associated with
+*                           location MESH(n) in direction n.
+*       Note that for radial regions INDXSR(n,i)=-1 does not
+*       exists.
+*
+*----------
+*
+      IMPLICIT         NONE
+*----
+*  Subroutine arguments
+*----
+      INTEGER          IPRINT,NDIM,MXMESH,MAXSUR,MAXREG
+      INTEGER          MESH(4)
+      DOUBLE PRECISION DMESH(-1:MXMESH,4)
+      INTEGER          NBSUR,NBREG,INDXSR(5,-MAXSUR:MAXREG)
+      DOUBLE PRECISION SURVOL(-MAXSUR:MAXREG)
+      DOUBLE PRECISION POSTRI(2,3,MXMESH*MXMESH,6)
+*----
+*  Local parameters
+*----
+      INTEGER          IOUT
+      CHARACTER        NAMSBR*6
+      PARAMETER       (IOUT=6,NAMSBR='NXTVHC')
+      INTEGER          MAXDIM
+      PARAMETER       (MAXDIM=4)
+      DOUBLE PRECISION DZERO
+      PARAMETER       (DZERO=0.0D0)
+*----
+*  Functions
+*----
+      INTEGER          NXTITA,ITYITA
+      DOUBLE PRECISION VOLINT
+*----
+*  Local variables
+*----
+      INTEGER          NX,NZ,NR,NRTP,NRP,NSTP,NSP,NRTPP,NRP1,
+     >                 IT,IZ,IR,ISECT,IOFS,IVSI,NBVHEX,NBSHEX
+      DOUBLE PRECISION POSANN(0:2),VOLZ
+      INTEGER          ILVI,ILVT,IDIR,IOFZ,IS1,ISURTF,ISURTI,
+     >                 IVOLF,IVOLI,NBREGT,
+     >                 NSKPF,NSKPI,NBSURT
+*----
+*  Prepare loop over Z-direction.
+*----
+      NX=MESH(1)
+      NZ=MESH(3)
+      NR=MESH(4)
+      NRTP=NX**2
+      NRP=6*NRTP
+      NSTP=2*NX-1
+      NSP=6*NSTP
+      NRP1=NR+1
+      NRTPP=NRP*NRP1
+*----
+*  Compute number of surfaces
+*  1- Surface parallel to cylinder axis
+*----
+      NBVHEX=NRP
+      NBREG=NRTPP
+      NBSUR=NSP
+      NBSHEX=NSP
+      IF(NDIM .EQ. 3) THEN
+        NBREG=NBREG*NZ
+        NBSUR=NBSHEX*NZ+2*NRTPP
+        NBVHEX=NBVHEX*NZ
+        NBSHEX=NBSHEX*NZ+2*NRP
+      ENDIF
+      NBREGT=NBVHEX
+      NBSURT=NBSHEX
+*----
+*  2- Surface normal to cylinder axis (if any)
+*----
+      NSKPI=0
+      NSKPF=0
+      IF(NDIM .NE. 3) THEN
+        NSKPI=NRP
+        NSKPF=(NR+1)*NSKPI
+      ENDIF
+      IF(IPRINT .GE. 100) THEN
+        WRITE(IOUT,6000) NAMSBR
+        WRITE(IOUT,6010) 'H',NX
+        WRITE(IOUT,6011) 'MESHH ='
+        WRITE(IOUT,6012) (DMESH(IT,1),IT=-1,2*NX)
+        IF(NDIM .EQ. 3) THEN
+          WRITE(IOUT,6010) 'Z',NZ
+          WRITE(IOUT,6011) 'MESHZ ='
+          WRITE(IOUT,6012) (DMESH(IZ,3),IZ=-1,NZ)
+        ENDIF
+          WRITE(IOUT,6010) 'R',NR
+        WRITE(IOUT,6011) 'MESHR ='
+        WRITE(IOUT,6012) (DMESH(IR,4),IR=-1,NR)
+      ENDIF
+      SURVOL(-MAXSUR:MAXREG)=DZERO
+*----
+*  Call NXTVHT to obtain the triangles volumes and external surfaces
+*----
+      CALL NXTVHT(IPRINT,NDIM  ,MXMESH,NBSHEX,NBVHEX,
+     >            MESH  ,DMESH ,NBSURT,NBREGT,
+     >            INDXSR(1,-NBSHEX),SURVOL(-NBSHEX))
+*----
+*  For 3-D case, displace hexagonal surfaces towards
+*  the end of SURVOL leaving space for the possible
+*  annular sub-surfaces (skip NR spaces for radial volumes).
+*----
+      ISURTF=NBSUR
+      ISURTI=NBSURT
+      IF(NDIM .EQ. 3) THEN
+        DO IZ=1,2
+          DO IT=1,NRP
+            SURVOL(-ISURTF)=SURVOL(-ISURTI)
+            SURVOL(-ISURTI)=DZERO
+            DO IDIR=1,5
+              INDXSR(IDIR,-ISURTF)=INDXSR(IDIR,-ISURTI)
+              INDXSR(IDIR,-ISURTI)=0
+            ENDDO
+            ISURTF=ISURTF-NRP1
+            ISURTI=ISURTI-1
+          ENDDO
+        ENDDO
+      ENDIF
+*----
+*  Displace triangular volumes towards the end of vector
+*  SURVOL leaving space for the possible annular sub-regions.
+*----
+      IVOLF=NBREG
+      DO IVOLI=NBREGT,1,-1
+        SURVOL(IVOLF)=SURVOL(IVOLI)
+        SURVOL(IVOLI)=DZERO
+        DO IDIR=1,5
+          INDXSR(IDIR,IVOLF)=INDXSR(IDIR,IVOLI)
+          INDXSR(IDIR,IVOLI)=0
+        ENDDO
+        IVOLF=IVOLF-NRP1
+      ENDDO
+      IF(IPRINT .GE. 100) THEN
+        WRITE(IOUT,6002) 'SurVolTri'
+        WRITE(IOUT,6005) (ILVT,(INDXSR(IDIR,ILVT),IDIR=1,5),
+     >                    SURVOL(ILVT),ILVT=-NBSUR,NBREG)
+        WRITE(IOUT,6003)
+      ENDIF
+*----
+*  Loop over radial regions for triangular region
+*  intersection.
+*----
+      POSANN(1)=0.0
+      POSANN(2)=0.0
+      DO IR=NR,1,-1
+        POSANN(0)=DMESH(IR,4)
+*----
+*  Loop over a single sector of triangle since centered annular
+*  region
+*----
+        ISECT=1
+        DO IT=1,NRTP
+          ITYITA=NXTITA(POSTRI(1,1,IT,ISECT),POSANN,VOLINT)
+          IF(ITYITA .NE. 0) THEN
+*            WRITE(IOUT,*) 'IR,ISECT,IT,ITYITA,VOLINT',
+*     >      IR,ISECT,IT,ITYITA,VOLINT
+            IF(NDIM .EQ. 3) THEN
+*----
+*  Remove contribution from top and bottom surfaces
+*----
+              IOFZ=NSP*NZ
+              DO IZ=1,2
+                IOFS=(IT-1)*NRP1+IR
+                DO IS1=1,6
+                  ISURTF=IOFZ+IOFS
+                  ISURTI=ISURTF+1
+                  SURVOL(-ISURTI)=SURVOL(-ISURTI)-VOLINT
+                  SURVOL(-ISURTF)=VOLINT
+                  INDXSR(1,-ISURTF)=INDXSR(1,-ISURTI)
+                  INDXSR(2,-ISURTF)=INDXSR(2,-ISURTI)
+                  INDXSR(3,-ISURTF)=INDXSR(3,-ISURTI)
+                  INDXSR(4,-ISURTF)=IR
+                  INDXSR(5,-ISURTF)=INDXSR(5,-ISURTI)
+                  IOFS=IOFS+NRP1*NRTP
+                ENDDO
+                IOFZ=IOFZ+NRTPP
+              ENDDO
+*----
+*  Remove contributions to volumes
+*----
+              IOFZ=0
+              DO IZ=1,NZ
+                IOFS=(IT-1)*NRP1+IR
+                VOLZ=VOLINT*(DMESH(IZ,3)-DMESH(IZ-1,3))
+                DO IS1=1,6
+                  ILVT=IOFZ+IOFS
+                  ILVI=ILVT+1
+                  SURVOL(ILVI)=SURVOL(ILVI)-VOLZ
+                  SURVOL(ILVT)=VOLZ
+                  INDXSR(1,ILVT)=INDXSR(1,ILVI)
+                  INDXSR(2,ILVT)=INDXSR(2,ILVI)
+                  INDXSR(3,ILVT)=INDXSR(3,ILVI)
+                  INDXSR(4,ILVT)=IR
+                  INDXSR(5,ILVT)=INDXSR(5,ILVI)
+                  IOFS=IOFS+NRP1*NRTP
+                ENDDO
+                IOFZ=IOFZ+NRTPP
+              ENDDO
+            ELSE
+              IOFS=(IT-1)*NRP1+IR
+              VOLZ=VOLINT
+              DO IS1=1,6
+                ILVT=IOFS
+                ILVI=ILVT+1
+                SURVOL(ILVI)=SURVOL(ILVI)-VOLZ
+                SURVOL(ILVT)=VOLZ
+                INDXSR(1,ILVT)=INDXSR(1,ILVI)
+                INDXSR(2,ILVT)=INDXSR(2,ILVI)
+                INDXSR(3,ILVT)=INDXSR(3,ILVI)
+                INDXSR(4,ILVT)=IR
+                INDXSR(5,ILVT)=INDXSR(5,ILVI)
+                IOFS=IOFS+NRP1*NRTP
+              ENDDO
+            ENDIF
+          ENDIF
+        ENDDO
+      ENDDO
+*----
+*  Print surfaces and volumes if required and return
+*----
+      IF(IPRINT .GE. 100) THEN
+        WRITE(IOUT,6002) 'SurVolTriAnn'
+        WRITE(IOUT,6005) (IVSI,(INDXSR(IR,IVSI),IR=1,5),SURVOL(IVSI),
+     >                     IVSI=-NBSUR,NBREG)
+        WRITE(IOUT,6003)
+        WRITE(IOUT,6001) NAMSBR
+      ENDIF
+      RETURN
+*----
+*  Output formats
+*----
+ 6000 FORMAT('(* Output from --',A6,'-- follows ')
+ 6001 FORMAT('   Output from --',A6,'-- completed *)')
+ 6002 FORMAT(A12,'={')
+ 6003 FORMAT('};')
+ 6005 FORMAT((6(I10,','),F20.10,:,','))
+ 6010 FORMAT(1X,'MESH DIMENSIONS IN ',A1,' =',I10)
+ 6011 FORMAT(1X,A7)
+ 6012 FORMAT(5F20.10)
+      END