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

@@ -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"

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,32 +49,26 @@ 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)
 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"),
 )

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),

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,

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],

src/Subzero.jl

@@ -69,34 +69,36 @@ 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:
 [[[x1, y1], [x2, y2], ..., [xn, yn], [x1, y1]], 
  [[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
 
 """
 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,
  V<:AbstractArray{T},
  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")

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
 
 """