PierceSurface Subroutine

private subroutine PierceSurface(ISurf, R1, RN, IPIERC, CPhit)

Arguments

Type	Intent		Name
integer,	intent(in)	::	ISurf
real(kind=r64),	intent(in)	::	R1(3)
real(kind=r64),	intent(in)	::	RN(3)
integer,	intent(out)	::	IPIERC
real(kind=r64),	intent(out)	::	CPhit(3)
Source Code

SUBROUTINE PierceSurface(ISurf, R1, RN, IPIERC, CPhit)


          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Fred Winkelmann
          !       DATE WRITTEN   July 1997
          !       MODIFIED       Sept 2003, FCW: modification of Daylighting routine DayltgPierceSurface
          !       RE-ENGINEERED  na


          ! PURPOSE OF THIS SUBROUTINE:
          ! Returns point CPhit that line through point R1 in direction of unit vector RN intersects
          ! the plan of surface ISurf. IPIERC = 1 if CPhit is inside the perimeter of ISurf. If not,
          ! IPIERC = 0. This routine works for convex and concave surfaces with 3 or more vertices.
          !
          ! METHODOLOGY EMPLOYED:


          ! REFERENCES:
          ! Based on DOE-2.1E subroutine DPIERC.


          ! USE STATEMENTS:na


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


          ! SUBROUTINE PARAMETER DEFINITIONS:na
          ! INTERFACE BLOCK SPECIFICATIONS:na
          ! DERIVED TYPE DEFINITIONS:na


          ! SUBROUTINE ARGUMENT DEFINITIONS:
  INTEGER, INTENT(IN)          :: ISurf      ! Surface index
  REAL(r64), INTENT(IN)        :: R1(3)      ! Point from which ray originates
  REAL(r64), INTENT(IN)        :: RN(3)      ! Unit vector along in direction of ray whose
                                             !  intersection with surface is to be determined
  INTEGER, INTENT(OUT)         :: IPIERC     ! =1 if line through point R1 in direction of unit vector
                                             !  RN intersects surface ISurf; =0 otherwise.
  REAL(r64), INTENT(OUT)       :: CPhit(3)   ! Point that ray along RN intersects plane of surface


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
  INTEGER   :: NV                       ! Number of vertices (3 or 4)
  REAL(r64) :: AXC(3)                   ! Cross product of A and C
  REAL(r64) :: SN(3)                    ! Vector normal to surface (SN = A1 X A2)
!unused  REAL(r64) :: AA(3)                    ! AA(I) = A(N,I)
!unused  REAL(r64) :: CC(3)                    ! CC(I) = C(N,I)
  REAL(r64) :: CCC(3)                   ! Vector from vertex 2 to CP
  REAL(r64) :: AAA(3)                   ! Vector from vertex 2 to vertex 1
  REAL(r64) :: BBB(3)                   ! Vector from vertex 2 to vertex 3
  INTEGER   :: N                        ! Vertex loop index
  REAL(r64) :: F1,F2                    ! Intermediate variables
  REAL(r64) :: SCALE                    ! Scale factor
  REAL(r64) :: DOTCB                    ! Dot product of vectors CCC and BBB
  REAL(r64) :: DOTCA                    ! Dot product of vectors CCC and AAA
!unused  REAL(r64) :: DOTAXCSN                 ! Dot product of vectors AXC and SN


  REAL(r64), ALLOCATABLE, SAVE, DIMENSION(:,:) :: V     ! Vertices of surfaces
  REAL(r64), ALLOCATABLE, SAVE, DIMENSION(:,:) :: A     ! Vertex-to-vertex vectors; A(1,i) is from vertex 1 to 2, etc.
  REAL(r64), ALLOCATABLE, SAVE, DIMENSION(:,:) :: C     ! Vectors from vertices to intersection point
  LOGICAL,SAVE :: FirstTime=.true.


          ! FLOW:
  IPIERC = 0


  ! Vertex vectors
  IF (FirstTime) THEN
    ALLOCATE(V(MaxVerticesPerSurface,3))
    V=0.0d0
    ALLOCATE(A(MaxVerticesPerSurface,3))
    A=0.0d0
    ALLOCATE(C(MaxVerticesPerSurface,3))
    C=0.0d0
    FirstTime=.false.
  ENDIF


  ! Set the first two V & A
  V(1,1) = Surface(ISurf)%Vertex(1)%X
  V(1,2) = Surface(ISurf)%Vertex(1)%Y
  V(1,3) = Surface(ISurf)%Vertex(1)%Z


  A(1,1) = Surface(ISurf)%Vertex(2)%X - V(1,1)
  A(1,2) = Surface(ISurf)%Vertex(2)%Y - V(1,2)
  A(1,3) = Surface(ISurf)%Vertex(2)%Z - V(1,3)


  V(2,1) = Surface(ISurf)%Vertex(2)%X
  V(2,2) = Surface(ISurf)%Vertex(2)%Y
  V(2,3) = Surface(ISurf)%Vertex(2)%Z


  A(2,1) = Surface(ISurf)%Vertex(3)%X - V(2,1)
  A(2,2) = Surface(ISurf)%Vertex(3)%Y - V(2,2)
  A(2,3) = Surface(ISurf)%Vertex(3)%Z - V(2,3)


  ! Vector normal to surface
  SN(1) = A(1,2) * A(2,3) - A(1,3) * A(2,2)
  SN(2) = A(1,3) * A(2,1) - A(1,1) * A(2,3)
  SN(3) = A(1,1) * A(2,2) - A(1,2) * A(2,1)


  ! Scale factor, the solution of SN.(CPhit-V2) = 0 and
  ! CPhit = R1 + SCALE*RN, where CPhit is the point that RN,
  ! when extended, intersects the plane of the surface.
  F2 = DOT_PRODUCT(SN, RN)
  IF (ABS(F2) < 0.01d0) RETURN   ! Skip surfaces that are parallel to RN
  F1 = SN(1)*(V(2,1)-R1(1)) + SN(2)*(V(2,2)-R1(2)) + SN(3)*(V(2,3)-R1(3))
  !F1 = DOT_PRODUCT(SN, V2 - R1)
  SCALE = F1 / F2
  IF (SCALE <= 0.0d0) RETURN  ! Skip surfaces that RN points away from
  CPhit = R1 + RN * SCALE  ! Point that RN intersects plane of surface


  ! Two cases: rectangle and non-rectangle; do rectangle
  ! first since most common shape and faster calculation
  IF (Surface(ISurf)%Shape == Rectangle .OR. Surface(ISurf)%Shape == RectangularDoorWindow .or.         &
      Surface(ISurf)%Shape == RectangularOverhang .OR. Surface(ISurf)%Shape == RectangularLeftFin .or.  &
      Surface(ISurf)%Shape == RectangularRightFin) THEN
    !
    ! Surface is rectangular
    !
    ! Vectors from vertex 2 to vertex 1 and vertex 2 to vertex 3


    ! Intersection point, CCC, is inside rectangle if
    ! 0 < CCC.BBB < BBB.BBB AND 0 < CCC.AAA < AAA.AAA


    !CCC = CPhit - V2  ! Vector from vertex 2 to CPhit
    CCC(1) = CPhit(1) - V(2,1)
    CCC(2) = CPhit(2) - V(2,2)
    CCC(3) = CPhit(3) - V(2,3)


    ! Set third V just for here
    V(3,1) = Surface(ISurf)%Vertex(3)%X
    V(3,2) = Surface(ISurf)%Vertex(3)%Y
    V(3,3) = Surface(ISurf)%Vertex(3)%Z


    !BBB = V3 - V2
    BBB(1) = V(3,1) - V(2,1)
    BBB(2) = V(3,2) - V(2,2)
    BBB(3) = V(3,3) - V(2,3)


    DOTCB = DOT_PRODUCT(CCC, BBB)
    IF (DOTCB < 0.0d0) RETURN
    IF (DOTCB > DOT_PRODUCT(BBB,BBB)) RETURN


    !AAA = V1 - V2
    AAA(1) = V(1,1) - V(2,1)
    AAA(2) = V(1,2) - V(2,2)
    AAA(3) = V(1,3) - V(2,3)


    DOTCA = DOT_PRODUCT(CCC, AAA)
    IF (DOTCA < 0.0d0) RETURN
    IF (DOTCA > DOT_PRODUCT(AAA,AAA)) RETURN
    ! Surface is intersected
    IPIERC = 1


 ELSE
    !
    ! Surface is not rectangular
    !


    !if (NV == 3) then
    !else ! NV=4
    !endif


    ! First two of V & A already set
    ! test first vertex:
    C(1,1) = CPhit(1) - V(1,1)
    C(1,2) = CPhit(2) - V(1,2)
    C(1,3) = CPhit(3) - V(1,3)
    AXC(1) = A(1,2) * C(1,3) - A(1,3) * C(1,2)
    AXC(2) = A(1,3) * C(1,1) - A(1,1) * C(1,3)
    AXC(3) = A(1,1) * C(1,2) - A(1,2) * C(1,1)
    IF (DOT_PRODUCT(AXC,SN) < 0.0d0) RETURN  ! If at least one dot product is negative, intersection outside of surface


    ! test second vertex:
    C(2,1) = CPhit(1) - V(2,1)
    C(2,2) = CPhit(2) - V(2,2)
    C(2,3) = CPhit(3) - V(2,3)
    AXC(1) = A(2,2) * C(2,3) - A(2,3) * C(2,2)
    AXC(2) = A(2,3) * C(2,1) - A(2,1) * C(2,3)
    AXC(3) = A(2,1) * C(2,2) - A(2,2) * C(2,1)
    IF (DOT_PRODUCT(AXC,SN) < 0.0d0) RETURN  ! If at least one dot product is negative, intersection outside of surface


    NV = Surface(ISurf)%Sides
    ! Since first two of V & A already set, start with 3.  (so if NV=3, this loop won't happen)
    DO N = 3,NV-1
      V(N,1) = Surface(ISurf)%Vertex(N)%X
      V(N,2) = Surface(ISurf)%Vertex(N)%Y
      V(N,3) = Surface(ISurf)%Vertex(N)%Z

      A(N,1) = Surface(ISurf)%Vertex(N+1)%X - V(N,1)
      A(N,2) = Surface(ISurf)%Vertex(N+1)%Y - V(N,2)
      A(N,3) = Surface(ISurf)%Vertex(N+1)%Z - V(N,3)


      C(N,1) = CPhit(1) - V(N,1)
      C(N,2) = CPhit(2) - V(N,2)
      C(N,3) = CPhit(3) - V(N,3)


      AXC(1) = A(N,2) * C(N,3) - A(N,3) * C(N,2)
      AXC(2) = A(N,3) * C(N,1) - A(N,1) * C(N,3)
      AXC(3) = A(N,1) * C(N,2) - A(N,2) * C(N,1)


      IF (DOT_PRODUCT(AXC,SN) < 0.0d0) RETURN  ! If at least one dot product is negative, intersection outside of surface


    END DO


    ! Last vertex (NV=3 or NV=4)
    V(NV,1) = Surface(ISurf)%Vertex(NV)%X
    V(NV,2) = Surface(ISurf)%Vertex(NV)%Y
    V(NV,3) = Surface(ISurf)%Vertex(NV)%Z

    A(NV,1) = V(1,1) - V(NV,1)
    A(NV,2) = V(1,2) - V(NV,2)
    A(NV,3) = V(1,3) - V(NV,3)

    C(NV,1) = CPhit(1) - V(NV,1)
    C(NV,2) = CPhit(2) - V(NV,2)
    C(NV,3) = CPhit(3) - V(NV,3)


    AXC(1) = A(NV,2) * C(NV,3) - A(NV,3) * C(NV,2)
    AXC(2) = A(NV,3) * C(NV,1) - A(NV,1) * C(NV,3)
    AXC(3) = A(NV,1) * C(NV,2) - A(NV,2) * C(NV,1)

    IF (DOT_PRODUCT(AXC,SN) < 0.0d0) RETURN  ! If at least one dot product is negative, intersection outside of surface


    IPIERC = 1      ! Surface is intersected
  END IF


  RETURN

END SUBROUTINE PierceSurface
All Procedures

PierceSurface Subroutine

Source Code

All Procedures

private subroutine PierceSurface(ISurf, R1, RN, IPIERC, CPhit)

Arguments

Calls

Called By

Source Code

Source Code

All Procedures
