From 3bb3aa94d3d445abf0abb37f9b576b48da17bb98 Mon Sep 17 00:00:00 2001 From: Skylar A Gering Date: Thu, 13 Jun 2024 12:55:06 -0400 Subject: [PATCH] Swap clipping from LG to GO (#92) * Swap cut polygons * Fix removing cut poly * Swap angle calculations * Remove scaling polygons function * Swap area function * Remove calc_point_poly_dist * Swap out points in poly * Remove sort regions function * Remove intersect_coords * Rewrite moment of inertia calculations * Finish separate_xy removal * Switch basic calls to use get_polygons * Add Extents * Remove LG from grid output writer * Debugging gridded data difference * Remove LG from output * Update poly_to_floes * Finish removing LG clipping calls from all but floe utils * Remove basic multipolygon functions * Remove remaining LibGEOS intersection calls * Remove bounding box utils * Fix multipolygon clipping calls * Abstract LG.Polygon calls out * Abstract rest of LG functions out * Update ridge raft test * Remove perturbation * Small changes for testing * Remove perturbation * Readd perturb vec for random seed * Remove unneeded print statements * Remove LibGEOS from Subzero usage * Update and run simulation checks * Update Makie compat * Small fixes --- Project.toml | 13 +- examples/moving_bounds.jl | 33 +- examples/shear_flow.jl | 2 +- examples/simple_strait.jl | 5 +- examples/test_run.jl | 2 +- src/Subzero.jl | 16 +- src/eulerian_OA_code.jl | 2 +- src/floe_utils.jl | 789 ++---------------- src/output.jl | 42 +- src/physical_processes/collisions.jl | 98 ++- src/physical_processes/coupling.jl | 25 +- src/physical_processes/fractures.jl | 55 +- src/physical_processes/ridge_raft.jl | 28 +- src/physical_processes/simplification.jl | 22 +- src/physical_processes/update_floe.jl | 13 +- src/physical_processes/welding.jl | 5 +- src/simulation_components/domain_and_grid.jl | 38 +- src/simulation_components/floe.jl | 206 +++-- src/tools/analyze_floe.jl | 4 +- src/tools/file_convert.jl | 5 +- test/runtests.jl | 8 +- test/test_floe.jl | 116 ++- test/test_floe_utils.jl | 335 +------- test/test_model.jl | 4 +- test/test_physical_processes/test_coupling.jl | 18 +- .../test_physical_processes/test_fractures.jl | 39 +- .../test_ridge_raft.jl | 49 +- .../test_simplification.jl | 8 +- .../test_update_floe.jl | 12 +- 29 files changed, 497 insertions(+), 1495 deletions(-) diff --git a/Project.toml b/Project.toml index 8afc3b0..815fbbb 100644 --- a/Project.toml +++ b/Project.toml @@ -8,17 +8,17 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" +GeometryOps = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" Measures = "442fdcdd-2543-5da2-b0f3-8c86c306513e" NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9" -PolygonInbounds = "e4521ec6-8c1d-418e-9da2-b3bc4ae105d6" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66" @@ -29,16 +29,17 @@ VoronoiCells = "e3e34ffb-84e9-5012-9490-92c94d0c60a4" [compat] BenchmarkTools = "1" +CairoMakie = "0.12" DataStructures = "0.18" +Extents = "0.1" GeometryBasics = "0.4" +GeometryOps = "0.1.8" Interpolations = "0.14, 0.15" JLD2 = "0.4" -LibGEOS = "0.8, 0.9" -Makie = "0.19, 0.20" +Makie = "0.19, 0.20, 0.21" Measures = "0.3" NCDatasets = "0.12, 0.13, 0.14" -NetCDF = "0.11" -PolygonInbounds = "0.2" +NetCDF = "0.11, 0.12" SplitApplyCombine = "1" StaticArrays = "1" Statistics = "1" diff --git a/examples/moving_bounds.jl b/examples/moving_bounds.jl index 3843bd5..4a81e6d 100644 --- a/examples/moving_bounds.jl +++ b/examples/moving_bounds.jl @@ -20,8 +20,8 @@ ocean = Ocean(grid, 0.0, 0.0, 0.0) atmos = Atmos(grid, 0.0, 0.0, -1.0) # Domain creation -nboundary = MovingBoundary(North, grid, 0.2, 0.0) -sboundary = MovingBoundary(South, grid, 0.2, 0.0) +nboundary = MovingBoundary(North, grid, 0.0, -0.1) +sboundary = MovingBoundary(South, grid, 0.0, 0.1) eboundary = PeriodicBoundary(East, grid) wboundary = PeriodicBoundary(West, grid) @@ -30,7 +30,7 @@ domain = Domain(nboundary, sboundary, eboundary, wboundary) # Floe creation floe_arr = initialize_floe_field( FT, - 75, + 500, [1.0], domain, hmean, @@ -49,7 +49,7 @@ consts = Constants(E = modulus, Cd_io = 0.0, f = 0.0, turnθ = 0.0) # Run simulation run_time!(simulation) = @time run!(simulation) -dir = "output/packed_compression" +dir = "output/moving_bounds" # Output setup initwriter = InitialStateOutputWriter(dir = dir, overwrite = true) @@ -57,24 +57,18 @@ floewriter = FloeOutputWriter(40, dir = dir, overwrite = true) writers = OutputWriters(initwriter, floewriter) # Simulation settings -# ridgeraft_settings = RidgeRaftSettings( -# ridge_raft_on = true, -# Δt = 150, -# domain_gain_probability = 0.5 -# ) +ridgeraft_settings = RidgeRaftSettings( + ridge_raft_on = true, + Δt = 150, + domain_gain_probability = 0.5 +) weld_settings = WeldSettings( weld_on = true, Δts = [150, 300, 600], Nxs = [2, 1, 1], Nys = [2, 2, 1], ) -weld_on::Bool = false - Δts::Vector{Int} = Vector{Int}() - Nxs::Vector{Int} = Vector{Int}() - Nys::Vector{Int} = Vector{Int}() - min_weld_area::FT = 1e6 - max_weld_area::FT = 2e10 - welding_coeff::FT = 150 + coupling_settings = CouplingSettings(two_way_coupling_on = true) @@ -87,12 +81,13 @@ simulation = Simulation( writers = writers, rng = Xoshiro(1), coupling_settings = coupling_settings, + ridgeraft_settings = ridgeraft_settings, ) run_time!(simulation) Subzero.plot_sim( - "output/packed_compression/floes.jld2", - "output/packed_compression/initial_state.jld2", + joinpath(dir, "floes.jld2"), + joinpath(dir, "initial_state.jld2"), 20, - "output/packed_compression/packed_compression.mp4", + joinpath(dir, "moving_bounds.mp4"), ) \ No newline at end of file diff --git a/examples/shear_flow.jl b/examples/shear_flow.jl index d9ece28..121a5b0 100644 --- a/examples/shear_flow.jl +++ b/examples/shear_flow.jl @@ -72,7 +72,7 @@ simulation = Simulation( model = model, consts = consts, Δt = Δt, - nΔt = 10000, + nΔt = 2000, verbose = true, writers = writers, rng = Xoshiro(1), diff --git a/examples/simple_strait.jl b/examples/simple_strait.jl index c6674e8..cb0e0be 100644 --- a/examples/simple_strait.jl +++ b/examples/simple_strait.jl @@ -46,7 +46,7 @@ floe_settings = FloeSettings( # Floe creation floe_arr = initialize_floe_field( FT, - 75, + 500, [0.7], domain, hmean, @@ -80,7 +80,8 @@ run_time!(simulation) = @time run!(simulation) initwriter = InitialStateOutputWriter(dir = dir, overwrite = true) floewriter = FloeOutputWriter(50, dir = dir, overwrite = true) -writers = OutputWriters(initwriter, floewriter) +gridwriter = GridOutputWriter(100, grid, (10, 10), dir = dir, overwrite = true) +writers = OutputWriters(initwriter, floewriter, gridwriter) simulation = Simulation( model = model, diff --git a/examples/test_run.jl b/examples/test_run.jl index 1c88d97..8494881 100644 --- a/examples/test_run.jl +++ b/examples/test_run.jl @@ -33,7 +33,7 @@ end function update_height(floes, i, new_height, consts) floes.height[i] = new_height floes.mass[i] = floes.area[i] * floes.height[i] * consts.ρi - floes.moment[i] = Subzero.calc_moment_inertia( + floes.moment[i] = Subzero._calc_moment_inertia( floes.coords[i], floes.centroid[i], floes.height[i], diff --git a/src/Subzero.jl b/src/Subzero.jl index cf377e4..ceaffbd 100644 --- a/src/Subzero.jl +++ b/src/Subzero.jl @@ -69,11 +69,11 @@ export plot_sim import Base.@kwdef # this is being exported as of version 1.9 -import LibGEOS as LG -using CairoMakie, DataStructures, Dates, GeometryBasics, Interpolations, JLD2, - LinearAlgebra, Logging, Makie, Measures, NCDatasets, NetCDF, - PolygonInbounds, Printf, Random, SplitApplyCombine, StaticArrays, - Statistics, StructArrays, VoronoiCells +import GeometryOps as GO +import GeometryOps.GeoInterface as GI +using CairoMakie, DataStructures, Dates, Extents, GeometryBasics, Interpolations, JLD2, + LinearAlgebra, Logging, Makie, Measures, NCDatasets, NetCDF, Printf, Random, + SplitApplyCombine, StaticArrays, Statistics, StructArrays, VoronoiCells """ Coordinates are vector of vector of vector of points of the form: @@ -81,7 +81,7 @@ Coordinates are vector of vector of vector of points of the form: [[w1, z1], [w2, z2], ..., [wn, zn], [w1, z1]], ...] where the xy coordinates are the exterior border of the floe and the wz coordinates, or any other following sets of coordinates, describe holes within the floe. - This form is for easy conversion to LibGEOS Polygons. + This form is for easy conversion to polygons. """ const PolyVec{T} = Vector{Vector{Vector{T}}} where T<:Real @@ -89,7 +89,7 @@ const PolyVec{T} = Vector{Vector{Vector{T}}} where T<:Real Coordinates are vector of vector of points of the form: [[x1, y1], [x2, y2], ..., [xn, yn], [x1, y1]] where the xy coordinates form a closed ring. PolyVec objects can be made out RingVec objects. -This form is for each conversion to LibGEOS LinearRings, which can also be made into Polygons. +This form is for each conversion to LinearRings, which can also be made into Polygons. """ const RingVec{T} = R where { T<:Real, @@ -97,6 +97,8 @@ const RingVec{T} = R where { R <: AbstractArray{V}, } +const Polys{T} = GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{Tuple{T, T}}, Nothing, Nothing}}, Nothing, Nothing} where T <: AbstractFloat + # Model include("simulation_components/floe.jl") include("floe_utils.jl") diff --git a/src/eulerian_OA_code.jl b/src/eulerian_OA_code.jl index 5f1ce45..a3535cf 100644 --- a/src/eulerian_OA_code.jl +++ b/src/eulerian_OA_code.jl @@ -48,7 +48,7 @@ Outputs: """ function cell_area_ratio(cell_poly, floe_poly) floe_in_cell = LG.intersection(floe_poly, cell_poly) - return LG.area(floe_in_cell)/LG.area(cell_poly) + return GO.area(floe_in_cell)/GO.area(cell_poly) end """ diff --git a/src/floe_utils.jl b/src/floe_utils.jl index d26d634..da4cbe5 100644 --- a/src/floe_utils.jl +++ b/src/floe_utils.jl @@ -31,128 +31,39 @@ function valid_polyvec!(coords) return coords end -""" - find_poly_centroid(poly) - -Syntactic sugar for using LibGEOS to find a polygon's centroid -Input: - poly -Output: - Vector{Float64} [x, y] where this represents the centroid in the xy plane -""" -find_poly_centroid(poly) = - LG.GeoInterface.coordinates(LG.centroid(poly))::Vector{Float64} - """ find_poly_coords(poly) -Syntactic sugar for using LibGEOS to find a polygon's coordinates +Syntactic sugar for to find a polygon's coordinates Input: - poly + poly Output: representing the floe's coordinates xy plane """ -find_poly_coords(poly::LG.Polygon) = - LG.GeoInterface.coordinates(poly)::PolyVec{Float64} - -""" - get_polygons(geom) - -Returns an empty polygon list as non-polygon element was provided -Inputs: - geom -Outputs: - -""" -function get_polygons(geom) - polys = - if geom isa LG.Polygon - LG.area(geom) > 0 ? [geom] : Vector{LG.Polygon}() - elseif geom isa LG.MultiPolygon - sub_geoms = LG.getGeometries(geom) - for i in reverse(eachindex(sub_geoms)) - if LG.area(sub_geoms[i]) == 0 - deleteat!(sub_geoms, i) - end - end - sub_geoms - else - Vector{LG.Polygon}() - end - return polys::Vector{LG.Polygon} -end - -function get_polygons(collection::LG.GeometryCollection) - polys = Vector{LG.Polygon}() - for geom in LG.getGeometries(collection) - append!(polys, get_polygons(geom)) - end - return polys -end +find_poly_coords(poly::Polys) = GI.coordinates(poly) """ intersect_polys(p1, p2) Intersect two geometries and return a list of polygons resulting. Inputs: - p1 - p2 + p1 + p2 Output: - Vector of LibGEOS Polygons + Vector of Polygons """ -intersect_polys(p1, p2) = get_polygons( - LG.intersection( - p1, - p2, - )::LG.Geometry, -)::Vector{LG.Polygon} - -""" - intersect_coords(c1, c2) - -Intersect geometries using their coordinates to return list of resulting -polygons. -Inputs: - c1 - c2 -Output: - Vector of LibGEOS Polygons -""" -intersect_coords(c1, c2) = intersect_polys( - LG.Polygon(c1), - LG.Polygon(c2), -)::Vector{LG.Polygon} - +intersect_polys(p1, p2; kwargs...) = GO.intersection(p1, p2; target = GI.PolygonTrait(), fix_multipoly = nothing) +diff_polys(p1, p2; kwargs...) = GO.difference(p1, p2; target = GI.PolygonTrait(), fix_multipoly = nothing) +union_polys(p1, p2; kwargs...) = GO.union(p1, p2; target = GI.PolygonTrait(), fix_multipoly = nothing) -""" - polyvec_extrema(coords) +simplify_poly(p, tol) = GO.simplify(p; tol = tol) -Finds extremal x and y values for given coordiantes, defining a tight bounding -box. -Inputs: - coords -Outputs: xmin, xmax, ymin, ymax -""" -function polyvec_extrema(coords::PolyVec{FT}) where FT - xmin = ymin = FT(Inf) - xmax = ymax = FT(-Inf) - for i in eachindex(coords[1]) - x, y = coords[1][i] - if x < xmin - xmin = x - end - if x > xmax - xmax = x - end - if y < ymin - ymin = y - end - if y > ymax - ymax = y - end - end - return xmin, xmax, ymin, ymax -end +make_polygon(coords::PolyVec) = GI.Polygon(GO.tuples(coords)) +make_polygon(tuple_coords) = GI.Polygon(tuple_coords) +make_polygon(ring::GI.LinearRing) = GI.Polygon([ring]) +make_multipolygon(coords::Vector{<:PolyVec}) = GI.MultiPolygon(GO.tuples(coords)) +make_multipolygon(tuple_coords) = GI.MultiPolygon(tuple_coords) +make_multipolygon(polys::Vector{<:GI.Polygon}) = GI.MultiPolygon(polys) """ deepcopy_floe(floe::LazyRow{Floe{FT}}) @@ -209,31 +120,6 @@ function deepcopy_floe(floe::LazyRow{Floe{FT}}) where {FT} return f end -""" - find_multipoly_coords(poly) - -Syntactic sugar for using LibGEOS to find a multipolygon's coordinates -Input: - poly -Output: - representing the floe's coordinates xy plane -""" -find_multipoly_coords(poly::LG.Polygon) = - [find_poly_coords(poly)] - - -""" - find_multipoly_coords(multipoly) - -Syntactic sugar for using LibGEOS to find a multi-polygon's coordinates -Input: - poly -Output: - representing the floe's coordinates xy plane -""" -find_multipoly_coords(multipoly::LG.MultiPolygon) = - LG.GeoInterface.coordinates(multipoly)::Vector{PolyVec{Float64}} - """ translate!(coords, Δx, Δy) @@ -273,21 +159,6 @@ function translate!(coords::PolyVec{FT}, Δx, Δy) where {FT<:AbstractFloat} return end -""" - scale(poly::LG.Polygon, factor) - -Scale given polygon with respect to the reference point (0,0). -Scaling factor is applied to both the x and y directions. -Inputs: - coords - factor -Output: -""" -function scale(poly::LG.Polygon, factor) - coords = find_poly_coords(poly) - return LG.Polygon(coords .* factor) -end - """ rotate_radians!(coords::PolyVec, α) @@ -335,37 +206,18 @@ function hashole(coords::PolyVec{FT}) where FT<:AbstractFloat end """ - hashole(poly::LG.Polygon) + hashole(poly::Polys) Determine if polygon has one or more holes Inputs: - poly LibGEOS polygon + poly polygon Outputs: true if there is a hole in the polygons, else false """ -function hashole(poly::LG.Polygon) - return LG.numInteriorRings(poly) > 0 +function hashole(poly::Polys) + return GI.nhole(poly) > 0 end -""" - hashole(multipoly::LG.MultiPolygon) - -Determine if any of multipolygon's internal polygons has holes -Inputs: - multipoly LibGEOS multipolygon -Outputs: - true if there is a hole in any of the polygons, else false -""" -function hashole(multipoly::LG.MultiPolygon) - poly_lst = LG.getGeometries(multipoly)::Vector{LG.Polygon} - for poly in poly_lst - if hashole(poly) - return true - end - end - return false -end - """ rmholes(coords::PolyVec{FT}) @@ -386,165 +238,73 @@ function rmholes!(coords::PolyVec{FT}) where {FT<:AbstractFloat} end """ - rmholes(poly::LG.Polygon) + rmholes(poly::Polys) Remove polygon's holes if they exist Inputs: - poly LibGEOS polygon + poly polygon Outputs: - LibGEOS polygon without any holes + polygon without any holes """ -function rmholes(poly::LG.Polygon) +function rmholes(poly::Polys) if hashole(poly) - return LG.Polygon(LG.exteriorRing(poly)) + return make_polygon(GI.getexterior(poly)) end return poly end -""" - rmholes(multipoly::LG.MultiPolygon) - -Remove holes from each polygon of a multipolygon if they exist -Inputs: - multipoly multipolygon coordinates -Outputs: - multipolygon without any holes -""" -function rmholes(multipoly::LG.MultiPolygon) - poly_lst = LG.getGeometries(multipoly)::Vector{LG.Polygon} - nohole_lst = LG.Polygon[] - for poly in poly_lst - push!(nohole_lst, rmholes(poly)) - end - return LG.MultiPolygon(nohole_lst) -end - -""" - sortregions(poly::LG.Polygon) - -Returns given polygon within a vector as it is the only region -Inputs: - poly LibGEOS polygon -Outputs: - single LibGEOS polygon within list -""" -function sortregions(poly::LG.Polygon) - return [poly] -end - -""" - sortregions(multipoly::LG.MultiPolygon) - -Sorts polygons within a multi-polygon by area in descending order -Inputs: - multipoly LibGEOS multipolygon -Outputs: - list of LibGEOS polygons sorted in descending - order by area -""" -function sortregions(multipoly::LG.MultiPolygon) - return sort!(LG.getGeometries(multipoly), by=LG.area, rev=true) -end - -""" - separate_xy(coords::PolyVec{T}) - -Pulls x and y coordinates from standard polygon vector coordinates into seperate -vectors. Only keeps external coordinates and disregards holes -Inputs: - coords polygon coordinates -Outputs: - x x coordinates - y y coordinates -""" -function separate_xy(coords::PolyVec{T}) where {T<:AbstractFloat} - x = first.(coords[1]) - y = last.(coords[1]) - return x, y -end +#= + _calc_moment_inertia(::Type{T} poly, cent, h; ρi = 920.0) -""" - calc_moment_inertia(coords::PolyVec{T}, h; rhoice = 920.0) +Calculate the mass moment of intertia from a polygon given the polygon, its centroid, +height, and the density of ice in the simulation. Answer will be of given type T. -Calculate the mass moment of intertia from polygon coordinates. -Inputs: - coords - h height of floe - rhoice Density of ice -Output: - mass moment of inertia -Note: - Assumes that first and last point within the coordinates are the same. - Will not give correct answer otherwise. +Note: Assumes that first and last point within the coordinates are the same and will not +produce correct answer otherwise. Based on paper: Marin, Joaquin."Computing columns, footings and gates through moments of area." Computers & Structures 18.2 (1984): 343-349. -""" -function calc_moment_inertia( - coords::PolyVec{T}, - centroid, - h; - ρi = 920.0 -) where {T<:AbstractFloat} - x, y = separate_xy(coords) - x .-= centroid[1] - y .-= centroid[2] - N = length(x) - wi = x[1:N-1] .* y[2:N] - x[2:N] .* y[1:N-1] - Ixx = 1/12 * sum(wi .* ((y[1:N-1] + y[2:N]).^2 - y[1:N-1] .* y[2:N])) - Iyy = 1/12 * sum(wi .* ((x[1:N-1] + x[2:N]).^2 - x[1:N-1] .* x[2:N])) - return abs(Ixx + Iyy)*h*ρi; -end - -""" - calc_moment_inertia(poly::LG.Polygon, h; rhoice = 920.0) - -Calculate the mass moment of intertia from a LibGEOS polygon object using above -coordinate-based moment of intertia function. -Inputs: - poly LibGEOS.Polygon - h height of floe - rhoice Density of ice -Output: - mass moment of inertia -""" -calc_moment_inertia(poly::LG.Polygon, h; ρi = 920.0) = - calc_moment_inertia( - find_poly_coords(poly), - find_poly_centroid(poly), - h, - ρi = ρi, - ) - -""" - polyedge(p1, p2) - -Outputs the coefficients of the line passing through p1 and p2. -The line is of the form w1x + w2y + w3 = 0. -Inputs: - p1 [x, y] point - p2 [x, y] point -Outputs: - Three-element vector for coefficents of line passing through p1 and p2 -Note: - See note on calc_poly_angles for credit for this function. -""" -function polyedge(p1::Vector{<:FT}, p2) where FT - x1 = p1[1] - y1 = p1[2] - x2 = p2[1] - y2 = p2[2] - w = if x1 == x2 - [-1/x1, 0, 1] - elseif y1 == y2 - [0, -1/y1, 1] - elseif x1 == y1 && x2 == y2 - [1, 1, 0] - else - v = (y1 - y2)/(x1*(y2 - y1) - y1*(x2 - x1) + eps(FT)) - [v, -v*(x2 - x1)/(y2 - y1), 1] +=# +function _calc_moment_inertia( + ::Type{T}, + poly, + cent, + height; + ρi = 920.0, +) where T + xc, yc = GO._tuple_point(cent, T) + Ixx, Iyy = zero(T), zero(T) + x1, y1 = zero(T), zero(T) + for (i, p2) in enumerate(GI.getpoint(GI.getexterior(poly))) + (x2, y2) = GO._tuple_point(p2, T) + x2, y2 = x2 - xc, y2 - yc + if i == 1 + x1, y1 = x2, y2 + continue end - return w + wi = (x1 - xc) * (y2 - yc) - (x2 - xc) * (y1 - yc) + Ixx += wi * (y1^2 + y1 * y2 + y2^2) + Iyy += wi * (x1^2 + x1 * x2 + x2^2) + x1, y1 = x2, y2 + end + Ixx *= 1/12 + Iyy *= 1/12 + return abs(Ixx + Iyy) * T(height) * T(ρi) +end + +# Find the length of the maximum radius of a given polygon +function _calc_max_radius(::Type{T}, poly, cent) where T + max_rad_sqrd = zero(T) + Δx, Δy = GO._tuple_point(cent, T) + for pt in GI.getpoint(GI.getexterior(poly)) + x, y = GO._tuple_point(pt, T) + x, y = x - Δx, y - Δy + rad_sqrd = x^2 + y^2 + if rad_sqrd > max_rad_sqrd + max_rad_sqrd = rad_sqrd + end + end + return sqrt(max_rad_sqrd) end """ @@ -588,230 +348,6 @@ function orient_coords(coords::RingVec) return new_coords end -""" - convex_angle_test(coords::RingVec{T}) - -Determine which angles in the polygon are convex, with the assumption that the -first angle is convex, no other vertex has a smaller x-coordinate, and the -vertices are assumed to be ordered in a clockwise sequence. The test is based on -the fact that every convex vertex is on the positive side of the line passing -through the two vertices immediately following each vertex being considered. -Inputs: - coords Vector of [x, y] vectors that make up the exterior - of a polygon -Outputs: - sgn One element for each [x,y] pair - if 1 then - the angle at that vertex is convex, if it is -1 then the angle is - concave. -""" -function convex_angle_test(coords::RingVec{T}) where T - L = 10^25 - # Extreme points used in following loop, apended by a 1 for dot product - top_left = [-L, -L, 1] - top_right = [-L, L, 1] - bottom_left = [L, -L, 1] - bottom_right = [L, L, 1] - sgn = [1] # First vertex is convex - - for k in collect(2:length(coords)-1) - p1 = coords[k-1] - p2 = coords[k] # Testing this point for concavity - p3 = coords[k+1] - # Coefficents of polygon edge passing through p1 and p2 - w = polyedge(p1, p2) - - #= Establish the positive side of the line w1x + w2y + w3 = 0. - The positive side of the line should be in the right side of the vector - (p2- p3).Δx and Δy give the direction of travel, establishing which of - the extreme points (see above) should be on the + side. If that point is - on the negative side of the line, then w is replaced by -w. =# - Δx = p2[1] - p1[1] - Δy = p2[2] - p1[2] - if Δx == Δy == 0 - throw(ArgumentError("Data into convextiy test is 0 or duplicated")) - end - vector_product = - if Δx <= 0 && Δy >= 0 # Bottom_right should be on + side. - dot(w, bottom_right) - elseif Δx <= 0 && Δy <=0 # Top_right should be on + side. - dot(w, top_right) - elseif Δx>=0 && Δy<=0 # Top_left should be on + side. - dot(w, top_left) - else # Bottom_left should be on + side. - dot(w, bottom_left) - end - w *= sign(vector_product) - - # For vertex at p2 to be convex, p3 has to be on + side of line - if (w[1]*p3[1] + w[2]*p3[2] + w[3]) < 0 - push!(sgn, -1) - else - push!(sgn, 1) - end - end - return sgn -end - -""" - calc_poly_angles(coords::PolyVec{T}) - -Computes internal polygon angles (in degrees) of an arbitrary simple polygon. -The program eliminates duplicate points, except that the first row must equal -the last, so that the polygon is closed. -Inputs: - coords coordinates from a polygon -Outputs: - Vector of polygon's interior angles in degrees - -Note - Translated into Julia from the following program (including helper - functions convex_angle_test and polyedge): - Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins - Digital Image Processing Using MATLAB, Prentice-Hall, 2004 - Revision: 1.6 Date: 2003/11/21 14:44:06 -Warning - Assumes polygon has clockwise winding order. Use orient_coords! to - update coordinates prior to use -""" -function calc_poly_angles(coords::PolyVec{FT}) where {FT<:AbstractFloat} - # Calculate needed vectors - pdiff = diff(coords[1]) - npoints = length(pdiff) - v1 = -pdiff - v2 = vcat(pdiff[2:end], pdiff[1:1]) - v1_dot_v2 = [sum(v1[i] .* v2[i]) for i in collect(1:npoints)] - mag_v1 = sqrt.([sum(v1[i].^2) for i in collect(1:npoints)]) - mag_v2 = sqrt.([sum(v2[i].^2) for i in collect(1:npoints)]) - # Protect against division by 0 caused by very close points - replace!(mag_v1, 0=>eps(FT)) - replace!(mag_v2, 0=>eps(FT)) - angles = real.( - acos.( - clamp!(v1_dot_v2 ./ mag_v1 ./ mag_v2, FT(-1), FT(1)) - ) * 180 / pi - ) - - #= The first angle computed was for the second vertex, and the last was for - the first vertex. Scroll one position down to make the last vertex be the - first. =# - sangles = circshift(angles, 1) - # Now determine if any vertices are concave and adjust angles accordingly. - sgn = convex_angle_test(coords[1]) - for i in eachindex(sangles) - sangles[i] = (sgn[i] == -1) ? (-sangles[i] + 360) : sangles[i] - end - return sangles -end - -""" - calc_point_poly_dist(xp::Vector{T},yp::Vector{T}, vec_poly::PolyVec{T}) - -Compute the distances from each one of a set of np points on a 2D plane to a -polygon. Distance from point j to an edge k is defined as a distance from this -point to a straight line passing through vertices v(k) and v(k+1), when the -projection of point j on this line falls INSIDE segment k; and to the closest of -v(k) or v(k+1) vertices, when the projection falls OUTSIDE segment k. -Inputs: - xp x-coordinates of points to find distance from vec_poly - yp y-coordiantes of points to find distance from vec_poly - vec_poly coordinates of polygon -Outputs: - List of distances from each point to the polygon. If - the point is inside of the polygon the value will be negative. This does not - take holes into consideration. - -Note - Translated into Julia from the following program: -p_poly_dist by Michael Yoshpe - last updated in 2006. -We mimic version 1 functionality with 4 inputs and 1 output. -Only needed code was translated. -""" -function calc_point_poly_dist( - xp::Vector{FT}, - yp::Vector{FT}, - vec_poly::PolyVec{FT} -) where {FT<:AbstractFloat} - @assert length(xp) == length(yp) - min_lst = if !isempty(xp) - # Vertices in polygon and given points - Pv = reduce(hcat, valid_polyvec!(vec_poly)[1])' - Pp = hcat(xp, yp) - np = length(xp) - nv = length(vec_poly[1]) - # Distances between all points and vertices in x and y - x_dist = repeat(Pv[:, 1], 1, np)' .- repeat(Pp[:, 1], 1, nv) - y_dist = repeat(Pv[:, 2], 1, np)' .- repeat(Pp[:, 2], 1, nv) - p2c_dist = hypot.(x_dist, y_dist) - # minimum distance to vertices - min_dist, min_idx = findmin(p2c_dist, dims = 2) - # Coordinates of consecutive vertices - V1 = Pv[1:end-1, :] - V2 = Pv[2:end, :] - Δv = V2 .- V1 - # Vector of distances between each pair of consecutive vertices - vds = hypot.(Δv[:, 1], Δv[:, 2]) - - if (cumsum(vds)[end-1] - vds[end]) < 10eps(FT) - throw(ArgumentError("Polygon vertices should not lie on a straight \ - line")) - end - - #= Each pair of consecutive vertices V1[j], V2[j] defines a rotated - coordinate system with origin at V1[j], and x axis along the vector - V2[j]-V1[j]. cθ and sθ rotate from original to rotated system =# - cθ = Δv[:, 1] ./ vds - sθ = Δv[:, 2] ./ vds - Cer = zeros(FT, 2, 2, nv-1) - Cer[1, 1, :] .= cθ - Cer[1, 2, :] .= sθ - Cer[2, 1, :] .= -sθ - Cer[2, 2, :] .= cθ - - # Build origin translation vector P1r in rotated frame by rotating V1 - V1r = hcat(cθ .* V1[:, 1] .+ sθ .* V1[:, 2], - -sθ .* V1[:, 1] .+ cθ .* V1[:, 2]) - - #= Ppr is a 3D array of size 2*np*(nv-1). Ppr(1,j,k) is an X coordinate - of point j in coordinate systems defined by segment k. Ppr(2,j,k) is its - Y coordinate. =# - Ppr = zeros(FT, 2, np, nv-1) - # Rotation and Translation - Ppr[1, :, :] .= Pp * Cer[1, :, :] .- - permutedims(repeat(V1r[:, 1], 1, 1, np), [2, 3, 1])[1, :, :] - Ppr[2, :, :] .= Pp * Cer[2, :, :] .- - permutedims(repeat(V1r[:, 2], 1, 1, np), [2, 3, 1])[1, :, :] - - # x and y coordinates of the projected (cross-over) points in original - # coordinate - r = Ppr[1, :, :] - cr = Ppr[2, :, :] - B = fill(convert(FT, Inf), np, nv-1) - #= For the projections that fall inside the segments, find the minimum - distances from points to their projections (note, that for some points - these might not exist) =# - for i in eachindex(r) - if r[i] > 0 && r[i] < vds[cld(i, np)] - B[i] = cr[i] - end - end - cr_min, cr_min_idx = findmin(abs.(B), dims = 2) - #= For projections that fall outside segments, closest point is a vertex - These points have a negative value if point is actually outside of - polygon =# - in_poly = inpoly2(Pp, Pv) - dmin = cr_min - for i in eachindex(dmin) - if isinf(dmin[i]) || (cr_min_idx[i] != min_idx[i] && cr_min[i] > min_dist[i]) - dmin[i] = min_dist[i] - end - if in_poly[i, 1] || in_poly[i, 2] - dmin[i] *= -1 - end - end - dmin[:, 1] # Turn array into a vector - else - FT[] - end - return min_lst -end - """ intersect_lines(l1, l2) @@ -1070,187 +606,4 @@ function find_shared_edges_midpoint(c1::PolyVec{FT}, c2; atol = 1e-1) where {FT} ) end return mid_x, mid_y -end - -""" - cut_polygon_coords(poly_coords::PolyVec, yp) - -Cut polygon through the line y = yp and return the polygon(s) coordinates below -the line -Inputs: - poly_coords polygon coordinates - yp value of line to split polygon through using line - y = yp -Outputs: - new_polygons List of coordinates of polygons below line - y = yp. -Note: - Code translated from MATLAB to Julia. Credit for initial code to Dominik - Brands (2010) and Jasper Menger (2009). Only needed pieces of function are - translated (horizonal cut). -""" -function cut_polygon_coords(poly_coords::PolyVec{<:FT}, yp) where {FT} - # Loop through each edge - coord1 = poly_coords[1][1:end-1] - coord2 = poly_coords[1][2:end] - for i in eachindex(coord1) - x1, y1 = coord1[i] - x2, y2 = coord2[i] - # If both edge endpoints are above cut line, remove edge - if y1 > yp && y2 > yp - coord1[i] = [NaN, NaN] - coord2[i] = [NaN, NaN] - # If start point is above cut line, move down to intersection point - elseif y1 > yp - coord1[i] = [(yp - y2)/(y1 - y2) * (x1 - x2) + x2, yp] - # If end point is above cut line, move down to intersection point - elseif y2 > yp - coord2[i] = [(yp - y1)/(y2 - y1) * (x2 - x1) + x1, yp] - end - end - # Add non-repeat points to coordinate list for new polygon - new_poly_coords = [coord1[1]] - for i in eachindex(coord1) - if !isequal(coord1[i], new_poly_coords[end]) - push!(new_poly_coords, coord1[i]) - end - if !isequal(coord2[i], new_poly_coords[end]) - push!(new_poly_coords, coord2[i]) - end - end - - new_polygons = Vector{PolyVec{FT}}() - # Multiple NaN's indicate new polygon if they seperate coordinates - nanidx_all = findall(c -> isnan(sum(c)), new_poly_coords) - # If no NaNs, just add coordiantes to list - if isempty(nanidx_all) - if new_poly_coords[1] != new_poly_coords[end] - push!(new_poly_coords, deepcopy(new_poly_coords[1])) - end - if length(new_poly_coords) > 3 - push!(new_polygons, [new_poly_coords]) - end - # Seperate out NaNs to seperate multiple polygons and add to list - else - if new_poly_coords[1] == new_poly_coords[end] - new_poly_coords = new_poly_coords[1:end-1] - end - # Shift so each polygon's vertices are together in a section - new_poly_coords = circshift(new_poly_coords, -nanidx_all[1] + 1) - nanidx_all .-= (nanidx_all[1] - 1) - # Determine start and stop point for each polygon's coordinates - start_poly = nanidx_all .+ 1 - end_poly = [nanidx_all[2:end] .- 1; length(new_poly_coords)] - for i in eachindex(start_poly) - if start_poly[i] <= end_poly[i] - poly = new_poly_coords[start_poly[i]:end_poly[i]] - if poly[1] != poly[end] - push!(poly, deepcopy(poly[1])) - end - if length(poly) > 3 - push!(new_polygons, [poly]) - end - end - end - end - - return new_polygons -end - -""" - split_polygon_hole(poly::LG.Polygon, ::Type{FT} = Float64) - -Splits polygon horizontally through first hole and return lists of polygons -created by split. -Inputs: - poly polygon to split - Float type to run simulation with -Outputs: - <(Vector{LibGEOS.Polyon}, (Vector{LibGEOS.Polyon}> list of polygons created - from split through first hole below line and polygons through first hole - above line. Note that if there is no hole, a list of the original polygon - and an empty list will be returned -""" -function split_polygon_hole(poly::LG.Polygon) - bottom_list = Vector{LG.Polygon}() - top_list = Vector{LG.Polygon}() - if hashole(poly) # Polygon has a hole - poly_coords = find_poly_coords(poly) - full_coords = [poly_coords[1]] - h1 = LG.Polygon([poly_coords[2]]) # First hole - h1_center = find_poly_centroid(h1) - poly_bottom = LG.MultiPolygon( - cut_polygon_coords(full_coords, h1_center[2]) - ) - # Adds in any other holes in poly - poly_bottom = LG.intersection(poly_bottom, poly) - poly_top = LG.difference(poly, poly_bottom) - bottom_list = LG.getGeometries(poly_bottom)::Vector{LG.Polygon} - top_list = LG.getGeometries(poly_top)::Vector{LG.Polygon} - else # No hole - bottom_list, top_list = Vector{LG.Polygon}([poly]), Vector{LG.Polygon}() - end - return bottom_list, top_list -end - -""" - points_in_poly(xy, coords::PolyVec{<:AbstractFloat}) - -Determines if the provided points are within the given polygon, including -checking that the points are not in any holes. -Inputs: - xy n-by-2 matrix of element where each row is a point and - the first column is the x-coordinates and the second is y-coordinates - coords coordinates of polygon, with the exterior - coordinates as the first element of the vector, any any hole coordinates - as subsequent entries. -Outputs: - in_idx vector of booleans the length of the given points xy - where an entry is true if the corresponding element in xy is within the - given polygon. -""" -function points_in_poly(xy, coords::PolyVec{<:AbstractFloat}) - in_idx = fill(false, length(xy[:, 1])) - if !isempty(xy) && !isempty(coords[1][1]) - # Loop over exterior coords and each hole - for i in eachindex(coords) - in_on = inpoly2(xy, reduce(hcat, coords[i])') - if i == 1 # Exterior outline of polygon - points must be within - in_idx = in_idx .|| (in_on[:, 1] .| in_on[:, 2]) - else # Holes in polygon - points can't be within - in_idx = in_idx .&& .!(in_on[:, 1] .| in_on[:, 2]) - end - end - end - return in_idx -end - -""" - points_in_poly(xy, multi_coords::Vector{<:PolyVec{<:AbstractFloat}}) - -Determines if the provided points are within the given multipolygon, including -checking that the points are not in any holes of any of the polygons. -Inputs: - xy n-by-2 matrix of element where each row is a point and - the first column is the x-coordinates and the second is y-coordinates - coords coordinates of the multi-polygon, - with each element of the vector being a PolyVec of coordinates for a - polygon and within each polygon the exterior coordinates as the first - element of the PolyVec, any any hole coordinates as subsequent entries. -Outputs: - in_idx vector of booleans the length of the given points xy - where an entry is true if the corresponding element in xy is within the - given polygon. -""" -function points_in_poly(xy, multi_coords::Vector{<:PolyVec{<:AbstractFloat}}) - # Check which of the points are within the domain coords - in_idx = fill(false, length(xy[:, 1])) - if !isempty(xy) && !isempty(multi_coords[1][1][1]) - # Loop over every polygon - for i in eachindex(multi_coords) - # See if the points are within current polygon - in_idx = in_idx .|| points_in_poly(xy, multi_coords[i]) - end - end - return in_idx end \ No newline at end of file diff --git a/src/output.jl b/src/output.jl index 65d3f41..4636cca 100644 --- a/src/output.jl +++ b/src/output.jl @@ -841,15 +841,14 @@ function calc_eulerian_data!(floes, topography, writer) pint = potential_interactions[i,j,:] # If there are any potential interactions if sum(pint) > 0 - cell_poly = LG.Polygon(rect_coords( - writer.xg[j], - writer.xg[j+1], - writer.yg[i], - writer.yg[i+1], - )) + cell_poly_list = [make_polygon(rect_coords(writer.xg[j], writer.xg[j+1], writer.yg[i], writer.yg[i+1]))] if length(topography) > 0 - topography_poly = LG.MultiPolygon(topography.coords) - cell_poly = LG.difference(cell_poly, topography_poly) + cell_poly_list = diff_polys(make_multipolygon(cell_poly_list), make_multipolygon(topography.coords)) + end + + if length(cell_poly_list) == 0 + writer.data[j, i, :] .= 0.0 + continue end floeidx = collect(1:length(floes))[pint .== 1] @@ -858,13 +857,12 @@ function calc_eulerian_data!(floes, topography, writer) pic -> partially in cell - only includes pieces of floes that are within grid bounds =# - pic_polys = [ - LG.intersection( - cell_poly, - LG.Polygon(floes.coords[idx]), - ) for idx in floeidx] - - pic_area = [LG.area(poly) for poly in pic_polys] + pic_area = zeros(length(floeidx)) + for (i, idx) in enumerate(floeidx) + floe_poly = make_polygon(floes.coords[idx]) + pic_area[i] = mapreduce(x -> sum(GO.area, Subzero.intersect_polys(floe_poly, x); init = 0.0), +, cell_poly_list; init = 0.0) + end + floeidx = floeidx[pic_area .> 0] pic_area = pic_area[pic_area .> 0] fic = floes[floeidx] @@ -873,11 +871,11 @@ function calc_eulerian_data!(floes, topography, writer) area_ratios = pic_area ./ fic_area area_tot = sum(pic_area) - mass_tot = sum(floes.mass[floeidx] .* area_ratios) + mass_tot = sum(fic_mass .* area_ratios) if mass_tot > 0 # mass and area ratios - ma_ratios = area_ratios .* (floes.mass[floeidx] ./ mass_tot) + ma_ratios = area_ratios .* (fic_mass ./ mass_tot) outputs = writer.outputs for k in eachindex(outputs) data = if outputs[k] == :u_grid @@ -889,7 +887,7 @@ function calc_eulerian_data!(floes, topography, writer) elseif outputs[k] == :dvdt_grid sum(floes.p_dvdt[floeidx] .* ma_ratios) elseif outputs[k] == :si_frac_grid - area_tot/LG.area(cell_poly) + area_tot/sum(GO.area, cell_poly_list; init = 0.0) elseif outputs[k] == :overarea_grid sum(floes.overarea[floeidx])/length(floeidx) elseif outputs[k] == :mass_grid @@ -917,13 +915,13 @@ function calc_eulerian_data!(floes, topography, writer) end stress elseif outputs[k] == :strain_ux_grid - sum([s[1, 1] for s in floes[floeidx].strain] .* ma_ratios) + sum([s[1, 1] for s in fic.strain] .* ma_ratios) elseif outputs[k] == :strain_vx_grid - sum([s[1, 2] for s in floes[floeidx].strain] .* ma_ratios) + sum([s[1, 2] for s in fic.strain] .* ma_ratios) elseif outputs[k] == :strain_uy_grid - sum([s[2, 1] for s in floes[floeidx].strain] .* ma_ratios) + sum([s[2, 1] for s in fic.strain] .* ma_ratios) elseif outputs[k] == :strain_vy_grid - sum([s[2, 2] for s in floes[floeidx].strain] .* ma_ratios) + sum([s[2, 2] for s in fic.strain] .* ma_ratios) end writer.data[j, i, k] = data end diff --git a/src/physical_processes/collisions.jl b/src/physical_processes/collisions.jl index cf2d1d7..a06ff5b 100644 --- a/src/physical_processes/collisions.jl +++ b/src/physical_processes/collisions.jl @@ -48,40 +48,20 @@ function calc_normal_force( Δx = FT(coords[1][idx2][1] - coords[1][idx1][1]) Δy = FT(coords[1][idx2][2] - coords[1][idx1][2]) Δl = sqrt(Δx^2 + Δy^2) - if Δl > 0.1 # should match our scale + if Δl > 0.1 # If overlap is large enough to consider force_dir .= [-Δy/Δl; Δx/Δl] end elseif m != 0 # Unusual number of intersection points - x, y = separate_xy(coords) - Δx = diff(x) - xmid = (x[2:end] .+ x[1:end-1]) ./ 2 - Δy = diff(y) - ymid = (y[2:end] .+ y[1:end-1]) ./ 2 - mag = sqrt.(Δx.^2 .+ Δy.^2) # vector magniture - uvec = [-Δy./mag Δx./mag] # unit vector - xt = xmid.+uvec[:, 1]./100 - yt = ymid+uvec[:, 2]./100 # should match our scale - in_idx = points_in_poly(hcat(xt, yt), coords) - uvec[in_idx, :] *= -1 - Fn = -force_factor * (mag * ones(FT, 1, 2)) .* uvec - dmin_lst = calc_point_poly_dist(xmid, ymid, c1) - on_idx = findall(d->abs(d)<1e-8, dmin_lst) - if 0 < length(on_idx) < length(dmin_lst) - Δl = mean(mag[on_idx]) - if Δl > 0.1 - Fn_tot = sum(Fn[on_idx, :], dims = 1) - force_dir .= Fn_tot'/sqrt(Fn_tot*Fn_tot')[1] - end - end + Δl = _many_intersect_normal_force!(FT, force_dir, region, make_polygon(c1), force_factor) end # Check if direction of the force desceases overlap, else negate direction if Δl > 0.1 c1new = translate(c1, force_dir[1], force_dir[2]) # Floe/boudary intersection after being moved in force direction - new_regions_list = intersect_coords(c1new, c2) + new_regions_list = intersect_polys(make_polygon(c1new), make_polygon(c2)) # See if the area of overlap has increased in corresponding region for new_region in new_regions_list - if LG.intersects(new_region, region) && LG.area(new_region)/area > 1 + if GO.intersects(new_region, region) && GO.area(new_region)/area > 1 force_dir .*= -1 end end @@ -89,6 +69,55 @@ function calc_normal_force( return force_dir * area * force_factor, Δl end +#= + _many_intersect_normal_force!(::Type{T}, force_dir, region, poly, force_factor) + +Calculate the force direction (`force_dir`) given more than two points of intersection +within the region between two polygons (the first of which is `poly`). +=# +function _many_intersect_normal_force!(::Type{T}, force_dir, region, poly, force_factor) where T + x1, y1 = zero(T), zero(T) + Δl, n_pts = zero(T), 0 + Fn_tot = (zero(T), zero(T)) + # Loop over each edge within the overlap region + for (i, p) in enumerate(GI.getpoint(GI.getexterior(region))) + x2, y2 = T(GI.x(p)), T(GI.y(p)) + if i == 1 + x1, y1 = x2, y2 + continue + end + # Find the edge midpoint and calculate distance to first polygon + xmid, ymid = 0.5 * (x2 + x1), 0.5 * (y2 + y1) + dist = abs(GO.signed_distance((xmid, ymid), poly, T)) + if dist < 1e-8 # Only consider region edge points on first polygon + Δx, Δy = x2 - x1, y2 - y1 + mag = sqrt(Δx^2 + Δy^2) + #= If force would push edge points closer to second polygon (past and out of the + overlap region), switch the force direction =# + xt = xmid + (-Δy / 100mag) + yt = ymid + (Δx / 100mag) + in_region = GO.coveredby((xt, yt), region) + f_sign = in_region ? 1 : -1 + # Calculate force from given edge and incorporate it into the total forces + Fn = (f_sign * force_factor) .* (-Δy, Δx) + Δl += mag + n_pts += 1 + Fn_tot = Fn_tot .+ Fn + end + x1, y1 = x2, y2 + end + # Take the average of the summed values + if 0 < n_pts < (GI.npoint(region) - 1) + Δl /= n_pts + if Δl > 0.1 # If overlap is large enough to consider, set new force direction + norm_Fn = sqrt(GI.x(Fn_tot)^2 + GI.y(Fn_tot)^2) + force_dir[1] = GI.x(Fn_tot) / norm_Fn + force_dir[2] = GI.y(Fn_tot) / norm_Fn + end + end + return Δl +end + """ calc_elastic_forces( floe1, @@ -104,7 +133,7 @@ regions created from floe collisions Inputs: c1 first floe's coordinates in collision c2 second floe's coordinates in collision - regions polygon regions of overlap during + regions polygon regions of overlap during collision region_areas area of each polygon in regions force_factor Spring constant equivalent for collisions @@ -147,7 +176,7 @@ function calc_elastic_forces( Δl_lst = zeros(FT, ncontact) for k in 1:ncontact if region_areas[k] != 0 - cx, cy = find_poly_centroid(regions[k])::Vector{Float64} + cx, cy = GO.centroid(regions[k]) fpoint[k, 1] = cx fpoint[k, 2] = cy normal_force, Δl = calc_normal_force( @@ -355,11 +384,11 @@ function floe_floe_interaction!( Δt, max_overlap::FT, ) where {FT<:AbstractFloat} - inter_regions = intersect_coords(ifloe.coords, jfloe.coords) + inter_regions = intersect_polys(make_polygon(ifloe.coords), make_polygon(jfloe.coords)) region_areas = Vector{FT}(undef, length(inter_regions)) total_area = FT(0) for i in eachindex(inter_regions) - a = FT(LG.area(inter_regions[i])) + a = FT(GO.area(inter_regions[i])) region_areas[i] = a total_area += a end @@ -434,10 +463,11 @@ function floe_domain_element_interaction!( Δt, max_overlap, ) - floe_poly = LG.Polygon(floe.coords) - bounds_poly = LG.Polygon(boundary.coords) + floe_poly = make_polygon(floe.coords) + bounds_poly = make_polygon(boundary.coords) # Check if the floe and boundary actually overlap - if LG.area(LG.intersection(floe_poly, bounds_poly)) > 0 + inter_area = sum(GO.area, intersect_polys(floe_poly, bounds_poly); init = 0.0) + if inter_area > 0 floe.status.tag = remove end return @@ -617,11 +647,11 @@ function floe_domain_element_interaction!( Δt, max_overlap::FT, ) where {FT} - inter_regions = intersect_coords(floe.coords, element.coords) + inter_regions = intersect_polys(make_polygon(floe.coords), make_polygon(element.coords)) region_areas = Vector{FT}(undef, length(inter_regions)) max_area = FT(0) for i in eachindex(inter_regions) - a = FT(LG.area(inter_regions[i])) + a = GO.area(inter_regions[i], FT) region_areas[i] = a if a > max_area max_area = a @@ -1047,7 +1077,7 @@ function ghosts_on_bounds!( ) nfloes = length(floes) nghosts = 1 - if !isempty(intersect_coords(floes.coords[elem_idx], boundary.coords)) + if !isempty(intersect_polys(make_polygon(floes.coords[elem_idx]), make_polygon(boundary.coords))) # ghosts of existing ghosts and original element for i in floes.ghosts[elem_idx] push!( diff --git a/src/physical_processes/coupling.jl b/src/physical_processes/coupling.jl index b892df5..fa204f3 100644 --- a/src/physical_processes/coupling.jl +++ b/src/physical_processes/coupling.jl @@ -183,9 +183,9 @@ function generate_subfloe_points( mc_y = zeros(FT, point_generator.npoints) mc_in = fill(false, point_generator.npoints) # Find bounding box - xmin, xmax, ymin, ymax = polyvec_extrema(coords) - Δx = xmax - xmin - Δy = ymax - ymin + poly = make_polygon(GO.tuples(coords)) + (xmin, xmax), (ymin, ymax) = GI.extent(poly) + Δx, Δy = xmax - xmin, ymax - ymin while err > point_generator.err if count > 10 err = 0.0 @@ -193,7 +193,7 @@ function generate_subfloe_points( else mc_x .= xmin .+ Δx * rand(rng, FT, point_generator.npoints) mc_y .= ymin .+ Δy * rand(rng, FT, point_generator.npoints) - mc_in .= points_in_poly(hcat(mc_x, mc_y), coords) + mc_in .= [GO.coveredby((mc_x[i], mc_y[i]), poly) for i in eachindex(mc_x)] err = abs(sum(mc_in)/point_generator.npoints * (Δx * Δy) - area)/area count += 1 end @@ -324,7 +324,8 @@ function generate_subfloe_points( end x_sub_floe = repeat(x_interior_points, n_ypoints) y_sub_floe = repeat(y_interior_points, inner = n_xpoints) - in_floe = points_in_poly(hcat(x_sub_floe, y_sub_floe), coords) + poly = make_polygon(GO.tuples(coords)) + in_floe = [GO.coveredby((x_sub_floe[i], y_sub_floe[i]), poly) for i in eachindex(x_sub_floe)] append!(xpoints, x_sub_floe[in_floe]) append!(ypoints, y_sub_floe[in_floe]) @@ -1653,21 +1654,17 @@ function calc_two_way_coupling!( domain.north, domain.east ) - cell_poly = LG.Polygon(cell_coords) + cell_poly = make_polygon(cell_coords) for i in eachindex(floe_locations.floeidx) floe_coords = translate( floes.coords[floe_locations.floeidx[i]], floe_locations.Δx[i], floe_locations.Δy[i], ) - floe_poly = LG.Polygon(floe_coords) - floe_area_in_cell = FT(sum( - LG.area.(intersect_polys(cell_poly, floe_poly)) - )) - # floe_area_in_cell = FT(LG.area(LG.intersection( - # cell_poly, - # floe_poly, - # ))) + floe_poly = make_polygon(floe_coords) + floe_area_in_cell = sum( + GO.area.(intersect_polys(cell_poly, floe_poly), FT) + ) if floe_area_in_cell > 0 # Add forces and area to ocean fields ocean.τx[cartidx] += (τocn.τx[i]/τocn.npoints[i]) * floe_area_in_cell diff --git a/src/physical_processes/fractures.jl b/src/physical_processes/fractures.jl index 9f3aa2b..412a023 100644 --- a/src/physical_processes/fractures.jl +++ b/src/physical_processes/fractures.jl @@ -302,10 +302,9 @@ function determine_fractures( ) # Determine if floe stresses are in or out of criteria allowable regions update_criteria!(criteria, floes) - σvals = combinedims(sort.(eigvals.(floes.stress)))' - in_idx = points_in_poly(σvals, criteria.vertices) + critical_poly = make_polygon(GO.tuples(criteria.vertices)) # If stresses are outside of criteria regions, we will fracture the floe - frac_idx = .!(in_idx) + frac_idx = [!GO.coveredby(eigvals(s), critical_poly) for s in floes.stress] frac_idx[floes.area .< min_floe_area] .= false return range(1, length(floes))[frac_idx] end @@ -332,30 +331,29 @@ Outputs: """ function deform_floe!( floe, - deformer_coords, + deformer_coords::PolyVec{FT}, deforming_forces, floe_settings, Δt, rng, -) - poly = LG.Polygon(floe.coords) - deformer_poly = LG.Polygon(deformer_coords) - overlap_region = sortregions(LG.intersection(poly, deformer_poly))[1] +) where FT + poly = make_polygon(floe.coords) + deformer_poly = make_polygon(deformer_coords) + overlap_regions =intersect_polys(poly, deformer_poly) + max_overlap_area, max_overlap_idx = findmax(GO.area, overlap_regions) + overlap_region = overlap_regions[max_overlap_idx] # If floe and the deformer floe have an overlap area - if LG.area(overlap_region) > 0 + if max_overlap_area > 0 # Determine displacement of deformer floe - rcent = find_poly_centroid(overlap_region) - dist = calc_point_poly_dist( - [rcent[1]], - [rcent[2]], - find_poly_coords(overlap_region), - ) + region_cent = GO.centroid(overlap_region) + dist = GO.signed_distance(region_cent,overlap_region, FT) force_fracs = deforming_forces ./ 2norm(deforming_forces) Δx, Δy = abs.(dist)[1] .* force_fracs # Temporarily move deformer floe to find new shape of floe - deformer_poly = LG.Polygon(translate(deformer_coords, Δx, Δy)) - new_floe_poly = sortregions(LG.difference(poly, deformer_poly))[1] - new_floe_area = LG.area(new_floe_poly) + deformer_poly = make_polygon(translate(deformer_coords, Δx, Δy)) + new_floes = diff_polys(poly, deformer_poly) + new_floe_area, new_floe_idx = findmax(GO.area, new_floes) + new_floe_poly = new_floes[new_floe_idx] # If didn't change floe area by more than 90% if new_floe_area > 0 && new_floe_area/floe.area > 0.9 # Update floe shape and conserve mass @@ -422,36 +420,29 @@ function split_floe( ) if !isempty(pieces) # Intersect voronoi tesselation pieces with floe - floe_poly = LG.Polygon(rmholes(floe.coords)) - pieces_polys = Vector{LG.Polygon}() - for p in pieces - append!( - pieces_polys, - get_polygons( - LG.intersection(LG.Polygon(p), floe_poly) - ), - ) - end + floe_poly = make_polygon(rmholes(floe.coords)) + pieces_polys = mapreduce(p -> intersect_polys(make_polygon(p), floe_poly), append!, pieces; init = Vector{Polys{FT}}()) # Conserve mass within pieces - pieces_areas = [LG.area(p) for p in pieces_polys] + pieces_areas = [GO.area(p) for p in pieces_polys] total_area = sum(pieces_areas) # Create floes out of each piece for i in eachindex(pieces_polys) if pieces_areas[i] > 0 mass = floe.mass * (pieces_areas[i]/total_area) height = mass / (floe_settings.ρi * pieces_areas[i]) - pieces_floes = poly_to_floes( + poly_to_floes!( FT, + new_floes, pieces_polys[i], height, - 0; # Δh - range of random height difference between floes + 0, # Δh - range of random height difference between floes + floe.rmax; floe_settings = floe_settings, rng = rng, u = floe.u, v = floe.v, ξ = floe.ξ, ) - append!(new_floes, pieces_floes) end end end diff --git a/src/physical_processes/ridge_raft.jl b/src/physical_processes/ridge_raft.jl index 134f9a2..c9dca04 100644 --- a/src/physical_processes/ridge_raft.jl +++ b/src/physical_processes/ridge_raft.jl @@ -88,13 +88,15 @@ function remove_floe_overlap!( simp_settings, rng, ) where {FT <: AbstractFloat} - # Find new floe shape and regions - new_floe_poly = LG.difference( - LG.Polygon(floes.coords[shrink_idx]), - LG.Polygon(grow_floe_coords), - ) - new_floe_poly = LG.simplify(new_floe_poly, simp_settings.tol) - total_area = LG.area(new_floe_poly) + # Find new floe shapes and regions + regions = diff_polys(make_polygon(floes.coords[shrink_idx]), make_polygon(grow_floe_coords)) + total_area = zero(FT) + nregions = 0 + for (i, region) in enumerate(regions) + regions[i] = simplify_poly(region, simp_settings.tol) + total_area += GO.area(region) + nregions += 1 + end floe_num = 0 # How many new floes have been created from the regions # Changes in area / volume transfer_area = floes.area[shrink_idx] - total_area @@ -102,9 +104,6 @@ function remove_floe_overlap!( # If the transfer area is sufficent, create new floes if transfer_area > ridgeraft_settings.min_overlap_frac * floes.area[shrink_idx] transfer_vol = floes.area[shrink_idx] * floes.height[shrink_idx] - # Find regions - regions = get_polygons(new_floe_poly) - nregions = length(regions) # Reset shrinking index to parent floe and determine floe shift parent_Δx = floes.centroid[shrink_parent_idx][1] - floes.centroid[shrink_idx][1] @@ -113,11 +112,10 @@ function remove_floe_overlap!( parent_centroid = floes.centroid[shrink_parent_idx] # Update existing floes/ghosts regions for region in regions - region_area = LG.area(region) + region_area = GO.area(region) new_coords = find_poly_coords(region)::PolyVec{FT} - xmin, xmax, ymin, ymax = polyvec_extrema(new_coords) - Δx = xmax - xmin - Δy = ymax - ymin + (xmin, xmax), (ymin, ymax) = GI.extent(region) + Δx, Δy = xmax - xmin, ymax - ymin # Region is big enought to be a floe and has okay aspect ratio if ( region_area > floe_settings.min_floe_area && @@ -130,7 +128,7 @@ function remove_floe_overlap!( parent_Δy, ) rmholes!(new_coords) # remove holes in floe - new_poly = LG.Polygon(new_coords) # parent floe new region poly + new_poly = make_polygon(new_coords) # parent floe new region poly new_vol = region_area * floes.height[shrink_idx] transfer_vol -= new_vol buffer_length = length(pieces_buffer) diff --git a/src/physical_processes/simplification.jl b/src/physical_processes/simplification.jl index 7aa7275..6c0044b 100644 --- a/src/physical_processes/simplification.jl +++ b/src/physical_processes/simplification.jl @@ -62,19 +62,17 @@ function smooth_floes!( Δt, rng, ) where {FT <: AbstractFloat} - topo_coords = topography.coords for i in eachindex(floes) if length(floes.coords[i][1]) > simp_settings.max_vertices - poly = LG.simplify(LG.Polygon(floes.coords[i]), simp_settings.tol) - if !isempty(topo_coords) - poly = LG.difference(poly, LG.MultiPolygon(topo_coords)) + poly_list = [simplify_poly(make_polygon(floes.coords[i]), simp_settings.tol)] + if !isempty(topography) + poly_list = diff_polys(make_multipolygon(poly_list), make_multipolygon(topography.coords); fix_multipoly = nothing) end - poly_list = get_polygons(rmholes(poly))::Vector{LG.Polygon} simp_poly = if length(poly_list) == 1 poly_list[1] else - areas = [LG.area(p) for p in poly_list] + areas = [GO.area(p) for p in poly_list] _, max_idx = findmax(areas) poly_list[max_idx] end @@ -108,9 +106,9 @@ function smooth_floes!( floes.status[i].tag = fuse push!(floes.status[i].fuse_idx, j) else - jpoly = LG.Polygon(floes.coords[j]) - intersect_area = LG.area(LG.intersection(simp_poly, jpoly)) - if intersect_area/LG.area(jpoly) > collision_settings.floe_floe_max_overlap + jpoly = make_polygon(floes.coords[j]) + intersect_area = sum(GO.area, intersect_polys(simp_poly, jpoly); init = 0.0) + if intersect_area/GO.area(jpoly) > collision_settings.floe_floe_max_overlap floes.status[i].tag = fuse push!(floes.status[i].fuse_idx, j) end @@ -156,9 +154,9 @@ function fuse_two_floes!( # Create new polygon if they fuse rmholes!(keep_floe.coords) rmholes!(remove_floe.coords) - poly1 = LG.Polygon(keep_floe.coords)::LG.Polygon - poly2 = LG.Polygon(remove_floe.coords)::LG.Polygon - new_poly_list = get_polygons(LG.union(poly1, poly2))::Vector{LG.Polygon} + poly1 = make_polygon(keep_floe.coords) + poly2 = make_polygon(remove_floe.coords) + new_poly_list = union_polys(poly1, poly2) if length(new_poly_list) == 1 # if they fused, they will make one polygon new_poly = rmholes(new_poly_list[1]) # mark smaller floe for removal diff --git a/src/physical_processes/update_floe.jl b/src/physical_processes/update_floe.jl index c8bbe49..dd5cf39 100644 --- a/src/physical_processes/update_floe.jl +++ b/src/physical_processes/update_floe.jl @@ -16,7 +16,7 @@ Updates existing floe shape and related physical properties based of the polygon defining the floe. Inputs: floe floe to update - new_poly polygon representing new outline of floe + new_poly polygon representing new outline of floe new_mass mass of floe floe_settings simulation's settings for making floes rng random number generator @@ -31,22 +31,23 @@ function replace_floe!( rng, ) where {FT} # Floe shape - floe.centroid = find_poly_centroid(new_poly) + floe.centroid = collect(GO.centroid(new_poly)) floe.coords = find_poly_coords(new_poly)::PolyVec{FT} floe.coords = [orient_coords(floe.coords[1])] if floe.coords[1][1] != floe.coords[1][end] push!(floe.coords, floe.coords[1][1]) end - floe.area = LG.area(new_poly) + floe.area = GO.area(new_poly) floe.height = new_mass/(floe.area * floe_settings.ρi) floe.mass = new_mass - floe.moment = calc_moment_inertia( - floe.coords, + floe.moment = _calc_moment_inertia( + FT, + new_poly, floe.centroid, floe.height; ρi = floe_settings.ρi, ) - floe.angles = calc_poly_angles(floe.coords) + floe.angles = GO.angles(make_polygon(floe.coords)) floe.α = FT(0) translate!(floe.coords, -floe.centroid[1], -floe.centroid[2]) floe.rmax = sqrt(maximum([sum(c.^2) for c in floe.coords[1]])) diff --git a/src/physical_processes/welding.jl b/src/physical_processes/welding.jl index fb3d1c2..4d8eb36 100644 --- a/src/physical_processes/welding.jl +++ b/src/physical_processes/welding.jl @@ -131,10 +131,7 @@ function timestep_welding!( ) ) # Find intersection area - inter_area = LG.area(LG.intersection( - LG.Polygon(floes.coords[i]), - LG.Polygon(floes.coords[j]) - )) + inter_area = FT(sum(GO.area, intersect_polys(make_polygon(floes.coords[i]), make_polygon(floes.coords[j])); init = 0.0)) # Probability two floes will weld weld_prob = weld_settings.welding_coeff * (inter_area / floes.area[i]) diff --git a/src/simulation_components/domain_and_grid.jl b/src/simulation_components/domain_and_grid.jl index dc5c28c..f2a4495 100644 --- a/src/simulation_components/domain_and_grid.jl +++ b/src/simulation_components/domain_and_grid.jl @@ -421,7 +421,7 @@ const NonPeriodicBoundary = Union{ } """ - TopographyE{FT}<:AbstractDomainElement{FT} + TopographyElement{FT}<:AbstractDomainElement{FT} Singular topographic element with coordinates field storing where the element is within the grid. These are used to create the desired topography within the @@ -443,7 +443,9 @@ struct TopographyElement{FT}<:AbstractDomainElement{FT} throw(ArgumentError("Topography element maximum radius must be \ positive and non-zero.")) end - new{FT}(valid_polyvec!(rmholes(coords)), centroid, rmax) + topo_poly = make_polygon(rmholes(coords)) + topo_poly = GO.ClosedRing()(topo_poly) + new{FT}(GI.coordinates(topo_poly), centroid, rmax) end end @@ -468,31 +470,31 @@ TopographyElement(args...) = TopographyElement{Float64}(args...) """ TopographyElement{FT}(poly) -Constructor for topographic element with LibGEOS Polygon +Constructor for topographic element with Polygon Inputs: - poly + poly Output: Topographic element of abstract float type FT """ -function TopographyElement{FT}(poly::LG.Polygon) where {FT <: AbstractFloat} +function TopographyElement{FT}(poly::Polys) where {FT <: AbstractFloat} topo = rmholes(poly) - centroid = find_poly_centroid(topo) + cx, cy = GO.centroid(topo) coords = find_poly_coords(topo) # Move coordinates to be centered at origin to calculate maximum radius translate!( coords, - -centroid[1], - -centroid[2], + -cx, + -cy, ) rmax = sqrt(maximum([sum(c.^2) for c in coords[1]])) translate!( coords, - centroid[1], - centroid[2], + cx, + cy, ) return TopographyElement{FT}( coords, - centroid, + [cx, cy], rmax, ) end @@ -509,11 +511,7 @@ Output: function TopographyElement{FT}( coords::PolyVec, ) where {FT <: AbstractFloat} - return TopographyElement{FT}( - LG.Polygon( - convert(PolyVec{Float64}, coords), - ), - ) + return TopographyElement{FT}(make_polygon(convert(PolyVec{Float64}, coords))) end """ @@ -544,10 +542,10 @@ function initialize_topography_field( ::Type{FT}, coords, ) where {FT <: AbstractFloat} - topo_arr = StructArray{TopographyElement{FT}}(undef, length(coords)) - for i in eachindex(coords) - c = Subzero.rmholes(coords[i]) - topo_arr[i] = TopographyElement{FT}(c) + topo_multipoly = GO.DiffIntersectingPolygons()(GI.MultiPolygon(coords)) + topo_arr = StructArray{TopographyElement{FT}}(undef, GI.npolygon(topo_multipoly)) + for (i, p) in enumerate(GI.getpolygon(topo_multipoly)) + topo_arr[i] = TopographyElement{FT}(p) end return topo_arr end diff --git a/src/simulation_components/floe.jl b/src/simulation_components/floe.jl index 9c944f8..0e6fb02 100644 --- a/src/simulation_components/floe.jl +++ b/src/simulation_components/floe.jl @@ -190,7 +190,7 @@ Base.:(:)(a::InteractionFields, b::InteractionFields) = Int(a):Int(b) """ Floe{FT}( - poly::LG.Polygon, + poly::Polys, hmean, Δh; floe_settings = FloeSettings(), @@ -198,9 +198,9 @@ Base.:(:)(a::InteractionFields, b::InteractionFields) = Int(a):Int(b) kwargs... ) -Constructor for floe with LibGEOS Polygon +Constructor for floe with a polygon Inputs: - poly + poly hmean mean height for floes Δh variability in height for floes floe_settings settings needed to initialize floe @@ -213,7 +213,7 @@ Output: points were able to be generated. """ function Floe{FT}( - poly::LG.Polygon, + poly::Polys, hmean, Δh; floe_settings = FloeSettings(), @@ -222,21 +222,18 @@ function Floe{FT}( ) where {FT <: AbstractFloat} floe = rmholes(poly) # Floe physical properties - centroid = find_poly_centroid(floe) + centroid = collect(GO.centroid(floe)) height = clamp( hmean + (-1)^rand(rng, 0:1) * rand(rng, FT) * Δh, floe_settings.min_floe_height, floe_settings.max_floe_height, ) - area_tot = LG.area(floe) + area_tot = GO.area(floe) mass = area_tot * height * floe_settings.ρi coords = find_poly_coords(floe) coords = [orient_coords(coords[1])] - moment = calc_moment_inertia( - coords, centroid, height; - ρi = floe_settings.ρi, - ) - angles = calc_poly_angles(coords) + moment = _calc_moment_inertia(FT, floe, centroid, height; ρi = floe_settings.ρi) + angles = GO.angles(make_polygon(coords)) translate!(coords, -centroid[1], -centroid[2]) rmax = sqrt(maximum([sum(c.^2) for c in coords[1]])) status = Status() @@ -306,12 +303,7 @@ Floe{FT}( kwargs..., ) where {FT <: AbstractFloat} = Floe{FT}( # Polygon convert is needed since LibGEOS only takes Float64 - LG.Polygon( - convert( - PolyVec{Float64}, - valid_polyvec!(rmholes(coords)), - ), - ), + make_polygon(convert(PolyVec{Float64}, valid_polyvec!(rmholes(coords)))), hmean, Δh; floe_settings = floe_settings, @@ -320,69 +312,72 @@ Floe{FT}( ) """ - poly_to_floes( + poly_to_floes!( ::Type{FT}, - floe_poly, + floes, + poly, hmean, - Δh; + Δh, + rmax; floe_settings, rng = Xoshiro(), kwargs... ) -Split a given polygon into regions and split around any holes before turning -each region with an area greater than the minimum floe area into a floe. +Split a given polygon around any holes before turning each region with an area greater than +the minimum floe area into a floe. Inputs: Type{FT} Type for grid's numberical fields - determines simulation run type - floe_poly polygon/ - multipolygon to turn into floes + floes vector of floes to add new floes to + poly polygons to turn into floes hmean average floe height Δh height range - floes will range in height from hmean - Δh to hmean + Δh + rmax maximum radius of floe (could be larger given context) floe_settings settings needed to initialize floe settings rng random number generator to generate random floe attributes - default uses Xoshiro256++ algorithm kwargs... Any additional keywords to pass to floe constructor -Output: - vector of floes making up input polygon(s) with area - above given minimum floe area. Floe's with holes split around holes. """ -function poly_to_floes( +function poly_to_floes!( ::Type{FT}, - floe_poly, + floes, + poly, hmean, - Δh; + Δh, + rmax; floe_settings = FloeSettings(min_floe_area = 0), rng = Xoshiro(), kwargs... ) where {FT <: AbstractFloat} - floes = StructArray{Floe{FT}}(undef, 0) - regions = LG.getGeometries(floe_poly)::Vector{LG.Polygon} - while !isempty(regions) - r = pop!(regions) - a = LG.area(r) - if a >= floe_settings.min_floe_area && a > 0 - if !hashole(r) - floe = Floe( - FT, - r::LG.Polygon, - hmean, - Δh; - floe_settings = floe_settings, - rng = rng, - kwargs... - ) - push!(floes, floe) - else - region_bottom, region_top = split_polygon_hole(r) - append!(regions, region_bottom) - append!(regions, region_top) + a = GO.area(poly) + if a >= floe_settings.min_floe_area && a > 0 + if !hashole(poly) + floe = Floe( + FT, + poly::Polys, + hmean, + Δh; + floe_settings = floe_settings, + rng = rng, + kwargs... + ) + push!(floes, floe) + return 1 + else + cx, cy = GO.centroid(GI.gethole(poly, 1), FT) + new_regions = GO.cut(poly, GI.Line([(cx - rmax, cy), (cx + rmax, cy)]), FT) + n = 0 + for r in new_regions + n += poly_to_floes!(FT, floes, r, hmean, Δh, rmax; + floe_settings = floe_settings, rng = rng, kwargs...) end + return n end end - return floes + return 0 end """ @@ -432,24 +427,22 @@ function initialize_floe_field( rng = Xoshiro(), ) where {FT <: AbstractFloat} floe_arr = StructArray{Floe{FT}}(undef, 0) - floe_polys = [LG.Polygon(valid_polyvec!(c)) for c in coords] + floe_polys = [make_polygon(valid_polyvec!(c)) for c in coords] # Remove overlaps with topography if !isempty(domain.topography) - topo_poly = LG.MultiPolygon(domain.topography.coords) - floe_polys = [LG.difference(f, topo_poly) for f in floe_polys] + floe_polys = diff_polys(make_multipolygon(floe_polys), make_multipolygon(domain.topography.coords)) end # Turn polygons into floes for p in floe_polys - append!( - floe_arr, - poly_to_floes( - FT, - p, - hmean, - Δh; - floe_settings = floe_settings, - rng = rng, - ), + poly_to_floes!( + FT, + floe_arr, + p, + hmean, + Δh, + domain.east.val - domain.west.val; + floe_settings = floe_settings, + rng = rng, ) end # Warn about floes with area less than minimum floe size @@ -518,7 +511,8 @@ function generate_voronoi_coords( ) where {FT <: AbstractFloat} xpoints = Vector{FT}() ypoints = Vector{FT}() - area_frac = LG.area(LG.MultiPolygon(domain_coords)) / reduce(*, scale_fac) + domain_poly = make_multipolygon(GO.tuples(domain_coords)) + area_frac = GO.area(domain_poly) / reduce(*, scale_fac) # Increase the number of points based on availible percent of bounding box npoints = ceil(Int, desired_points / area_frac) current_points = 0 @@ -526,13 +520,11 @@ function generate_voronoi_coords( while current_points < desired_points && tries <= max_tries x = rand(rng, FT, npoints) y = rand(rng, FT, npoints) - # Scaled and translated points - st_xy = hcat( - scale_fac[1] * x .+ trans_vec[1], - scale_fac[2] * y .+ trans_vec[2] - ) - # Check which of the points are within the domain coords - in_idx = points_in_poly(st_xy, domain_coords) + # Check which of the scaled and translated points are within the domain coords + in_idx = [GO.coveredby( + (scale_fac[1] * x[i] .+ trans_vec[1], scale_fac[2] * y[i] .+ trans_vec[2]), + domain_poly + ) for i in eachindex(x)] current_points += sum(in_idx) tries += 1 append!(xpoints, x[in_idx]) @@ -566,7 +558,7 @@ function generate_voronoi_coords( ] tcoords[i] = [valid_ringvec!([ Vector(c) .* scale_fac .+ - trans_vec .+ perturb_vec + trans_vec #.+ perturb_vec for c in tess_cells[i] ])] end @@ -628,17 +620,16 @@ function initialize_floe_field( rng = Xoshiro(), ) where {FT <: AbstractFloat} floe_arr = StructArray{Floe{FT}}(undef, 0) + nfloes_added = 0 # Availible space in domain - open_water = LG.Polygon(floe_bounds) - open_water = LG.intersection(open_water, LG.Polygon(rect_coords(domain.west.val, domain.east.val, domain.south.val, domain.north.val))) + domain_poly = make_polygon(rect_coords(domain.west.val, domain.east.val, domain.south.val, domain.north.val)) + open_water = intersect_polys(make_polygon(floe_bounds), domain_poly) if !isempty(domain.topography) - open_water = LG.difference( - open_water, - LG.MultiPolygon(domain.topography.coords) - ) + open_water = diff_polys(make_multipolygon(open_water), make_multipolygon(domain.topography.coords)) end - (bounds_xmin, bounds_xmax), (bounds_ymin, bounds_ymax) = GeometryBasics.GeoInterface.extent(open_water) - open_water_area = LG.area(open_water) + open_water_mp = make_multipolygon(open_water) + (bounds_xmin, bounds_xmax), (bounds_ymin, bounds_ymax) = GI.extent(open_water_mp) + open_water_area = GO.area(open_water_mp, FT) # Split domain into cells with given concentrations nrows, ncols = size(concentrations[:, :]) @@ -654,19 +645,15 @@ function initialize_floe_field( if c > 0 c = c > 1 ? 1 : c # Grid cell bounds - xmin = bounds_xmin + collen * (j - 1) # TODO: change this + xmin = bounds_xmin + collen * (j - 1) ymin = bounds_ymin + rowlen * (i - 1) - cell_bounds = rect_coords( - xmin, - xmin + collen, - ymin, - ymin + rowlen, - ) trans_vec = [xmin, ymin] # Open water in cell - open_cell = LG.intersection(LG.Polygon(cell_bounds), open_water) - open_coords = find_multipoly_coords(open_cell) - open_area = LG.area(open_cell) + cell_init = make_polygon(rect_coords(xmin, xmin + collen, ymin, ymin + rowlen)) + open_cell = intersect_polys(cell_init, open_water_mp) + open_cell_mpoly = make_multipolygon(open_cell) + open_coords = [GI.coordinates(c) for c in open_cell] + open_area = sum(GO.area, open_cell; init = 0.0) # Generate coords with voronoi tesselation and make into floes ncells = ceil(Int, nfloes * open_area / open_water_area / c) floe_coords = generate_voronoi_coords( @@ -678,25 +665,28 @@ function initialize_floe_field( ncells, ) nfloes = length(floe_coords) - if nfloes > 0 - floes_area = FT(0.0) + if !isempty(floe_coords) + floe_poly_list = [make_polygon(c) for c in floe_coords] + nfloes = length(floe_poly_list) floe_idx = shuffle(rng, range(1, nfloes)) + floes_area = FT(0.0) while !isempty(floe_idx) && floes_area/open_area <= c idx = pop!(floe_idx) - floe_poly = LG.intersection( - LG.Polygon(floe_coords[idx]), - open_cell - ) - floes = poly_to_floes( - FT, - floe_poly, - hmean, - Δh; - floe_settings = floe_settings, - rng = rng, - ) - append!(floe_arr, floes) - floes_area += sum(floes.area) + poly_pieces_list = intersect_polys(floe_poly_list[idx], open_cell_mpoly) + for piece in poly_pieces_list + n_new_floes = poly_to_floes!( + FT, + floe_arr, + piece, + hmean, + Δh, + domain.east.val - domain.west.val; + floe_settings = floe_settings, + rng = rng, + ) + floes_area += sum(Iterators.drop(floe_arr.area, nfloes_added)) + nfloes_added += n_new_floes + end end end end diff --git a/src/tools/analyze_floe.jl b/src/tools/analyze_floe.jl index 6609923..62f724a 100644 --- a/src/tools/analyze_floe.jl +++ b/src/tools/analyze_floe.jl @@ -54,6 +54,6 @@ close(f) # plotting voronoi w/ labels scatter(xp, yp, markersize = 3, label = "centroids") annotate!([(xp[n] + 0.2, yp[n] + 0.03, Plots.text(id[n], 4)) for n in 1:159]) -verts = [Subzero.separate_xy(c) for c in coords] -plot!(first.(verts), last.(verts), seriestype = [:shape]) +xs, ys = first.(c[1]), last.(c[1]) +plot!(xs, ys, seriestype = [:shape]) savefig("output/sim/labeled_voronoi.png") \ No newline at end of file diff --git a/src/tools/file_convert.jl b/src/tools/file_convert.jl index 5f9a6b7..de6adc4 100644 --- a/src/tools/file_convert.jl +++ b/src/tools/file_convert.jl @@ -53,9 +53,8 @@ end # JLD2 to MATLAB function julfloe2matfloe(floes, Δg, out_fn) - coords = Subzero.separate_xy.(floes.coords) - xcoords = first.(coords) - ycoords = last.(coords) + xcoords = first.(floes.coords[1]) + ycoords = last.(floes.coords[1]) reshaped_x = Vector{Matrix{Float64}}() reshaped_y = Vector{Matrix{Float64}}() for i in range(1, length(xcoords)) diff --git a/test/runtests.jl b/test/runtests.jl index 46edb52..cb2542c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,7 @@ -using DataStructures, GeometryBasics, JLD2, Logging, NCDatasets, - PolygonInbounds, Random, SplitApplyCombine, Statistics, StructArrays, - Subzero, VoronoiCells -import LibGEOS as LG +using DataStructures, GeometryBasics, JLD2, Logging, NCDatasets, Random, SplitApplyCombine, + Statistics, StructArrays, Subzero, VoronoiCells +import GeometryOps as GO +import GeometryOps.GeoInterface as GI using Test @testset "Subzero.jl" begin diff --git a/test/test_floe.jl b/test/test_floe.jl index 5cdf4c7..471ce89 100644 --- a/test/test_floe.jl +++ b/test/test_floe.jl @@ -4,9 +4,9 @@ file = jldopen("inputs/floe_shapes.jld2", "r") floe_coords = file["floe_vertices"][1:end] close(file) - poly1 = LG.Polygon(Subzero.valid_polyvec!(floe_coords[1])) - centroid1 = LG.GeoInterface.coordinates(LG.centroid(poly1)) - area = LG.area(poly1) + poly1 = Subzero.make_polygon(Subzero.valid_polyvec!(floe_coords[1])) + centroid1 = GO.centroid(poly1) + area = GO.area(poly1) # Test InteractionFields enum interactions = range(1, 7)' @@ -25,7 +25,7 @@ @test typeof(floe_from_coords) <: Floe @test floe_from_coords.u == 0.2 @test 0.49 <= floe_from_coords.height <= 0.51 - @test floe_from_coords.centroid == centroid1 + @test floe_from_coords.centroid == collect(centroid1) @test floe_from_coords.area == area @test floe_from_coords.status.tag == Subzero.active @@ -41,69 +41,53 @@ @test floe_from_poly.u == 0.0 @test floe_from_poly.v == -0.2 @test 0.49 <= floe_from_poly.height <= 0.51 - @test floe_from_poly.centroid == centroid1 + @test floe_from_poly.centroid == collect(centroid1) @test floe_from_poly.area == area @test floe_from_poly.status.tag == Subzero.active # Test poly_to_floes rect_poly = [[[0.0, 0.0], [0.0, 5.0], [10.0, 5.0], [10.0, 0.0], [0.0, 0.0]]] - + rmax_rect = 2sqrt(5^2 + 2.5^2) c_poly_hole = [[[0.0, 0.0], [0.0, 10.0], [10.0, 10.0], [10.0, 0.0], [4.0, 0.0], [4.0, 6.0], [2.0, 6.0], [2.0, 0.0], [0.0, 0.0]], [[6.0, 4.0], [6.0, 6.0], [7.0, 6.0], [7.0, 4.0], [6.0, 4.0]]] - + rmax_cpoly = 2sqrt(5^2 + 5^2) # Test polygon with no holes - floe_arr = Subzero.poly_to_floes( + floe_arr = StructArray{Floe{FT}}(undef, 0) + n_new = Subzero.poly_to_floes!( FT, - LG.Polygon(rect_poly), + floe_arr, + Subzero.make_polygon(rect_poly), 0.25, 0.01, + rmax_rect; ) - @test length(floe_arr) == 1 - @test typeof(floe_arr[1]) <: Floe + @test n_new == 1 && length(floe_arr) == 1 @test !Subzero.hashole(floe_arr.coords[1]) # Test with polygon below minimum floe area - floe_arr = Subzero.poly_to_floes( + n_new = Subzero.poly_to_floes!( FT, - LG.Polygon(rect_poly), + floe_arr, + Subzero.make_polygon(rect_poly), 0.25, - 0.01; + 0.01, + rmax_rect; floe_settings = FloeSettings(min_floe_area = 55), ) - @test isempty(floe_arr) + @test n_new == 0 && length(floe_arr) == 1 # Test with polygon with a hole that is split into 3 polyons - floe_arr = Subzero.poly_to_floes( + n_new = Subzero.poly_to_floes!( FT, - LG.Polygon(c_poly_hole), + floe_arr, + Subzero.make_polygon(c_poly_hole), 0.25, 0.01, + rmax_cpoly, ) - @test length(floe_arr) == 3 + @test n_new == 3 && length(floe_arr) == 4 @test !any(Subzero.hashole.(floe_arr.coords)) - @test typeof(floe_arr) <: StructArray{<:Floe} - - # Test with multipolygon - floe_arr = Subzero.poly_to_floes( - FT, - LG.MultiPolygon([c_poly_hole, rect_poly]), - 0.25, - 0.01 - ) - @test length(floe_arr) == 4 - @test typeof(floe_arr) <: StructArray{<:Floe} - - # Test with multipolygon with minimum area - floe_arr = Subzero.poly_to_floes( - FT, - LG.MultiPolygon([c_poly_hole, rect_poly]), - 0.25, - 0.01; - floe_settings = FloeSettings(min_floe_area = 30), - ) - @test length(floe_arr) == 2 - @test typeof(floe_arr) <: StructArray{<:Floe} # Test initialize_floe_field from coord list grid = RegRectilinearGrid( @@ -126,7 +110,7 @@ wbound, StructVector([TopographyElement(t) for t in [island, topo1]]) ) - topo_polys = LG.MultiPolygon([island, topo1]) + topo_polys = Subzero.make_multipolygon([island, topo1]) # From file without topography -> floes below recommended area floe_arr = (@test_logs (:warn, "Some user input floe areas are less than the suggested minimum floe area.") initialize_floe_field( @@ -176,10 +160,8 @@ rng = Xoshiro(0) ) @test typeof(floe_arr) <: StructArray{<:Floe} - @test all([LG.area( - LG.intersection(LG.Polygon(c), topo_polys) - ) for c in floe_arr.coords] .< 1e-6) - + @test all([sum(GO.area, Subzero.intersect_polys(Subzero.make_polygon(c), topo_polys); init = 0.0) for c in floe_arr.coords] .< 1e-6) + # Test generate_voronoi_coords - general case domain_coords = [[[[1, 2], [1.5, 3.5], [1, 5], [2.5, 5], [2.5, 2], [1, 2]]]] bounding_box = [[[1, 2], [1, 5], [2.5, 5], [2.5, 2], [1, 2]]] @@ -192,16 +174,15 @@ 10, max_tries = 20, # 20 tries makes it very likely to reach 10 polygons ) - bounding_poly = LG.Polygon(bounding_box) + bounding_poly = Subzero.make_polygon(bounding_box) @test length(voronoi_coords) == 10 for c in voronoi_coords - fpoly = LG.Polygon(c) + fpoly = Subzero.make_polygon(c) @test isapprox( - LG.area(LG.intersection(fpoly, bounding_poly)), - LG.area(fpoly), + sum(GO.area, Subzero.intersect_polys(fpoly, bounding_poly); init = 0.0), + GO.area(fpoly), atol = 1e-3, ) - @test LG.isValid(fpoly) end # Test warning and no points generated @@ -228,14 +209,16 @@ rng = Xoshiro(1) ) @test isapprox( - sum(floe_arr.area)/(1.6e5*1.6e5 - LG.area(topo_polys)), + sum(floe_arr.area)/(1.6e5*1.6e5 - GO.area(topo_polys)), 0.5, atol = 1e-1 ) @test all(floe_arr.area .> 1e4) - @test all([LG.area( - LG.intersection(LG.Polygon(c), topo_polys) - ) for c in floe_arr.coords] .< 1e-6) + for (i, c) in enumerate(floe_arr.coords) + Subzero.intersect_polys(Subzero.make_polygon(c), topo_polys) + end + @test all([sum(GO.area, Subzero.intersect_polys(Subzero.make_polygon(c), topo_polys); init = 0.0) for c in floe_arr.coords] .< 1e-6) + nfloes = length(floe_arr) @test all(floe_arr.id .== range(1, nfloes)) @@ -251,25 +234,26 @@ rng = rng ) nfloes = length(floe_arr) - floe_polys = [LG.Polygon(f) for f in floe_arr.coords] + floe_polys = [Subzero.make_polygon(f) for f in floe_arr.coords] first_cell = [[[-8e4, -8e4], [-8e4, 0], [0, 0], [0, -8e4], [-8e4, -8e4]]] for j in 1:2 for i in 1:2 - cell = LG.Polygon( - Subzero.translate(first_cell, 8e4*(j-1), 8e4*(i-1)) - ) - open_cell_area = LG.area(LG.difference(cell, topo_polys)) + cell = Subzero.make_polygon(Subzero.translate(first_cell, 8e4*(j-1), 8e4*(i-1))) + cell_without_topos = Subzero.diff_polys(cell, topo_polys) + open_cell_area = sum(GO.area, cell_without_topos; init = 0.0) c = concentrations[i, j] - floes_in_cell = [LG.intersection(p, cell) for p in floe_polys] - @test c - 100eps() <= sum(LG.area.(floes_in_cell))/open_cell_area < 1 + eps() + floes_in_cell_area = 0 + n_polys = 0 + for floe in floe_polys + for mask in cell_without_topos + floes_in_cell_area += GO.area(Subzero.intersect_polys(floe, mask)) + end + end + @test c - 100eps() <= floes_in_cell_area/open_cell_area < 1 + eps() end end - @test all([LG.area( - LG.intersection(p, topo_polys) - ) for p in floe_polys] .< 1e-3) - @test all([LG.isValid(p) for p in floe_polys]) + @test all([sum(GO.area, Subzero.intersect_polys(p, topo_polys); init = 0.0) for p in floe_polys] .< 1e-3) @test all(floe_arr.id .== range(1, nfloes)) - @test typeof(initialize_floe_field( Float32, 25, diff --git a/test/test_floe_utils.jl b/test/test_floe_utils.jl index edaa2d4..c9097b2 100644 --- a/test/test_floe_utils.jl +++ b/test/test_floe_utils.jl @@ -2,12 +2,10 @@ ext = [[0.0, 1.0], [0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]] hole1 = [[0.2, 0.3], [0.2, 0.2], [0.3, 0.2], [0.3, 0.3], [0.2, 0.3]] hole2 = [[0.5, 0.6], [0.5, 0.5], [0.6, 0.5], [0.6, 0.6], [0.5, 0.6]] - poly_nohole = LG.Polygon([ext]) - poly_hole1 = LG.Polygon([ext, hole1]) - poly_hole2 = LG.Polygon([ext, hole1, hole2]) - multipoly_hole1 = LG.MultiPolygon([poly_nohole, poly_hole1]) - multipoly_hole2 = LG.MultiPolygon([poly_nohole, poly_hole1, poly_hole2]) - + poly_nohole = Subzero.make_polygon([ext]) + poly_hole1 = Subzero.make_polygon([ext, hole1]) + poly_hole2 = Subzero.make_polygon([ext, hole1, hole2]) + # Test validating/correcting RingVecs and PolyVecs @test Subzero.valid_ringvec!(ext) == ext invalid_ext = [[0.0, 1.0], [0.0, 0.0], [1.0, 0.0], [1.0, 1.0]] @@ -41,7 +39,6 @@ @test !Subzero.hashole(poly_nohole) @test Subzero.hashole(poly_hole1) @test Subzero.hashole(poly_hole2) - @test Subzero.hashole(multipoly_hole1) # Test removing holes from polygons and multipolygons @test [ext] == Subzero.rmholes([ext]) @@ -49,16 +46,9 @@ copy_holes = [ext, hole1] Subzero.rmholes!(copy_holes) @test copy_holes == [ext] - @test LG.equals(Subzero.rmholes(poly_nohole), poly_nohole) - @test LG.equals(Subzero.rmholes(poly_hole1), poly_nohole) - @test LG.equals(Subzero.rmholes(poly_hole2), poly_nohole) - @test !Subzero.hashole(Subzero.rmholes(multipoly_hole1)) - - # Test sorting regions from polygons and multipolygons - @test Subzero.sortregions(poly_nohole) == [poly_nohole] - poly_lst = Subzero.sortregions(multipoly_hole2) - @test LG.equals(poly_lst[1], poly_nohole) - @test LG.equals(poly_lst[3], poly_hole2) + @test GO.equals(Subzero.rmholes(poly_nohole), poly_nohole) + @test GO.equals(Subzero.rmholes(poly_hole1), poly_nohole) + @test GO.equals(Subzero.rmholes(poly_hole2), poly_nohole) # Test translating coordinates and polygons @test Subzero.translate([ext], 0.0, 0.0) == [ext] @@ -74,25 +64,12 @@ Subzero.translate!(test_trans, 1.5, -1.5) @test test_trans == [[[-0.5, 0.5], [-0.5, -0.5], [0.5, -0.5], [0.5, 0.5]]] - # Test scaling polygons - centered_coords = [[[-1.0, -1.0], [1.0, -1.0], [1.0, 1.0], - [-1.0, 1.0], [-1.0, -1.0]]] - centered_poly = LG.Polygon(centered_coords) - @test LG.equals(Subzero.scale(centered_poly, 1), centered_poly) - @test LG.equals(Subzero.scale(centered_poly, 2.0), - LG.Polygon(centered_coords .* 2.0)) - - # Test seperating x and y coordiantes from PolyVec form - x_ext, y_ext = Subzero.separate_xy([ext]) - @test x_ext == [0.0, 0.0, 1.0, 1.0, 0.0] - @test y_ext == [1.0, 0.0, 0.0, 1.0, 1.0] - # Test moment of intertia calculations - compared to values output my MATLAB - poly_moment = Subzero.calc_moment_inertia([ext], [0.5, 0.5], 0.25) + poly_moment = Subzero._calc_moment_inertia(Float64, Subzero.make_polygon([ext]), [0.5, 0.5], 0.25) @test isapprox(poly_moment, 38.333, atol = 0.001) - @test Subzero.calc_moment_inertia(poly_nohole, 0.25) == poly_moment - tri_coords = [[[0, 1], [0, 0], [1, 0], [0, 1]]] .* 6.67 - tri_moment = Subzero.calc_moment_inertia(LG.Polygon(tri_coords), 0.5) + @test Subzero._calc_moment_inertia(Float64, poly_nohole, GO.centroid(poly_nohole), 0.25) == poly_moment + tri_poly = Subzero.make_polygon([[[0, 1], [0, 0], [1, 0], [0, 1]]] .* 6.67) + tri_moment = Subzero._calc_moment_inertia(Float64, tri_poly, GO.centroid(tri_poly), 0.5) @test isapprox(tri_moment, 50581.145, atol = 0.001) # Test orient_coords @@ -127,63 +104,6 @@ [6.5e4, 4.5e4], ] - # Test polygon angles - some basic shapes and then compared to MATLAB - rect_coords = [[[1.0, 1.0], [1.0, 2.0], [2.0, 2.0], [2.0, 1.0]]] - @test Subzero.calc_poly_angles( - [Subzero.orient_coords(rect_coords[1])] - ) == [90.0, 90.0, 90.0, 90.0] - tri_coords = [[[0.0, 0.0], [0.0, 4.0], [3.0, 0.0]]] - @test prod(isapprox.( - Subzero.calc_poly_angles([Subzero.orient_coords(tri_coords[1])]), - [90.0, 36.8699, 53.1301], - atol = 0.001, - )) - concave_tri_coords = [[[-3.0, -2.0], [0.0,0.0], [5.0, 0.0]]] - @test prod(isapprox.( - Subzero.calc_poly_angles([Subzero.orient_coords(concave_tri_coords[1])]), - [19.6538, 146.3099, 14.0362], - atol = 0.001, - )) - # generate list of random polygons - polygon_lst = voronoicells( - rand(10), - rand(10), - Rectangle(Point2(0.0, 0.0), Point2(1.0, 1.0)), - ).Cells - for poly in polygon_lst - @test isapprox( - sum(Subzero.calc_poly_angles( - [Subzero.orient_coords(Vector{Vector{Float64}}(poly))] - )), - 180 * (length(poly) - 2), - atol = 1e-3, - ) - end - - # Test calc_point_poly_dist - some basic shapes and compared to MATLAB - xpoints = [0.0, 1.1, 1.1, 2.0] - ypoints = [0.0, 1.0, 1.1, 2.0] - @test prod(isapprox( - Subzero.calc_point_poly_dist(xpoints, ypoints, rect_coords), - [sqrt(2); 0; -0.1; 0.0], - atol = 0.001), - ) - @test_throws AssertionError Subzero.calc_point_poly_dist( - xpoints, - ypoints, - [[[1.0, 1.0], [1.0, 2.0]]], - ) - @test Subzero.calc_point_poly_dist( - Float64[], - Float64[], - rect_coords, - ) == Float64[] - @test_throws AssertionError Subzero.calc_point_poly_dist( - Float64[], - [1.0, 2.0], - rect_coords, - ) - # ------------------------- Test intersection of lines --------------------- l1 = [[[0.0, 0.0], [2.5, 0.0], [5.0, 0.0]]] l2 = [[[2.0, -3.0], [3.0, 0.0], [4.0, 3.0]]] @@ -302,239 +222,6 @@ offset_shared_v[2], offset_shared_v[1], ) == (7.5, 20) - # -------------- Test cutting polygon through horizontal line -------------- - # Cut a hexagon through the line y = -1 - poly_coords = [[ - [2.0, -3.0], - [0.0, 0.0], - [2.0, 2.0], - [6.0, 2.0], - [8.0, 0.0], - [6.0, -3.0], - [2.0, -3.0,], - ]] - poly = LG.Polygon(Subzero.cut_polygon_coords(poly_coords, -1)[1]) - @test LG.isValid(poly) - cut_coords = [[ - [2.0, -3.0], - [2-4/3, -1.0], - [22/3, -1.0], - [6.0, -3.0], - [2.0, -3.0], - ]] - @test Set(cut_coords[1]) == Set(LG.GeoInterface.coordinates(poly)[1]) - @test isapprox( - LG.GeoInterface.coordinates(LG.centroid(poly)), - LG.GeoInterface.coordinates(LG.centroid(LG.Polygon(cut_coords))), - atol = 1e-6, - ) - - # Cut a c-shaped polygon through the line y = 5, creating two polygons - poly_coords = [[ - [0.0, 0.0], - [0.0, 10.0], - [10.0, 10.0], - [10.0, 0.0], - [4.0, 0.0], - [4.0, 6.0], - [2.0, 6.0], - [2.0, 0.0], - [0.0, 0.0], - ]] - poly1_coords, poly2_coords = Subzero.cut_polygon_coords(poly_coords, 5) - poly1 = LG.Polygon(poly1_coords) - poly2 = LG.Polygon(poly2_coords) - @test LG.isValid(poly1) - @test LG.isValid(poly2) - @test Set([ - [4.0, 0.0], - [4.0, 5.0], - [10.0, 5.0], - [10.0, 0.0], - [4.0, 0.0], - ]) == Set(LG.GeoInterface.coordinates(poly1)[1]) - @test Set([ - [0.0, 0.0], - [0.0, 5.0], - [2.0, 5.0], - [2.0, 0.0], - [0.0, 0.0], - ]) == Set(LG.GeoInterface.coordinates(poly2)[1]) - - # Cut a triangle through the line y = -10 so only tip is left -> No polygon - poly_coords = [[[0.0, 0.0], [10.0, 0.0], [5.0, -10.0], [0.0, 0.0]]] - poly = Subzero.cut_polygon_coords(poly_coords, -10) - @test isempty(poly) - - # Cut a rectangle through line y = 0 so only bottom edge is left -> No polygon - poly_coords = [[ - [0.0, 0.0], - [0.0, 5.0], - [10.0, 5.0], - [10.0, 0.0], - [0.0, 0.0], - ]] - poly = Subzero.cut_polygon_coords(poly_coords, 0) - @test isempty(poly) - - # Cut a rectangle through the line y = -10 so no polygon is left - poly = Subzero.cut_polygon_coords(poly_coords, -10) - @test isempty(poly) - - # ------------------ Test polygon splitting through hole ------------------ - # Polygon with no holes -> returns original polygon and an empty list - poly_coords = [[ - [0.0, 0.0], - [0.0, 5.0], - [10.0, 5.0], - [10.0, 0.0], - [0.0, 0.0], - ]] - poly_bottom, poly_top = Subzero.split_polygon_hole(LG.Polygon(poly_coords)) - @test length(poly_bottom) == 1 - @test length(poly_top) == 0 - @test LG.area(LG.difference(poly_bottom[1], LG.Polygon(poly_coords))) == 0 - # Polygon with one hole -> creates two polygons - poly_coords = [ - [ - [0.0, 0.0], - [0.0, 5.0], - [10.0, 5.0], - [10.0, 0.0], - [0.0, 0.0], - ], - [ # hole! - [3.0, 2.0], - [3.0, 4.0], - [6.0, 4.0], - [6.0, 2.0], - [3.0, 2.0], - ], - ] - poly_bottom, poly_top = Subzero.split_polygon_hole(LG.Polygon(poly_coords)) - poly_bottom_coords = [[ - [0.0, 0.0], - [0.0, 3.0], - [3.0, 3.0], - [3.0, 2.0], - [6.0, 2.0], - [6.0, 3.0], - [10.0, 3.0], - [10.0, 0.0], - [0.0, 0.0], - ]] - poly_top_coords = [[ - [0.0, 3.0], - [0.0, 5.0], - [10.0, 5.0], - [10.0, 3.0], - [6.0, 3.0], - [6.0, 4.0], - [3.0, 4.0], - [3.0, 3.0], - [0.0, 3.0], - ]] - @test length(poly_bottom) == length(poly_top) == 1 - @test LG.area(LG.difference( - poly_bottom[1], - LG.Polygon(poly_bottom_coords), - )) == 0 - @test LG.area(LG.difference(poly_top[1], LG.Polygon(poly_top_coords))) == 0 - # Polygon with one holes -> creates three polygons - poly_coords = [ - [ - [0.0, 0.0], - [0.0, 10.0], - [10.0, 10.0], - [10.0, 0.0], - [4.0, 0.0], - [4.0, 6.0], - [2.0, 6.0], - [2.0, 0.0], - [0.0, 0.0], - ], - [ - [6.0, 4.0], - [6.0, 6.0], - [7.0, 6.0], - [7.0, 4.0], - [6.0, 4.0], - ], - ] - poly_bottom, poly_top = Subzero.split_polygon_hole(LG.Polygon(poly_coords)) - poly_bottom1_coords = [[ - [4.0, 0.0], - [4.0, 5.0], - [6.0, 5.0], - [6.0, 4.0], - [7.0, 4.0], - [7.0, 5.0], - [10.0, 5.0], - [10.0, 0.0], - [4.0, 0.0], - ]] - poly_bottom2_coords = [[ - [0.0, 0.0], - [0.0, 5.0], - [2.0, 5.0], - [2.0, 0.0], - [0.0, 0.0], - ]] - poly_top_coords = [[ - [0.0, 5.0], - [0.0, 10.0], - [10.0, 10.0], - [10.0, 5.0], - [7.0, 5.0], - [7.0, 6.0], - [6.0, 6.0], - [6.0, 5.0], - [4.0, 5.0], - [4.0, 6.0], - [2.0, 6.0], - [2.0, 5.0], - [0.0, 5.0], - ]] - @test length(poly_top) == 1 - @test length(poly_bottom) == 2 - @test LG.area(LG.difference( - poly_bottom[1], - LG.Polygon(poly_bottom1_coords), - )) == 0 - @test LG.area(LG.difference( - poly_bottom[2], - LG.Polygon(poly_bottom2_coords), - )) == 0 - @test LG.area(LG.difference( - poly_top[1], - LG.Polygon(poly_top_coords), - )) == 0 - - # Test points_in_poly for polygons - points = [-5 5; 2.5 2.5; 2 7; 5 6; 5.5 5.5; 9 2] - # No holes - coords = [[[0.0, 10.0], [0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0]]] - @test all(Subzero.points_in_poly(points, coords) .== - [false, true, true, true, true, true]) - # two holes - coords = [ - [[0.0, 10.0], [0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0]], - [[2.0, 3], [2, 2], [3, 2], [3, 3], [2, 3]], - [[5.0, 6], [5, 5], [6, 5], [6, 6], [5, 6]], - ] - @test all(Subzero.points_in_poly(points, coords) .== - [false, false, true, false, false, true]) - - # Test points_in_poly for multi-polygons - multi_coords = [coords, [[[12.0, 5], [15, 10], [18, 5], [12, 5]]]] - points = vcat(points, [15 7; 11 10]) - @test all(Subzero.points_in_poly(points, multi_coords) .== - [false, false, true, false, false, true, true, false]) - - # No coords / no points - @test all(Subzero.points_in_poly(points, [[[Vector{Float64}()]]]) .== - [false, false, false, false, false, false, false, false]) - @test all(Subzero.points_in_poly([], multi_coords) .== Vector{Bool}()) # ------------------ Test rotating coordinates ------------------ og_coords = [ diff --git a/test/test_model.jl b/test/test_model.jl index f7cde87..fa06aea 100644 --- a/test/test_model.jl +++ b/test/test_model.jl @@ -214,7 +214,7 @@ @testset "Topography" begin coords = [[[0.0, 1.0], [0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]] - poly = LG.Polygon(coords) + poly = Subzero.make_polygon(coords) # Polygon Constructor topo1 = Subzero.TopographyElement(poly) @test topo1.coords == coords @@ -244,7 +244,7 @@ # Create field of topography coords_w_hole = [ - [[0.0, 10.0], [0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0]], + [[0.5, 10.0], [0.5, 0.0], [10.0, 0.0], [10.0, 10.0], [0.5, 10.0]], [[2.0, 8.0], [2.0, 4.0], [8.0, 4.0], [8.0, 8.0], [2.0, 8.0]] ] topo_field_64 = initialize_topography_field( diff --git a/test/test_physical_processes/test_coupling.jl b/test/test_physical_processes/test_coupling.jl index 8dc469c..4cac171 100644 --- a/test/test_physical_processes/test_coupling.jl +++ b/test/test_physical_processes/test_coupling.jl @@ -6,16 +6,17 @@ file = jldopen("inputs/floe_shapes.jld2", "r") floe_coords = file["floe_vertices"][1:end] close(file) - poly1 = LG.Polygon(Subzero.valid_polyvec!(floe_coords[1])) - centroid1 = LG.GeoInterface.coordinates(LG.centroid(poly1)) + poly1 = Subzero.make_polygon(Subzero.valid_polyvec!(floe_coords[1])) + centroid1 = GO.centroid(poly1) origin_coords = Subzero.translate( floe_coords[1], -centroid1[1], -centroid1[2], ) - xo, yo = Subzero.separate_xy(origin_coords) + origin_poly = Subzero.make_polygon(origin_coords) + xo, yo = first.(origin_coords[1]), last.(origin_coords[1]) rmax = sqrt(maximum([sum(xo[i]^2 + yo[i]^2) for i in eachindex(xo)])) - area = LG.area(poly1) + area = GO.area(poly1) mc_x, mc_y, status = Subzero.generate_subfloe_points( MonteCarloPointsGenerator(), origin_coords, @@ -25,8 +26,7 @@ Xoshiro(1) ) @test length(mc_x) == length(mc_y) && length(mc_x) > 0 - in_on = inpoly2(hcat(mc_x, mc_y), hcat(xo, yo)) - mc_in = in_on[:, 1] .| in_on[:, 2] + mc_in = [GO.coveredby((mc_x[i], mc_y[i]), origin_poly) for i in eachindex(mc_x)] @test all(mc_in) xmin, xmax = extrema(xo) ymin, ymax = extrema(yo) @@ -281,9 +281,9 @@ periodic_bound, periodic_bound, ) - cell_poly = LG.Polygon(cell) - @test LG.area(cell_poly)::Float64 == 8 - @test LG.GeoInterface.coordinates(cell_poly) == + cell_poly = Subzero.make_polygon(cell) + @test GO.area(cell_poly)::Float64 == 8 + @test GI.coordinates(cell_poly) == [[[-9, -2], [-9, 2], [-7, 2], [-7, -2], [-9, -2]]] @test Subzero.center_cell_coords( 1, diff --git a/test/test_physical_processes/test_fractures.jl b/test/test_physical_processes/test_fractures.jl index 17a24fe..52b6171 100644 --- a/test/test_physical_processes/test_fractures.jl +++ b/test/test_physical_processes/test_fractures.jl @@ -11,14 +11,14 @@ ) isa HiblerYieldCurve # Test calculate_hibler hibler_verts = Subzero.calculate_hibler(0.5, 5e5, -1) - hibler_poly = LG.Polygon(hibler_verts) - @test isapprox(LG.area(hibler_poly), 49054437859.374, atol = -1e3) + hibler_poly = Subzero.make_polygon(hibler_verts) + @test isapprox(GO.area(hibler_poly), 49054437859.374, atol = -1e3) @test all(isapprox.( - Subzero.find_poly_centroid(hibler_poly), - [-1.25e5, -1.25e5], + GO.centroid(hibler_poly), + (-1.25e5, -1.25e5), atol = 1e-3 )) - x_verts, y_verts = Subzero.separate_xy(hibler_verts) + x_verts, y_verts = first.(hibler_verts[1]), last.(hibler_verts[1]) @test all(isapprox.( extrema(x_verts), [-264743.588, 14727.999], @@ -30,14 +30,14 @@ atol = 1e-3 )) hibler_verts = Subzero.calculate_hibler(0.25, 2.25e5, 20.0) - hibler_poly = LG.Polygon(hibler_verts) - @test isapprox(LG.area(hibler_poly), 2483380916.630, atol = -1e3) + hibler_poly = Subzero.make_polygon(hibler_verts) + @test isapprox(GO.area(hibler_poly), 2483380916.630, atol = -1e3) @test all(isapprox.( - Subzero.find_poly_centroid(hibler_poly), - [-28125, -28125], + GO.centroid(hibler_poly), + (-28125, -28125), atol = 1e-3 )) - x_verts, y_verts = Subzero.separate_xy(hibler_verts) + x_verts, y_verts = first.(hibler_verts[1]), last.(hibler_verts[1]) @test all(isapprox.( extrema(x_verts), [-59567.307, 3313.799], @@ -185,10 +185,7 @@ floe1_copy = deepcopy(floes[1]) colliding_coords = no_frac_floe.coords deforming_forces = frac_deform_floe.interactions[xforce:yforce] - init_overlap = LG.area(LG.intersection( - LG.Polygon(floe1_copy.coords), - LG.Polygon(colliding_coords), - )) + init_overlap = sum(GO.area, Subzero.intersect_polys(Subzero.make_polygon(floe1_copy.coords), Subzero.make_polygon(colliding_coords)); init = 0.0) Subzero.deform_floe!( floe1_copy, colliding_coords, @@ -197,10 +194,8 @@ 10, Xoshiro(1), ) - @test init_overlap > LG.area(LG.intersection( - LG.Polygon(floe1_copy.coords), # These coords have changed - LG.Polygon(colliding_coords), - )) + post_deform_overlap = sum(GO.area, Subzero.intersect_polys(Subzero.make_polygon(floe1_copy.coords), Subzero.make_polygon(colliding_coords)); init = 0.0) + @test init_overlap > post_deform_overlap @test all(isapprox.( floe1_copy.centroid, @@ -223,11 +218,11 @@ 10, ) # Test that the pieces all fit within original floe - og_floe_poly = LG.Polygon(floes.coords[1]) - new_floes_polys = LG.MultiPolygon(new_floes.coords) + og_floe_poly = Subzero.make_polygon(floes.coords[1]) + new_floes_polys = Subzero.make_multipolygon(new_floes.coords) @test isapprox( - LG.area(LG.intersection(new_floes_polys, og_floe_poly)), - LG.area(og_floe_poly), + sum(GO.area, Subzero.intersect_polys(new_floes_polys, og_floe_poly); init = 0.0), + GO.area(og_floe_poly), atol = 1e-6, ) # Conserve mass diff --git a/test/test_physical_processes/test_ridge_raft.jl b/test/test_physical_processes/test_ridge_raft.jl index 16b64c3..42e0da9 100644 --- a/test/test_physical_processes/test_ridge_raft.jl +++ b/test/test_physical_processes/test_ridge_raft.jl @@ -1,10 +1,10 @@ -using LibGEOS @testset "Ridging and Rafting" begin function update_height(floes, i, new_height, floe_settings) floes.height[i] = new_height floes.mass[i] = floes.area[i] * floes.height[i] * floe_settings.ρi - floes.moment[i] = Subzero.calc_moment_inertia( - floes.coords[i], + floes.moment[i] = Subzero._calc_moment_inertia( + Float64, + Subzero.make_polygon(floes.coords[i]), floes.centroid[i], floes.height[i], ρi = floe_settings.ρi, @@ -123,10 +123,8 @@ using LibGEOS p_x_momentum_init, p_y_momentum_init, ) if floe1_subsume || floe2_subsume - @test LibGEOS.area(LibGEOS.intersection( - LibGEOS.Polygon(floes.coords[1]), - LibGEOS.Polygon(floes.coords[2]) - )) == 0 # floes DO NOT overlap anymore! + inter_polys = Subzero.intersect_polys(Subzero.make_polygon(floes.coords[1]), Subzero.make_polygon(floes.coords[2])) + @test sum(GO.area, inter_polys; init = 0.0) == 0 # floes DO NOT overlap anymore! @test floe1_subsume ? (mass1 < floes.mass[1] && mass2 > floes.mass[2]) : (mass1 > floes.mass[1] && mass2 < floes.mass[2]) @@ -183,14 +181,10 @@ using LibGEOS @test area1 - bounds_overlap_area == floes.area[1] @test area2 - topo_overlap_area == floes.area[2] @test cent1 != floes.centroid[1] && cent2 != floes.centroid[2] - @test LibGEOS.area(LibGEOS.intersection( - LibGEOS.Polygon(floes.coords[1]), - boundary_poly, - )) == 0 - @test LibGEOS.area(LibGEOS.intersection( - LibGEOS.Polygon(floes.coords[2]), - topo_poly, - )) == 0 + inter_polys = Subzero.intersect_polys(Subzero.make_polygon(floes.coords[1]), boundary_poly) + @test sum(GO.area, inter_polys; init = 0.0) == 0 + inter_polys = Subzero.intersect_polys(Subzero.make_polygon(floes.coords[2]), topo_poly) + @test sum(GO.area, inter_polys; init = 0.0) == 0 else @test total_mass == sum(floes.mass) @test isapprox(h1, floes.height[1]) && isapprox(h2, floes.height[2]) @@ -277,11 +271,11 @@ using LibGEOS PeriodicBoundary(East, grid), PeriodicBoundary(West, grid), ) - boundary_poly = LibGEOS.union( - LibGEOS.MultiPolygon([collision_domain.north.coords, collision_domain.south.coords]), - LibGEOS.MultiPolygon([collision_domain.east.coords, collision_domain.west.coords]), - ) - topo_poly = LibGEOS.Polygon(topo_coords) + boundary_poly = GO.UnionIntersectingPolygons()(Subzero.make_multipolygon([ + collision_domain.north.coords, collision_domain.south.coords, + collision_domain.east.coords, collision_domain.west.coords])) + + topo_poly = Subzero.make_polygon(topo_coords) consts = Constants() floe_settings = FloeSettings(min_floe_area = 1e7) coupling_settings = CouplingSettings() @@ -418,14 +412,11 @@ using LibGEOS floes_base = setup_floes_with_inters(coords, collision_domain, consts, collision_settings, lock, [0.0, 0.5e4], [0.0, 0.5e4], ) - bounds_overlap_area = LibGEOS.area(LibGEOS.intersection( - LibGEOS.Polygon(floes_base.coords[1]), - boundary_poly, - )) - topo_overlap_area = LibGEOS.area(LibGEOS.intersection( - LibGEOS.Polygon(floes_base.coords[2]), - topo_poly, - )) + bounds_polys = Subzero.intersect_polys(Subzero.make_polygon(floes_base.coords[1]), boundary_poly) + bounds_overlap_area = sum(GO.area, bounds_polys; init = 0.0) + topos_polys = Subzero.intersect_polys(Subzero.make_polygon(floes_base.coords[2]), topo_poly) + topo_overlap_area = sum(GO.area, topos_polys; init = 0.0) + # Ridging with domain floes = deepcopy(floes_base) ridge_settings = Subzero.RidgeRaftSettings( @@ -617,7 +608,7 @@ using LibGEOS # Make sure floe ridged onto domain bounary and broke @test h1 < floes.height[1] @test h1 < pieces_list.height[1] - @test pieces_list.height[1] > floes.height[1] + @test isapprox(pieces_list.height[1], floes.height[1]) @test cent1 != floes.centroid[1] # Make sure IDs are correct @test max_id == 3 diff --git a/test/test_physical_processes/test_simplification.jl b/test/test_physical_processes/test_simplification.jl index 5791427..033e828 100644 --- a/test/test_physical_processes/test_simplification.jl +++ b/test/test_physical_processes/test_simplification.jl @@ -464,11 +464,9 @@ # Test mass is conserved @test total_mass == sum(floe_set2.mass) # Test first floe was cut by topography and only larger piece was kept - @test LG.area(LG.intersection( - LG.Polygon(floe_set2.coords[1]), - LG.Polygon(open_domain_with_topo.topography.coords[1]), - )) == 0 - @test LG.area(LG.Polygon(floe_set2.coords[1])) > 2og_f1_area/3 + @test sum(GO.area, Subzero.intersect_polys(Subzero.make_polygon(floe_set2.coords[1]), Subzero.make_polygon(open_domain_with_topo.topography.coords[1])); init = 0.0) == 0 + + @test GO.area(Subzero.make_polygon(floe_set2.coords[1])) > 2og_f1_area/3 # Test that both floes are tagged for fusion @test floe_set2.status[1].fuse_idx == [2] @test floe_set2.status[2].fuse_idx == [1] diff --git a/test/test_physical_processes/test_update_floe.jl b/test/test_physical_processes/test_update_floe.jl index df3afda..a9017de 100644 --- a/test/test_physical_processes/test_update_floe.jl +++ b/test/test_physical_processes/test_update_floe.jl @@ -76,7 +76,7 @@ [10.0, 10.0], [0.0, 0.0], ]] - tri_poly = LG.Polygon(triangle_coords) # this is a triangle + tri_poly = Subzero.make_polygon(triangle_coords) # this is a triangle Subzero.replace_floe!( f1, tri_poly, @@ -84,9 +84,9 @@ FloeSettings(), Xoshiro(1) ) - @test f1.centroid == Subzero.find_poly_centroid(tri_poly) + @test all(f1.centroid .== GO.centroid(tri_poly)) @test f1.coords == Subzero.find_poly_coords(tri_poly) - @test f1.area == LG.area(tri_poly) + @test f1.area == GO.area(tri_poly) @test f1.mass == mass1 @test f1.height * f1.area * 920.0 == f1.mass @test f1.α == 0 @@ -266,12 +266,10 @@ mass1 = sqr_floe.mass moment1 = sqr_floe.moment x1, y1 = sqr_floe.centroid + tri_rect_poly = Subzero.union_polys(Subzero.make_polygon(square_coords), Subzero.make_polygon(triangle_coords))[1] Subzero.replace_floe!( sqr_floe, - LG.union( - LG.Polygon(square_coords), - LG.Polygon(triangle_coords) - ), + tri_rect_poly, sqr_floe.mass + tri_floe.mass, FloeSettings(), Xoshiro(1)