Skip to content

Commit

Permalink
Adding interior_combine() to Fortran implementation.
Browse files Browse the repository at this point in the history
  • Loading branch information
@dhermes
dhermes committed Nov 21, 2017
1 parent aeb1853 commit 95c8ac3
Show file tree
Hide file tree
Showing 2 changed files with 247 additions and 7 deletions.
128 changes: 125 additions & 3 deletions src/bezier/surface_intersection.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,11 +25,11 @@ module surface_intersection
evaluate_barycentric, jacobian_both, subdivide_nodes, compute_edge_nodes
implicit none
private &
LocateCandidate, MAX_LOCATE_SUBDIVISIONS, LOCATE_EPS, &
LocateCandidate, MAX_LOCATE_SUBDIVISIONS, LOCATE_EPS, MAX_EDGES, &
newton_refine_solve, split_candidate, allocate_candidates, &
update_candidates, ignored_edge_corner, ignored_double_corner, &
ignored_corner, classify_tangent_intersection, no_intersections, &
remove_node
remove_node, finalize_segment
public &
Intersection, CurvedPolygonSegment, &
IntersectionClassification_FIRST, IntersectionClassification_SECOND, &
Expand All @@ -40,7 +40,7 @@ module surface_intersection
SurfaceContained_FIRST, SurfaceContained_SECOND, newton_refine, &
locate_point, classify_intersection, add_st_vals, &
surfaces_intersection_points, get_next, to_front, add_segment, &
surfaces_intersect
interior_combine, surfaces_intersect

! NOTE: This (for now) is not meant to be C-interoperable.
type :: Intersection
Expand Down Expand Up @@ -68,6 +68,7 @@ module surface_intersection
! NOTE: These values are also defined in equivalent Python source.
integer(c_int), parameter :: MAX_LOCATE_SUBDIVISIONS = 20
real(c_double), parameter :: LOCATE_EPS = 0.5_dp**47
integer(c_int), parameter :: MAX_EDGES = 10
! Values of IntersectionClassification enum:
integer(c_int), parameter :: IntersectionClassification_FIRST = 0
integer(c_int), parameter :: IntersectionClassification_SECOND = 1
Expand Down Expand Up @@ -1140,6 +1141,127 @@ subroutine add_segment(curr_node, next_node, index, segments)

end subroutine add_segment

subroutine finalize_segment(segment_ends, index, count)

integer(c_int), allocatable, intent(inout) :: segment_ends(:)
integer(c_int), intent(inout) :: index
integer(c_int), intent(in) :: count
! Variables outside of signature.
integer(c_int) :: curr_size
integer(c_int), allocatable :: segment_ends_swap(:)

if (allocated(segment_ends)) then
curr_size = size(segment_ends)
if (curr_size < index) then
allocate(segment_ends_swap(index))
segment_ends_swap(:curr_size) = segment_ends(:curr_size)
call move_alloc(segment_ends_swap, segment_ends)
end if
else
allocate(segment_ends(index))
end if

segment_ends(index) = count
index = index + 1

end subroutine finalize_segment

subroutine interior_combine( &
num_intersections, intersections, &
segment_ends, segments, status)

! NOTE: This subroutine is not meant to be part of the interface for this
! module, but it is (for now) public, so that it can be tested.
! NOTE: This assumes, but does not check, that ``num_intersections > 0``.
! NOTE: If ``segments`` is already allocated / with data, the data will be
! overwritten in place.

! Possible error states:
! * Status_SUCCESS: On success.
! * Status_UNKNOWN: If a curved polygon requires more than ``MAX_EDGES``
! sides. (This could be due to either a particular
! complex intersection or a programming error which
! causes ``at_start`` to never be true.)

integer(c_int), intent(in) :: num_intersections
type(Intersection), intent(in) :: intersections(num_intersections)
integer(c_int), allocatable, intent(out) :: segment_ends(:)
type(CurvedPolygonSegment), allocatable, intent(inout) :: segments(:)
integer(c_int), intent(out) :: status
! Variables outside of signature.
integer(c_int) :: unused(num_intersections)
integer(c_int) :: segment_index, segment_ends_index
integer(c_int) :: remaining, i, start
type(Intersection) :: curr_node, next_node
logical(c_bool) :: at_start

! In this case, we know we'll have at least ``num_intersections`` segments
! (and may have more due to corners).
if (allocated(segments)) then
if (size(segments) < num_intersections) then
deallocate(segments)
allocate(segments(num_intersections))
end if
else
allocate(segments(num_intersections))
end if

status = Status_SUCCESS

! Set all of the unused indices.
unused = [ (i, i = 1, num_intersections) ]
remaining = num_intersections
segment_index = 1
segment_ends_index = 1
do while (remaining > 0)
! "Pop" off intersection from the end to start a curved polygon.
start = unused(remaining)
remaining = remaining - 1

curr_node = intersections(start)
at_start = .FALSE.
! Now move around the edge segments until the curved polygon's
! entire boundary has been traversed.
edge_loop: do i = 1, MAX_EDGES
! Follow along the current edge to the next node.
call get_next( &
num_intersections, intersections, unused, remaining, &
start, curr_node, next_node, at_start)
call add_segment(curr_node, next_node, segment_index, segments)
! After adding the segment, check if we are back where we started.
if (at_start) then
exit edge_loop
end if

! Now the ``next_node`` becomes the current node, but we may
! need to "rotate" it to the next edge if on a corner.
call to_front( &
num_intersections, intersections, unused, remaining, &
start, next_node, curr_node, at_start)
! After "rotating" to the next edge, check if we are back where
! we started.
if (at_start) then
exit edge_loop
end if
end do edge_loop

if (at_start) then
! NOTE: We pass ``segment_index - 1`` because it has been incremented
! one past the last written index.
call finalize_segment( &
segment_ends, segment_ends_index, segment_index - 1)
else
! If the loop terminated without reaching the start node, then
! we have encountered an error.
! LCOV_EXCL_START
status = Status_UNKNOWN
return
! LCOV_EXCL_STOP
end if
end do

end subroutine interior_combine

subroutine surfaces_intersect( &
num_nodes1, nodes1, degree1, &
num_nodes2, nodes2, degree2, &
Expand Down
126 changes: 122 additions & 4 deletions tests/fortran/test_surface_intersection.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,15 +27,16 @@ module test_surface_intersection
SurfaceContained_FIRST, SurfaceContained_SECOND, newton_refine, &
locate_point, classify_intersection, add_st_vals, &
surfaces_intersection_points, get_next, to_front, add_segment, &
surfaces_intersect
interior_combine, surfaces_intersect
use types, only: dp
use unit_test_helpers, only: print_status
implicit none
private &
test_newton_refine, test_locate_point, test_classify_intersection, &
test_add_st_vals, test_surfaces_intersection_points, &
intersection_check, intersection_equal, test_get_next, test_to_front, &
test_add_segment, test_surfaces_intersect
test_add_segment, test_interior_combine, segment_check, &
test_surfaces_intersect
public surface_intersection_all_tests

contains
Expand All @@ -51,6 +52,7 @@ subroutine surface_intersection_all_tests(success)
call test_get_next(success)
call test_to_front(success)
call test_add_segment(success)
call test_interior_combine(success)
call test_surfaces_intersect(success)

end subroutine surface_intersection_all_tests
Expand Down Expand Up @@ -839,7 +841,7 @@ function intersection_check( &
intersection_, s, t, index_first, index_second, &
interior_curve) result(same)

type(intersection), intent(in) :: intersection_
type(Intersection), intent(in) :: intersection_
real(c_double), intent(in) :: s, t
integer(c_int), intent(in) :: index_first, index_second
integer(c_int), intent(in) :: interior_curve
Expand All @@ -856,7 +858,7 @@ end function intersection_check

function intersection_equal(intersection1, intersection2) result(same)

type(intersection), intent(in) :: intersection1, intersection2
type(Intersection), intent(in) :: intersection1, intersection2
logical(c_bool) :: same

same = intersection_check( &
Expand Down Expand Up @@ -1179,6 +1181,122 @@ subroutine test_add_segment(success)

end subroutine test_add_segment

subroutine test_interior_combine(success)
logical(c_bool), intent(inout) :: success
! Variables outside of signature.
logical :: case_success
integer :: case_id
character(16) :: name
integer(c_int) :: first, second
type(Intersection) :: intersections(4)
integer(c_int), allocatable :: segment_ends(:)
type(CurvedPolygonSegment), allocatable :: segments(:)
integer(c_int) :: status

case_id = 1
name = "interior_combine"

first = IntersectionClassification_FIRST
second = IntersectionClassification_SECOND

! CASE 1: ``segments`` and ``segment_ends`` are not allocated,
! from 1L-6L (ID: 23).
intersections(1) = Intersection(0.25_dp, 0.75_dp, 1, 3, second)
intersections(2) = Intersection(0.75_dp, 0.25_dp, 3, 1, first)

case_success = ( &
.NOT. allocated(segments) .AND. &
.NOT. allocated(segment_ends))
call interior_combine( &
2, intersections(:2), segment_ends, segments, status)
case_success = ( &
case_success .AND. &
allocated(segment_ends) .AND. &
all(segment_ends == [4]) .AND. &
allocated(segments) .AND. &
size(segments) == 4 .AND. &
segment_check(segments(1), 0.75_dp, 1.0_dp, 3) .AND. &
segment_check(segments(2), 0.0_dp, 0.25_dp, 1) .AND. &
segment_check(segments(3), 0.75_dp, 1.0_dp, 6) .AND. &
segment_check(segments(4), 0.0_dp, 0.25_dp, 4) .AND. &
status == Status_SUCCESS)
call print_status(name, case_id, case_success, success)

! CASE 2: ``segments`` and ``segment_ends`` must be re-allocated;
! two disjoint intersections, from 13L-38Q (ID: 39).
intersections(1) = Intersection(0.625_dp, 0.25_dp, 1, 1, first)
intersections(2) = Intersection(0.375_dp, 0.75_dp, 1, 1, second)
intersections(3) = Intersection(0.3125_dp, 0.25_dp, 1, 2, first)
intersections(4) = Intersection(0.6875_dp, 0.75_dp, 1, 3, second)

case_success = ( &
allocated(segment_ends) .AND. &
size(segment_ends) == 1 .AND. &
allocated(segments) .AND. &
size(segments) == 4)
call interior_combine( &
4, intersections, segment_ends, segments, status)
case_success = ( &
case_success .AND. &
allocated(segment_ends) .AND. &
all(segment_ends == [3, 6]) .AND. &
allocated(segments) .AND. &
size(segments) == 6 .AND. &
segment_check(segments(1), 0.75_dp, 1.0_dp, 6) .AND. &
segment_check(segments(2), 0.0_dp, 0.25_dp, 4) .AND. &
segment_check(segments(3), 0.625_dp, 0.6875_dp, 1) .AND. &
segment_check(segments(4), 0.3125_dp, 0.375_dp, 1) .AND. &
segment_check(segments(5), 0.75_dp, 1.0_dp, 4) .AND. &
segment_check(segments(6), 0.0_dp, 0.25_dp, 5) .AND. &
status == Status_SUCCESS)
call print_status(name, case_id, case_success, success)

! CASE 3: ``segments`` has been allocated but is not large enough;
! intersection terminates after node is "rotated" via
! ``to_front()``, from 1L-28Q (ID: 25) with rounded parameters
! from ``s = 0.40587727318528016`` and
! ``t == 0.59412272681471978``.
intersections(1) = Intersection(0.375_dp, 0.5625_dp, 2, 2, second)
intersections(2) = Intersection(0.1875_dp, 0.0_dp, 1, 3, first)
! Note that the corner is the last intersection, which will be used
! as ``start``.
deallocate(segments)
allocate(segments(1)) ! Too small.
case_success = ( &
allocated(segment_ends) .AND. &
size(segment_ends) == 2 .AND. &
allocated(segments) .AND. &
size(segments) == 1)
call interior_combine( &
2, intersections(:2), segment_ends, segments, status)
case_success = ( &
case_success .AND. &
allocated(segment_ends) .AND. &
all(segment_ends == [3]) .AND. &
allocated(segments) .AND. &
size(segments) == 3 .AND. &
segment_check(segments(1), 0.1875_dp, 1.0_dp, 1) .AND. &
segment_check(segments(2), 0.0_dp, 0.375_dp, 2) .AND. &
segment_check(segments(3), 0.5625_dp, 1.0_dp, 5) .AND. &
status == Status_SUCCESS)
call print_status(name, case_id, case_success, success)

end subroutine test_interior_combine

function segment_check(segment, start, end_, edge_index) result(same)

type(CurvedPolygonSegment), intent(in) :: segment
real(c_double), intent(in) :: start, end_
integer(c_int), intent(in) :: edge_index
logical(c_bool) :: same

same = ( &
segment%start == start .AND. &
segment%end_ == end_ .AND. &
segment%edge_index == edge_index)

end function segment_check

subroutine test_surfaces_intersect(success)
logical(c_bool), intent(inout) :: success
! Variables outside of signature.
Expand Down

0 comments on commit 95c8ac3

Please sign in to comment.