Initial commit from Polytechnique Montreal

1 files changed, 277 insertions, 0 deletions

diff --git a/Dragon/src/NXTVCA.f b/Dragon/src/NXTVCA.f
new file mode 100644
index 0000000..449ee79
--- /dev/null
+++ b/Dragon/src/NXTVCA.f
@@ -0,0 +1,277 @@
+*DECK NXTVCA
+      SUBROUTINE NXTVCA(IPRINT,NDIM  ,IDIRC ,MXMESH,MAXSUR,MAXREG,
+     >                  MESH  ,DMESH ,NBSUR ,NBREG ,INDXSR,SURVOL)
+*
+*----------
+*
+*Purpose:
+* Compute the volume and area associated with each region
+* or surface for a Cartesian 1-D, 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.
+* IDIRC   the direction of the first axis of a Cartesian geometry
+*         assuming the axis are in a cyclic rotation.
+* 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.
+*
+*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.
+*
+*Reference:
+*  G. Marleau,
+*  New Geometries Processing in DRAGON: The NXT: Module,
+*  Report IGE-260, Polytechnique Montreal,
+*  Montreal, 2005.
+*
+*Comments:
+*  1- Contents of IDIRC:
+*         IDIRC      axes in 1-D   axes in 2-D   axes in 3-D
+*         1          x             (x,y)         (x,y,z)
+*         2          y             (y,z)         (y,z,x)
+*         3          z             (z,x)         (z,x,y)
+*  2- Contents of the DMESH array:
+*     mesh in $X$  is x(i)=DMESH(i,1) for i=0,MESH(1);
+*     mesh in $Y$  is y(j)=DMESH(j,2) for j=0,MESH(2);
+*     mesh in $Z$  is z(k)=DMESH(k,3) for k=0,MESH(3);
+*     if(IDIRC = 1) then
+*     ->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)
+*     else if(IDIRC = 2) then
+*     ->annular regions in the $Y-Z$ plane
+*       centre of cylinder in (y,z)=(DMESH(-1,2),DMESH(-1,3))
+*       radius of shells      r(l)=DMESH(l,4), l=1,MESH(4)
+*     else if(IDIRC = 3) then
+*     ->annular regions in the $Z-X$ plane
+*       centre of cylinder in (z,x)=(DMESH(-1,3),DMESH(-1,1))
+*       radius of shells      r(l)=DMESH(l,4), l=1,MESH(4)
+*     endif
+*  3- Contents of the INDXSR array:
+*     For i>0
+*       INDXSR(1,i)= ix is the $X$ location of region i
+*       INDXSR(2,i)= iy is the $Y$ location of region i
+*       INDXSR(3,i)= iz is the $Z$ location of region i
+*       INDXSR(4,i)= ir =0 is the $R$ location of region i.
+*       INDXSR(5,i)= not used.
+*     For i<0
+*       INDXSR(1,i)= ix is the $X$ location of surface i
+*       INDXSR(2,i)= iy is the $Y$ 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)= not used.
+*       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,IDIRC,MXMESH,MAXSUR,MAXREG
+      INTEGER          MESH(4)
+      DOUBLE PRECISION DMESH(-1:MXMESH,4)
+      INTEGER          NBSUR,NBREG,INDXSR(5,-MAXSUR:MAXREG)
+      DOUBLE PRECISION SURVOL(-MAXSUR:MAXREG)
+*----
+*  Local parameters
+*----
+      INTEGER          IOUT
+      CHARACTER        NAMSBR*6
+      PARAMETER       (IOUT=6,NAMSBR='NXTVCA')
+      INTEGER          MAXDIM
+      PARAMETER       (MAXDIM=4)
+      DOUBLE PRECISION DZERO,DONE
+      PARAMETER       (DZERO=0.0D0,DONE=1.0D0)
+*----
+*  Local variables
+*----
+      INTEGER          IDIR(MAXDIM),NM(MAXDIM),IDM(MAXDIM)
+      INTEGER          ID,IDG,IDGP1,IDGP2,IDGPP,IX,IY,IZ,IXYZ,IR,
+     >                 ISUR,ISUR2,IVOL,IDIRCX
+      DOUBLE PRECISION DX,DY,DZ
+*----
+*  Data
+*----
+      CHARACTER        CDIR(MAXDIM)*1
+      SAVE             CDIR
+      DATA             CDIR /'X','Y','Z','R'/
+*----
+*  Prepare loops over spatial directions as a function
+*  of IDIRC and NDIM.
+*----
+      SURVOL(-MAXSUR:MAXREG)=DZERO
+      IDG=0
+      IDGP1=0
+      IDGP2=0
+      IR=0
+      NBREG=1
+      DO 100 ID=1,NDIM
+        IDG=MOD(IDIRC+ID-2,3)+1
+        IDIR(ID)=IDG
+        NM(IDG)=MESH(IDG)
+        NBREG=NBREG*NM(IDG)
+        IDM(IDG)=1
+ 100  CONTINUE
+      DO 101 ID=NDIM+1,3
+        IDG=MOD(IDIRC+ID-2,3)+1
+        IDIR(ID)=IDG
+        NM(IDG)=1
+        IDM(IDG)=0
+ 101  CONTINUE
+      IF(MAXREG .LT. NBREG) CALL XABORT(NAMSBR//
+     >': Insufficient space to store region volumes')
+*----
+*  number of surfaces
+*----
+      NBSUR=0
+      DO 102 ID=1,NDIM
+        IDG=MOD(IDIRC+ID-2,3)+1
+        NBSUR=NBSUR+2*NBREG/NM(IDG)
+ 102  CONTINUE
+      IF(MAXSUR .LT. NBSUR) CALL XABORT(NAMSBR//
+     >': Insufficient space to store surface areas')
+*----
+*  Print mesh if required
+*----
+      IF(IPRINT .GE. 100) THEN
+        WRITE(IOUT,6000) NAMSBR
+        DO 600 ID=1,NDIM
+          IDG=IDIR(ID)
+          WRITE(IOUT,6002) 'MESH'//CDIR(IDG)
+          WRITE(IOUT,6006) (DMESH(IXYZ,IDG),IXYZ=0,NM(IDG))
+          WRITE(IOUT,6003)
+ 600    CONTINUE
+      ENDIF
+*----
+*  Compute surface area
+*  1- Loop over directions (110)
+*----
+      ISUR=0
+      IDIRCX=1
+      DO 110 ID=1,NDIM
+        IDG=MOD(IDIRCX+ID-2,3)+1
+        IF(IDM(IDG) .EQ. 1) THEN
+          IDGP1=MOD(IDIRCX+ID-1,3)+1
+          IDGP2=MOD(IDIRCX+ID,3)+1
+*----
+*  2- Loop over second normal direction (111)
+*----
+          ISUR2=ISUR+NM(IDGP2)*NM(IDGP1)
+          DO 111 IY=1,NM(IDGP2)
+            DY=DONE
+            IF(IDM(IDGP2) .EQ. 1)
+     >      DY=DMESH(IY,IDGP2)-DMESH(IY-1,IDGP2)
+*----
+*  3- Loop over first normal direction (112)
+*----
+            DO 112 IX=1,NM(IDGP1)
+              DX=DONE
+              IF(IDM(IDGP1) .EQ. 1)
+     >        DX=DMESH(IX,IDGP1)-DMESH(IX-1,IDGP1)
+              ISUR=ISUR+1
+              SURVOL(-ISUR)=DX*DY
+              INDXSR(IDGP1,-ISUR)=IX
+              INDXSR(IDGP2,-ISUR)=IY
+              INDXSR(IDG,-ISUR)=-1
+              INDXSR(4,-ISUR)=IR
+              ISUR2=ISUR2+1
+              SURVOL(-ISUR2)=DX*DY
+              INDXSR(IDGP1,-ISUR2)=IX
+              INDXSR(IDGP2,-ISUR2)=IY
+              INDXSR(IDG,-ISUR2)=-2
+              INDXSR(4,-ISUR2)=IR
+ 112        CONTINUE
+ 111      CONTINUE
+          ISUR=ISUR2
+        ENDIF
+ 110  CONTINUE
+*----
+*  Computes regional volumes
+*  1- Loop on $Z$ (120)
+*----
+      IR=0
+      IVOL=0
+      SURVOL(IVOL)=DZERO
+      INDXSR(IDGP1,IVOL)=0
+      INDXSR(IDGP2,IVOL)=0
+      INDXSR(IDG,IVOL)=0
+      INDXSR(4,IVOL)=0
+      IDGPP=IDIR(3)
+      IDGP2=IDIR(2)
+      IDGP1=IDIR(1)
+      IDGPP=3
+      IDGP2=2
+      IDGP1=1
+      DO 120 IZ=1,NM(IDGPP)
+        DZ=DONE
+        IF(IDM(IDGPP) .EQ. 1)
+     >  DZ=DMESH(IZ,IDGPP)-DMESH(IZ-1,IDGPP)
+*----
+*  2- Loop on $Y$ (121)
+*----
+        DO 121 IY=1,NM(IDGP2)
+          DY=DONE
+          IF(IDM(IDGP2) .EQ. 1)
+     >    DY=DMESH(IY,IDGP2)-DMESH(IY-1,IDGP2)
+*----
+*  3- Loop on $X$ (122)
+*----
+          DO 122 IX=1,NM(IDGP1)
+            DX=DONE
+            IF(IDM(IDGP1) .EQ. 1)
+     >      DX=DMESH(IX,IDGP1)-DMESH(IX-1,IDGP1)
+            IVOL=IVOL+1
+            SURVOL(IVOL)=DX*DY*DZ
+            INDXSR(IDGP1,IVOL)=IX
+            INDXSR(IDGP2,IVOL)=IY
+            INDXSR(IDGPP,IVOL)=IZ
+            INDXSR(4,IVOL)=IR
+ 122      CONTINUE
+ 121    CONTINUE
+ 120  CONTINUE
+*----
+*  Print volumes if required
+*----
+      IF(IPRINT .GE. 100) THEN
+        WRITE(IOUT,6002) 'SurVol'
+        WRITE(IOUT,6005) (IVOL,(INDXSR(IR,IVOL),IR=1,4),SURVOL(IVOL),
+     >                     IVOL=-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((5(I10,','),F20.10,:,','))
+ 6006 FORMAT(4(F20.10,:,','))
+      END