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.
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 CalcSurfaceCentroid
! SUBROUTINE INFORMATION:
! AUTHOR B. Griffith
! DATE WRITTEN Feb. 2004
! MODIFIED na
! RE-ENGINEERED na
! PURPOSE OF THIS SUBROUTINE:
! compute centroid of all the surfaces in the main
! surface structure. Store the vertex coordinates of
! the centroid in the 'SURFACE' derived type.
! METHODOLOGY EMPLOYED:
! The centroid of triangle is easily computed by averaging the vertices
! The centroid of a quadrilateral is computed by area weighting the centroids
! of two triangles.
! (Algorithm would need to be changed for higher order
! polygons with more than four sides).
! REFERENCES:
! na
! USE STATEMENTS:
USE Vectors
USE General, ONLY: RoundSigDigits
IMPLICIT NONE ! Enforce explicit typing of all variables in this routine
! SUBROUTINE ARGUMENT DEFINITIONS:
! na
! SUBROUTINE PARAMETER DEFINITIONS:
! na
! INTERFACE BLOCK SPECIFICATIONS
! na
! DERIVED TYPE DEFINITIONS
! na
! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
Type(vector), Dimension(3) :: Triangle1 !working struct for a 3-sided surface
Type(vector), Dimension(3) :: Triangle2 !working struct for a 3-sided surface
INTEGER :: thisSurf ! working variable for do loop
REAL(r64) :: Tri1Area ! working variable for denominator
REAL(r64) :: Tri2Area ! working variable for denominator
REAL(r64) :: TotalArea ! working variable for denominator
REAL(r64) :: Xcm ! temporary X coord for centriod
REAL(r64) :: Ycm ! temporary Y coord for centriod
REAL(r64) :: Zcm ! temporary Z coord for centriod
REAL(r64) :: XcmTri1 ! temporary X coord for centriod of triangle 1
REAL(r64) :: YcmTri1 ! temporary Y coord for centriod of triangle 1
REAL(r64) :: ZcmTri1 ! temporary Z coord for centriod of triangle 1
REAL(r64) :: XcmTri2 ! temporary X coord for centriod of triangle 2
REAL(r64) :: YcmTri2 ! temporary Y coord for centriod of triangle 2
REAL(r64) :: ZcmTri2 ! temporary Z coord for centriod of triangle 2
type(vector) :: VecAvg ! Average (calc for multisided polygons (>4 sides))
integer :: vert
integer :: negZcount ! for warning error in surface centroids
negZcount=0
! loop through all the surfaces
Do thisSurf=1, TotSurfaces
! IF (Surface(thisSurf)%CLASS == 'INTMASS') CYCLE
IF (Surface(thisSurf)%CLASS == SurfaceClass_IntMass) CYCLE
!re-init
Xcm = 0.0d0
Ycm = 0.0d0
Zcm = 0.0d0
SELECT CASE (surface(thissurf)%sides) !is this a 3- or 4-sided surface
CASE (3) !3-sided polygon
! centriod is simple average
Xcm = sum(surface(thisSurf)%vertex%x) / 3.0d0
Ycm = sum(surface(thisSurf)%vertex%y) / 3.0d0
Zcm = sum(surface(thisSurf)%vertex%z) / 3.0d0
CASE (4) !4-sided polygon
! re-init
Triangle1%X = 0.0d0
Triangle1%Y = 0.0d0
Triangle1%Z = 0.0d0
Triangle2%X = 0.0d0
Triangle2%Y = 0.0d0
Triangle2%Z = 0.0d0
XcmTri1 = 0.0d0
YcmTri1 = 0.0d0
ZcmTri1 = 0.0d0
XcmTri2 = 0.0d0
YcmTri2 = 0.0d0
ZcmTri2 = 0.0d0
Tri1Area = 0.0d0
Tri2Area = 0.0d0
! split into 2 3-sided polygons (Triangle 1 and Triangle 2)
Triangle1%X = surface(thisSurf)%vertex( (/1,2,3/) )%X
Triangle1%Y = surface(thisSurf)%vertex( (/1,2,3/) )%Y
Triangle1%Z = surface(thisSurf)%vertex( (/1,2,3/) )%Z
Triangle2%X = surface(thisSurf)%vertex( (/1,3,4/) )%X
Triangle2%Y = surface(thisSurf)%vertex( (/1,3,4/) )%Y
Triangle2%Z = surface(thisSurf)%vertex( (/1,3,4/) )%Z
! get area of triangles.
Tri1Area = AreaPolygon(3,Triangle1)
Tri2Area = AreaPolygon(3,Triangle2)
! get total Area of quad.
TotalArea = surface(thisSurf)%grossarea
IF (TotalArea <= 0.0d0 ) then
!catch a problem....
CALL ShowWarningError('CalcSurfaceCentroid: zero area surface, for surface='//TRIM(surface(thisSurf)%Name))
cycle
endif
! get centroid of Triangle 1
XcmTri1 = sum(Triangle1%x) / 3.0d0
YcmTri1 = sum(Triangle1%y) / 3.0d0
ZcmTri1 = sum(Triangle1%z) / 3.0d0
! get centroid of Triangle 2
XcmTri2 = sum(Triangle2%x) / 3.0d0
YcmTri2 = sum(Triangle2%y) / 3.0d0
ZcmTri2 = sum(Triangle2%z) / 3.0d0
! find area weighted combination of the two centroids.
Xcm = (XcmTri1*Tri1Area + XcmTri2*Tri2Area) / TotalArea
Ycm = (YcmTri1*Tri1Area + YcmTri2*Tri2Area) / TotalArea
Zcm = (ZcmTri1*Tri1Area + ZcmTri2*Tri2Area) / TotalArea
CASE (5:) !multi-sided polygon
! (Maybe triangulate? For now, use old "z" average method")
! and X and Y -- straight average
! X1=MINVAL(Surface(thisSurf)%Vertex(1:Surface(thisSurf)%Sides)%X)
! X2=MAXVAL(Surface(thisSurf)%Vertex(1:Surface(thisSurf)%Sides)%X)
! Y1=MINVAL(Surface(thisSurf)%Vertex(1:Surface(thisSurf)%Sides)%Y)
! Y2=MAXVAL(Surface(thisSurf)%Vertex(1:Surface(thisSurf)%Sides)%Y)
! Z1=MINVAL(Surface(thisSurf)%Vertex(1:Surface(thisSurf)%Sides)%Z)
! Z2=MAXVAL(Surface(thisSurf)%Vertex(1:Surface(thisSurf)%Sides)%Z)
! Xcm=(X1+X2)/2.0d0
! Ycm=(Y1+Y2)/2.0d0
! Zcm=(Z1+Z2)/2.0d0
! Calc centroid as average of surfaces
VecAvg=vector(0.0d0,0.0d0,0.0d0)
do vert=1,Surface(thisSurf)%Sides
VecAvg=VecAvg + Surface(thisSurf)%Vertex(vert)
enddo
VecAvg=VecAvg/REAL(Surface(thisSurf)%Sides,r64)
Xcm=VecAvg%x
Ycm=VecAvg%y
Zcm=VecAvg%z
CASE DEFAULT
IF (Surface(thisSurf)%Name /= Blank) THEN
CALL ShowWarningError('CalcSurfaceCentroid: caught problem with # of sides, for surface='//TRIM(surface(thisSurf)%Name))
CALL ShowContinueError('... number of sides must be >= 3, this surface # sides='// &
TRIM(RoundSigDigits(surface(thisSurf)%Sides)))
ELSE
CALL ShowWarningError('CalcSurfaceCentroid: caught problem with # of sides, for surface=#'// &
TRIM(RoundSigDigits(thisSurf)))
CALL ShowContinueError('...surface name is blank. Examine surfaces -- '// &
'this may be a problem with ill-formed interzone surfaces.')
CALL ShowContinueError('... number of sides must be >= 3, this surface # sides='// &
TRIM(RoundSigDigits(surface(thisSurf)%Sides)))
ENDIF
END SELECT
! store result in the surface structure in DataSurfaces
Surface(thissurf)%centroid%x = Xcm
Surface(thissurf)%centroid%y = Ycm
Surface(thissurf)%centroid%z = Zcm
if (Zcm < 0.0d0) then
if (Surface(thisSurf)%ExtWind .or. Surface(thisSurf)%ExtBoundCond == ExternalEnvironment) negZcount=negZcount+1
endif
ENDDO !loop through surfaces
if (negZcount > 0) then
CALL ShowWarningError('CalcSurfaceCentroid: '//TRIM(RoundSigDigits(negZcount))// &
' Surfaces have the Z coordinate < 0.')
CALL ShowContinueError('...in any calculations, Wind Speed will be 0.0 for these surfaces.')
CALL ShowContinueError('...in any calculations, Outside temperatures will be the '// &
'outside temperature + '//TRIM(RoundSigDigits(WeatherFileTempModCoeff,3))//' for these surfaces.')
CALL ShowContinueError('...that is, these surfaces will have conditions as though at ground level.')
endif
RETURN
END SUBROUTINE CalcSurfaceCentroid