Extensions to Graph, Line, and Point classes; bug fix in hex_switch_orientation (ec997c89) · Commits · Futility / Futility

src/Geom_Graph.f90

+44 −30

Original line number	Diff line number	Diff line
		@@ -539,12 +539,14 @@ ENDFUNCTION isMinimumCycle_graphType
		!> @brief Inserts a new vertex into a graph object
		!> @param this the graph to modify
		!> @param point the point containing coordinates of the new vertex to add
		!> @param idx the index the point was inserted at; optional
		!>
		SUBROUTINE insertVertex_point(this,point)
		SUBROUTINE insertVertex_point(this,point,idx)
		CLASS(GraphType),INTENT(INOUT) :: this
		TYPE(PointType),INTENT(IN) :: point
		INTEGER(SIK),INTENT(OUT),OPTIONAL :: idx

		CALL this%insertVertex(point%coord)
		CALL this%insertVertex(point%coord,idx)

		ENDSUBROUTINE insertVertex_point
		!
		@@ -552,15 +554,19 @@ ENDSUBROUTINE insertVertex_point
		!> @brief Inserts an array of new vertexex into a graph object
		!> @param this the graph to modify
		!> @param coord the coordinates of the new vertex to add
		!> @param idx the indexes the points were inserted at; optional
		!>
		SUBROUTINE insertVertex_coords_array(this,coords)
		SUBROUTINE insertVertex_coords_array(this,coords,idx)
		CLASS(GraphType),INTENT(INOUT) :: this
		REAL(SRK),INTENT(IN) :: coords(:,:)
		INTEGER(SIK),ALLOCATABLE,INTENT(OUT),OPTIONAL :: idx(:)
		!
		INTEGER(SIK) :: i
		INTEGER(SIK) :: i,ind

		IF(PRESENT(idx)) ALLOCATE(idx(SIZE(coords,DIM=2)))
		DO i = 1,SIZE(coords,DIM=2)
		CALL this%insertVertex_coords(coords(:,i))
		CALL this%insertVertex_coords(coords(:,i),ind)
		IF(PRESENT(idx)) idx(i) = ind
		ENDDO !i

		ENDSUBROUTINE insertVertex_coords_array
		@@ -569,13 +575,17 @@ ENDSUBROUTINE insertVertex_coords_array
		!> @brief Inserts a new vertex into a graph object
		!> @param this the graph to modify
		!> @param coord the coordinates of the new vertex to add
		!> @param idx the index the point was inserted at; optional
		!>
		SUBROUTINE insertVertex_coords(this,coord)
		SUBROUTINE insertVertex_coords(this,coord,idx)
		CLASS(GraphType),INTENT(INOUT) :: this
		REAL(SRK),INTENT(IN) :: coord(2)
		INTEGER(SIK) :: i,j,n,k,idx
		INTEGER(SIK),INTENT(OUT),OPTIONAL :: idx
		!
		INTEGER(SIK) :: i,j,n,k,ind
		INTEGER(SIK),ALLOCATABLE :: tmpE(:,:)
		REAL(SRK),ALLOCATABLE :: tmpVertices(:,:),tmpQE(:,:,:)

		IF(ALLOCATED(this%vertices)) THEN
		n=SIZE(this%vertices,DIM=2)
		CALL dmallocA(tmpVertices,2,n+1)
		@@ -585,10 +595,10 @@ ENDSUBROUTINE insertVertex_coords_array
		IF(coord(1) .APPROXLE. this%vertices(1,i)) THEN
		IF(coord(1) .APPROXEQA. this%vertices(1,i)) THEN
		IF(coord(2) .APPROXEQA. this%vertices(2,i)) THEN
		idx=-1 !Duplicate vertex
		ind=-1 !Duplicate vertex
		EXIT
		ELSEIF(coord(2) < this%vertices(2,i)) THEN
		idx=i !Before i
		ind=i !Before i
		ELSE
		!After i
		DO j=i+1,n
		@@ -600,54 +610,54 @@ ENDSUBROUTINE insertVertex_coords_array
		!Search on y through sequence of same x
		DO k=i+1,j
		IF(coord(2) .APPROXEQA. this%vertices(2,k)) THEN
		idx=-1
		ind=-1
		EXIT
		ELSEIF(coord(2) < this%vertices(2,k)) THEN
		idx=k
		ind=k
		EXIT
		ENDIF
		ENDDO
		IF(k == j+1) idx=k
		IF(k == j+1) ind=k
		ENDIF
		ELSE
		idx=i !Before i
		ind=i !Before i
		ENDIF
		EXIT
		ENDIF
		ENDDO
		IF(j /= 0) i=j
		IF(i == n+1) idx=n+1 !Last point
		IF(idx > 0) THEN
		IF(idx > 1) tmpVertices(:,1:idx-1)=this%vertices(:,1:idx-1)
		tmpVertices(:,idx)=coord
		IF(idx <= n) tmpVertices(:,idx+1:n+1)=this%vertices(:,idx:n)
		IF(i == n+1) ind=n+1 !Last point
		IF(ind > 0) THEN
		IF(ind > 1) tmpVertices(:,1:ind-1)=this%vertices(:,1:ind-1)
		tmpVertices(:,ind)=coord
		IF(ind <= n) tmpVertices(:,ind+1:n+1)=this%vertices(:,ind:n)
		CALL demallocA(this%vertices)
		CALL MOVE_ALLOC(tmpVertices,this%vertices)
		!Expand Edge Matrices
		CALL dmallocA(tmpE,n+1,n+1)
		CALL dmallocA(tmpQE,3,n+1,n+1)
		DO j=1,idx-1
		DO i=1,idx-1
		DO j=1,ind-1
		DO i=1,ind-1
		tmpE(i,j)=this%edgeMatrix(i,j)
		tmpE(j,i)=this%edgeMatrix(j,i)
		tmpQE(:,i,j)=this%quadEdges(:,i,j)
		tmpQE(:,j,i)=this%quadEdges(:,j,i)
		ENDDO
		DO i=idx+1,n+1
		DO i=ind+1,n+1
		tmpE(i,j)=this%edgeMatrix(i-1,j)
		tmpE(j,i)=this%edgeMatrix(j,i-1)
		tmpQE(:,i,j)=this%quadEdges(:,i-1,j)
		tmpQE(:,j,i)=this%quadEdges(:,j,i-1)
		ENDDO
		ENDDO
		DO j=idx+1,n+1
		DO i=1,idx-1
		DO j=ind+1,n+1
		DO i=1,ind-1
		tmpE(i,j)=this%edgeMatrix(i,j-1)
		tmpE(j,i)=this%edgeMatrix(j-1,i)
		tmpQE(:,i,j)=this%quadEdges(:,i,j-1)
		tmpQE(:,j,i)=this%quadEdges(:,j-1,i)
		ENDDO
		DO i=idx+1,n+1
		DO i=ind+1,n+1
		tmpE(i,j)=this%edgeMatrix(i-1,j-1)
		tmpE(j,i)=this%edgeMatrix(j-1,i-1)
		tmpQE(:,i,j)=this%quadEdges(:,i-1,j-1)
		@@ -660,11 +670,15 @@ ENDSUBROUTINE insertVertex_coords_array
		CALL MOVE_ALLOC(tmpQE,this%quadEdges)
		ENDIF
		ELSE
		ind = 1
		CALL dmallocA(this%vertices,2,1)
		this%vertices(:,1)=coord
		CALL dmallocA(this%edgeMatrix,1,1)
		CALL dmallocA(this%quadEdges,3,1,1)
		ENDIF

		IF(PRESENT(idx)) idx = ind

		ENDSUBROUTINE insertVertex_coords
		!
		!-------------------------------------------------------------------------------

src/Geom_Hex.f90

+1 −1

Original line number	Diff line number	Diff line
		@@ -333,7 +333,7 @@ FUNCTION hex_switch_orientation(orientation) RESULT(reversed)
		INTEGER(SIK) :: reversed

		REQUIRE(ANY(orientation == HEX_ORIENTATIONS))
		reversed = 2_SIK - orientation
		reversed = 1_SIK - orientation

		ENDFUNCTION hex_switch_orientation
		!

src/Geom_Line.f90

+104 −0

Original line number	Diff line number	Diff line
		@@ -15,6 +15,8 @@
		!> line, its midpoint, or an intersection with another line.
		!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
		MODULE Geom_Line
		#include "Futility_DBC.h"
		USE Futility_DBC
		USE IntrType
		USE Geom_Points

		@@ -43,6 +45,9 @@ TYPE :: LineType
		!> @copybrief GeomLine::clear_LineType
		!> @copydetails GeomLine::clear_LineType
		PROCEDURE,PASS :: clear => clear_LineType
		!> @copybrief GeomLine::isSet
		!> @copydetails GeomLine::isSet
		PROCEDURE,PASS :: isSet
		!> @copybrief GeomLine::dim_LineType
		!> @copydetails GeomLine::dim_LineType
		PROCEDURE,PASS :: getDim => dim_LineType
		@@ -68,6 +73,18 @@ TYPE :: LineType
		!> @copybrief GeomLine::pointIsRight_LineType
		!> @copydetails GeomLine::pointIsRight_LineType
		PROCEDURE,PASS :: pointIsRight => pointIsRight_LineType
		!> @copybrief GeomLine::slope
		!> @copydetails GeomLine::slope
		PROCEDURE,PASS :: slope
		!> @copybrief GeomLine::isVertical
		!> @copydetails GeomLine::isVertical
		PROCEDURE,PASS :: isVertical
		!> @copybrief GeomLine::isHorizontal
		!> @copydetails GeomLine::isHorizontal
		PROCEDURE,PASS :: isHorizontal
		!> @copybrief GeomLine::isPerpendicular
		!> @copydetails GeomLine::isPerpendicular
		PROCEDURE,PASS :: isPerpendicular
		ENDTYPE LineType

		!> @brief Generic interface for 'is equal to' operator (==)
		@@ -118,6 +135,19 @@ ELEMENTAL SUBROUTINE clear_LineType(line)
		ENDSUBROUTINE clear_LineType
		!
		!-------------------------------------------------------------------------------
		!> @brief Indicates if line has been set or not
		!> @param this the line to query
		!> @returns bool logical indicating if @c this has been set (true) or not (false)
		!>
		ELEMENTAL FUNCTION isSet(this) RESULT(bool)
		CLASS(LineType),INTENT(IN) :: this
		LOGICAL(SBK) :: bool

		bool = (this%p1%dim > 0 .AND. this%p2%dim > 0)

		ENDFUNCTION isSet
		!
		!-------------------------------------------------------------------------------
		!> @brief Returns the dimension of the LineType @e line
		!> @param line line segment of type @c LineType
		!> @returns @c n the dimension of LineType @c line
		@@ -543,4 +573,78 @@ ELEMENTAL FUNCTION approxequal_LineType(l0,l1) RESULT(bool)
		IF(l0%p1 .APPROXEQA. l1%p1) bool=l0%p2 .APPROXEQA. l1%p2
		ENDFUNCTION approxequal_LineType
		!
		!-------------------------------------------------------------------------------
		!> @brief Calculates the slope of a line
		!> @param this the line to query
		!> @returns m the resulting slope
		!>
		!> The client is reponsible for catching infinite slope for verticle lines.
		!>
		IMPURE ELEMENTAL FUNCTION slope(this) RESULT(m)
		CLASS(LineType),INTENT(IN) :: this
		REAL(SRK) :: m

		REQUIRE(this%isSet())

		m = (this%p2%coord(2) - this%p1%coord(2)) / (this%p2%coord(1) - this%p1%coord(1))

		ENDFUNCTION slope
		!
		!-------------------------------------------------------------------------------
		!> @brief Determines if a line is vertical
		!> @param this the line to query
		!> @returns vert the resulting slope
		!>
		IMPURE ELEMENTAL FUNCTION isVertical(this) RESULT(vert)
		CLASS(LineType),INTENT(IN) :: this
		LOGICAL(SBK) :: vert

		REQUIRE(this%isSet())

		vert = (this%p1%coord(1) .APPROXEQA. this%p2%coord(1))

		ENDFUNCTION isVertical
		!
		!-------------------------------------------------------------------------------
		!> @brief Determines if a line is horizontal
		!> @param this the line to query
		!> @returns horiz the resulting slope
		!>
		IMPURE ELEMENTAL FUNCTION isHorizontal(this) RESULT(horiz)
		CLASS(LineType),INTENT(IN) :: this
		LOGICAL(SBK) :: horiz

		REQUIRE(this%isSet())

		horiz = (this%p1%coord(2) .APPROXEQA. this%p2%coord(2))

		ENDFUNCTION isHorizontal
		!
		!-------------------------------------------------------------------------------
		!> @brief Determines whether a line is perpendicular to line1 line
		!> @param line1 the line
		!> @param line2 the other line to test
		!> @returns perp logical indicating if @c line2 is perpendicular
		!> to @c line1 (true) or not (false)
		!>
		IMPURE ELEMENTAL FUNCTION isPerpendicular(line1,line2) RESULT(perp)
		CLASS(LineType),INTENT(IN) :: line1
		TYPE(LineType),INTENT(IN) :: line2
		LOGICAL(SBK) :: perp
		!
		LOGICAL(SBK) :: l1vert,l2vert

		REQUIRE(line1%isSet())
		REQUIRE(line2%isSet())

		IF(line1%isVertical()) THEN
		perp = line2%isHorizontal()
		ELSEIF(line2%isVertical()) THEN
		perp = line1%isHorizontal()
		ELSE
		perp = (line1%slope() .APPROXEQA. -1.0_SRK/line2%slope())
		ENDIF

		ENDFUNCTION isPerpendicular
		!
		ENDMODULE Geom_Line

src/Geom_Points.f90

+18 −2

Original line number	Diff line number	Diff line
		@@ -54,6 +54,9 @@ TYPE :: PointType
		!> @copybrief GeomPoints::clear_PointType
		!> @copydetails GeomPoints::clear_PointType
		PROCEDURE,PASS :: clear => clear_PointType
		!> @copybrief GeomLine::isSet
		!> @copydetails GeomLine::isSet
		PROCEDURE,PASS :: isSet
		!> @copybrief GeomPoints::RotateQtrClockwise_PointType
		!> @copydetails GeomPoints::RotateQtrClockwise_PointType
		PROCEDURE,PASS :: RotateQtrClockwise => RotateQtrClockwise_PointType
		@@ -261,6 +264,19 @@ ELEMENTAL SUBROUTINE clear_PointType(p)
		ENDSUBROUTINE clear_PointType
		!
		!-------------------------------------------------------------------------------
		!> @brief Indicates if point has been set or not
		!> @param this the point to query
		!> @returns bool logical indicating if @c this has been set (true) or not (false)
		!>
		ELEMENTAL FUNCTION isSet(this) RESULT(bool)
		CLASS(PointType),INTENT(IN) :: this
		LOGICAL(SBK) :: bool

		bool = (this%dim > 0)

		ENDFUNCTION isSet
		!
		!-------------------------------------------------------------------------------
		!> @brief Routine rotates the x-y coordinates of a point type variable @nrot number
		!> of clockwise quarter rotations.
		!> @param p the point type variable to rotate

Admin message