CHKSBS Subroutine

private subroutine CHKSBS(HTS, GRSNR, SBSNR)

Arguments

Type	Intent		Name
integer,	intent(in)	::	HTS
integer,	intent(in)	::	GRSNR
integer,	intent(in)	::	SBSNR
Calls

Help
Called By

Help
Source Code

CHKSBS
Source Code

SUBROUTINE CHKSBS(HTS,GRSNR,SBSNR)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Legacy Code
          !       DATE WRITTEN
          !       MODIFIED       na
          !       RE-ENGINEERED  Lawrie, Oct 2000

          ! PURPOSE OF THIS SUBROUTINE:
          ! Checks that a subsurface is completely
          ! enclosed by its base surface.

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! BLAST/IBLAST code, original author George Walton

          ! 3D Planar Polygons
          ! In 3D applications, one sometimes wants to test a point and polygon that are in the same plane.
          ! For example, one may have the intersection point of a ray with the plane of a polyhedron's face,
          ! and want to test if it is inside the face.  Or one may want to know if the base of a 3D perpendicular
          ! dropped from a point is inside a planar polygon.

          ! 3D inclusion is easily determined by projecting the point and polygon into 2D.  To do this, one simply
          ! ignores one of the 3D coordinates and uses the other two.  To optimally select the coordinate to ignore,
          ! compute a normal vector to the plane, and select the coordinate with the largest absolute value [Snyder & Barr, 1987].
          ! This gives the projection of the polygon with maximum area, and results in robust computations.
          ! John M. Snyder & Alan H. Barr, "Ray Tracing Complex Models Containing Surface Tessellations",
          ! Computer Graphics 21(4), 119-126 (1987) [also in the Proceedings of SIGGRAPH 1987]
          !--- using adapted routine from Triangulation code -- EnergyPlus.


          ! USE STATEMENTS:
          ! na

  IMPLICIT NONE    ! Enforce explicit typing of all variables in this routine

          ! SUBROUTINE ARGUMENT DEFINITIONS:
  INTEGER, INTENT(IN) :: HTS      ! Heat transfer surface number of the general receiving surf
  INTEGER, INTENT(IN) :: GRSNR    ! Surface number of general receiving surface
  INTEGER, INTENT(IN) :: SBSNR    ! Surface number of subsurface

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! MSG - for error message
  CHARACTER(len=8), PARAMETER, DIMENSION(4) :: MSG = (/'misses  ','        ','within  ','overlaps'/)
  CHARACTER(len=*), PARAMETER :: AFormat='(1X,A)'
  CHARACTER(len=*), PARAMETER :: A4Format='(1X,A,A,A,A)'
  CHARACTER(len=*), PARAMETER :: A2I2Format='(1X,A,I5,A,I5)'
  CHARACTER(len=*), PARAMETER :: VFormat='(1X,A,I5,A,"(",2(F15.2,","),F15.2,")")'

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
  INTEGER N        ! Loop Control
  INTEGER NVT      ! Number of vertices
  REAL(r64), ALLOCATABLE, DIMENSION(:), SAVE :: XVT      ! X Vertices of
  REAL(r64), ALLOCATABLE, DIMENSION(:), SAVE :: YVT      ! Y vertices of
  REAL(r64), ALLOCATABLE, DIMENSION(:), SAVE :: ZVT      ! Z vertices of

  INTEGER NS1      ! Number of the figure being overlapped
  INTEGER NS2      ! Number of the figure doing overlapping
  INTEGER NS3      ! Location to place results of overlap

  LOGICAL, SAVE :: OneTimeFlag=.true.
  logical :: inside

  LOGICAL Out
  REAL(r64) X1,Y1,Z1,X2,Y2,Z2,BX,BY,BZ,BMAX !,SX,SY,SZ
!  INTEGER M

  IF (OneTimeFlag) THEN
    ALLOCATE(XVT(MaxVerticesPerSurface+1))
    ALLOCATE(YVT(MaxVerticesPerSurface+1))
    ALLOCATE(ZVT(MaxVerticesPerSurface+1))
    XVT=0.0d0
    YVT=0.0d0
    ZVT=0.0d0
    OneTimeFlag=.false.
  ENDIF

  NS1    = 1
  NS2    = 2
  NS3    = 3
  HCT(1) = 0.0d0
  HCT(2) = 0.0d0

          ! Put coordinates of base surface into clockwise sequence on the x'-y' plane.

  XVT=0.0d0
  YVT=0.0d0
  ZVT=0.0d0
  XVS=0.0d0
  YVS=0.0d0
  CALL CTRANS(GRSNR,HTS,NVT,XVT,YVT,ZVT)
  DO N = 1, NVT
    XVS(N) = XVT(NVT+1-N)
    YVS(N) = YVT(NVT+1-N)
  END DO

  CALL HTRANS1(NS2,NVT)

          ! Put coordinates of the subsurface into clockwise sequence.

  NVS=Surface(SBSNR)%Sides
  DO N = 1, NVS
    XVS(N) = ShadeV(SBSNR)%XV(NVS+1-N)
    YVS(N) = ShadeV(SBSNR)%YV(NVS+1-N)
  END DO
  CALL HTRANS1(NS1,NVS)

          ! Determine the overlap condition.

  CALL DeterminePolygonOverlap(NS1,NS2,NS3)

          ! Print error condition if necessary.

  IF (OverlapStatus /= FirstSurfWithinSecond) THEN
      OUT=.FALSE.
!C                            COMPUTE COMPONENTS OF VECTOR
!C                            NORMAL TO BASE SURFACE.
      X1 = Surface(GRSNR)%Vertex(1)%X-Surface(GRSNR)%Vertex(2)%X !XV(1,GRSNR)-XV(2,GRSNR)
      Y1 = Surface(GRSNR)%Vertex(1)%Y-Surface(GRSNR)%Vertex(2)%Y !YV(1,GRSNR)-YV(2,GRSNR)
      Z1 = Surface(GRSNR)%Vertex(1)%Z-Surface(GRSNR)%Vertex(2)%Z !ZV(1,GRSNR)-ZV(2,GRSNR)
      X2 = Surface(GRSNR)%Vertex(3)%X-Surface(GRSNR)%Vertex(2)%X !XV(3,GRSNR)-XV(2,GRSNR)
      Y2 = Surface(GRSNR)%Vertex(3)%Y-Surface(GRSNR)%Vertex(2)%Y !YV(3,GRSNR)-YV(2,GRSNR)
      Z2 = Surface(GRSNR)%Vertex(3)%Z-Surface(GRSNR)%Vertex(2)%Z !ZV(3,GRSNR)-ZV(2,GRSNR)
      BX = Y1*Z2-Y2*Z1
      BY = Z1*X2-Z2*X1
      BZ = X1*Y2-X2*Y1
!C                            FIND LARGEST COMPONENT.
      BMAX=MAX(ABS(BX),ABS(BY),ABS(BZ))
!C
      IF(ABS(BX).EQ.BMAX) THEN
!        write(outputfiledebug,*) ' looking bx-bmax',bmax
        DO N=1,Surface(SBSNR)%Sides !NV(SBSNR)
          inside=polygon_contains_point(Surface(GRSNR)%Sides,Surface(GRSNR)%Vertex,Surface(SBSNR)%Vertex(N),.true.,.false.,.false.)
          IF (.not. inside) THEN
            OUT=.true.
!            do m=1,surface(grsnr)%sides
!            write(outputfiledebug,*) 'grsnr,side=',m,surface(grsnr)%vertex
!            write(outputfiledebug,*) 'point outside=',surface(sbsnr)%vertex(n)
!            enddo
!            EXIT
          ENDIF
!          Y1 = Surface(GRSNR)%Vertex(Surface(GRSNR)%Sides)%Y-Surface(SBSNR)%Vertex(N)%Y !YV(NV(GRSNR),GRSNR)-YV(N,SBSNR)
!          Z1 = Surface(GRSNR)%Vertex(Surface(GRSNR)%Sides)%Z-Surface(SBSNR)%Vertex(N)%Z !ZV(NV(GRSNR),GRSNR)-ZV(N,SBSNR)
!          DO M=1,Surface(GRSNR)%Sides !NV(GRSNR)
!            Y2 = Y1
!            Z2 = Z1
!            Y1 = Surface(GRSNR)%Vertex(M)%Y-Surface(SBSNR)%Vertex(N)%Y !YV(M,GRSNR)-YV(N,SBSNR)
!            Z1 = Surface(GRSNR)%Vertex(M)%Z-Surface(SBSNR)%Vertex(N)%Z !ZV(M,GRSNR)-ZV(N,SBSNR)
!            SX = Y1*Z2-Y2*Z1
!            IF(SX*BX.LT.-1.0E-6) THEN
!              OUT=.TRUE.
!              write(outputfiledebug,*) 'sx*bx=',sx*bx
!              write(outputfiledebug,*) 'grsnr=',surface(grsnr)%vertex(m)
!              write(outputfiledebug,*) 'sbsnr=',surface(sbsnr)%vertex(n)
!            endif
!          ENDDO
!          IF (OUT) EXIT
        ENDDO
      ELSE IF(ABS(BY).EQ.BMAX) THEN
!        write(outputfiledebug,*) ' looking by-bmax',bmax
        DO N=1,Surface(SBSNR)%Sides !NV(SBSNR)
          inside=polygon_contains_point(Surface(GRSNR)%Sides,Surface(GRSNR)%Vertex,Surface(SBSNR)%Vertex(N),.false.,.true.,.false.)
          IF (.not. inside) THEN
            OUT=.true.
!            do m=1,surface(grsnr)%sides
!            write(outputfiledebug,*) 'grsnr,side=',m,surface(grsnr)%vertex
!            write(outputfiledebug,*) 'point outside=',surface(sbsnr)%vertex(n)
!            enddo
!            EXIT
          ENDIF
!          Z1 = Surface(GRSNR)%Vertex(Surface(GRSNR)%Sides)%Z-Surface(SBSNR)%Vertex(N)%Z !ZV(NV(GRSNR),GRSNR)-ZV(N,SBSNR)
!          X1 = Surface(GRSNR)%Vertex(Surface(GRSNR)%Sides)%X-Surface(SBSNR)%Vertex(N)%X !XV(NV(GRSNR),GRSNR)-XV(N,SBSNR)
!          DO M=1,Surface(GRSNR)%Sides !NV(GRSNR)
!            Z2 = Z1
!            X2 = X1
!            Z1 = Surface(GRSNR)%Vertex(M)%Z-Surface(SBSNR)%Vertex(N)%Z !ZV(M,GRSNR)-ZV(N,SBSNR)
!            X1 = Surface(GRSNR)%Vertex(M)%X-Surface(SBSNR)%Vertex(N)%X !XV(M,GRSNR)-XV(N,SBSNR)
!            SY = Z1*X2-Z2*X1
!            IF(SY*BY.LT.-1.0E-6) THEN
!              OUT=.TRUE.
!              write(outputfiledebug,*) 'sy*by=',sy*by
!              write(outputfiledebug,*) 'grsnr=',surface(grsnr)%vertex(m)
!              write(outputfiledebug,*) 'sbsnr=',surface(sbsnr)%vertex(n)
!            ENDIF
!          ENDDO
!          IF (OUT) EXIT
        ENDDO
      ELSE
!        write(outputfiledebug,*) ' looking bz-bmax',bmax
        DO N=1,Surface(SBSNR)%Sides !NV(SBSNR)
          inside=polygon_contains_point(Surface(GRSNR)%Sides,Surface(GRSNR)%Vertex,Surface(SBSNR)%Vertex(N),.false.,.false.,.true.)
          IF (.not. inside) THEN
            OUT=.true.
!            do m=1,surface(grsnr)%sides
!            write(outputfiledebug,*) 'grsnr,side=',m,surface(grsnr)%vertex
!            write(outputfiledebug,*) 'point outside=',surface(sbsnr)%vertex(n)
!            enddo
!            EXIT
          ENDIF
!          X1 = Surface(GRSNR)%Vertex(Surface(GRSNR)%Sides)%X-Surface(SBSNR)%Vertex(N)%X !XV(NV(GRSNR),GRSNR)-XV(N,SBSNR)
!          Y1 = Surface(GRSNR)%Vertex(Surface(GRSNR)%Sides)%Y-Surface(SBSNR)%Vertex(N)%Y !YV(NV(GRSNR),GRSNR)-YV(N,SBSNR)
!          DO M=1,Surface(GRSNR)%Sides !NV(GRSNR)
!            X2 = X1
!            Y2 = Y1
!            X1 = Surface(GRSNR)%Vertex(M)%X-Surface(SBSNR)%Vertex(N)%X !XV(M,GRSNR)-XV(N,SBSNR)
!            Y1 = Surface(GRSNR)%Vertex(M)%Y-Surface(SBSNR)%Vertex(N)%Y !YV(M,GRSNR)-YV(N,SBSNR)
!            SZ = X1*Y2-X2*Y1
!            IF(SZ*BZ.LT.-1.0E-6) THEN
!              OUT=.TRUE.
!              write(outputfiledebug,*) 'sz*bz=',sz*bz
!              write(outputfiledebug,*) 'grsnr=',surface(grsnr)%vertex(m)
!              write(outputfiledebug,*) 'sbsnr=',surface(sbsnr)%vertex(n)
!            ENDIF
!          ENDDO
!          IF (OUT) EXIT
        ENDDO
      END IF
!    CALL ShowWarningError('Base surface does not surround subsurface (CHKSBS), Overlap Status='//  &
!                           TRIM(cOverLapStatus(OverLapStatus)))
!    CALL ShowContinueError('Surface "'//TRIM(Surface(GRSNR)%Name)//'" '//TRIM(MSG(OverlapStatus))//  &
!                     ' SubSurface "'//TRIM(Surface(SBSNR)%Name)//'"')
!    IF (FirstSurroundError) THEN
!      CALL ShowWarningError('Base Surface does not surround subsurface errors occuring...'//  &
!                     'Check that the SurfaceGeometry object is expressing the proper starting corner and '//  &
!                     'direction [CounterClockwise/Clockwise]')
!      FirstSurroundError=.false.
!    ENDIF
IF (Out) THEN
    NumBaseSubSurround=NumBaseSubSurround+1
    IF (NumBaseSubSurround > 1) THEN
      ALLOCATE(TempSurfErrorTracking(NumBaseSubSurround))
      TempSurfErrorTracking(1:NumBaseSubSurround-1)=TrackBaseSubSurround
      DEALLOCATE(TrackBaseSubSurround)
      ALLOCATE(TrackBaseSubSurround(NumBaseSubSurround))
      TrackBaseSubSurround=TempSurfErrorTracking
      DEALLOCATE(TempSurfErrorTracking)
    ELSE
      ALLOCATE(TrackBaseSubSurround(NumBaseSubSurround))
    ENDIF
    TrackBaseSubSurround(NumBaseSubSurround)%SurfIndex1=GRSNR
    TrackBaseSubSurround(NumBaseSubSurround)%SurfIndex2=SBSNR
    TrackBaseSubSurround(NumBaseSubSurround)%MiscIndex=OverLapStatus
!    CALL ShowRecurringWarningErrorAtEnd('Base surface does not surround subsurface (CHKSBS), Overlap Status='//  &
!                       TRIM(cOverLapStatus(OverLapStatus)), &
!                       TrackBaseSubSurround(GRSNR)%ErrIndex1)
!    CALL ShowRecurringContinueErrorAtEnd('Surface "'//TRIM(Surface(GRSNR)%Name)//'" '//TRIM(MSG(OverlapStatus))//  &
!                       ' SubSurface "'//TRIM(Surface(SBSNR)%Name)//'"',  &
!                      TrackBaseSubSurround(SBSNR)%ErrIndex2)

    write(OutputFileShading,AFormat) '==== Base does not Surround subsurface details ===='
    write(OutputFileShading,A4Format) 'Surface=',TRIM(Surface(GRSNR)%Name),' ',TRIM(cOverLapStatus(OverLapStatus))
    write(OutputFileShading,A2I2Format) 'Surface#=',GRSNR,' NSides=',Surface(GRSNR)%Sides
    DO N=1,Surface(GRSNR)%Sides
      write(OutputFileShading,VFormat) 'Vertex ',N,'=',Surface(GRSNR)%Vertex(N)
    ENDDO
    write(OutputFileShading,A4Format) 'SubSurface=',TRIM(Surface(SBSNR)%Name)
    write(OutputFileShading,A2I2Format) 'Surface#=',SBSNR,' NSides=',Surface(SBSNR)%Sides
    DO N=1,Surface(SBSNR)%Sides
      write(OutputFileShading,VFormat) 'Vertex ',N,'=',Surface(SBSNR)%Vertex(N)
    ENDDO
    write(OutputFileShading,AFormat) '================================'
  END IF
ENDIF
  RETURN

contains

function polygon_contains_point ( nsides, polygon_3d, point_3d, ignorex, ignorey, ignorez) result(inside)

          ! Function information:
          !       Author         Linda Lawrie
          !       Date written   October 2005
          !       Modified       na
          !       Re-engineered  na

          ! Purpose of this function:
          ! Determine if a point is inside a simple 2d polygon.  For a simple polygon (one whose
          ! boundary never crosses itself).  The polygon does not need to be convex.

          ! Methodology employed:
          ! <Description>

          ! References:
          ! M Shimrat, Position of Point Relative to Polygon, ACM Algorithm 112,
          ! Communications of the ACM, Volume 5, Number 8, page 434, August 1962.

          ! Use statements:
  Use DataVectorTypes

  Implicit none ! Enforce explicit typing of all variables in this routine

          ! Function argument definitions:
  integer :: nsides      ! number of sides (vertices)
  type(vector) :: polygon_3d(nsides)  ! points of polygon
  type(vector) :: point_3d  ! point to be tested
  logical :: ignorex
  logical :: ignorey
  logical :: ignorez
  logical :: inside  ! return value, true=inside, false = not inside

          ! Function parameter definitions:
  REAL(r64), parameter :: point_tolerance=.00001d0

          ! Interface block specifications:
          ! na

          ! Derived type definitions:
          ! na

          ! Function local variable declarations:
  integer i
  integer ip1
  type(vector_2d) polygon(nsides)
  type(vector_2d) point

  inside = .false.
  if (ignorex) then
    polygon%x=polygon_3d%y
    polygon%y=polygon_3d%z
    point%x=point_3d%y
    point%y=point_3d%z
  elseif (ignorey) then
    polygon%x=polygon_3d%x
    polygon%y=polygon_3d%z
    point%x=point_3d%x
    point%y=point_3d%z
  elseif (ignorez) then
    polygon%x=polygon_3d%x
    polygon%y=polygon_3d%y
    point%x=point_3d%x
    point%y=point_3d%y
  endif

  do i = 1, nsides

    if ( i < nsides ) then
      ip1 = i + 1
    else
      ip1 = 1
    end if

    if ( ( polygon(i)%y   <  point%y .and. point%y <= polygon(ip1)%y   ) .or. &
         ( point%y <= polygon(i)%y   .and. polygon(ip1)%y   < point%y ) ) then
      if ( ( point%x - polygon(i)%x ) - ( point%y - polygon(i)%y ) &
         * ( polygon(ip1)%x - polygon(i)%x ) / ( polygon(ip1)%y - polygon(i)%y ) < 0 ) then
        inside = .not. inside
      end if
    end if

  end do

  return

end function polygon_contains_point

END SUBROUTINE CHKSBS
All Procedures
