path: root/Dragon/src/NXTPHC.f

*DECK NXTPHC
      SUBROUTINE NXTPHC(IPRINT,NDIM  ,MXMESH,MAXSUR,MAXREG,
     >                  MESH  ,DAMESH,NPIN  ,ITPIN ,DPIN  ,
     >                  NBSUR ,NBREG ,INDXSR,SURVOL,POSTRI)
*
*----------
*
*Purpose:
* Remove from the volumes or surfaces
* associated with an annular/hexagonal 2-D or 3-D geometry
* the volumes or surfaces of the overlapping pins.
*
*Copyright:
* Copyright (C) 2010 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$).
* DAMESH  spatial description of the Cartesian geometry.
* NPIN    number of pins to superimpose on geometry.
* ITPIN   type of pin.
* DPIN    pin location and dimensions.
* NBSUR   number of surfaces in the geometry.
* NBREG   final number of non void regions in the geometry.
* 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: input/output
* INDXSR  local indexing of regions.
* SURVOL  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 the DAMESH array:
*     hexagonal mesh is DAMESH(i,1) for i=0,MESH(1);
*     mesh in $Z$  is z(k)=DAMESH(k,3) for k=0,MESH(3);
*  2- Contents of the DPIN array for pin IPIN:
*     -> annular pin
*       ->annular regions in the $X-Y$ plane
*         centre (x,y,z)=(DPIN(0,IPIN)*COS(DPIN(-1,IPIN))
*                         DPIN(0,IPIN)*SIN(DPIN(-1,IPIN)),0.0D0)
*         outer pin radius      r=DPIN(4,IPIN)
*         pin height            dz(iz)=DPIN(3,IPIN)
*  3- 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 DAMESH(-1:MXMESH,4)
      INTEGER          NPIN,ITPIN(3,NPIN)
      DOUBLE PRECISION DPIN(-1:4,NPIN)
      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='NXTPHC')
      INTEGER          MAXDIM
      PARAMETER       (MAXDIM=4)
      DOUBLE PRECISION DCUTOF
      PARAMETER       (DCUTOF=1.0D-8)
      DOUBLE PRECISION DZERO,DONE,DTWO,DHALF,DSQ3O2
      PARAMETER       (DZERO=0.0D0,DONE=1.0D0,DTWO=2.0D0,
     >                 DHALF=0.5D0,DSQ3O2=0.86602540378444D0)
*----
*  Functions
*----
      DOUBLE PRECISION XDRCST,PI
      INTEGER          NXTIAA,NXTITA,NXTPRA
      INTEGER          ITYITP,ITYITA,ITYTAP
      DOUBLE PRECISION VOLITP,VOLITA,VOLTAP
*----
*  Local variables
*----
      INTEGER          NX,NZ,NANN,NRTP,NRP,NSTP,NSP,NRP1,NRTPP,NRTPS,
     >                 ISBOT,ISTOP,IX,IZ,IR,ISECT,IPIN,ISUR,IVOL,ILOCT
      DOUBLE PRECISION ZB,ZT,PPRMIN,PPRMAX,PPPMIN,PPPMAX,DPP
      DOUBLE PRECISION VOLPIN,VOLIAO,VOLIAI
      DOUBLE PRECISION POSPIN(0:2),POSANN(0:2)
      INTEGER          NFACES
      INTEGER          IA
*----
*  Allocatable arrays
*----
      INTEGER, ALLOCATABLE, DIMENSION(:) :: ITYIAP
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: VOLIAP
*----
*  Get dimensioning information
*----
      ALLOCATE(ITYIAP(MXMESH),VOLIAP(MXMESH))
      NFACES=3
      PI=XDRCST('Pi',' ')
      NX=MESH(1)
      NZ=MESH(3)
      NANN=MESH(4)
      NRTP=NX**2
      NRP=6*NRTP
      NSTP=2*NX-1
      NSP=6*NSTP
      NRP1=NANN+1
      NRTPP=NRP*NRP1
      NRTPS=NRTP*NRP1
      ZB=0.0D0
      ZT=0.0D0
      IF(NDIM .EQ. 3) THEN
        ISBOT=-NSP*NZ
        ISTOP=ISBOT-NRTPP
      ELSE
        ISBOT=0
        ISTOP=0
      ENDIF
      POSANN(1)=0.0D0
      POSANN(2)=0.0D0
*----
*  Print mesh if required
*----
      IF(IPRINT .GE. 100) THEN
        WRITE(IOUT,6000) NAMSBR
        WRITE(IOUT,6010) 'H',NX
        WRITE(IOUT,6011) 'MESHH ='
        WRITE(IOUT,6012) (DAMESH(IX,1),IX=-1,2*NX)
        IF(NZ .GT. 0) THEN
          WRITE(IOUT,6010) 'Z',NZ
          WRITE(IOUT,6011) 'MESHZ ='
          WRITE(IOUT,6012) (DAMESH(IZ,3),IZ=-1,NZ)
        ENDIF
          WRITE(IOUT,6010) 'R',NANN
        WRITE(IOUT,6011) 'MESHR ='
        WRITE(IOUT,6012) (DAMESH(IR,4),IR=-1,NANN)
*----
* Pin description
*----
        DO IPIN=1,NPIN
          WRITE(IOUT,6017) 'PinTyp',IPIN,ITPIN(3,IPIN)
          IF(ITPIN(3,IPIN) .EQ. 3) THEN
            WRITE(IOUT,6019) 'PinCXY',IPIN,
     >                     DPIN(0,IPIN)*COS(DPIN(-1,IPIN)),
     >                     DPIN(0,IPIN)*SIN(DPIN(-1,IPIN))
            WRITE(IOUT,6018) 'PinRad',IPIN,DPIN(4,IPIN)
          ENDIF
          IF(NDIM .EQ. 3) THEN
            WRITE(IOUT,6019) 'PinPoZ',IPIN,
     >                    DAMESH(-1,3)-DPIN(3,IPIN)/DTWO,
     >                    DAMESH(-1,3)+DPIN(3,IPIN)/DTWO
          ENDIF
        ENDDO
      ENDIF
*----
*  Loop over pins
*----
      DO IPIN=1,NPIN
*----
*  For 3-D problem,
*  Find pin bottom (ZB) and top (ZT) z location.
*----
        IF(NDIM .EQ. 3) THEN
          ZB=DAMESH(-1,3)-DPIN(3,IPIN)/DTWO
          ZT=ZB+DPIN(3,IPIN)
        ENDIF
*----
*  Annular pin properties
*----
        POSPIN(0)=DPIN(4,IPIN)
        POSPIN(1)=DPIN(0,IPIN)*COS(DPIN(-1,IPIN))
        POSPIN(2)=DPIN(0,IPIN)*SIN(DPIN(-1,IPIN))
        VOLPIN=PI*POSPIN(0)*POSPIN(0)
*----
*  Determine pin-annular regions intersections
*----
        DO IA=1,NANN
          POSANN(0)=DAMESH(IA,4)
          ITYIAP(IA)=0
          VOLIAP(IA)=0.0D0
*----
*  Find 2-D annular region/annular pin intersection
*----
          ITYIAP(IA)=NXTIAA(POSANN,POSPIN,VOLIAP(IA))
        ENDDO
*        IF(NANN .GE. 0) THEN
*          WRITE(IOUT,*) 'Pin-Annular intersection for pin =',IPIN
*          WRITE(IOUT,'(2I10,F20.10)')
*     >     (IA,ITYIAP(IA),VOLIAP(IA),IA=1,NANN)
*        ENDIF
*----
*  1- Loop over sectors
*----
        DO ISECT=1,6
*----
*  Loop over region in sector
*----
          DO IR=1,NRTP
            ITYITP=NXTITA(POSTRI(1,1,IR,ISECT),POSPIN,VOLITP)
*            WRITE(IOUT,'(A25,3I10,F20.10)') 'Pin-Triangle intersection',
*     >      ISECT,IR,ITYITP,VOLITP
            ILOCT=(ISECT-1)*NRTPS+(IR-1)*NRP1
            IF(ITYITP .NE. 0) THEN
              VOLIAO=DZERO
              VOLIAI=DZERO
              DO IA=1,NANN
*                IF(IPRINT .GE. 100) THEN
*                  WRITE(IOUT,6022) IPIN,ISECT,IR,IA,ITYIAP(IA)
*                ENDIF
                POSANN(0)=DAMESH(IA,4)
                VOLIAI=VOLIAO
                IF(ITYIAP(IA) .NE. 0) THEN
                  ITYITA=NXTITA(POSTRI(1,1,IR,ISECT),POSANN,VOLITA)
*                  WRITE(IOUT,*) 'NXTITA ',ITYITA,VOLITA
                  IF(ITYITA .EQ. -1) THEN
*----
*  Partial Cartesian/annular region intersection
*  Examine Rectangle/pin intersection
*  Note: ITYITA=0 already considered above
*----
                    IF(ITYITP .EQ. -1) THEN
*----
*  Partial Cartesian/pin intersection
*  Examine Annular/pin intersection
*  Note: ITYIAP=0 already considered above
*----
                      IF(ITYIAP(IA) .EQ. -1) THEN
*----
*  Partial Annular/pin intersection
*  Find intersection volume of three regions
*----
                        ITYTAP=NXTPRA(NFACES,POSTRI(1,1,IR,ISECT),
     >                                POSANN,POSPIN,VOLTAP)
                        VOLIAO=VOLTAP
                        VOLIAI=VOLIAO-VOLIAI
                      ELSE IF(ITYIAP(IA) .EQ. 1) THEN
*----
*  Annular region in pin
*  Volume is given by Rectangle/annular intersection
*----
                        VOLIAO=VOLITA
                        VOLIAI=VOLIAO-VOLIAI
                      ELSE IF(ITYIAP(IA) .EQ. 2) THEN
*----
*  Annular region contains pin
*  Volume is given by Annular/pin intersection
*----
                        VOLIAO=VOLITP
                        VOLIAI=VOLIAO-VOLIAI
                      ENDIF
                    ELSE IF(ITYITP .EQ. 1) THEN
*----
*  Cartesian region in pin
*  Volume is given by Rectangle/annular intersection
*----
                      VOLIAO=VOLITA
                      VOLIAI=VOLIAO-VOLIAI
                    ELSE IF(ITYITP .EQ. 2) THEN
*----
*  Cartesian region contains pin
*  Volume is given by Annular/pin intersection
*----
                      VOLIAO=VOLIAP(IA)
                      VOLIAI=VOLIAO-VOLIAI
                    ENDIF
                  ELSE IF(ITYITA .EQ. 0) THEN
*----
*  No Cartesian/annular region intersection
*  go to next annular region.
*----
                    GO TO 125
                  ELSE IF(ITYITA .EQ. 1) THEN
*----
*  Cartesian region in annular region
*  Volume is given by Rectangle/pin intersection
*----
                    VOLIAO=VOLITP
                    VOLIAI=VOLIAO-VOLIAI
                  ELSE IF(ITYITA .EQ. 2) THEN
*----
*  Cartesian region contains annular region
*  Volume is given by Annular/pin intersection
*----
                    VOLIAO=VOLIAP(IA)
                    VOLIAI=VOLIAO-VOLIAI
                  ENDIF
                ENDIF
 125            CONTINUE
*----
*  There is an intersection possible between the pin and
*  the triangle.
*  1- Look for botton and top surface in 3-D
*----
*                IF(IPRINT .GE. 100) THEN
*                  WRITE(IOUT,6023) VOLIAI,VOLIAO
*                ENDIF
                IF(NDIM .EQ. 3) THEN
                  IF(ZB .LE. DAMESH(0,3) .AND.
     >               ZT .GE. DAMESH(0,3)) THEN
*----
*  Remove area contribution from bottom surface
*----
                    ISUR=ISBOT-ILOCT-IA
                    SURVOL(ISUR)=SURVOL(ISUR)-VOLIAI
                  ENDIF
                  IF(ZB .LE. DAMESH(MESH(3),3) .AND.
     >               ZT .GE. DAMESH(MESH(3),3)) THEN
*----
*  Remove area contribution from top surface
*----
                    ISUR=ISTOP-ILOCT-IA
                    SURVOL(ISUR)=SURVOL(ISUR)-VOLIAI
                  ENDIF
                  PPPMIN=ZB
                  PPPMAX=PPPMIN+DPIN(3,IPIN)
                  DO IZ=1,MESH(3)
                    PPRMIN=MAX(DAMESH(IZ-1,3),PPPMIN)
                    PPRMAX=MIN(DAMESH(IZ,3),PPPMAX)
                    IF(PPRMIN .LT. PPRMAX) THEN
                      DPP=VOLIAI*(PPRMAX-PPRMIN)
                      IVOL=(IZ-1)*NRTPP+ILOCT+IA
                      SURVOL(IVOL)=SURVOL(IVOL)-DPP
                    ENDIF
                  ENDDO
                ELSE
*                  WRITE(IOUT,'(A12,4I10,3F20.10)') 'Volume id',
*     >ISECT,IR,IA,ILOCT+IA,SURVOL(IVOL),VOLIAI,
*     >SURVOL(IVOL)-VOLIAI
                  IVOL=ILOCT+IA
                  SURVOL(IVOL)=SURVOL(IVOL)-VOLIAI
                ENDIF
*----
*  If pin all extracted, go to next pin
*----
                IF(VOLPIN .EQ. VOLIAO) GO TO 115
*                IF(VOLPIN .LE. DZERO) GO TO 115
              ENDDO
*----
*  Use DELV to correct volumes and surface area for
*  regions inside a rectangle but outside annular ring
*----
*              write(6,*) VOLIAI,VOLIAO,VOLITP,VOLITP-VOLIAI
              VOLIAI=VOLIAO
              VOLIAO=VOLITP
              VOLIAI=VOLIAO-VOLIAI
              IA=NRP1
              IF(NDIM .EQ. 3) THEN
                IF(ZB .LE. DAMESH(0,3) .AND.
     >             ZT .GE. DAMESH(0,3)) THEN
*----
*  Remove area contribution from bottom surface
*----
                  ISUR=ISBOT-ILOCT-IA
                  SURVOL(ISUR)=SURVOL(ISUR)-VOLIAI
                ENDIF
                IF(ZB .LE. DAMESH(NZ,3) .AND.
     >             ZT .GE. DAMESH(NZ,3)) THEN
*----
*  Remove area contribution from top surface
*----
                  ISUR=ISTOP-ILOCT-IA
                  SURVOL(ISUR)=SURVOL(ISUR)-VOLIAI
                ENDIF
                PPPMIN=ZB
                PPPMAX=PPPMIN+DPIN(3,IPIN)
                DO IZ=1,NZ
                  PPRMIN=MAX(DAMESH(IZ-1,3),PPPMIN)
                  PPRMAX=MIN(DAMESH(IZ,3),PPPMAX)
                  IF(PPRMIN .LT. PPRMAX) THEN
                    DPP=VOLIAI*(PPRMAX-PPRMIN)
                    IVOL=(IZ-1)*NRTPP+ILOCT+IA
                    SURVOL(IVOL)=SURVOL(IVOL)-DPP
                  ENDIF
                ENDDO
              ELSE
                IVOL=ILOCT+IA
                SURVOL(IVOL)=SURVOL(IVOL)-VOLIAI
              ENDIF
            ENDIF  
          ENDDO
        ENDDO
 115    CONTINUE
      ENDDO
*----
*  Test for negative surface area and volumes
*----
      DO ISUR=1,NBSUR
        IF(SURVOL(-ISUR) .LT. -DCUTOF) THEN
          WRITE(IOUT,9000) NAMSBR,-ISUR
          WRITE(IOUT,9002) (INDXSR(IR,-ISUR),IR=1,5),SURVOL(-ISUR)
          CALL XABORT(NAMSBR//
     >    ': Region with negative surface area detected')
        ELSE IF(SURVOL(-ISUR) .LT. DCUTOF) THEN
          SURVOL(-ISUR)=DZERO
        ENDIF
      ENDDO
      DO  IVOL=1,NBREG
        IF(SURVOL(IVOL) .LT. -DCUTOF) THEN
          WRITE(IOUT,9001) NAMSBR,IVOL
          WRITE(IOUT,9002) (INDXSR(IR,IVOL),IR=1,5),SURVOL(IVOL)
          CALL XABORT(NAMSBR//
     >    ': Region with negative volume detected')
        ELSE IF(SURVOL(IVOL) .LT. DCUTOF) THEN
          SURVOL(IVOL)=DZERO
        ENDIF
      ENDDO
*----
*  Print volumes if required
*----
      IF(IPRINT .GE. 100) THEN
        WRITE(IOUT,6002) 'SurVol'
        WRITE(IOUT,6005) (IVOL,(INDXSR(IR,IVOL),IR=1,5),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((6(I10,','),F20.10,:,','))
 6010 FORMAT(1X,'MESH DIMENSIONS IN ',A1,' =',I10)
 6011 FORMAT(1X,A7)
 6012 FORMAT(5F20.10)
 6017 FORMAT(A6,I4.4,'=',I10,';')
 6018 FORMAT(A6,I4.4,'=',F20.10,';')
 6019 FORMAT(A6,I4.4,'={',F20.10,',',F20.10,'};')
* 6022 FORMAT('Pin =',I10,' Sector = ',I10,' Region = ',I10,
*     >       ' Annulus =',I10,' Intersection type',I10)
* 6023 FORMAT(20X,' Volume of intersection= ',F20.10,
*     >       ' Volume remaining= ',F20.10)
*----
*  Error and Warning formats
*----
 9000 FORMAT('**** ERROR in -- ',A6,'-- found'/
     >       '     Area of region ',I5,' is negative')
 9001 FORMAT('**** ERROR in -- ',A6,'-- found'/
     >       '     Volume of region ',I5,' is negative')
 9002 FORMAT(5I10,F20.10)
      END