Commit eb482f2b authored by aarograh's avatar aarograh
Browse files

Updates to Geom_Hex 

Adds hexagon_vertexes, hex_2D_to_1D, and hex_1D_to_2D.
Also fixes a few minor bugs in other routines and
adds a unit test for all current routines.
parent bc7c1988

src/Geom.f90

+3 −0
Original line number Diff line number Diff line
@@ -60,6 +60,9 @@ PUBLIC :: bound_angle 
PUBLIC :: hex_side_length_from_pitch

PUBLIC :: hex_array_element_centroid

PUBLIC :: rotate_2d_vector

PUBLIC :: hexagon_vertexes

PUBLIC :: hex_2D_to_1D

PUBLIC :: hex_1D_to_2D



INTERFACE newGeom

  !> @copybrief Geom::newGeom_line

src/Geom_Hex.f90

+118 −28
Original line number Diff line number Diff line
@@ -27,6 +27,9 @@ PUBLIC :: bound_angle 
PUBLIC :: hex_side_length_from_pitch

PUBLIC :: hex_array_element_centroid

PUBLIC :: rotate_2d_vector

PUBLIC :: hexagon_vertexes

PUBLIC :: hex_2D_to_1D

PUBLIC :: hex_1D_to_2D



!> Enumerations for hexagonal unit orientations

INTEGER(SIK),PARAMETER,PUBLIC :: HEX_POINTY_TOP=0
@@ -116,11 +119,12 @@ ENDFUNCTION hex_elements_from_rings 
!> are performed to ensure the final result is in the range.  The client code

!> is responsible for ensuring this.

!>

ELEMENTAL SUBROUTINE bound_angle(a, lower, upper, rad)

  REAL(SRK),INTENT(INOUT) :: a

ELEMENTAL FUNCTION bound_angle(a, lower, upper, rad) RESULT(angle)

  REAL(SRK),INTENT(IN) :: a

  REAL(SRK),INTENT(IN) :: lower

  REAL(SRK),INTENT(IN) :: upper

  LOGICAL(SBK),INTENT(IN),OPTIONAL :: rad

  REAL(SRK) :: angle

  !

  REAL(SRK) :: increment
@@ -128,14 +132,15 @@ ELEMENTAL SUBROUTINE bound_angle(a, lower, upper, rad) 
  IF(PRESENT(rad)) THEN

    IF(rad) increment = TWOPI

  ENDIF

  DO WHILE(a > upper)

    a = a - increment

  angle = a

  DO WHILE(angle > upper)

    angle = angle - increment

  ENDDO

  DO WHILE(a <= lower)

    a = a + increment

  DO WHILE(angle <= lower)

    angle = angle + increment

  ENDDO



ENDSUBROUTINE bound_angle

ENDFUNCTION bound_angle

!

!-------------------------------------------------------------------------------

!> @brief Calculates the hexagon side length from the pitch
@@ -146,7 +151,7 @@ FUNCTION hex_side_length_from_pitch(pitch) RESULT(length) 
  REAL(SRK),INTENT(IN) :: pitch

  REAL(SRK) :: length



  length = pitch / SQRT(3.0)

  length = pitch / SQRT(3.0_SRK)



ENDFUNCTION hex_side_length_from_pitch

!
@@ -164,7 +169,7 @@ ENDFUNCTION hex_side_length_from_pitch 
!> the origin (0.0, 0.0) is assumed to be at the center of the hexagon at the

!> center of the array.

!>

FUNCTION hex_array_element_centroid(pitch, row, col, rings, orientation) RESULT(coords)

PURE FUNCTION hex_array_element_centroid(pitch, row, col, rings, orientation) RESULT(coords)

  REAL(SRK),INTENT(IN) :: pitch

  INTEGER(SIK),INTENT(IN) :: row

  INTEGER(SIK),INTENT(IN) :: col
@@ -175,22 +180,14 @@ FUNCTION hex_array_element_centroid(pitch, row, col, rings, orientation) RESULT( 
  INTEGER(SIK) :: mid

  REAL(SRK) :: vec(2)



  REQUIRE(pitch > 0.0)

  REQUIRE(rings > 0)

  REQUIRE(row > 0)

  REQUIRE(row <= hex_rows(rings))

  REQUIRE(col > 0)

  REQUIRE(col <= hex_columns(rings, row))

  REQUIRE(ANY(orientation == HEX_ORIENTATIONS))



  mid = hex_mid(rings)



  ! Calculate x,y for array of flat-topped hexagons first

  vec(1) = (-(mid - 1) + ABS(mid - row) * 0.5 + col) * pitch

  ! Calculate x,y for array of pointy-topped hexagons first

  vec(1) = (-(mid - 1) + ABS(mid - row) * HALF + (col - 1)) * pitch

  vec(2) = (mid - row) * pitch * COS(SIXTHPI)



  ! Now rotate 30 degrees for pointy-topped

  IF(orientation == HEX_POINTY_TOP) THEN

  IF(orientation == HEX_FLAT_TOP) THEN

    coords = rotate_2d_vector(vec, SIXTHPI, rad=.TRUE.)

  ELSE

    coords = vec
@@ -208,28 +205,121 @@ ENDFUNCTION hex_array_element_centroid 
!> The @c vec and @c angle must both be degrees (@c rad=false) or both be radians

!> (@c rad=true).

!>

FUNCTION rotate_2d_vector(vec, angle, rad) RESULT(rot_vec)

PURE FUNCTION rotate_2d_vector(vec, angle, rad) RESULT(rot_vec)

  REAL(SRK),INTENT(IN) :: vec(2)

  REAL(SRK),INTENT(IN) :: angle

  LOGICAL(SBK),INTENT(IN),OPTIONAL :: rad

  REAL(SRK) :: rot_vec(2)

  !

  LOGICAL(SBK) :: lrad

  REAL(SRK) :: a, lower,upper

  REAL(SRK) :: a



  lrad = .FALSE.

  IF(PRESENT(rad)) lrad = rad

  a = angle

  lower = ZERO

  IF(lrad) THEN

    upper = TWOPI

  ELSE

    upper = 360.0

  IF(.NOT.lrad) THEN

    a = a * degrees2rad

  ENDIF

  CALL bound_angle(a, lower, upper, RAD=lrad)

  a = bound_angle(a, ZERO, TWOPI, RAD=.TRUE.)

  rot_vec(1) = COS(a)*vec(1) - SIN(a)*vec(2)

  rot_vec(2) = SIN(a)*vec(1) + COS(a)*vec(2)



ENDFUNCTION rotate_2d_vector

!

!-------------------------------------------------------------------------------

!> @brief Returns the list of points defining a hexagon

!> @param apothem the apothem of the hexagon

!> @param orientation the orientation of the hexagon

!> @param centroid the centroid of the hexagon; optional, default to (0.0, 0.0)

!>

!> The points will be ordered starting at the positive x-axis and progressing

!> clockwise.

!>

FUNCTION hexagon_vertexes(apothem, orientation, centroid) RESULT(hexagon)

  REAL(SRK),INTENT(IN) :: apothem

  INTEGER(SIK),INTENT(IN) :: orientation

  REAL(SRK),INTENT(IN),OPTIONAL :: centroid(2)

  REAL(SRK) :: hexagon(2,6)

  !

  REAL(SRK) :: side_length,c(2)



  REQUIRE(apothem > 0)

  REQUIRE(ANY(orientation == HEX_ORIENTATIONS))



  IF(PRESENT(centroid)) c = centroid



  side_length = hex_side_length_from_pitch(TWO * apothem)

  IF(orientation == HEX_FLAT_TOP) THEN

    hexagon = RESHAPE(SOURCE=[c(1) + side_length, c(2), &

                              c(1) + side_length * HALF, c(2) - apothem, &

                              c(1) - side_length * HALF, c(2) - apothem, &

                              c(1) - side_length, c(2), &

                              c(1) - side_length * HALF, c(2) + apothem, &

                              c(1) + side_length * HALF, c(2) + apothem], &

                      SHAPE=[2,6])

  ELSE

    hexagon = RESHAPE(SOURCE=[c(1) + apothem, c(2) - side_length * HALF, &

                              c(1), c(2) - side_length, &

                              c(1) - apothem, c(2) - side_length * HALF, &

                              c(1) - apothem, c(2) + side_length * HALF, &

                              c(1), c(2) + side_length, &

                              c(1) + apothem, c(2) + side_length * HALF], &

                      SHAPE=[2,6])

  ENDIF



ENDFUNCTION hexagon_vertexes

!

!-------------------------------------------------------------------------------

!> @brief Converts 2D hex array indexes to 1D index

!> @param rings the number of rings in the array

!> @param row the row index

!> @param col the column index

!> @returns index

!>

!> The row index begins with the top row, but the 1D index begins with the bottom row.

!> The column index and 1D index both move left to right down the row.

!>

ELEMENTAL FUNCTION hex_2D_to_1D(rings, row, col) RESULT(index)

  INTEGER(SIK),INTENT(IN) :: rings

  INTEGER(SIK),INTENT(IN) :: row

  INTEGER(SIK),INTENT(IN) :: col

  INTEGER(SIK) :: index

  !

  INTEGER(SIK) :: i



  index = 0

  DO i = hex_rows(rings),row+1,-1

    index = index + hex_columns(rings, i)

  ENDDO !i

  index = index + col



ENDFUNCTION hex_2D_to_1D

!

!-------------------------------------------------------------------------------

!> @brief Converts 1D hex array index to 2D indexes

!> @param rings the number of rings in the array

!> @param index the 1D hex array index

!> @returns indexes the 2D hex array indexes (col, row)

!>

!> The row index begins with the top row, but the 1D index begins with the bottom row.

!> The column index and 1D index both move left to right down the row.

!>

PURE FUNCTION hex_1D_to_2D(rings, index) RESULT(indexes)

  INTEGER(SIK),INTENT(IN) :: rings

  INTEGER(SIK),INTENT(IN) :: index

  INTEGER(SIK) :: indexes(2)

  !

  INTEGER(SIK) :: i



  indexes(1) = 1

  i = index

  DO WHILE(i > hex_columns(rings, indexes(1)))

    i = i - hex_columns(rings, indexes(1))

    indexes(1) = indexes(1) + 1

  ENDDO

  indexes(2) = i

  indexes(1) = hex_rows(rings) - indexes(1) + 1 !invert to match top->bottom row indexing



ENDFUNCTION hex_1D_to_2D

!

ENDMODULE Geom_Hex
 
 No newline at end of file

unit_tests/CMakeLists.txt

+1 −0
Original line number Diff line number Diff line
@@ -37,6 +37,7 @@ SET(UNIT_TEST_NAMES 
    testGeom_Box

    testGeom_Graph

    testGeom_Poly

    testGeom_Hex

    testGeom

    testSorting

    testMeshTransfer

unit_tests/testGeom_Hex/CMakeLists.txt

0 → 100644
+10 −0
Original line number Diff line number Diff line 
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!

#                          Futility Development Group                          !

#                             All rights reserved.                             !

#                                                                              !

# Futility is a jointly-maintained, open-source project between the University !

# of Michigan and Oak Ridge National Laboratory.  The copyright and license    !

# can be found in LICENSE.txt in the head directory of this repository.        !

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!

INCLUDE(Futility_CreateUnitTest)

Futility_CreateUnitTest(testGeom_Hex)

unit_tests/testGeom_Hex/testGeom_Hex.f90

0 → 100644
+420 −0

File added.

Preview size limit exceeded, changes collapsed.