diff options
Diffstat (limited to 'Dragon/src/NXTVCA.f')
| -rw-r--r-- | Dragon/src/NXTVCA.f | 277 |
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 |