Nodes of different colours represent the following:
Solid arrows point from a parent (sub)module to the submodule which is descended from it. Dashed arrows point from a module being used to the module or program unit using it. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | SurfNum | |||
| integer, | intent(in) | :: | NSides |
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
SUBROUTINE CheckConvexity(SurfNum,NSides)
! SUBROUTINE INFORMATION:
! AUTHOR Tyler Hoyt
! DATE WRITTEN December 2010
! MODIFIED CR8752 - incorrect note of non-convex polygons
! RE-ENGINEERED na
! PURPOSE OF THIS SUBROUTINE: This subroutine verifies the convexity of a
! surface that is exposed to the sun in the case that full shading calculations
! are required. The calculation conveniently detects collinear points as well,
! and returns a list of indices that are collinear within the plane of the surface.
! METHODOLOGY EMPLOYED: First the surface is determined to have dimension 2 in
! either the xy, yz, or xz plane. That plane is selected to do the testing.
! Vectors representing the edges of the polygon and the perpendicular dot product
! between adjacent edges are computed. This allows the turning angle to be determined.
! If the turning angle is greater than pi/2, it turns to the right, and if it is
! less than pi/2, it turns left. The direction of the turn is stored, and if it
! changes as the edges are iterated the surface is not convex. Meanwhile it stores
! the indices of vertices that are collinear and are later removed.
! REFERENCES:
! http://mathworld.wolfram.com/ConvexPolygon.html
! USE STATEMENTS:
USE General, ONLY: RoundSigDigits
USE DataErrorTracking
IMPLICIT NONE ! Enforce explicit typing of all variables in this routine
! SUBROUTINE ARGUMENT DEFINITIONS:
INTEGER, INTENT(IN) :: SurfNum ! Current surface number
INTEGER, INTENT(IN) :: NSides ! Number of sides to figure
! SUBROUTINE PARAMETER DEFINITIONS:
REAL(r64),PARAMETER :: TurnThreshold = 0.000001d0 ! Sensitivity of convexity test, in radians
CHARACTER(len=*), PARAMETER :: ErrFmt="(' (',F8.3,',',F8.3,',',F8.3,')')"
! INTERFACE BLOCK SPECIFICATIONS
! na
! DERIVED TYPE DEFINITIONS
! na
! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
INTEGER :: N ! Loop index
INTEGER :: Np1 ! Loop index
INTEGER :: Np2 ! Loop index
REAL(r64) :: Det ! Determinant for picking projection plane
REAL(r64) :: DotProd ! Dot product for determining angle
REAL(r64) :: Theta ! Angle between edge vectors
REAL(r64) :: LastTheta ! Angle between edge vectors
REAL(r64) :: V1len ! Edge vector length
REAL(r64) :: V2len ! Edge vector length
LOGICAL :: SignFlag ! Direction of edge turn : .true. is right, .false. is left
LOGICAL :: PrevSignFlag ! Container for the sign of the previous iteration's edge turn
LOGICAL :: FirstTimeFlag ! Flag indicating first iteration
REAL(r64), ALLOCATABLE, SAVE, DIMENSION(:) :: X ! containers for x,y,z vertices of the surface
REAL(r64), ALLOCATABLE, SAVE, DIMENSION(:) :: Y
REAL(r64), ALLOCATABLE, SAVE, DIMENSION(:) :: Z
REAL(r64), ALLOCATABLE, SAVE, DIMENSION(:) :: A ! containers for convexity test
REAL(r64), ALLOCATABLE, SAVE, DIMENSION(:) :: B
INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: SurfCollinearVerts ! Array containing indices of collinear vertices
INTEGER, SAVE :: VertSize ! size of X,Y,Z,A,B arrays
REAL(r64) :: cosarg
INTEGER :: M ! Array index for SurfCollinearVerts container
INTEGER :: J ! Loop index
INTEGER :: K ! Loop index
INTEGER :: Ind ! Location of surface vertex to be removed
LOGICAL, SAVE :: firstTime=.true.
REAL(r64), SAVE :: ACosZero ! set on firstTime
LOGICAL :: SurfCollinearWarning
CHARACTER(len=132) :: ErrLineOut=Blank ! Character string for producing error messages
if (firstTime) then
ACosZero=ACOS(0.0d0)
ALLOCATE(X(MaxVerticesPerSurface+2))
ALLOCATE(Y(MaxVerticesPerSurface+2))
ALLOCATE(Z(MaxVerticesPerSurface+2))
ALLOCATE(A(MaxVerticesPerSurface+2))
ALLOCATE(B(MaxVerticesPerSurface+2))
ALLOCATE(SurfCollinearVerts(MaxVerticesPerSurface))
VertSize=MaxVerticesPerSurface
firstTime=.false.
endif
if (NSides > VertSize) then
DEALLOCATE(X)
DEALLOCATE(Y)
DEALLOCATE(Z)
DEALLOCATE(A)
DEALLOCATE(B)
DEALLOCATE(SurfCollinearVerts)
ALLOCATE(X(MaxVerticesPerSurface+2))
ALLOCATE(Y(MaxVerticesPerSurface+2))
ALLOCATE(Z(MaxVerticesPerSurface+2))
ALLOCATE(A(MaxVerticesPerSurface+2))
ALLOCATE(B(MaxVerticesPerSurface+2))
ALLOCATE(SurfCollinearVerts(MaxVerticesPerSurface))
VertSize=MaxVerticesPerSurface
endif
DO N = 1, NSides
X(N) = SurfaceTmp(SurfNum)%Vertex(N)%X
Y(N) = SurfaceTmp(SurfNum)%Vertex(N)%Y
Z(N) = SurfaceTmp(SurfNum)%Vertex(N)%Z
END DO
X(Nsides+1) = SurfaceTmp(SurfNum)%Vertex(1)%X
Y(Nsides+1) = SurfaceTmp(SurfNum)%Vertex(1)%Y
Z(Nsides+1) = SurfaceTmp(SurfNum)%Vertex(1)%Z
X(Nsides+2) = SurfaceTmp(SurfNum)%Vertex(2)%X
Y(Nsides+2) = SurfaceTmp(SurfNum)%Vertex(2)%Y
Z(Nsides+2) = SurfaceTmp(SurfNum)%Vertex(2)%Z
! Determine a suitable plane in which to do the tests
Det = 0.0d0
DO N = 1, NSides
Det = Det + X(N)*Y(N+1) - X(N+1)*Y(N)
END DO
IF (ABS(Det) > 1.d-4) THEN
A = X
B = Y
ELSE
Det = 0.0d0
DO N = 1, NSides
Det = Det + X(N)*Z(N+1) - X(N+1)*Z(N)
END DO
IF (ABS(Det) > 1.d-4) THEN
A = X
B = Z
ELSE
Det = 0.0d0
DO N = 1, NSides
Det = Det + Y(N)*Z(N+1) - Y(N+1)*Z(N)
END DO
IF (ABS(Det) > 1.d-4) THEN
A = Y
B = Z
ELSE
! This condition should not be reached if the surfaces are guaranteed to be planar already
CALL ShowSevereError('CheckConvexity: Surface="'//TRIM(SurfaceTmp(SurfNum)%Name)// &
'" is non-planar.')
CALL ShowContinueError('Coincident Vertices will be removed as possible.')
DO N=1,SurfaceTmp(SurfNum)%Sides
WRITE(ErrLineOut,ErrFmt) SurfaceTmp(SurfNum)%Vertex(N)
CALL ShowContinueError(ErrLineOut)
ENDDO
END IF
END IF
END IF
M = 0
FirstTimeFlag = .true.
SurfCollinearWarning=.false.
DO N = 1, NSides ! perform convexity test in the plane determined above.
DotProd = (A(N+1)-A(N))*(B(N+2)-B(N+1))-(B(N+1)-B(N))*(A(N+2)-A(N+1))
V1len = SQRT((A(N+1)-A(N))**2+(B(N+1)-B(N))**2)
V2len = SQRT((A(N+2)-A(N+1))**2+(B(N+2)-B(N+1))**2)
IF (V1Len <= 1.d-8 .or. V2Len <= 1.d-8) CYCLE
cosarg=DotProd/(V1len * V2len)
if (cosarg < -1.0d0) then
cosarg=-1.0d0
elseif (cosarg > 1.0d0) then
cosarg=1.0d0
endif
Theta = ACOS(cosarg)
IF (Theta < (ACosZero-TurnThreshold)) THEN
SignFlag = .true.
ELSE
IF (Theta > (ACosZero+TurnThreshold)) THEN
SignFlag = .false.
ELSE ! Store the index of the collinear vertex for removal
IF (.not. SurfCollinearWarning) THEN
IF (DisplayExtraWarnings) THEN
CALL ShowWarningError('CheckConvexity: Surface="'//TRIM(SurfaceTmp(SurfNum)%Name)// &
'", Collinear points have been removed.')
ENDIF
SurfCollinearWarning=.true.
ENDIF
TotalCoincidentVertices=TotalCoincidentVertices+1
M = M + 1
SurfCollinearVerts(M) = N + 1
CYCLE
END IF
END IF
IF (FirstTimeFlag) THEN
PrevSignFlag = SignFlag
FirstTimeFlag = .false.
CYCLE
END IF
IF (SignFlag .neqv. PrevSignFlag) THEN
IF (SolarDistribution /= MinimalShadowing .and. SurfaceTmp(SurfNum)%ExtSolar) THEN
IF (DisplayExtraWarnings) THEN
CALL ShowWarningError('CheckConvexity: Zone="'//trim(Zone(SurfaceTmp(SurfNum)%Zone)%Name)// &
'", Surface="'//TRIM(SurfaceTmp(SurfNum)%Name)//'" is non-convex.')
Np1=N+1
IF (Np1 > NSides) Np1=Np1-NSides
Np2=N+2
IF (Np2 > NSides) Np2=Np2-NSides
CALL ShowContinueError('...vertex '//trim(RoundSigDigits(N))//' to vertex '//trim(RoundSigDigits(Np1))// &
' to vertex '//trim(RoundSigDigits(Np2)))
CALL ShowContinueError('...vertex '//trim(RoundSigDigits(N))//'=['// &
trim(RoundSigDigits(X(N),2))//','// &
trim(RoundSigDigits(Y(N),2))//','// &
trim(RoundSigDigits(Z(N),2))//']')
CALL ShowContinueError('...vertex '//trim(RoundSigDigits(Np1))//'=['// &
trim(RoundSigDigits(X(N+1),2))//','// &
trim(RoundSigDigits(Y(N+1),2))//','// &
trim(RoundSigDigits(Z(N+1),2))//']')
CALL ShowContinueError('...vertex '//trim(RoundSigDigits(Np2))//'=['// &
trim(RoundSigDigits(X(N+2),2))//','// &
trim(RoundSigDigits(Y(N+2),2))//','// &
trim(RoundSigDigits(Z(N+2),2))//']')
! CALL ShowContinueError('...theta angle=['//trim(RoundSigDigits(Theta,6))//']')
! CALL ShowContinueError('...last theta angle=['//trim(RoundSigDigits(LastTheta,6))//']')
ENDIF
SurfaceTmp(SurfNum)%IsConvex=.false.
EXIT
END IF
END IF
PrevSignFlag = SignFlag
LastTheta=Theta
END DO
! must check to make sure don't remove NSides below 3
IF (M > 0) THEN ! Remove the collinear points determined above
IF (NSides-M > 2) THEN
SurfaceTmp(SurfNum)%Sides = NSides - M
ELSE ! too many
IF (DisplayExtraWarnings) THEN
CALL ShowWarningError('CheckConvexity: Surface="'//TRIM(SurfaceTmp(SurfNum)%Name)//'" has ['// &
trim(RoundSigDigits(M))//'] collinear points.')
CALL ShowContinueError('...too many to remove all. Will leave the surface with 3 sides. '// &
'But this is now a degenerate surface')
ENDIF
TotalDegenerateSurfaces=TotalDegenerateSurfaces+1
SurfaceTmp(SurfNum)%Sides = MAX(NSides-M,3)
M=NSides-SurfaceTmp(SurfNum)%Sides
ENDIF
DO J = 1, M
Ind = SurfCollinearVerts(J)
DO K = Ind, NSides - J
SurfaceTmp(SurfNum)%Vertex(K-J+1)%X = SurfaceTmp(SurfNum)%Vertex(K-J+2)%X
SurfaceTmp(SurfNum)%Vertex(K-J+1)%Y = SurfaceTmp(SurfNum)%Vertex(K-J+2)%Y
SurfaceTmp(SurfNum)%Vertex(K-J+1)%Z = SurfaceTmp(SurfNum)%Vertex(K-J+2)%Z
END DO
END DO
END IF
END SUBROUTINE CheckConvexity