From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Dragon/src/NXTVHT.f | 449 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 449 insertions(+) create mode 100644 Dragon/src/NXTVHT.f (limited to 'Dragon/src/NXTVHT.f') diff --git a/Dragon/src/NXTVHT.f b/Dragon/src/NXTVHT.f new file mode 100644 index 0000000..ac8ba4d --- /dev/null +++ b/Dragon/src/NXTVHT.f @@ -0,0 +1,449 @@ +*DECK NXTVHT + SUBROUTINE NXTVHT(IPRINT,NDIM ,MXMESH,MAXSUR,MAXREG, + > MESH ,DMESH ,NBSUR ,NBREG ,INDXSR,SURVOL) +* +*---------- +* +*Purpose: +* Compute the volume and area associated with each region +* or surface for an 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. +* +*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) +*---- +* Local parameters +*---- + INTEGER IOUT + CHARACTER NAMSBR*6 + PARAMETER (IOUT=6,NAMSBR='NXTVHT') + INTEGER MAXDIM + PARAMETER (MAXDIM=4) + DOUBLE PRECISION DZERO,DONE,DHALF,DSQ3O2 + PARAMETER (DZERO=0.0D0,DONE=1.0D0, + > DHALF=0.5D0,DSQ3O2=0.86602540378444D0) +*---- +* Local variables +*---- + INTEGER NX,NZ,NRTP,NRP,NSTP,NSP,IX,IZ,ISECT,IR,IS, + > IOFS,IVSI,JVSI,IVSIR,NVSIR,ISUVW,IZ1 + DOUBLE PRECISION SIDEU,SIDEUP,SIDER,SIDERP,SIDEL,SIDELP, + > AREAT,AREAL,AREAR,DZ +*---- +* Prepare loop over Z-direction. +*---- + NX=MESH(1) + NZ=MESH(3) + NRTP=NX**2 + NRP=6*NRTP + NSTP=2*NX-1 + NSP=6*NSTP + IF(NDIM .EQ. 3) THEN + NBREG=NRP*NZ + NBSUR=NSP*NZ+2*NRP + ELSE + NBREG=NRP + NBSUR=NSP + ENDIF + IF(IPRINT .GE. 100) THEN + WRITE(IOUT,6000) NAMSBR + WRITE(IOUT,6010) NX,NZ + WRITE(IOUT,6011) 'MESHH =' + WRITE(IOUT,6012) (DMESH(IX,1),IX=-1,2*NX) + IF(NDIM .EQ. 3) THEN + WRITE(IOUT,6011) 'MESHZ =' + WRITE(IOUT,6012) (DMESH(IZ,3),IZ=-1,NZ) + ENDIF + WRITE(IOUT,6002) 'SurVol init' + WRITE(IOUT,6005) (IVSI,(INDXSR(IR,IVSI),IR=1,5),SURVOL(IVSI), + > IVSI=-NBSUR,NBREG) + WRITE(IOUT,6003) + ENDIF + SURVOL(-MAXSUR:MAXREG)=DZERO +*---- +* All the triangles are equilateral and have the same surface area +* except possibly those in the last crown because this last crown +* may be thinner than the other +* 1) Find first the area of generic triangle +*---- + SIDEU=DMESH(NX+1,1)-DMESH(NX,1) + SIDEUP=SIDEU/DSQ3O2 + AREAT=DHALF*SIDEU*SIDEUP + SURVOL(1:NRP)=AREAT +*---- +* Process last crown +* for right triangle remove volume past last cell boundary +* for left triangle volume extends only to cell boundary +* Add surfaces +*---- + SIDEL=DMESH(1,1)-DMESH(0,1) + SIDELP=SIDEL/DSQ3O2 + SIDER=SIDEU-SIDEL + SIDERP=SIDER/DSQ3O2 + AREAL=DHALF*SIDEL*SIDELP + AREAR=AREAT-DHALF*SIDER*SIDERP +*---- +* Finish filling first plan with AREA +*---- + IS=0 + IOFS=NRTP-NSTP + DO ISECT=1,6 +*---- +* right triangles +*---- + DO IX=1,NX-1 + IR=IX+IOFS + IS=IS-1 + SURVOL(IR)=AREAR + SURVOL(IS)=SIDERP + ENDDO +*---- +* Left triangles +*---- + DO IX=NX,NSTP + IR=IX+IOFS + IS=IS-1 + SURVOL(IR)=AREAL + SURVOL(IS)=SIDELP + ENDDO + IOFS=IOFS+NRTP + ENDDO + IF(NDIM .EQ. 3) THEN +*---- +* Hexagonal Surfaces +*---- + DO IZ=NZ,1,-1 + DZ=DMESH(IZ,3)-DMESH(IZ-1,3) + IX=-NSP*(IZ-1) + DO IS=-1,-NSP,-1 + SURVOL(IX+IS)=DZ*SURVOL(IS) + ENDDO + ENDDO +*---- +* Fill bottom face and top face +*---- + IS=-NSP*NZ + DO IR=1,NRP + IS=IS-1 + SURVOL(IS)=SURVOL(IR) + SURVOL(IS-NRP)=SURVOL(IR) + ENDDO +*---- +* Fill all planes with volumes +*---- + DO IZ=NZ,1,-1 + DZ=DMESH(IZ,3)-DMESH(IZ-1,3) + IX=NRP*(IZ-1) +*---- +* Volumes +*---- + DO IR=1,NRP + SURVOL(IX+IR)=DZ*SURVOL(IR) + ENDDO + ENDDO + ENDIF +*---- +* Fill INDXSR to identify regions and volumes +* Process first sector (-pi/6 to pi/6) for volume +* First plane or 2-D +*---- + IF(NDIM .EQ. 3) THEN + IZ1=1 + ELSE + IZ1=0 + ENDIF + IVSI=0 + DO IX=1,NX +*---- +* First line of triangle in the crown +*---- + DO IR=1,IX-1 + IVSI=IVSI+1 + INDXSR(1,IVSI)=NX+IX + INDXSR(5,IVSI)=NX+IX-IR + INDXSR(3,IVSI)=IZ1 + INDXSR(2,IVSI)=NX+IR + ENDDO +*---- +* Second line of triangle in the crown +*---- + DO IR=1,IX + IVSI=IVSI+1 + INDXSR(1,IVSI)=NX+IX + INDXSR(5,IVSI)=NX+IX-IR+1 + INDXSR(3,IVSI)=IZ1 + INDXSR(2,IVSI)=NX+IR + ENDDO + ENDDO +*---- +* Complete for sectors 2-6 on first plane +*---- + NVSIR=IVSI + IVSIR=NVSIR + ISUVW=2*NX+1 + DO ISECT=2,6 + IF(ISECT .EQ. 6) THEN + DO JVSI=1,NVSIR + IVSIR=IVSIR+1 + INDXSR(1,IVSIR)=INDXSR(2,JVSI) + INDXSR(2,IVSIR)=ISUVW-INDXSR(5,JVSI) + INDXSR(3,IVSIR)=IZ1 + INDXSR(5,IVSIR)=INDXSR(1,JVSI) + ENDDO + ELSE IF (ISECT.EQ.5) THEN + DO JVSI=1,NVSIR + IVSIR=IVSIR+1 + INDXSR(1,IVSIR)=ISUVW-INDXSR(5,JVSI) + INDXSR(2,IVSIR)=ISUVW-INDXSR(1,JVSI) + INDXSR(3,IVSIR)=IZ1 + INDXSR(5,IVSIR)=INDXSR(2,JVSI) + ENDDO + ELSE IF (ISECT.EQ.4) THEN + DO JVSI=1,NVSIR + IVSIR=IVSIR+1 + INDXSR(1,IVSIR)=ISUVW-INDXSR(1,JVSI) + INDXSR(2,IVSIR)=ISUVW-INDXSR(2,JVSI) + INDXSR(3,IVSIR)=IZ1 + INDXSR(5,IVSIR)=ISUVW-INDXSR(5,JVSI) + ENDDO + ELSE IF (ISECT.EQ.3) THEN + DO JVSI=1,NVSIR + IVSIR=IVSIR+1 + INDXSR(1,IVSIR)=ISUVW-INDXSR(2,JVSI) + INDXSR(2,IVSIR)=INDXSR(5,JVSI) + INDXSR(3,IVSIR)=IZ1 + INDXSR(5,IVSIR)=ISUVW-INDXSR(1,JVSI) + ENDDO + ELSE IF (ISECT.EQ.2) THEN + DO JVSI=1,NVSIR + IVSIR=IVSIR+1 + INDXSR(1,IVSIR)=INDXSR(5,JVSI) + INDXSR(2,IVSIR)=INDXSR(1,JVSI) + INDXSR(3,IVSIR)=IZ1 + INDXSR(5,IVSIR)=ISUVW-INDXSR(2,JVSI) + ENDDO + ENDIF + ENDDO +*---- +* Process other planes in 3-D +*---- + IVSI=NRP + DO IZ=2,NZ + DO IR=1,NRP + IVSI=IVSI+1 + INDXSR(1,IVSI)=INDXSR(1,IR) + INDXSR(2,IVSI)=INDXSR(2,IR) + INDXSR(3,IVSI)=IZ + INDXSR(5,IVSI)=INDXSR(5,IR) + ENDDO + ENDDO +*---- +* Process surfaces +* First plane or 2-D +*---- + IVSI=0 + IOFS=3-6*NX + IX=NX +*---- +* Sector 1 and 4 +* First line of triangle on the last crown +*---- + DO IR=1,IX-1 + IVSI=IVSI-1 + INDXSR(1,IVSI+IOFS)=-1 + INDXSR(5,IVSI+IOFS)=IR+1 + INDXSR(3,IVSI+IOFS)=IZ1 + INDXSR(2,IVSI+IOFS)=NX+1-IR + INDXSR(1,IVSI)=-2 + INDXSR(5,IVSI)=NX+IX-IR + INDXSR(3,IVSI)=IZ1 + INDXSR(2,IVSI)=NX+IR + ENDDO +*---- +* Second line of triangle on the last crown +*---- + DO IR=1,IX + IVSI=IVSI-1 + INDXSR(1,IVSI+IOFS)=-1 + INDXSR(5,IVSI+IOFS)=IR + INDXSR(3,IVSI+IOFS)=IZ1 + INDXSR(2,IVSI+IOFS)=NX+1-IR + INDXSR(1,IVSI)=-2 + INDXSR(5,IVSI)=NX+IX-IR+1 + INDXSR(3,IVSI)=IZ1 + INDXSR(2,IVSI)=NX+IR + ENDDO +*---- +* Sector 2 and 5 +* First line of triangle on the last crown +*---- + DO IR=1,IX-1 + IVSI=IVSI-1 + INDXSR(1,IVSI+IOFS)=IR+1 + INDXSR(5,IVSI+IOFS)=NX+IR + INDXSR(3,IVSI+IOFS)=IZ1 + INDXSR(2,IVSI+IOFS)=-1 + INDXSR(1,IVSI)=2*NX-IR + INDXSR(5,IVSI)=NX+1-IR + INDXSR(3,IVSI)=IZ1 + INDXSR(2,IVSI)=-2 + ENDDO +*---- +* Second line of triangle on the last crown +*---- + DO IR=1,IX + IVSI=IVSI-1 + INDXSR(1,IVSI+IOFS)=IR + INDXSR(5,IVSI+IOFS)=NX+IR + INDXSR(3,IVSI+IOFS)=IZ1 + INDXSR(2,IVSI+IOFS)=-1 + INDXSR(1,IVSI)=2*NX+1-IR + INDXSR(5,IVSI)=NX+1-IR + INDXSR(3,IVSI)=IZ1 + INDXSR(2,IVSI)=-2 + ENDDO +*---- +* Sector 3 and 6 +* First line of triangle on the last crown +*---- + DO IR=1,IX-1 + IVSI=IVSI-1 + INDXSR(1,IVSI)=NX+1-IR + INDXSR(5,IVSI)=-1 + INDXSR(3,IVSI)=IZ1 + INDXSR(2,IVSI)=2*NX-IR + INDXSR(1,IVSI+IOFS)=NX+IR + INDXSR(5,IVSI+IOFS)=-2 + INDXSR(3,IVSI+IOFS)=IZ1 + INDXSR(2,IVSI+IOFS)=IR+1 + ENDDO +*---- +* Second line of triangle on the last crown +*---- + DO IR=1,IX + IVSI=IVSI-1 + INDXSR(1,IVSI)=NX+1-IR + INDXSR(5,IVSI)=-1 + INDXSR(3,IVSI)=IZ1 + INDXSR(2,IVSI)=2*NX+1-IR + INDXSR(1,IVSI+IOFS)=NX+IR + INDXSR(5,IVSI+IOFS)=-2 + INDXSR(3,IVSI+IOFS)=IZ1 + INDXSR(2,IVSI+IOFS)=IR + ENDDO + IF(NDIM .EQ. 3) THEN +*---- +* Process other planes in 3-D +*---- + IVSI=IVSI-6*NX+3 + DO IZ=2,NZ + DO IR=-1,-NSP,-1 + IVSI=IVSI-1 + INDXSR(1,IVSI)=INDXSR(1,IR) + INDXSR(2,IVSI)=INDXSR(2,IR) + INDXSR(3,IVSI)=IZ + INDXSR(5,IVSI)=INDXSR(5,IR) + ENDDO + ENDDO +*---- +* Process bottom and top faces if required +*---- + DO IR=1,NRP + IVSI=IVSI-1 + IVSIR=IVSI-NRP + INDXSR(1,IVSI)=INDXSR(1,IR) + INDXSR(2,IVSI)=INDXSR(2,IR) + INDXSR(3,IVSI)=-1 + INDXSR(5,IVSI)=INDXSR(5,IR) + INDXSR(1,IVSIR)=INDXSR(1,IR) + INDXSR(2,IVSIR)=INDXSR(2,IR) + INDXSR(3,IVSIR)=-2 + INDXSR(5,IVSIR)=INDXSR(5,IR) + ENDDO + ENDIF +*---- +* Print surfaces and volumes if required and return +*---- + IF(IPRINT .GE. 100) THEN + WRITE(IOUT,6002) 'SurVol' + 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 =',2I10) + 6011 FORMAT(1X,A7) + 6012 FORMAT(5F20.10) + END -- cgit v1.2.3