From bd6ff6103930b9f41cc3204f3ddc19c64d8ed830 Mon Sep 17 00:00:00 2001 From: Skylar A Gering Date: Tue, 25 Jun 2024 15:59:32 -0400 Subject: [PATCH] Add polygon field to structs (#96) * 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 * Add polygon field to struct * Swap vertex finder code to use polygon * Finish merge * Update generate_subfloe_points to work with polygon * Add poly field to topography and boundaries * Clean up collision tests * Use translate polygon helper function * Remove unneeded make_polygon calls * Remove rmholes function * Update polygon rotation and movement * Debug angles and optimize * Add FT types to clipping * Update grid cell polygon process * Remove orient_coords function * Remove remaining un-needed make_polygon calls * Clean up Project.toml * Small PR fixes --- Project.toml | 8 +- examples/shear_flow.jl | 5 +- examples/simple_strait.jl | 2 +- src/Subzero.jl | 11 +- src/floe_utils.jl | 212 ++--- src/output.jl | 28 +- src/physical_processes/collisions.jl | 181 ++-- src/physical_processes/coupling.jl | 88 +- src/physical_processes/fractures.jl | 115 +-- src/physical_processes/process_settings.jl | 2 +- src/physical_processes/ridge_raft.jl | 46 +- src/physical_processes/simplification.jl | 37 +- src/physical_processes/update_floe.jl | 65 +- src/physical_processes/welding.jl | 6 +- src/simulation_components/domain_and_grid.jl | 48 +- src/simulation_components/floe.jl | 69 +- test/runtests.jl | 2 +- test/test_floe.jl | 95 +-- test/test_floe_utils.jl | 69 +- test/test_model.jl | 9 +- .../test_collisions.jl | 798 ++++++------------ test/test_physical_processes/test_coupling.jl | 101 +-- .../test_physical_processes/test_fractures.jl | 31 +- .../test_process_settings.jl | 2 +- .../test_ridge_raft.jl | 6 +- .../test_simplification.jl | 4 +- test/test_physical_processes/test_welding.jl | 2 +- 27 files changed, 720 insertions(+), 1322 deletions(-) diff --git a/Project.toml b/Project.toml index 815fbbb..ded0781 100644 --- a/Project.toml +++ b/Project.toml @@ -4,12 +4,11 @@ authors = ["Skylar Gering"] version = "0.1.0" [deps] -BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" 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" @@ -21,6 +20,7 @@ NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc" SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" @@ -28,12 +28,10 @@ StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" 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" +GeometryOps = "0.1.10" Interpolations = "0.14, 0.15" JLD2 = "0.4" Makie = "0.19, 0.20, 0.21" diff --git a/examples/shear_flow.jl b/examples/shear_flow.jl index 121a5b0..fdd9270 100644 --- a/examples/shear_flow.jl +++ b/examples/shear_flow.jl @@ -44,7 +44,7 @@ floe_settings = FloeSettings( # Floe creation floe_arr = initialize_floe_field( FT, - 50, + 500, [0.8], domain, hmean, @@ -72,14 +72,13 @@ simulation = Simulation( model = model, consts = consts, Δt = Δt, - nΔt = 2000, + nΔt = 5000, verbose = true, writers = writers, rng = Xoshiro(1), coupling_settings = coupling_settings, ) run_time!(simulation) - plot_sim( "output/shear_flow/floes.jld2", "output/shear_flow/initial_state.jld2", diff --git a/examples/simple_strait.jl b/examples/simple_strait.jl index cb0e0be..9bf500e 100644 --- a/examples/simple_strait.jl +++ b/examples/simple_strait.jl @@ -87,7 +87,7 @@ simulation = Simulation( model = model, consts = consts, Δt = Δt, - nΔt = 10000, + nΔt = 2500, verbose = true, floe_settings = floe_settings, coupling_settings = coupling_settings, diff --git a/src/Subzero.jl b/src/Subzero.jl index ceaffbd..860d332 100644 --- a/src/Subzero.jl +++ b/src/Subzero.jl @@ -71,9 +71,11 @@ export import Base.@kwdef # this is being exported as of version 1.9 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 +import StaticArrays as SA +using CairoMakie, CoordinateTransformations, DataStructures, Dates, Extents, + Interpolations, JLD2, LinearAlgebra, Logging, Makie, Measures, NCDatasets, NetCDF, + Printf, Random, Rotations, SplitApplyCombine, Statistics, StructArrays, + VoronoiCells """ Coordinates are vector of vector of vector of points of the form: @@ -99,6 +101,9 @@ const RingVec{T} = R where { const Polys{T} = GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{Tuple{T, T}}, Nothing, Nothing}}, Nothing, Nothing} where T <: AbstractFloat +Base.convert(::Type{Polys{Float32}}, p::Polys{Float64}) = GO.tuples(p, Float32) +Base.convert(::Type{Polys{Float64}}, p::Polys{Float32}) = GO.tuples(p, Float64) + # Model include("simulation_components/floe.jl") include("floe_utils.jl") diff --git a/src/floe_utils.jl b/src/floe_utils.jl index da4cbe5..3ff4e3d 100644 --- a/src/floe_utils.jl +++ b/src/floe_utils.jl @@ -40,7 +40,7 @@ Input: Output: representing the floe's coordinates xy plane """ -find_poly_coords(poly::Polys) = GI.coordinates(poly) +find_poly_coords(poly) = GI.coordinates(poly) """ intersect_polys(p1, p2) @@ -52,12 +52,46 @@ Inputs: Output: Vector of Polygons """ -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) - +intersect_polys(p1, p2, ::Type{FT} = Float64; kwargs...) where FT = GO.intersection(p1, p2, FT; target = GI.PolygonTrait(), fix_multipoly = nothing) +diff_polys(p1, p2, ::Type{FT} = Float64; kwargs...) where FT = GO.difference(p1, p2, FT; target = GI.PolygonTrait(), fix_multipoly = nothing) +union_polys(p1, p2, ::Type{FT} = Float64; kwargs...) where FT = GO.union(p1, p2, FT; target = GI.PolygonTrait(), fix_multipoly = nothing) simplify_poly(p, tol) = GO.simplify(p; tol = tol) +function _translate_poly(::Type{FT}, p, Δx, Δy) where FT + t = CoordinateTransformations.Translation(Δx, Δy) + # TODO: can remove the tuples call after GO SVPoint PR + return GO.tuples(GO.transform(t, p), FT) +end + +function _translate_floe!(::Type{FT}, floe, Δx, Δy) where FT + translate!(floe.coords, Δx, Δy) + floe.centroid[1] += Δx + floe.centroid[2] += Δy + floe.poly = _translate_poly(FT, floe.poly, Δx, Δy) + return +end + +function _move_poly(::Type{FT}, poly, Δx, Δy, Δα, cx = zero(FT), cy = zero(FT)) where FT + rot = CoordinateTransformations.LinearMap(Rotations.Angle2d(Δα)) + cent_rot = CoordinateTransformations.recenter(rot, (cx, cy)) + trans = CoordinateTransformations.Translation(Δx, Δy) + # TODO: can remove the tuples call after GO SVPoint PR + return GO.tuples(GO.transform(trans ∘ cent_rot, poly), FT) +end + +function _move_floe!(::Type{FT}, floe, Δx, Δy, Δα) where FT + cx, cy = floe.centroid + # Move coordinates and centroid + translate!(floe.coords, -cx, -cy) + rotate_radians!(floe.coords, Δα) + floe.centroid[1] += Δx + floe.centroid[2] += Δy + translate!(floe.coords, cx + Δx, cy + Δy) + # Move Polygon + floe.poly = _move_poly(FT, floe.poly, Δx, Δy, Δα, cx, cy) + return +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]) @@ -65,6 +99,14 @@ 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) +get_floe(floes::StructArray, i::Int) = LazyRow(floes, i) + +function _make_bounding_box_polygon(::Type{FT}, xmin, xmax, ymin, ymax) where FT + points = ((xmin, ymin), (xmin, ymax), (xmax, ymax), (xmax, ymin), (xmin, ymin)) + ring = GI.LinearRing(SA.SVector{5, Tuple{FT, FT}}(points)) + return GI.Polygon(SA.SVector(ring)) +end + """ deepcopy_floe(floe::LazyRow{Floe{FT}}) @@ -76,9 +118,11 @@ Outputs: they share values, but not referance. """ function deepcopy_floe(floe::LazyRow{Floe{FT}}) where {FT} + poly = GO.tuples(floe.poly, FT) f = Floe{FT}( + poly = poly, centroid = copy(floe.centroid), - coords = translate(floe.coords, 0, 0), + coords = find_poly_coords(poly), height = floe.height, area = floe.area, mass = floe.mass, @@ -180,18 +224,6 @@ function rotate_radians!(coords::PolyVec, α) return end -""" - rotate_degrees!(coords::PolyVec, α) - -Rotate a polygon's coordinates by α degrees around the origin. -Inputs: - coords polygon coordinates - α degrees to rotate the coordinates -Outputs: - Updates coordinates in place -""" -rotate_degrees!(coords::PolyVec, α) = rotate_radians!(coords, α * π/180) - """ hashole(coords::PolyVec{FT}) @@ -218,18 +250,6 @@ function hashole(poly::Polys) return GI.nhole(poly) > 0 end -""" - rmholes(coords::PolyVec{FT}) - -Remove polygon coordinates's holes if they exist -Inputs: - coords polygon coordinates -Outputs: - polygonc coordinates after removing holes -""" -function rmholes(coords::PolyVec{FT}) where {FT<:AbstractFloat} - return [coords[1]] -end function rmholes!(coords::PolyVec{FT}) where {FT<:AbstractFloat} if length(coords) > 1 @@ -237,20 +257,9 @@ function rmholes!(coords::PolyVec{FT}) where {FT<:AbstractFloat} end end -""" - rmholes(poly::Polys) - -Remove polygon's holes if they exist -Inputs: - poly polygon -Outputs: - polygon without any holes -""" -function rmholes(poly::Polys) - if hashole(poly) - return make_polygon(GI.getexterior(poly)) - end - return poly +function rmholes!(poly::Polys) + deleteat!(poly.geom, 2:GI.nring(poly)) + return end #= @@ -293,7 +302,7 @@ function _calc_moment_inertia( end # Find the length of the maximum radius of a given polygon -function _calc_max_radius(::Type{T}, poly, cent) where T +function calc_max_radius(poly, cent, ::Type{T}) where T max_rad_sqrd = zero(T) Δx, Δy = GO._tuple_point(cent, T) for pt in GI.getpoint(GI.getexterior(poly)) @@ -307,102 +316,6 @@ function _calc_max_radius(::Type{T}, poly, cent) where T return sqrt(max_rad_sqrd) end -""" - orient_coords(coords) - -Take given coordinates and make it so that the first point has the smallest -x-coordiante and so that the coordinates are ordered in a clockwise sequence. -Duplicates vertices will be removed and the coordiantes will be closed (first -and last point are the same). - -Input: - coords vector of points [x, y] -Output: - coords oriented clockwise with smallest x-coordinate first -""" -function orient_coords(coords::RingVec) - # extreem_idx is point with smallest x-value - if tie, choose lowest y-value - extreem_idx = 1 - for i in eachindex(coords) - ipoint = coords[i] - epoint = coords[extreem_idx] - if ipoint[1] < epoint[1] - extreem_idx = i - elseif ipoint[1] == epoint[1] && ipoint[2] < epoint[2] - extreem_idx = i - end - end - # extreem point must be first point in list - new_coords = similar(coords) - circshift!(new_coords, coords, -extreem_idx + 1) - valid_ringvec!(new_coords) - - # if coords are counterclockwise, switch to clockwise - orient_matrix = hcat( - ones(3), - vcat(new_coords[1]', new_coords[2]', new_coords[end-1]') # extreem/adjacent points - ) - if det(orient_matrix) > 0 - reverse!(new_coords) - end - return new_coords -end - -""" - intersect_lines(l1, l2) - -Finds the intersection points of two curves l1 and l2. The curves l1, l2 can be -either closed or open. In this version, l1 and l2 must be distinct. If no -intersections are found, the returned P is empty. -Inputs: - l1 line/polygon coordinates - l2 line/polygon coordinates -Outputs: - Set of points that are at the intersection of the - two line segments. -""" -function intersect_lines(l1::PolyVec{FT}, l2) where {FT} - points = Vector{Tuple{FT, FT}}() - for i in 1:(length(l1[1]) - 1) - # First line represented as a1x + b1y = c1 - x11 = l1[1][i][1] - x12 = l1[1][i+1][1] - y11 = l1[1][i][2] - y12 = l1[1][i+1][2] - a1 = y12 - y11 - b1 = x11 - x12 - c1 = a1 * x11 + b1 * y11 - for j in 1:(length(l2[1]) - 1) - # Second line represented as a2x + b2y = c2 - x21 = l2[1][j][1] - x22 = l2[1][j+1][1] - y21 = l2[1][j][2] - y22 = l2[1][j+1][2] - a2 = y22 - y21 - b2 = x21 - x22 - c2 = a2 * x21 + b2 * y21 - - determinant = a1 * b2 - a2 * b1 - # Find place there two lines cross - # Note that lines extend beyond given line segments - if determinant != 0 - x = (b2*c1 - b1*c2)/determinant - y = (a1*c2 - a2*c1)/determinant - p = (x, y) - # Make sure intersection is on given line segments - if min(x11, x12) <= x <= max(x11, x12) && - min(x21, x22) <= x <= max(x21, x22) && - min(y11, y12) <= y <= max(y11, y12) && - min(y21, y22) <= y <= max(y21, y22) && - !(p in points) - push!(points, p) - end - end - end - end - return points -end - """ which_vertices_match_points(ipoints, coords, atol) @@ -410,7 +323,7 @@ Find which vertices in coords match given points Inputs: points points to match to vertices within polygon - coords polygon coordinates + region polygon atol distance vertex can be away from target point before being classified as different points Output: @@ -419,24 +332,17 @@ Note: If last coordinate is a repeat of first coordinate, last coordinate index is NOT recorded. """ -function which_vertices_match_points( - points, - coords::PolyVec{FT}, - atol = 1, -) where {FT} +function which_vertices_match_points(points, region::Polys{FT}, atol = 1) where FT idxs = Vector{Int}() npoints = length(points) if points[1] == points[end] npoints -= 1 end - @views for i in 1:npoints # find which vertex matches point + for i in 1:npoints # find which vertex matches point min_dist = FT(Inf) min_vert = 1 - for j in eachindex(coords[1]) - dist = sqrt( - (coords[1][j][1] - points[i][1])^2 + - (coords[1][j][2] - points[i][2])^2, - ) + for (j, pt) in enumerate(GI.getpoint(GI.getexterior(region))) + dist = sqrt(GO.distance(pt, points[i], FT)) if dist < min_dist min_dist = dist min_vert = j diff --git a/src/output.jl b/src/output.jl index 4636cca..67f9d42 100644 --- a/src/output.jl +++ b/src/output.jl @@ -764,24 +764,6 @@ function auto_extension(filename, ext) return filename end -""" - rect_coords(xmin, xmax, ymin, ymax) - -PolyVec coordinates of a rectangle given minimum and maximum x and y coordinates -Inputs: - xmin minimum x coordinate of rectangle - xmax maximum x coordinate of rectangle - ymin minimum y coordiante of rectangle - ymax maximum y coordiante of rectangle -Output: - PolyVect coordinates for edges of rectangle with given minimums and maximums -""" -function rect_coords(xmin, xmax, ymin, ymax) -return [[[xmin, ymin], [xmin, ymax], - [xmax, ymax], [xmax, ymin], - [xmin, ymin]]] -end - """ grids_from_lines(xlines, ylines) @@ -809,7 +791,7 @@ Inputs: Output: Floe data averaged on eularian grid provided and saved in writer.data field """ -function calc_eulerian_data!(floes, topography, writer) +function calc_eulerian_data!(floes::FLT, topography, writer) where {FT <: AbstractFloat, FLT <: StructArray{<:Floe{FT}}} # Calculate/collect needed values Δx = writer.xg[2] - writer.xg[1] Δy = writer.yg[2] - writer.yg[1] @@ -841,9 +823,9 @@ function calc_eulerian_data!(floes, topography, writer) pint = potential_interactions[i,j,:] # If there are any potential interactions if sum(pint) > 0 - cell_poly_list = [make_polygon(rect_coords(writer.xg[j], writer.xg[j+1], writer.yg[i], writer.yg[i+1]))] + cell_poly_list = [_make_bounding_box_polygon(FT, writer.xg[j], writer.xg[j+1], writer.yg[i], writer.yg[i+1])] if length(topography) > 0 - cell_poly_list = diff_polys(make_multipolygon(cell_poly_list), make_multipolygon(topography.coords)) + cell_poly_list = diff_polys(make_multipolygon(cell_poly_list), make_multipolygon(topography.poly), FT) end if length(cell_poly_list) == 0 @@ -859,8 +841,8 @@ function calc_eulerian_data!(floes, topography, writer) =# 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) + floe_poly = floes.poly[idx] + pic_area[i] = mapreduce(x -> sum(GO.area, Subzero.intersect_polys(floe_poly, x, FT); init = 0.0), +, cell_poly_list; init = 0.0) end floeidx = floeidx[pic_area .> 0] diff --git a/src/physical_processes/collisions.jl b/src/physical_processes/collisions.jl index a06ff5b..49028c1 100644 --- a/src/physical_processes/collisions.jl +++ b/src/physical_processes/collisions.jl @@ -12,12 +12,11 @@ Functions needed for collisions between floes, boundaries, and topography force_factor, ) -Calculate normal force for collision between polygons with coordinates c1 and c2 -given an overlapping region, the area of that region, their intersection points -in the region, and a force factor. +Calculate normal force for collision between polygons p1 and p2 given an overlapping region, +the area of that region, their intersection points in the region, and a force factor. Inputs: - c1 first polygon coordinates - c2 second polygon coordinates + p1 first polygon + p2 second polygon region coordiantes for one region of intersection between the polygons area area of region @@ -29,36 +28,37 @@ Outputs: Δl mean length of distance between intersection points """ function calc_normal_force( - c1, - c2, + p1, + p2, region, area, ipoints, force_factor::FT, ) where {FT<:AbstractFloat} force_dir = zeros(FT, 2) - coords = find_poly_coords(region) # Identify which region coordinates are the intersection points (ipoints) - p = which_vertices_match_points(ipoints, coords) + p = which_vertices_match_points(ipoints, region) m = length(p) # Calculate force direction Δl = FT(0) if m == 2 # Two intersection points idx1, idx2 = p - Δx = FT(coords[1][idx2][1] - coords[1][idx1][1]) - Δy = FT(coords[1][idx2][2] - coords[1][idx1][2]) + pt1 = GI.getpoint(GI.getexterior(region), idx1) + pt2 = GI.getpoint(GI.getexterior(region), idx2) + Δx = FT(GI.x(pt2) - GI.x(pt1)) + Δy = FT(GI.y(pt2) - GI.y(pt1)) Δl = sqrt(Δx^2 + Δy^2) 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 - Δl = _many_intersect_normal_force!(FT, force_dir, region, make_polygon(c1), force_factor) + Δl = _many_intersect_normal_force!(FT, force_dir, region, p1, 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]) + p1new = _translate_poly(FT, p1, force_dir[1], force_dir[2])::Polys{FT} # Floe/boudary intersection after being moved in force direction - new_regions_list = intersect_polys(make_polygon(c1new), make_polygon(c2)) + new_regions_list = intersect_polys(p1new, p2, FT) # See if the area of overlap has increased in corresponding region for new_region in new_regions_list if GO.intersects(new_region, region) && GO.area(new_region)/area > 1 @@ -120,8 +120,8 @@ end """ calc_elastic_forces( - floe1, - floe2, + p1, + p2, regions, region_areas, force_factor, @@ -131,8 +131,8 @@ end Calculate normal forces, the point the force is applied, and the overlap area of regions created from floe collisions Inputs: - c1 first floe's coordinates in collision - c2 second floe's coordinates in collision + p1 first floe's polygon in collision + p2 second floe's polygon in collision regions polygon regions of overlap during collision region_areas area of each polygon in regions @@ -147,18 +147,17 @@ Outputs: Δl mean length of distance between intersection points """ function calc_elastic_forces( - c1, - c2, + p1, + p2, regions, region_areas, force_factor::FT, ) where {FT<:AbstractFloat} - ipoints = intersect_lines(c1, c2) # Intersection points + ipoints = GO.intersection_points(p1, p2) ncontact = 0 if !isempty(ipoints) && length(ipoints) >= 2 # Find overlapping regions greater than minumum area - n1 = length(c1[1]) - 1 - n2 = length(c2[1]) - 1 + n1, n2 = GI.npoint(p1) - 1, GI.npoint(p2) - 1 min_area = min(n1, n2) * 100 / 1.75 for i in reverse(eachindex(region_areas)) if region_areas[i] < min_area @@ -179,14 +178,7 @@ function calc_elastic_forces( cx, cy = GO.centroid(regions[k]) fpoint[k, 1] = cx fpoint[k, 2] = cy - normal_force, Δl = calc_normal_force( - c1, - c2, - regions[k], - region_areas[k], - ipoints, - force_factor - ) + normal_force, Δl = calc_normal_force(p1, p2, regions[k], region_areas[k], ipoints, force_factor) force[k, 1] = normal_force[1] force[k, 2] = normal_force[2] Δl_lst[k] = Δl @@ -212,7 +204,7 @@ Outputs: v v velocity at point (x, y) assuming it is on given floe """ function get_velocity( - floe::Union{LazyRow{Floe{FT}}, Floe{FT}}, + floe::FloeType{FT}, x::FT, y::FT, ) where {FT} @@ -384,7 +376,7 @@ function floe_floe_interaction!( Δt, max_overlap::FT, ) where {FT<:AbstractFloat} - inter_regions = intersect_polys(make_polygon(ifloe.coords), make_polygon(jfloe.coords)) + inter_regions = intersect_polys(ifloe.poly, jfloe.poly, FT) region_areas = Vector{FT}(undef, length(inter_regions)) total_area = FT(0) for i in eachindex(inter_regions) @@ -410,8 +402,8 @@ function floe_floe_interaction!( end # Calculate normal forces, force points, and overlap areas normal_forces, fpoints, overlaps, Δl = calc_elastic_forces( - ifloe.coords, - jfloe.coords, + ifloe.poly, + jfloe.poly, inter_regions, region_areas, force_factor @@ -456,17 +448,15 @@ Output: remove. Else, nothing is changed. """ function floe_domain_element_interaction!( - floe, + floe::FloeType{FT}, boundary::OpenBoundary, element_idx, consts, Δt, max_overlap, -) - floe_poly = make_polygon(floe.coords) - bounds_poly = make_polygon(boundary.coords) +) where FT # Check if the floe and boundary actually overlap - inter_area = sum(GO.area, intersect_polys(floe_poly, bounds_poly); init = 0.0) + inter_area = sum(GO.area, intersect_polys(floe.poly, boundary.poly, FT); init = 0.0) if inter_area > 0 floe.status.tag = remove end @@ -647,7 +637,7 @@ function floe_domain_element_interaction!( Δt, max_overlap::FT, ) where {FT} - inter_regions = intersect_polys(make_polygon(floe.coords), make_polygon(element.coords)) + inter_regions = intersect_polys(floe.poly, element.poly, FT) region_areas = Vector{FT}(undef, length(inter_regions)) max_area = FT(0) for i in eachindex(inter_regions) @@ -667,8 +657,8 @@ function floe_domain_element_interaction!( force_factor = consts.E * floe.height / sqrt(floe.area) # Calculate normal forces, force points, and overlap areas normal_forces, fpoints, overlaps, Δl = calc_elastic_forces( - floe.coords, - element.coords, + floe.poly, + element.poly, inter_regions, region_areas, force_factor, @@ -724,7 +714,8 @@ function update_boundary!( ) where {D <: Union{North, South}, FT <: AbstractFloat} Δd = boundary.v * Δt boundary.val += Δd - translate!(boundary.coords, FT(0), Δd) + translate!(boundary.coords, zero(FT), Δd) + _translate_poly(FT, boundary.poly, zero(FT), Δd) end """ update_boundary!(boundary, Δt) @@ -744,7 +735,8 @@ function update_boundary!( ) where {D <: Union{East, West}, FT <: AbstractFloat} Δd = boundary.u * Δt boundary.val += Δd - translate!(boundary.coords, Δd, FT(0)) + translate!(boundary.coords, Δd, zero(FT)) + _translate_poly(FT, boundary.poly, Δd, zero(FT)) end """ @@ -965,9 +957,9 @@ function timestep_collisions!( (ghost_id_pair[1] == g1) ⊻ (ghost_id_pair[2] == g2) ) floe_floe_interaction!( - LazyRow(floes, i), + get_floe(floes, i), i, - LazyRow(floes, j), + get_floe(floes, j), j, consts, Δt, @@ -977,7 +969,7 @@ function timestep_collisions!( end end floe_domain_interaction!( - LazyRow(floes, i), + get_floe(floes, i), domain, consts, Δt, @@ -1005,7 +997,7 @@ function timestep_collisions!( if j <= length(floes) && j > i jidx = Int(j) # add matching interaction (j with i) - add_interactions!(1, LazyRow(floes, jidx), FT(i), + add_interactions!(1, get_floe(floes, jidx), FT(i), i_inters[inter_idx:inter_idx, xforce:yforce], i_inters[inter_idx:inter_idx, xpoint:ypoint], i_inters[inter_idx:inter_idx, overlap:overlap] @@ -1028,7 +1020,7 @@ function timestep_collisions!( floes.interactions[g][1:gnp, ypoint] .-= floes.centroid[g][2] - floes.centroid[i][2] # add interactions - add_interactions!(gnp, LazyRow(floes, i), i, + add_interactions!(gnp, get_floe(floes, i), i, floes.interactions[g][1:gnp, xforce:yforce], floes.interactions[g][1:gnp, xpoint:ypoint], floes.interactions[g][1:gnp, overlap] @@ -1040,7 +1032,7 @@ function timestep_collisions!( end end # calculate interactions torque and total forces / torque - calc_torque!(LazyRow(floes, i)) + calc_torque!(get_floe(floes, i)) floes.collision_force[i][1] += sum( @view floes.interactions[i][1:floes.num_inters[i], xforce] ) @@ -1070,36 +1062,29 @@ Outputs: are added to the floe list """ function ghosts_on_bounds!( - floes, + floes::FLT, elem_idx, boundary, trans_vec, -) +) where {FT <: AbstractFloat, FLT <: StructArray{<:Floe{FT}}} nfloes = length(floes) nghosts = 1 - if !isempty(intersect_polys(make_polygon(floes.coords[elem_idx]), make_polygon(boundary.coords))) + if !isempty(intersect_polys(floes.poly[elem_idx], boundary.poly, FT)) # ghosts of existing ghosts and original element for i in floes.ghosts[elem_idx] - push!( - floes, - deepcopy_floe(LazyRow(floes, i)), - ) + push!(floes, deepcopy_floe(get_floe(floes, i))) nghosts += 1 end - push!( - floes, - deepcopy_floe(LazyRow(floes, elem_idx)), - ) + push!(floes, deepcopy_floe(get_floe(floes, elem_idx))) for i in (nfloes + 1):(nfloes + nghosts) - translate!(floes.coords[i], trans_vec[1], trans_vec[2]) - floes.centroid[i] .+= trans_vec + _translate_floe!(FT, get_floe(floes, i), trans_vec...) end end return end """ - find_ghosts( + find_ghosts!( floes, elem_idx, ebound::PeriodicBoundary, @@ -1120,7 +1105,7 @@ Outputs: the domain, else it is swapped with one of its ghost's which has a centroid within the domain. """ -function find_ghosts( +function find_ghosts!( floes, elem_idx, ebound::PeriodicBoundary{East, FT}, @@ -1130,40 +1115,27 @@ function find_ghosts( nfloes = length(floes) # passing through western boundary if (floes.centroid[elem_idx][1] - floes.rmax[elem_idx] < wbound.val) - ghosts_on_bounds!( - floes, - elem_idx, - wbound, - [Lx, FT(0)], - ) + ghosts_on_bounds!(floes, elem_idx, wbound, [Lx, FT(0)]) # passing through eastern boundary elseif (floes.centroid[elem_idx][1] + floes.rmax[elem_idx] > ebound.val) - ghosts_on_bounds!( - floes, - elem_idx, - ebound, - [-Lx, FT(0)], - ) + ghosts_on_bounds!(floes, elem_idx, ebound, [-Lx, FT(0)]) end # if element centroid isn't in domain's east/west direction, swap with ghost - if length(floes) > nfloes + new_nfloes = length(floes) + if new_nfloes > nfloes if floes.centroid[elem_idx][1] < wbound.val - translate!(floes.coords[elem_idx], Lx, FT(0)) - floes.centroid[elem_idx][1] += Lx - translate!(floes.coords[end], -Lx, FT(0)) - floes.centroid[end][1] -= Lx + _translate_floe!(FT, get_floe(floes, elem_idx), Lx, zero(FT)) + _translate_floe!(FT, get_floe(floes, new_nfloes), -Lx, zero(FT)) elseif ebound.val < floes.centroid[elem_idx][1] - translate!(floes.coords[elem_idx], -Lx, FT(0)) - floes.centroid[elem_idx][1] -= Lx - translate!(floes.coords[end], Lx, FT(0)) - floes.centroid[end][1] += Lx + _translate_floe!(FT, get_floe(floes, elem_idx), -Lx, zero(FT)) + _translate_floe!(FT, get_floe(floes, new_nfloes), Lx, zero(FT)) end end return end """ - find_ghosts( + find_ghosts!( floes, elem_idx, nbound::PeriodicBoundary{North, <:AbstractFloat}, @@ -1184,7 +1156,7 @@ Outputs: the domain, else it is swapped with one of its ghost's which has a centroid within the domain. """ -function find_ghosts( +function find_ghosts!( floes, elem_idx, nbound::PeriodicBoundary{North, FT}, @@ -1194,33 +1166,20 @@ function find_ghosts( nfloes = length(floes) # passing through southern boundary if (floes.centroid[elem_idx][2] - floes.rmax[elem_idx] < sbound.val) - ghosts_on_bounds!( - floes, - elem_idx, - sbound, - [FT(0), Ly], - ) + ghosts_on_bounds!(floes, elem_idx, sbound, [FT(0), Ly]) # passing through northern boundary elseif (floes.centroid[elem_idx][2] + floes.rmax[elem_idx] > nbound.val) - ghosts_on_bounds!( - floes, - elem_idx, - nbound, - [FT(0), -Ly], - ) + ghosts_on_bounds!(floes, elem_idx, nbound, [FT(0), -Ly]) end # if element centroid isn't in domain's north/south direction, swap with ghost - if length(floes) > nfloes + new_nfloes = length(floes) + if new_nfloes > nfloes if floes.centroid[elem_idx][2] < sbound.val - translate!(floes.coords[elem_idx], FT(0), Ly) - floes.centroid[elem_idx][2] += Ly - translate!(floes.coords[end], FT(0), -Ly) - floes.centroid[end][2] -= Ly + _translate_floe!(FT, get_floe(floes, elem_idx), zero(FT), Ly) + _translate_floe!(FT, get_floe(floes, new_nfloes), zero(FT), -Ly) elseif nbound.val < floes.centroid[elem_idx][2] - translate!(floes.coords[elem_idx], FT(0), -Ly) - floes.centroid[elem_idx][2] -= Ly - translate!(floes.coords[end], FT(0), Ly) - floes.centroid[end][2] += Ly + _translate_floe!(FT, get_floe(floes, elem_idx), zero(FT), -Ly) + _translate_floe!(FT, get_floe(floes, new_nfloes), zero(FT), Ly) end end return @@ -1248,7 +1207,7 @@ function add_floe_ghosts!( for i in eachindex(floes) # the floe is active in the simulation and a parent floe if floes.status[i].tag == active && floes.ghost_id[i] == 0 - find_ghosts( + find_ghosts!( floes, i, max_boundary, diff --git a/src/physical_processes/coupling.jl b/src/physical_processes/coupling.jl index fa204f3..7e1ffbe 100644 --- a/src/physical_processes/coupling.jl +++ b/src/physical_processes/coupling.jl @@ -146,8 +146,8 @@ SubGridPointsGenerator{FT}( """ generate_subfloe_points( point_generator - coords, - rmax, + poly, + centroid, area, status, rng @@ -157,8 +157,8 @@ Generate monte carlo points centered on the origin within the floe according to parameters defined in the point_generator argument. Inputs: point_generator monte carlo point generator - coords PolyVec of floe coords centered on origin - rmax floe's maximum radius + poly Polygon representing floe shape + centroid floe's centroid area floe's area status floe status (i.e. active, fuse in simulation) rng random number generator to generate monte @@ -171,8 +171,8 @@ Ouputs: """ function generate_subfloe_points( point_generator::MonteCarloPointsGenerator{FT}, - coords, - rmax, + poly, + centroid, area, status, rng @@ -183,7 +183,7 @@ function generate_subfloe_points( mc_y = zeros(FT, point_generator.npoints) mc_in = fill(false, point_generator.npoints) # Find bounding box - poly = make_polygon(GO.tuples(coords)) + poly = _translate_poly(FT, poly, -GI.x(centroid), -GI.y(centroid))::Polys{FT} (xmin, xmax), (ymin, ymax) = GI.extent(poly) Δx, Δy = xmax - xmin, ymax - ymin while err > point_generator.err @@ -210,8 +210,8 @@ end """ generate_subfloe_points( point_generator, - coords, - rmax, + poly, + centroid, area, status, rng @@ -221,8 +221,8 @@ Generate evenly spaced points within given floe coordinates to be used for coupling. If only one point falls within the floe, return the floe's centroid. Inputs: point_generator sub-grid point generator - coords PolyVec of floe coords centered on origin - rmax floe's maximum radius + poly Polygon representing floe shape + centroid floe's centroid area floe's area status floe status (i.e. active, fuse in simulation) rng random number generator is not used in @@ -234,44 +234,32 @@ Ouputs: """ function generate_subfloe_points( point_generator::SubGridPointsGenerator{FT}, - coords, - rmax, + poly, + centroid, area, status, rng ) where {FT <: AbstractFloat} - xmax = coords[1][1][1] - xmin = coords[1][1][1] - ymax = coords[1][1][2] - ymin = coords[1][1][2] - nverts = length(coords[1]) - + poly = _translate_poly(FT, poly, -GI.x(centroid), -GI.y(centroid))::Polys{FT} + (xmin, xmax), (ymin, ymax) = GI.extent(poly) xpoints = Vector{FT}() ypoints = Vector{FT}() - # Add points along edges - for i in 1:(nverts - 1) + local x1, y1 + for (i, p) in enumerate(GI.getpoint(GI.getexterior(poly))) # Determine points on edges - x1, y1 = coords[1][i] - x2, y2 = coords[1][i+1] - Δx = x2 - x1 - Δy = y2 - y1 - l = sqrt((Δx)^2 + (Δy)^2) - # Find maximum and minimum x and y points - if x1 > xmax - xmax = x1 - elseif x1 < xmin - xmin = x1 - end - if y1 > ymax - ymax = y1 - elseif y1 < ymin - ymin = y1 + x2, y2 = GI.x(p), GI.y(p) + if i == 1 + x1, y1 = x2, y2 + continue end + Δx, Δy = x2 - x1, y2 - y1 + l = sqrt((Δx)^2 + (Δy)^2) # Add current vertex push!(xpoints, x1) push!(ypoints, y1) # If distance between i and i+1 vertex is less than 2 sub-grid cells # but greater than one, add a point inbetween those two vertices + x2_unshifted, y2_unshifted = x2, y2 if l <= 2point_generator.Δg if l > point_generator.Δg push!(xpoints, x1 + Δx/2) @@ -298,13 +286,14 @@ function generate_subfloe_points( append!(xpoints, range(x1, x2, length = n_edge_points)) append!(ypoints, range(y1, y2, length = n_edge_points)) end + x1, y1 = x2_unshifted, y2_unshifted end # Add points in the interior of the floe n_xpoints = ceil(Int, (xmax - xmin) / point_generator.Δg) n_ypoints = ceil(Int, (ymax - ymin) / point_generator.Δg) x_interior_points = if n_xpoints < 3 n_xpoints = 1 - FT(0):FT(0) # coords are centered at the origin + FT(0):FT(0) # polygon is centered at the origin else range( xmin + point_generator.Δg/2, @@ -324,9 +313,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) - poly = make_polygon(GO.tuples(coords)) + # Check which points are within the polygon and add to list 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]) return xpoints, ypoints, status @@ -1126,12 +1114,13 @@ Output: outside of grid, could return a line at the edge of the boundary. """ function center_cell_coords( + ::Type{FT}, xidx::Int, yidx::Int, grid::RegRectilinearGrid, ns_bound, ew_bound, -) +) where FT xmin = (xidx - 1.5) * grid.Δx + grid.x0 xmax = xmin + grid.Δx ymin = (yidx - 1.5) * grid.Δy + grid.y0 @@ -1147,9 +1136,7 @@ function center_cell_coords( ns_bound, ew_bound, ) - return [[[xmin, ymin], [xmin, ymax], - [xmax, ymax], [xmax, ymin], - [xmin, ymin]]] + return _make_bounding_box_polygon(FT, xmin, xmax, ymin, ymax) end """ @@ -1511,7 +1498,7 @@ function calc_one_way_coupling!( for i in eachindex(floes) # Monte carlo point cartesian coordinates and grid cell indices npoints = calc_subfloe_values!( - LazyRow(floes, i), + get_floe(floes, i), grid, domain, cart_vals, @@ -1639,7 +1626,7 @@ function calc_two_way_coupling!( ) where {FT} # Determine force from floe on each grid cell it is in cell_area = grid.Δx * grid.Δy - Threads.@threads for cartidx in CartesianIndices(ocean.scells) + Threads.@threads for cartidx in CartesianIndices(ocean.scells) ocean.τx[cartidx] = FT(0) ocean.τy[cartidx] = FT(0) ocean.si_frac[cartidx] = FT(0) @@ -1647,21 +1634,20 @@ function calc_two_way_coupling!( floe_locations = grid.floe_locations[cartidx] if !isempty(floe_locations.floeidx) # Coordinates of grid cell - cell_coords = center_cell_coords( + cell_poly = center_cell_coords( + FT, cartidx[1], cartidx[2], grid, domain.north, domain.east ) - cell_poly = make_polygon(cell_coords) for i in eachindex(floe_locations.floeidx) - floe_coords = translate( - floes.coords[floe_locations.floeidx[i]], + floe_poly = _translate_poly(FT, + floes.poly[floe_locations.floeidx[i]], floe_locations.Δx[i], floe_locations.Δy[i], - ) - floe_poly = make_polygon(floe_coords) + )::Polys{FT} floe_area_in_cell = sum( GO.area.(intersect_polys(cell_poly, floe_poly), FT) ) diff --git a/src/physical_processes/fractures.jl b/src/physical_processes/fractures.jl index 412a023..4c9b88a 100644 --- a/src/physical_processes/fractures.jl +++ b/src/physical_processes/fractures.jl @@ -45,22 +45,7 @@ Note: mutable struct HiblerYieldCurve{FT<:AbstractFloat}<:AbstractFractureCriteria pstar::FT c::FT - vertices::PolyVec{FT} - - function HiblerYieldCurve{FT}( - pstar::Real, - c::Real, - vertices::PolyVec - ) where {FT<:AbstractFloat} - try - valid_polyvec!(vertices) - catch - throw(ArgumentError("The given vertices for the HiblerYieldCurve \ - can't be made into a valid polygon and thus the initial yield \ - curve can't be created.")) - end - new{FT}(pstar, c, vertices) - end + poly::Polys{FT} end """ @@ -82,7 +67,7 @@ correct constructor will be called with all other arguments. HiblerYieldCurve(args...) = HiblerYieldCurve{Float64}(args...) """ - calculate_hibler(floes, pstar, c) + _calculate_hibler(FT, floes, pstar, c) Calculate Hibler's Elliptical Yield Curve as described in his 1979 paper "A Dynamic Thermodynamic Sea Ice Model". @@ -97,18 +82,17 @@ Note: ice thickness and compactness. c is determined to that 10% open water reduces the strength substantially and pstar is considered a free parameter. """ -function calculate_hibler(mean_height, pstar, c) +function _calculate_hibler(::Type{FT}, mean_height, pstar, c) where FT compactness = 1 # Could be a user input with future development p = pstar*mean_height*exp(-c*(1-compactness)) - t = range(0, 2π, length = 100) + α_range = range(zero(FT), FT(2π), length = 100) a = p*sqrt(2)/2 b = a/2 - x = a*cos.(t) - y = b*sin.(t) - vertices = [splitdims([x'; y'])] - rotate_degrees!(vertices, 45) - translate!(vertices, -p/2, -p/2) - return valid_polyvec!(vertices) + ring_coords = [(a*cos(α), b*sin(α)) for α in α_range] + ring_coords[end] = ring_coords[1] # make sure first and last element are exactly the same + # TODO: eventually make with SVectors! + poly = GI.Polygon([ring_coords]) + return _move_poly(FT, poly, -p/2, -p/2, π/4) end """ @@ -121,18 +105,18 @@ Inputs: pstar used to tune ellipse for optimal fracturing c used to tune ellipse for optimal fracturing Outputs: - HiblerYieldCurve struct with vertices determined using the calculate_hibler + HiblerYieldCurve struct with vertices determined using the _calculate_hibler function. """ HiblerYieldCurve{FT}( - floes, + floes::StructArray{<:Floe{FT}}, pstar = 2.25e5, c = 20.0, -) where {FT <: AbstractFloat}= +) where {FT <: AbstractFloat} = HiblerYieldCurve{FT}( pstar, c, - calculate_hibler(mean(floes.height), pstar, c), + _calculate_hibler(FT, mean(floes.height), pstar, c), ) """ @@ -150,20 +134,7 @@ Note: Applied Physics 42.21 (2009): 214017. """ struct MohrsCone{FT<:AbstractFloat}<:AbstractFractureCriteria - vertices::PolyVec - - function MohrsCone{FT}( - vertices::PolyVec - ) where {FT <: AbstractFloat} - try - valid_polyvec!(vertices) - catch - throw(ArgumentError("The given vertices for the Mohr's Cone can't \ - be made into a valid polygon and thus the initial yield \ - curve can't be created.")) - end - new{FT}(vertices) - end + poly::Polys{FT} end """ @@ -185,7 +156,7 @@ constructor will be called with all other arguments. MohrsCone(args...) = MohrsCone{Float64}(args...) """ - calculate_mohrs(σ1, σ2, σ11, σ22) + _calculate_mohrs(FT, σ1, σ2, σ11, σ22) Creates PolyVec from vertex values for Mohr's Cone (triangle in 2D) Inputs: @@ -198,15 +169,15 @@ Inputs: Output: Mohr's Cone vertices (triangle since we are in 2D) in principal stress space """ -calculate_mohrs(σ1, σ2, σ11, σ22) = valid_polyvec!([[ - [σ1, σ2], - [σ11, σ22], - [σ22, σ11], - [σ1, σ2], -]]) +function _calculate_mohrs(::Type{FT}, σ1, σ2, σ11, σ22) where FT + # TODO: eventually make with SVectors! + points = [(σ1, σ2), (σ11, σ22), (σ22, σ11), (σ1, σ2)] + return GI.Polygon([points]) +end """ - calculate_mohrs( + _calculate_mohrs( + FT, q, σc, σ11; @@ -232,23 +203,24 @@ Note: Applied Physics 42.21 (2009): 214017. Equations taken from original version of Subzero written in MATLAB """ -function calculate_mohrs( +function _calculate_mohrs( + ::Type{FT}, q = 5.2, σc = 2.5e5, σ11 = -3.375e4, -) +) where FT σ1 = ((1/q) + 1) * σc / ((1/q) - q) σ2 = q * σ1 + σc σ22 = q * σ11 + σc - return calculate_mohrs(-σ1, -σ2, -σ11, -σ22) + return _calculate_mohrs(FT, -σ1, -σ2, -σ11, -σ22) end """ MohrsCone{FT}(val::AbstractFloat, args...) -Calculate Mohr's Cone vertices given calculate_mohrs arguments. +Calculate Mohr's Cone vertices given _calculate_mohrs arguments. """ -MohrsCone{FT}(args...) where FT = MohrsCone{FT}(calculate_mohrs(args...)) +MohrsCone{FT}(args...) where FT = MohrsCone{FT}(_calculate_mohrs(FT, args...)) """ update_criteria!(criteria::HiblerYieldCurve, floes) @@ -261,8 +233,9 @@ Inputs: Outputs: None. Updates the criteria's vertices field to update new criteria. """ -function update_criteria!(criteria::HiblerYieldCurve, floes) - criteria.vertices = calculate_hibler( +function update_criteria!(criteria::HiblerYieldCurve{FT}, floes) where FT + criteria.poly = _calculate_hibler( + FT, mean(floes.height), criteria.pstar, criteria.c @@ -302,9 +275,8 @@ function determine_fractures( ) # Determine if floe stresses are in or out of criteria allowable regions update_criteria!(criteria, floes) - critical_poly = make_polygon(GO.tuples(criteria.vertices)) # If stresses are outside of criteria regions, we will fracture the floe - frac_idx = [!GO.coveredby(eigvals(s), critical_poly) for s in floes.stress] + frac_idx = [!GO.coveredby(eigvals(s), criteria.poly) for s in floes.stress] frac_idx[floes.area .< min_floe_area] .= false return range(1, length(floes))[frac_idx] end @@ -312,7 +284,7 @@ end """ deform_floe!( floe, - deformer_coords, + deformer_poly, deforming_forces, ) @@ -331,15 +303,14 @@ Outputs: """ function deform_floe!( floe, - deformer_coords::PolyVec{FT}, + deformer_poly::Polys{FT}, deforming_forces, floe_settings, Δt, rng, ) where FT - poly = make_polygon(floe.coords) - deformer_poly = make_polygon(deformer_coords) - overlap_regions =intersect_polys(poly, deformer_poly) + poly = floe.poly + overlap_regions = intersect_polys(poly, deformer_poly, FT) 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 @@ -350,8 +321,8 @@ function deform_floe!( 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 = make_polygon(translate(deformer_coords, Δx, Δy)) - new_floes = diff_polys(poly, deformer_poly) + deformer_poly = _translate_poly(FT, deformer_poly, Δx, Δy) + new_floes = diff_polys(poly, deformer_poly, FT) 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% @@ -420,8 +391,8 @@ function split_floe( ) if !isempty(pieces) # Intersect voronoi tesselation pieces with floe - floe_poly = make_polygon(rmholes(floe.coords)) - pieces_polys = mapreduce(p -> intersect_polys(make_polygon(p), floe_poly), append!, pieces; init = Vector{Polys{FT}}()) + rmholes!(floe.poly) + pieces_polys = mapreduce(p -> intersect_polys(make_polygon(p), floe.poly, FT), append!, pieces; init = Vector{Polys{FT}}()) # Conserve mass within pieces pieces_areas = [GO.area(p) for p in pieces_polys] total_area = sum(pieces_areas) @@ -519,8 +490,8 @@ function fracture_floes!( deforming_floe_idx = Int(deforming_inter[floeidx]) if deforming_floe_idx <= length(floes) deform_floe!( - LazyRow(floes, frac_idx[i]), - floes.coords[deforming_floe_idx], + get_floe(floes, frac_idx[i]), + floes.poly[deforming_floe_idx], deforming_inter[xforce:yforce], floe_settings, Δt, @@ -531,7 +502,7 @@ function fracture_floes!( end # Split flie into pieces new_floes = split_floe( - LazyRow(floes, frac_idx[i]), + get_floe(floes, frac_idx[i]), rng, fracture_settings, floe_settings, diff --git a/src/physical_processes/process_settings.jl b/src/physical_processes/process_settings.jl index 8453e6a..7c9bfaf 100644 --- a/src/physical_processes/process_settings.jl +++ b/src/physical_processes/process_settings.jl @@ -19,7 +19,7 @@ Settings needed to create floes within the model. """ @kwdef struct FloeSettings{ FT <: AbstractFloat, - GT <: AbstractSubFloePointsGenerator, + GT <: AbstractSubFloePointsGenerator{FT}, } ρi::FT = 920.0 min_floe_area::FT = 1e6 diff --git a/src/physical_processes/ridge_raft.jl b/src/physical_processes/ridge_raft.jl index c9dca04..37a0e83 100644 --- a/src/physical_processes/ridge_raft.jl +++ b/src/physical_processes/ridge_raft.jl @@ -45,7 +45,7 @@ end remove_floe_overlap!( floes, shrink_idx, - grow_floe_coords, + grow_floe_poly, pieces_buffer, max_floe_id, broken, @@ -58,7 +58,7 @@ Removes area/volume of overlap from floe that loses area during ridging/rafting Inputs: floes list of floes shrink_idx index of floe that loses area - grow_floe_coords coordinate of floe/domain that subsumes area + grow_floe_poly polygon of floe/domain that subsumes area pieces_buffer list of new floe pieces caused by breakage of floes max_floe_id maximum floe ID before this ridging/rafting @@ -79,7 +79,7 @@ function remove_floe_overlap!( floes::StructArray{<:Floe{FT}}, shrink_idx, shrink_parent_idx, - grow_floe_coords, + grow_floe_poly, pieces_buffer, max_floe_id, broken, @@ -89,7 +89,7 @@ function remove_floe_overlap!( rng, ) where {FT <: AbstractFloat} # Find new floe shapes and regions - regions = diff_polys(make_polygon(floes.coords[shrink_idx]), make_polygon(grow_floe_coords)) + regions = diff_polys(floes.poly[shrink_idx], grow_floe_poly, FT) total_area = zero(FT) nregions = 0 for (i, region) in enumerate(regions) @@ -113,7 +113,6 @@ function remove_floe_overlap!( # Update existing floes/ghosts regions for region in regions region_area = GO.area(region) - new_coords = find_poly_coords(region)::PolyVec{FT} (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 @@ -122,20 +121,15 @@ function remove_floe_overlap!( (Δx > Δy ? Δy/Δx : Δx/Δy) > floe_settings.min_aspect_ratio ) floe_num += 1 - translate!( # shift region coords to parent floe location - new_coords, - parent_Δx, - parent_Δy, - ) - rmholes!(new_coords) # remove holes in floe - new_poly = make_polygon(new_coords) # parent floe new region poly + new_poly = _translate_poly(FT, region, parent_Δx, parent_Δy)::Polys{FT} + rmholes!(new_poly) # remove holes in floe new_vol = region_area * floes.height[shrink_idx] transfer_vol -= new_vol buffer_length = length(pieces_buffer) # If this is the first region created, replace original floe if floe_num == 1 replace_floe!( # replace parent floe - LazyRow(floes, shrink_parent_idx), + get_floe(floes, shrink_parent_idx), new_poly, new_vol * floe_settings.ρi, floe_settings, @@ -149,16 +143,14 @@ function remove_floe_overlap!( g_Δy = floes.centroid[gidx][2] - parent_centroid[2] # replace ghost floe replace_floe!( - LazyRow(floes, gidx), + get_floe(floes, gidx), new_poly, new_vol * floe_settings.ρi, floe_settings, rng, ) # shift ghost floe - translate!(floes.coords[gidx], g_Δx, g_Δy) - floes.centroid[gidx][1] += g_Δx - floes.centroid[gidx][1] += g_Δy + _translate_floe!(FT, get_floe(floes, gidx), g_Δx, g_Δy) end else # if floe breaks, mark floe and ghosts as broken @@ -179,11 +171,11 @@ function remove_floe_overlap!( else # >1 region, so floe must break and add pieces to buffer push!( pieces_buffer, - deepcopy_floe(LazyRow(floes, shrink_parent_idx)) + deepcopy_floe(get_floe(floes, shrink_parent_idx)) ) buffer_length += 1 replace_floe!( - LazyRow(pieces_buffer, buffer_length), + get_floe(pieces_buffer, buffer_length), new_poly, new_vol * floe_settings.ρi, floe_settings, @@ -298,7 +290,7 @@ function floe_floe_ridge!( floes, lose_mass_idx, lose_parent_idx, - floes.coords[gain_mass_idx], + floes.poly[gain_mass_idx], pieces_buffer, max_floe_id, broken, @@ -319,8 +311,8 @@ function floe_floe_ridge!( if nregions < 1 conserve_momentum_change_floe_shape!( mg, Ig, xg, yg, Δt, - LazyRow(floes, gain_mass_idx), - LazyRow(floes, lose_mass_idx), + get_floe(floes, gain_mass_idx), + get_floe(floes, lose_mass_idx), ) elseif nregions == 1 conserve_momentum_transfer_mass!(floes, @@ -412,7 +404,7 @@ function floe_domain_ridge!( floes, idx, parent_idx, - domain_element.coords, + domain_element.poly, pieces_buffer, max_floe_id, broken, @@ -456,7 +448,7 @@ function floe_domain_ridge!( x_tmp, y_tmp, Δt, - LazyRow(floes, idx), + get_floe(floes, idx), ) end if !broken[idx] @@ -544,7 +536,7 @@ function floe_floe_raft!( floes, lose_mass_idx, lose_parent_idx, - floes.coords[gain_mass_idx], + floes.poly[gain_mass_idx], pieces_buffer, max_floe_id, broken, @@ -565,8 +557,8 @@ function floe_floe_raft!( if nregions == 0 conserve_momentum_change_floe_shape!( mg, Ig, xg, yg, Δt, - LazyRow(floes, gain_mass_idx), - LazyRow(floes, lose_mass_idx), + get_floe(floes, gain_mass_idx), + get_floe(floes, lose_mass_idx), ) elseif nregions == 1 conserve_momentum_transfer_mass!(floes, diff --git a/src/physical_processes/simplification.jl b/src/physical_processes/simplification.jl index 6c0044b..379dff2 100644 --- a/src/physical_processes/simplification.jl +++ b/src/physical_processes/simplification.jl @@ -63,10 +63,10 @@ function smooth_floes!( rng, ) where {FT <: AbstractFloat} for i in eachindex(floes) - if length(floes.coords[i][1]) > simp_settings.max_vertices - poly_list = [simplify_poly(make_polygon(floes.coords[i]), simp_settings.tol)] + if GI.npoint(GI.getexterior(floes.poly[i])) > simp_settings.max_vertices + poly_list = [simplify_poly(floes.poly[i], simp_settings.tol)] if !isempty(topography) - poly_list = diff_polys(make_multipolygon(poly_list), make_multipolygon(topography.coords); fix_multipoly = nothing) + poly_list = diff_polys(make_multipolygon(poly_list), make_multipolygon(topography.poly), FT) end simp_poly = if length(poly_list) == 1 @@ -79,7 +79,7 @@ function smooth_floes!( x_tmp, y_tmp = floes.centroid[i] moment_tmp = floes.moment[i] replace_floe!( - LazyRow(floes, i), + get_floe(floes, i), simp_poly, floes.mass[i], floe_settings, @@ -92,7 +92,7 @@ function smooth_floes!( x_tmp, y_tmp, Δt, - LazyRow(floes, i), + get_floe(floes, i), ) # Mark interactions for fusion for j in eachindex(floes) @@ -106,8 +106,8 @@ function smooth_floes!( floes.status[i].tag = fuse push!(floes.status[i].fuse_idx, j) else - jpoly = make_polygon(floes.coords[j]) - intersect_area = sum(GO.area, intersect_polys(simp_poly, jpoly); init = 0.0) + jpoly = floes.poly[j] + intersect_area = sum(GO.area, intersect_polys(simp_poly, jpoly, FT); 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) @@ -144,21 +144,22 @@ Outputs: Note that the smaller floe's ID is NOT updated! """ function fuse_two_floes!( - keep_floe, + keep_floe::FloeType{FT}, remove_floe, Δt, floe_settings, prefuse_max_floe_id, rng, -) +) where FT # Create new polygon if they fuse - rmholes!(keep_floe.coords) - rmholes!(remove_floe.coords) - poly1 = make_polygon(keep_floe.coords) - poly2 = make_polygon(remove_floe.coords) - new_poly_list = union_polys(poly1, poly2) + rmholes!(keep_floe.poly) + rmholes!(remove_floe.poly) + poly1 = keep_floe.poly + poly2 = remove_floe.poly + new_poly_list = union_polys(poly1, poly2, FT) if length(new_poly_list) == 1 # if they fused, they will make one polygon - new_poly = rmholes(new_poly_list[1]) + new_poly = new_poly_list[1] + rmholes!(new_poly) # mark smaller floe for removal remove_floe.status.tag = remove # record as value will change with replace @@ -242,8 +243,8 @@ function fuse_floes!( keep_idx, remove_idx = floes.area[i] < floes.area[j] ? (j, i) : (i, j) fuse_two_floes!( - LazyRow(floes, keep_idx), - LazyRow(floes, remove_idx), + get_floe(floes, keep_idx), + get_floe(floes, remove_idx), Δt, floe_settings, prefuse_max_floe_id, @@ -293,7 +294,7 @@ function remove_floes!( ) # Dissolve small/thin floes and add mass to ocean dissolve_floe!( - LazyRow(floes, i), + get_floe(floes, i), grid, domain, dissolved, diff --git a/src/physical_processes/update_floe.jl b/src/physical_processes/update_floe.jl index dd5cf39..6b349b6 100644 --- a/src/physical_processes/update_floe.jl +++ b/src/physical_processes/update_floe.jl @@ -31,9 +31,10 @@ function replace_floe!( rng, ) where {FT} # Floe shape + floe.poly = new_poly floe.centroid = collect(GO.centroid(new_poly)) + # TODO: can't replace until we remove coords from floe entirely 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 @@ -47,20 +48,18 @@ function replace_floe!( floe.height; ρi = floe_settings.ρi, ) - floe.angles = GO.angles(make_polygon(floe.coords)) + floe.angles = GO.angles(floe.poly, FT) 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]])) + floe.rmax = calc_max_radius(floe.poly, floe.centroid, FT) # Floe monte carlo points x_subfloe_points, y_subfloe_points, status = generate_subfloe_points( floe_settings.subfloe_point_generator, - floe.coords, - floe.rmax, + floe.poly, + floe.centroid, floe.area, floe.status, rng, ) - translate!(floe.coords, floe.centroid[1], floe.centroid[2]) floe.x_subfloe_points = x_subfloe_points floe.y_subfloe_points = y_subfloe_points # Floe status / identification @@ -211,7 +210,6 @@ function update_new_rotation_conserve!( Δt, ) # Find radius of each polygon to shared midpoint - #x, y = find_shared_edges_midpoint(floe1.coords, floe2.coords) rad1 = sqrt( (floe1.centroid[1] - x)^2 + (floe1.centroid[2] - y)^2 @@ -392,7 +390,7 @@ Inputs: Outputs: Caculates stress on floe at current timestep from interactions """ -function calc_stress!(floe::Union{LazyRow{Floe{FT}}, Floe{FT}}) where {FT} +function calc_stress!(floe::FloeType{FT}) where {FT} # Stress calcultions xi, yi = floe.centroid inters = floe.interactions @@ -423,31 +421,33 @@ Inputs: Outputs: strain 2x2 matrix for floe strain """ -function calc_strain!(floe::Union{LazyRow{Floe{FT}}, Floe{FT}}) where {FT} +function calc_strain!(floe::FloeType{FT}) where {FT} + fill!(floe.strain, zero(FT)) # coordinates of floe centered at centroid - translate!(floe.coords, -floe.centroid[1], -floe.centroid[2]) - fill!(floe.strain, FT(0)) - for i in 1:(length(floe.coords[1]) - 1) - xdiff = floe.coords[1][i + 1][1] - floe.coords[1][i][1] - ydiff = floe.coords[1][i + 1][2] - floe.coords[1][i][2] - rad1 = sqrt(floe.coords[1][i][1]^2 + floe.coords[1][i][2]^2) - θ1 = atan(floe.coords[1][i][2], floe.coords[1][i][1]) - rad2 = sqrt(floe.coords[1][i + 1][1]^2 + floe.coords[1][i + 1][2]^2) - θ2 = atan(floe.coords[1][i + 1][2], floe.coords[1][i + 1][1]) + trans_poly = _translate_poly(FT, floe.poly, -floe.centroid[1], -floe.centroid[2])::Polys{FT} + local x1, y1 + for (i, p2) in enumerate(GI.getpoint(GI.getexterior(trans_poly))) + x2, y2 = GO._tuple_point(p2, FT) + if i == 1 + x1, y1 = x2, y2 + continue + end + xdiff, ydiff = x2 - x1, y2 - y1 + rad1, rad2 = sqrt(x1^2 + y1^2), sqrt(x2^2 + y2^2) + θ1, θ2 = atan(y1, x1), atan(y2, x2) u1 = floe.u - floe.ξ * rad1 * sin(θ1) u2 = floe.u - floe.ξ * rad2 * sin(θ2) v1 = floe.u + floe.ξ * rad1 * cos(θ1) v2 = floe.u + floe.ξ * rad2 * cos(θ2) - udiff = u2 - u1 - vdiff = v2 - v1 + udiff, vdiff = u2 - u1, v2 - v1 floe.strain[1, 1] += udiff * ydiff floe.strain[1, 2] += udiff * xdiff + vdiff * ydiff floe.strain[2, 2] += vdiff * xdiff + x1, y1 = x2, y2 end floe.strain[1, 2] *= FT(0.5) floe.strain[2, 1] = floe.strain[1, 2] floe.strain ./= 2floe.area - translate!(floe.coords, floe.centroid[1], floe.centroid[2]) return end @@ -466,17 +466,17 @@ Output: None. Floe's fields are updated with values. """ function timestep_floe_properties!( - floes, + floes::StructArray{<:Floe{FT}}, tstep, Δt, floe_settings, -) +) where FT Threads.@threads for i in eachindex(floes) cforce = floes.collision_force[i] ctrq = floes.collision_trq[i] # Update stress if floes.num_inters[i] > 0 - calc_stress!(LazyRow(floes, i)) + calc_stress!(get_floe(floes, i)) end # Ensure no extreem values due to model instability if floes.height[i] > floe_settings.max_floe_height @@ -505,18 +505,7 @@ function timestep_floe_properties!( Δα = 1.5Δt*floes.ξ[i] - 0.5Δt*floes.p_dαdt[i] floes.α[i] += Δα - translate!( - floes.coords[i], - -floes.centroid[i][1], - -floes.centroid[i][2], - ) - rotate_radians!(floes.coords[i], Δα) - floes.centroid[i] .+= [Δx, Δy] - translate!( - floes.coords[i], - floes.centroid[i][1], - floes.centroid[i][2], - ) + _move_floe!(FT, get_floe(floes, i), Δx, Δy, Δα) floes.p_dxdt[i] = floes.u[i] floes.p_dydt[i] = floes.v[i] floes.p_dαdt[i] = floes.ξ[i] @@ -556,7 +545,7 @@ function timestep_floe_properties!( floes.p_dξdt[i] = dξdt # Update strain - calc_strain!(LazyRow(floes, i)) + calc_strain!(get_floe(floes, i)) end return end \ No newline at end of file diff --git a/src/physical_processes/welding.jl b/src/physical_processes/welding.jl index 4d8eb36..34efaa3 100644 --- a/src/physical_processes/welding.jl +++ b/src/physical_processes/welding.jl @@ -131,7 +131,7 @@ function timestep_welding!( ) ) # Find intersection area - inter_area = FT(sum(GO.area, intersect_polys(make_polygon(floes.coords[i]), make_polygon(floes.coords[j])); init = 0.0)) + inter_area = FT(sum(GO.area, intersect_polys(floes.poly[i], floes.poly[j], FT); init = 0.0)) # Probability two floes will weld weld_prob = weld_settings.welding_coeff * (inter_area / floes.area[i]) @@ -159,8 +159,8 @@ function timestep_welding!( new_area > weld_settings.max_weld_area && break # Weld floe i and j, replacing floe i with welded floe fuse_two_floes!( - LazyRow(floes, i), - LazyRow(floes, j), + get_floe(floes, i), + get_floe(floes, j), Δt, floe_settings, max_floe_id, diff --git a/src/simulation_components/domain_and_grid.jl b/src/simulation_components/domain_and_grid.jl index f2a4495..4abba07 100644 --- a/src/simulation_components/domain_and_grid.jl +++ b/src/simulation_components/domain_and_grid.jl @@ -181,6 +181,7 @@ A sub-type of AbstractBoundary that allows a floe to pass out of the domain edge without any effects on the floe. """ struct OpenBoundary{D, FT}<:AbstractBoundary{D, FT} + poly::Polys{FT} coords::PolyVec{FT} val::FT end @@ -222,6 +223,7 @@ function OpenBoundary{D, FT}( ) where {D <: AbstractDirection, FT <: AbstractFloat} val, coords = boundary_coords(grid, D) OpenBoundary{D, FT}( + make_polygon(convert(PolyVec{FT}, coords)), coords, val, ) @@ -235,6 +237,7 @@ the opposite side of the domain when it crosses the boundary, bringing the floe back into the domain. """ struct PeriodicBoundary{D, FT}<:AbstractBoundary{D, FT} + poly::Polys{FT} coords::PolyVec{FT} val::FT end @@ -276,6 +279,7 @@ function PeriodicBoundary{D, FT}( ) where {D <: AbstractDirection, FT <: AbstractFloat} val, coords = boundary_coords(grid, D) PeriodicBoundary{D, FT}( + make_polygon(convert(PolyVec{FT}, coords)), coords, val, ) @@ -289,6 +293,7 @@ having the floe collide with the boundary. The boundary acts as an immovable, unbreakable ice floe in the collision. """ struct CollisionBoundary{D, FT}<:AbstractBoundary{D, FT} + poly::Polys{FT} coords::PolyVec{FT} val::FT end @@ -332,6 +337,7 @@ function CollisionBoundary{D, FT}( ) where {D <: AbstractDirection, FT <: AbstractFloat} val, coords = boundary_coords(grid, D) CollisionBoundary{D, FT}( + make_polygon(convert(PolyVec{FT}, coords)), coords, val, ) @@ -353,6 +359,7 @@ the the combination of u and v velocities. The same, but opposite is true for th and south walls. """ mutable struct MovingBoundary{D, FT}<:AbstractBoundary{D, FT} + poly::Polys{FT} coords::PolyVec{FT} val::FT u::FT @@ -402,6 +409,7 @@ function MovingBoundary{D, FT}( ) where {D <: AbstractDirection, FT <: AbstractFloat} val, coords = boundary_coords(grid, D) MovingBoundary{D, FT}( + make_polygon(convert(PolyVec{FT}, coords)), coords, val, u, @@ -430,12 +438,14 @@ in that they will not move or break due to floe interactions, but they will affect floes. """ struct TopographyElement{FT}<:AbstractDomainElement{FT} + poly::Polys{FT} coords::PolyVec{FT} centroid::Vector{FT} rmax::FT function TopographyElement{FT}( - coords, + poly, + coords, # TODO: Can remove when not used anymore in other parts of the code centroid, rmax, ) where {FT <: AbstractFloat} @@ -443,9 +453,8 @@ struct TopographyElement{FT}<:AbstractDomainElement{FT} throw(ArgumentError("Topography element maximum radius must be \ positive and non-zero.")) end - topo_poly = make_polygon(rmholes(coords)) - topo_poly = GO.ClosedRing()(topo_poly) - new{FT}(GI.coordinates(topo_poly), centroid, rmax) + poly = GO.ClosedRing()(poly) + new{FT}(poly, find_poly_coords(poly), centroid, rmax) end end @@ -477,24 +486,14 @@ Output: Topographic element of abstract float type FT """ function TopographyElement{FT}(poly::Polys) where {FT <: AbstractFloat} - topo = rmholes(poly) - cx, cy = GO.centroid(topo) - coords = find_poly_coords(topo) - # Move coordinates to be centered at origin to calculate maximum radius - translate!( - coords, - -cx, - -cy, - ) - rmax = sqrt(maximum([sum(c.^2) for c in coords[1]])) - translate!( - coords, - cx, - cy, - ) + rmholes!(poly) + centroid = collect(GO.centroid(poly)) # TODO: Remove collect once type is changed + coords = find_poly_coords(poly) + rmax = calc_max_radius(poly, centroid, FT) return TopographyElement{FT}( + poly, coords, - [cx, cy], + centroid, rmax, ) end @@ -508,11 +507,8 @@ Inputs: Output: Topographic element of abstract float type FT """ -function TopographyElement{FT}( - coords::PolyVec, -) where {FT <: AbstractFloat} - return TopographyElement{FT}(make_polygon(convert(PolyVec{Float64}, coords))) -end +TopographyElement{FT}(coords::PolyVec) where {FT <: AbstractFloat} = + TopographyElement{FT}(make_polygon(convert(PolyVec{FT}, coords))) """ initialize_topography_field(args...) @@ -657,7 +653,7 @@ function get_domain_element(domain, idx) return domain.west else topo_idx = -(idx + 4) - return LazyRow(domain.topography, topo_idx) + return get_floe(domain.topography, topo_idx) end end diff --git a/src/simulation_components/floe.jl b/src/simulation_components/floe.jl index 0e6fb02..e28dd40 100644 --- a/src/simulation_components/floe.jl +++ b/src/simulation_components/floe.jl @@ -98,6 +98,7 @@ Singular sea ice floe with fields describing current state. """ @kwdef mutable struct Floe{FT<:AbstractFloat} # Physical Properties ------------------------------------------------- + poly::Polys{FT} # polygon that represents the floe's shape centroid::Vector{FT} # center of mass of floe (might not be in floe!) coords::PolyVec{FT} # floe coordinates height::FT # floe height (m) @@ -149,6 +150,8 @@ Singular sea ice floe with fields describing current state. p_dαdt::FT = 0.0 # previous timestep angular-velocity end +const FloeType{FT} = Union{LazyRow{Floe{FT}}, Floe{FT}} where FT + """ Floe(::Type{FT}, args...; kwargs...) @@ -220,7 +223,8 @@ function Floe{FT}( rng = Xoshiro(), kwargs... ) where {FT <: AbstractFloat} - floe = rmholes(poly) + floe = GO.tuples(poly, FT) + rmholes!(floe) # Floe physical properties centroid = collect(GO.centroid(floe)) height = clamp( @@ -231,27 +235,25 @@ function Floe{FT}( 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(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]])) + angles = GO.angles(floe, FT) + rmax = calc_max_radius(floe, centroid, FT) status = Status() # Generate Monte Carlo Points x_subfloe_points, y_subfloe_points, status = generate_subfloe_points( floe_settings.subfloe_point_generator, - coords, - rmax, + floe, + centroid, area_tot, status, rng, ) - translate!(coords, centroid[1], centroid[2]) # Generate Stress History stress_history = StressCircularBuffer{FT}(floe_settings.nhistory) fill!(stress_history, zeros(FT, 2, 2)) return Floe{FT}(; + poly = GO.tuples(floe, FT), centroid = centroid, coords = coords, height = height, @@ -294,16 +296,18 @@ Output: forcings start at 0 and floe's status is "active" as long as monte carlo points were able to be generated. """ -Floe{FT}( +function Floe{FT}( coords::PolyVec, hmean, Δh; floe_settings = FloeSettings(), rng = Xoshiro(), kwargs..., -) where {FT <: AbstractFloat} = - Floe{FT}( # Polygon convert is needed since LibGEOS only takes Float64 - make_polygon(convert(PolyVec{Float64}, valid_polyvec!(rmholes(coords)))), +) where {FT <: AbstractFloat} + valid_polyvec!(coords) + rmholes!(coords) + return Floe{FT}( + make_polygon(coords), hmean, Δh; floe_settings = floe_settings, @@ -311,6 +315,8 @@ Floe{FT}( kwargs..., ) +end + """ poly_to_floes!( ::Type{FT}, @@ -425,12 +431,13 @@ function initialize_floe_field( Δh; floe_settings = FloeSettings(min_floe_area = 0.0), rng = Xoshiro(), + supress_warnings = false, ) where {FT <: AbstractFloat} floe_arr = StructArray{Floe{FT}}(undef, 0) floe_polys = [make_polygon(valid_polyvec!(c)) for c in coords] # Remove overlaps with topography if !isempty(domain.topography) - floe_polys = diff_polys(make_multipolygon(floe_polys), make_multipolygon(domain.topography.coords)) + floe_polys = diff_polys(make_multipolygon(floe_polys), make_multipolygon(domain.topography.poly), FT) end # Turn polygons into floes for p in floe_polys @@ -452,7 +459,7 @@ function initialize_floe_field( 4 * (domain.east.val - domain.west.val) * (domain.north.val - domain.south.val) / 1e4 ) - if any(floe_arr.area .< min_floe_area) + if any(floe_arr.area .< min_floe_area) && !supress_warnings @warn "Some user input floe areas are less than the suggested minimum \ floe area." end @@ -460,7 +467,7 @@ function initialize_floe_field( if !all( domain.west.val .< first.(floe_arr.centroid) .< domain.east.val) && !all(domain.south.val .< last.(floe_arr.centroid) .< domain.north.val - ) + ) && !supress_warnings @warn "Some floe centroids are out of the domain." end # Initialize floe IDs @@ -509,8 +516,8 @@ function generate_voronoi_coords( min_to_warn::Int; max_tries::Int = 10, ) where {FT <: AbstractFloat} - xpoints = Vector{FT}() - ypoints = Vector{FT}() + xpoints = Vector{Float64}() + ypoints = Vector{Float64}() 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 @@ -518,8 +525,8 @@ function generate_voronoi_coords( current_points = 0 tries = 0 while current_points < desired_points && tries <= max_tries - x = rand(rng, FT, npoints) - y = rand(rng, FT, npoints) + x = rand(rng, npoints) + y = rand(rng, npoints) # 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]), @@ -552,13 +559,8 @@ function generate_voronoi_coords( # Scale and translate voronoi coordinates tcoords = Vector{PolyVec{FT}}(undef, length(tess_cells)) for i in eachindex(tess_cells) - perturb_vec = [ - (-1)^rand(rng, 0:1) * rand(rng, FT)*1e-10, - (-1)^rand(rng, 0:1) * rand(rng, FT)*1e-10, - ] tcoords[i] = [valid_ringvec!([ - Vector(c) .* scale_fac .+ - trans_vec #.+ perturb_vec + Vector(c) .* scale_fac .+ trans_vec for c in tess_cells[i] ])] end @@ -615,17 +617,17 @@ function initialize_floe_field( domain, hmean, Δh; - floe_bounds = rect_coords(domain.west.val, domain.east.val, domain.south.val, domain.north.val), - floe_settings = FloeSettings(min_floe_area = 0.0), + floe_bounds = _make_bounding_box_polygon(FT, domain.west.val, domain.east.val, domain.south.val, domain.north.val), + floe_settings = FloeSettings(FT, min_floe_area = 0), rng = Xoshiro(), ) where {FT <: AbstractFloat} floe_arr = StructArray{Floe{FT}}(undef, 0) nfloes_added = 0 # Availible space in domain - 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) + domain_poly = _make_bounding_box_polygon(FT, domain.west.val, domain.east.val, domain.south.val, domain.north.val) + open_water = intersect_polys(floe_bounds, domain_poly, FT) if !isempty(domain.topography) - open_water = diff_polys(make_multipolygon(open_water), make_multipolygon(domain.topography.coords)) + open_water = diff_polys(make_multipolygon(open_water), make_multipolygon(domain.topography.poly), FT) end open_water_mp = make_multipolygon(open_water) (bounds_xmin, bounds_xmax), (bounds_ymin, bounds_ymax) = GI.extent(open_water_mp) @@ -637,7 +639,6 @@ function initialize_floe_field( Ly = bounds_ymax - bounds_ymin rowlen = Ly / nrows collen = Lx / ncols - # Loop over cells for j in range(1, ncols) for i in range(1, nrows) @@ -649,10 +650,10 @@ function initialize_floe_field( ymin = bounds_ymin + rowlen * (i - 1) trans_vec = [xmin, ymin] # Open water in cell - cell_init = make_polygon(rect_coords(xmin, xmin + collen, ymin, ymin + rowlen)) - open_cell = intersect_polys(cell_init, open_water_mp) + cell_init = _make_bounding_box_polygon(FT, xmin, xmin + collen, ymin, ymin + rowlen) + open_cell = intersect_polys(cell_init, open_water_mp, FT) open_cell_mpoly = make_multipolygon(open_cell) - open_coords = [GI.coordinates(c) for c in open_cell] + open_coords = [find_poly_coords(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) diff --git a/test/runtests.jl b/test/runtests.jl index cb2542c..875198a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,4 @@ -using DataStructures, GeometryBasics, JLD2, Logging, NCDatasets, Random, SplitApplyCombine, +using DataStructures, JLD2, Logging, NCDatasets, Random, SplitApplyCombine, Statistics, StructArrays, Subzero, VoronoiCells import GeometryOps as GO import GeometryOps.GeoInterface as GI diff --git a/test/test_floe.jl b/test/test_floe.jl index 471ce89..1f1877f 100644 --- a/test/test_floe.jl +++ b/test/test_floe.jl @@ -1,12 +1,22 @@ @testset "Floe" begin + # Test Setup FT = Float64 - # test generating monte carlo points + Lx = 8e4 + Ly = Lx + Δgrid = 1e4 + hmean = 0.5 + Δh = 0.01 + rng = Xoshiro(1) + fs_small_min_area = FloeSettings(min_floe_area = 55) + + # Test generating monte carlo points file = jldopen("inputs/floe_shapes.jld2", "r") floe_coords = file["floe_vertices"][1:end] close(file) + poly1 = Subzero.make_polygon(Subzero.valid_polyvec!(floe_coords[1])) centroid1 = GO.centroid(poly1) - area = GO.area(poly1) + area1 = GO.area(poly1) # Test InteractionFields enum interactions = range(1, 7)' @@ -15,87 +25,51 @@ @test interactions[overlap] == 7 # Test with coords inputs - floe_from_coords = Floe( - floe_coords[1], - 0.5, - 0.01, - u = 0.2, - rng = Xoshiro(1), - ) + floe_from_coords = Floe(floe_coords[1], hmean, Δh; u = 0.2, rng = rng) @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 == collect(centroid1) - @test floe_from_coords.area == area + @test floe_from_coords.area == area1 @test floe_from_coords.status.tag == Subzero.active # Test with polygon input - floe_from_poly = Floe( - poly1, - 0.5, - 0.01, - v = -0.2, - rng = Xoshiro(1), - ) + floe_from_poly = Floe(poly1, hmean, Δh; v = -0.2, rng = Xoshiro(1)) @test typeof(floe_from_poly) <: Floe @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 == collect(centroid1) - @test floe_from_poly.area == area + @test floe_from_poly.area == area1 @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]]] + rect_coords = [[[0.0, 0.0], [0.0, 5.0], [10.0, 5.0], [10.0, 0.0], [0.0, 0.0]]] + rect_poly = Subzero.make_polygon(rect_coords) 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]]] + c_hole_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]], + ] + c_hole_poly = Subzero.make_polygon(c_hole_coords) rmax_cpoly = 2sqrt(5^2 + 5^2) # Test polygon with no holes floe_arr = StructArray{Floe{FT}}(undef, 0) - n_new = Subzero.poly_to_floes!( - FT, - floe_arr, - Subzero.make_polygon(rect_poly), - 0.25, - 0.01, - rmax_rect; - ) + n_new = Subzero.poly_to_floes!(FT, floe_arr, rect_poly, hmean, Δh, rmax_rect) @test n_new == 1 && length(floe_arr) == 1 @test !Subzero.hashole(floe_arr.coords[1]) # Test with polygon below minimum floe area - n_new = Subzero.poly_to_floes!( - FT, - floe_arr, - Subzero.make_polygon(rect_poly), - 0.25, - 0.01, - rmax_rect; - floe_settings = FloeSettings(min_floe_area = 55), - ) + n_new = Subzero.poly_to_floes!(FT, floe_arr, rect_poly, hmean, Δh, rmax_rect; floe_settings = fs_small_min_area) @test n_new == 0 && length(floe_arr) == 1 # Test with polygon with a hole that is split into 3 polyons - n_new = Subzero.poly_to_floes!( - FT, - floe_arr, - Subzero.make_polygon(c_poly_hole), - 0.25, - 0.01, - rmax_cpoly, - ) + n_new = Subzero.poly_to_floes!(FT, floe_arr, c_hole_poly, hmean, Δh, rmax_cpoly) @test n_new == 3 && length(floe_arr) == 4 @test !any(Subzero.hashole.(floe_arr.coords)) # Test initialize_floe_field from coord list - grid = RegRectilinearGrid( - (-8e4, 8e4), - (-8e4, 8e4), - 1e4, - 1e4, - ) + grid = RegRectilinearGrid((-Lx, Lx), (-Ly, Ly), Δgrid, Δgrid) nbound = CollisionBoundary(North, grid) sbound = CollisionBoundary(South, grid) ebound = CollisionBoundary(East, grid) @@ -126,12 +100,7 @@ @test all(floe_arr.id .== range(1, nfloes)) # From file with small domain -> floes outside of domain - small_grid = RegRectilinearGrid( - (-5e4, 5e4), - (-5e4, 5e4), - 1e4, - 1e4, - ) + small_grid = RegRectilinearGrid((-Lx/2, Lx/2), (-Ly/2, Ly), Δgrid, Δgrid) small_domain_no_topo = Domain( CollisionBoundary(North, small_grid), CollisionBoundary(South, small_grid), @@ -142,8 +111,8 @@ FT, floe_coords, small_domain_no_topo, - 0.5, - 0.1; + hmean, + Δh; floe_settings = FloeSettings(min_floe_area = 1e5), )) @test typeof(floe_arr) <: StructArray{<:Floe} @@ -154,8 +123,8 @@ FT, floe_coords, domain_with_topo, - 0.5, - 0.1; + hmean, + Δh; floe_settings = FloeSettings(min_floe_area = 10), rng = Xoshiro(0) ) diff --git a/test/test_floe_utils.jl b/test/test_floe_utils.jl index c9097b2..4dd65c9 100644 --- a/test/test_floe_utils.jl +++ b/test/test_floe_utils.jl @@ -40,15 +40,13 @@ @test Subzero.hashole(poly_hole1) @test Subzero.hashole(poly_hole2) - # Test removing holes from polygons and multipolygons - @test [ext] == Subzero.rmholes([ext]) - @test [ext] == Subzero.rmholes([ext, hole1]) + # Test removing holes from polygons copy_holes = [ext, hole1] + poly_copy_holes = Subzero.make_polygon(copy_holes) Subzero.rmholes!(copy_holes) @test copy_holes == [ext] - @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) + Subzero.rmholes!(poly_copy_holes) + @test GI.nhole(poly_copy_holes) == 0 # Test translating coordinates and polygons @test Subzero.translate([ext], 0.0, 0.0) == [ext] @@ -72,57 +70,6 @@ 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 - c1 = [ - [9.75e4, 7e4], - [9.75e4, 5e4], - [9.75e4, 5e4], - [10.05e4, 5e4], - [10.05e4, 7e4], - ] - c1_new = Subzero.orient_coords(c1) - @test c1_new == [ - [97500.0, 50000.0], - [97500.0, 70000.0], - [100500.0, 70000.0], - [100500.0, 50000.0], - [97500.0, 50000.0], - ] - c2 = [ - [6.5e4, 6.5e4], - [8.5e4, 6.5e4], - [8.5e4, 4.5e4], - [8.5e4, 4.5e4], - [6.5e4, 4.5e4], - ] - c2_new = Subzero.orient_coords(c2) - @test c2_new == [ - [6.5e4, 4.5e4], - [6.5e4, 6.5e4], - [8.5e4, 6.5e4], - [8.5e4, 4.5e4], - [6.5e4, 4.5e4], - ] - - # ------------------------- 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]]] - @test issetequal( - Subzero.intersect_lines(l1, l2), - Set([(3.0, 0.0)]), - ) - l1 = [[[0., -1], [1, 1], [2, -1], [3, 1]]] - l2 = [[[0., 0], [1, 0], [3, 0]]] - @test issetequal( - Subzero.intersect_lines(l1, l2), - Set([(0.5, -0.0), (1.5, 0), (2.5, -0.0)]), - ) - l2 = [[[10., 10]]] - @test issetequal( - Subzero.intersect_lines(l1, l2), - Set{Tuple{Float64, Float64}}(), - ) - # ------------------------- Test finding shared points --------------------- two_shared_v = [ [[[0.0, 0.0], [0.0, 20.0], [20.0, 20.0], [20.0, 0.0], [0.0, 0.0]]], @@ -130,7 +77,7 @@ ] @test Subzero.which_vertices_match_points( two_shared_v[1][1], - two_shared_v[2], + Subzero.make_polygon(two_shared_v[2]), ) == [1, 2] @test Subzero.which_points_on_edges( @@ -144,7 +91,7 @@ ] @test Subzero.which_vertices_match_points( three_shared_v[1][1], - three_shared_v[2], + Subzero.make_polygon(three_shared_v[2]), ) == [3, 4, 5] @test Subzero.which_points_on_edges( @@ -164,7 +111,7 @@ ] @test Subzero.which_vertices_match_points( four_shared_v[1][1], - four_shared_v[2], + Subzero.make_polygon(four_shared_v[2]), ) == [1, 2, 5, 6] offset_shared_v = [ @@ -186,7 +133,7 @@ ] @test Subzero.which_vertices_match_points( triange_shared_v[1][1], - triange_shared_v[2], + Subzero.make_polygon(triange_shared_v[2]), ) == [1, 2, 3] @test Subzero.which_points_on_edges( diff --git a/test/test_model.jl b/test/test_model.jl index fa06aea..befa77a 100644 --- a/test/test_model.jl +++ b/test/test_model.jl @@ -230,6 +230,7 @@ @test topo2.rmax == sqrt(0.5) # Basic constructor topo3 = TopographyElement( + Subzero.make_polygon(coords), coords, [0.5, 0.5], sqrt(0.5), @@ -237,6 +238,7 @@ @test topo3.coords == coords # check when radius is less than or equal to 0 @test_throws ArgumentError TopographyElement( + Subzero.make_polygon(coords), coords, [0.5, 0.5], -sqrt(0.5), @@ -306,11 +308,13 @@ Subzero.PeriodicBoundary(West, g), ) # domain with north < south + p_placeholder = GI.Polygon([[(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (0.0, 0.0)]]) @test_throws ArgumentError Subzero.Domain( b1, Subzero.OpenBoundary( South, - Subzero.PolyVec{Float64}(undef, 0), + p_placeholder, + GI.coordinates(p_placeholder), 6e5, ), b2, @@ -323,7 +327,8 @@ b2, Subzero.OpenBoundary( West, - Subzero.PolyVec{Float64}(undef, 0), + p_placeholder, + GI.coordinates(p_placeholder), 6e5, ), ) diff --git a/test/test_physical_processes/test_collisions.jl b/test/test_physical_processes/test_collisions.jl index c18dad7..a2e7567 100644 --- a/test/test_physical_processes/test_collisions.jl +++ b/test/test_physical_processes/test_collisions.jl @@ -1,300 +1,193 @@ @testset "Collisions" begin FT = Float64 + Δt = 10 + Δh = 0.0 + consts = Constants() + collision_settings = CollisionSettings() + spinlock = Threads.SpinLock() + supress_warnings = true + # Grid setup + Lx = 1e5 + Ly = Lx + grid = RegRectilinearGrid((-Lx, Lx), (-Ly, Ly), 1e4, 1e4) + + # Boundary setup + p_n_bound = PeriodicBoundary(North, grid) + p_s_bound = PeriodicBoundary(South, grid) + p_e_bound = PeriodicBoundary(East, grid) + p_w_bound = PeriodicBoundary(West, grid) + + c_n_bound = CollisionBoundary(North, grid) + c_s_bound = CollisionBoundary(South, grid) + c_e_bound = CollisionBoundary(East, grid) + c_w_bound = CollisionBoundary(West, grid) + + o_n_bound = OpenBoundary(North, grid) + o_s_bound = OpenBoundary(South, grid) + o_e_bound = OpenBoundary(East, grid) + o_w_bound = OpenBoundary(West, grid) + + topos = initialize_topography_field(FT, [[[[1e4, 0.0], [0.0, 1e4], [1e4, 2e4], [2e4, 1e4], [1e4, 0.0]]]]) + # topos = StructArray([TopographyElement([[[1e4, 0.0], [0.0, 1e4], [1e4, 2e4], [2e4, 1e4], [1e4, 0.0]]])]) + + topo_domain = Domain(p_n_bound, p_s_bound, c_e_bound, o_w_bound, topos) + collision_domain = Domain(c_n_bound, c_s_bound, c_e_bound, c_w_bound) + open_domain = Domain(o_n_bound, o_s_bound, o_e_bound, o_w_bound) + ew_periodic_domain = Domain(o_n_bound, o_s_bound, p_e_bound, p_w_bound) + ns_periodic_domain = Domain(p_n_bound, p_s_bound, o_e_bound, o_w_bound) + double_periodic_domain = Domain(p_n_bound, p_s_bound, p_e_bound, p_w_bound) + @testset "Floe-Floe Interactions" begin - Δt = 10 - max_overlap = 0.55 hmean = 0.25 - Δh = 0.0 - tri = Floe( - [[[0.0, 0.0], [1e4, 3e4], [2e4, 0], [0.0, 0.0]]], - hmean, - Δh, - ) - tri.u = 0.1 - rect = Floe( - [[ - [0.0, 2.5e4], - [0.0, 2.9e4], - [2e4, 2.9e4], - [2e4, 2.5e4], - [0.0, 2.5e4], - ]], - hmean, - Δh, - ) - rect.v = -0.1 - cfloe = Floe( - [[ - [0.5e4, 2.7e4], - [0.5e4, 3.5e4], - [1.5e4, 3.5e4], - [1.5e4, 2.7e4], - [1.25e4, 2.7e4], - [1.25e4, 3e4], - [1e4, 3e4], - [1e4, 2.7e4], - [0.5e4, 2.7e4], - ]], - hmean, - Δh, - ) - cfloe.u = 0.3 - consts = Constants() - # Triange tip intersected with a rectangle - Subzero.floe_floe_interaction!(tri, 1, rect, 2, consts, Δt, max_overlap) + max_overlap = 0.55 + # Floe coordinates + tri_coord = [[[0.0, 0.0], [1e4, 3e4], [2e4, 0], [0.0, 0.0]]] + corner_rect_coord = [[[0.0, 2.5e4], [0.0, 2.9e4], [2e4, 2.9e4], [2e4, 2.5e4], [0.0, 2.5e4]]] + small_shift_corner_rect_coord = Subzero.translate(corner_rect_coord, 0.5e4, 0.0) + big_shift_corner_rect_coord = Subzero.translate(corner_rect_coord, 1.9999999e4, 0.0) + middle_rect_coord = [[ [1.8e4, 2.7e4], [1.8e4, 2.8e4], [2.1e4, 2.8e4], [2.1e4, 2.7e4], [1.8e4, 2.7e4]]] + cshape_coord = [[[0.5e4, 2.7e4], [0.5e4, 3.5e4], [1.5e4, 3.5e4], [1.5e4, 2.7e4], [1.25e4, 2.7e4], [1.25e4, 3e4], [1e4, 3e4], [1e4, 2.7e4], [0.5e4, 2.7e4]]] + + # Test interactions between triange tip intersecting with a rectangle + tri, corner_rect = Floe(tri_coord, hmean, Δh), Floe(corner_rect_coord, hmean, Δh) + tri.u, corner_rect.v = 0.1, -0.1 + Subzero.floe_floe_interaction!(tri, 1, corner_rect, 2, consts, Δt, max_overlap) @test isapprox(tri.interactions[1, xforce], -64613382.47, atol = 1e-2) @test isapprox(tri.interactions[1, yforce], -521498991.51, atol = 1e-2) @test isapprox(tri.interactions[1, xpoint], 10000.00, atol = 1e-2) @test isapprox(tri.interactions[1, ypoint], 26555.55, atol = 1e-2) @test isapprox(tri.interactions[1, overlap], 8000000, atol = 1e-2) - @test tri.status.tag != Subzero.fuse && rect.status.tag != Subzero.fuse + @test tri.status.tag != Subzero.fuse && corner_rect.status.tag != Subzero.fuse @test isempty(tri.status.fuse_idx) Subzero.calc_torque!(tri) @test isapprox(tri.interactions[1, torque], 1069710443203.99, atol = 1e-2) - tri.interactions = zeros(0, 7) + # Sideways C intersected with rectangle so there are two areas of overlap - Subzero.floe_floe_interaction!(cfloe, 1, rect, 2, consts, Δt, max_overlap) - @test isapprox(cfloe.interactions[1, xforce], -163013665.41, atol = 1e-2) - @test isapprox(cfloe.interactions[2, xforce], -81506832.70, atol = 1e-2) - @test isapprox(cfloe.interactions[1, yforce], 804819565.60, atol = 1e-2) - @test isapprox(cfloe.interactions[2, yforce], 402409782.80, atol = 1e-2) - @test isapprox(cfloe.interactions[1, xpoint], 7500.00, atol = 1e-2) - @test isapprox(cfloe.interactions[2, xpoint], 13750.00, atol = 1e-2) - @test isapprox(cfloe.interactions[1, ypoint], 28000.00, atol = 1e-2) - @test isapprox(cfloe.interactions[2, ypoint], 28000.00, atol = 1e-2) - @test isapprox(cfloe.interactions[1, overlap], 10000000, atol = 1e-2) - @test isapprox(cfloe.interactions[2, overlap], 5000000, atol = 1e-2) - @test tri.status.tag != Subzero.fuse && rect.status.tag != Subzero.fuse - Subzero.calc_torque!(cfloe) - @test isapprox(cfloe.interactions[1, torque], -2439177121266.03, atol = 1e-2) - @test isapprox(cfloe.interactions[2, torque], 1295472581868.05, atol = 1e-2) - cfloe.interactions = zeros(0, 7) + cshape_floe, corner_rect = Floe(cshape_coord, hmean, Δh), Floe(corner_rect_coord, hmean, Δh) + cshape_floe.u, corner_rect.v = 0.3, -0.1 + Subzero.floe_floe_interaction!(cshape_floe, 1, corner_rect, 2, consts, Δt, max_overlap) + @test isapprox(cshape_floe.interactions[1, xforce], -163013665.41, atol = 1e-2) + @test isapprox(cshape_floe.interactions[2, xforce], -81506832.70, atol = 1e-2) + @test isapprox(cshape_floe.interactions[1, yforce], 804819565.60, atol = 1e-2) + @test isapprox(cshape_floe.interactions[2, yforce], 402409782.80, atol = 1e-2) + @test isapprox(cshape_floe.interactions[1, xpoint], 7500.00, atol = 1e-2) + @test isapprox(cshape_floe.interactions[2, xpoint], 13750.00, atol = 1e-2) + @test isapprox(cshape_floe.interactions[1, ypoint], 28000.00, atol = 1e-2) + @test isapprox(cshape_floe.interactions[2, ypoint], 28000.00, atol = 1e-2) + @test isapprox(cshape_floe.interactions[1, overlap], 10000000, atol = 1e-2) + @test isapprox(cshape_floe.interactions[2, overlap], 5000000, atol = 1e-2) + @test tri.status.tag != Subzero.fuse && corner_rect.status.tag != Subzero.fuse + Subzero.calc_torque!(cshape_floe) + @test isapprox(cshape_floe.interactions[1, torque], -2439177121266.03, atol = 1e-2) + @test isapprox(cshape_floe.interactions[2, torque], 1295472581868.05, atol = 1e-2) + # Floes overlapping more than 55% - rectangle and shifted rectangle - shift_rect = deepcopy(rect) - shift_rect.coords = Subzero.translate(shift_rect.coords, 0.5e4, 0.0) - Subzero.floe_floe_interaction!(rect, 1, shift_rect, 2, consts, Δt, max_overlap) - @test rect.status.tag == Subzero.fuse - @test rect.status.fuse_idx == [2] - @test isempty(rect.interactions) - empty!(rect.status.fuse_idx) - rect.status.tag == Subzero.active - # Floe 2 (small rectangle) is overlapping with rectangle by more than 55% - small_rect = Floe( - [[ - [1.8e4, 2.7e4], - [1.8e4, 2.8e4], - [2.1e4, 2.8e4], - [2.1e4, 2.7e4], - [1.8e4, 2.7e4], - ]], - hmean, - Δh, - ) - Subzero.floe_floe_interaction!( - rect, - 1, - small_rect, - 2, - consts, - Δt, - max_overlap, - ) - @test rect.status.tag == Subzero.fuse - @test rect.status.fuse_idx == [2] + corner_rect, small_shift_corner_rect = Floe(corner_rect_coord, hmean, Δh), Floe(small_shift_corner_rect_coord, hmean, Δh) + corner_rect.v, small_shift_corner_rect.v = -0.1, -0.1 + Subzero.floe_floe_interaction!(corner_rect, 1, small_shift_corner_rect, 2, consts, Δt, max_overlap) + @test corner_rect.status.tag == Subzero.fuse + @test corner_rect.status.fuse_idx == [2] + @test isempty(corner_rect.interactions) + + # Two floes (original and original with a small shift) overlapping by more than 55% + corner_rect, middle_rect = Floe(corner_rect_coord, hmean, Δh), Floe(middle_rect_coord, hmean, Δh) + corner_rect.v = -0.1 + Subzero.floe_floe_interaction!(corner_rect, 1, middle_rect, 2, consts, Δt, max_overlap) + @test corner_rect.status.tag == Subzero.fuse + @test corner_rect.status.fuse_idx == [2] - # Overlapping barely floes - such a small overlap that forces are not calculated - shift_rect = deepcopy(rect) - shift_rect.coords = Subzero.translate(shift_rect.coords, 1.9999999e4, 0.0) - Subzero.floe_floe_interaction!(shift_rect, 1, rect, 2, consts, Δt, max_overlap) - @test isempty(shift_rect.interactions) + # Two floes (original and original with a big shift) have such a small overlap that forces are not calculated + corner_rect, big_shift_corner_rect = Floe(corner_rect_coord, hmean, Δh), Floe(big_shift_corner_rect_coord, hmean, Δh) + corner_rect.v, big_shift_corner_rect.v = -0.1, -0.1 + Subzero.floe_floe_interaction!(big_shift_corner_rect, 1, corner_rect, 2, consts, Δt, max_overlap) + @test isempty(big_shift_corner_rect.interactions) end - @testset "Floe Boundary Interactions" begin - Δt = 10 - max_overlap = 0.75 - Lx = 1e5 - Ly = Lx hmean = 0.25 - Δh = 0.0 - grid = RegRectilinearGrid( - (-Lx, Lx), - (-Ly, Ly), - 1e4, - 1e4, - ) - nboundary = PeriodicBoundary(North, grid) - sboundary = PeriodicBoundary(South, grid) - eboundary = CollisionBoundary(East, grid) - wboundary = OpenBoundary(West, grid) - topo_elem = TopographyElement( - [[[1e4, 0.0], [0.0, 1e4], [1e4, 2e4], [2e4, 1e4], [1e4, 0.0]]], - ) - domain = Domain( - nboundary, - sboundary, - eboundary, - wboundary, - StructArray([topo_elem]), - ) + max_overlap = 0.75 + ## Floe coordinates + # Floe overlapping with northern periodic boundary + north_coords = [[[5e4, 9.75e4], [5e4, 10.05e4], [7e4, 10.05e4], [7e4, 9.75e4], [5e4, 9.75e4]]] # Diagonal floe barely overlaping with eastern collision boundary - efloe_small = Floe( - [[ - [9.5e4, 0.0], - [9e4, 0.5e4], - [10e4, 2.5e4], - [10.05e4, 2e4], - [9.5e4, 0.0], - ]], - hmean, - Δh, - ) - efloe_small.u = 0.5 - efloe_small.v = 0.25 + east_small_coords = [[[9.5e4, 0.0], [9e4, 0.5e4], [10e4, 2.5e4], [10.05e4, 2e4], [9.5e4, 0.0]]] # Floe overlapping with eastern collision boundary by more than 75% to trigger overlap condition - efloe_large = Floe( - [[ - [9e4, -7e4], - [9e4, -5e4], - [1.4e5, -5e4], - [1.4e5, -7e4], - [9e4, -7e4], - ]], - hmean, - Δh, - ) - efloe_large.u = 0.1 - efloe_large.v = -0.35 - # Floe overlapping with boundary with more than one region - cfloe = Floe( - [[ - [9.5e4, 7e4], - [9.5e4, 9e4], - [1.05e5, 9e4], - [1.05e5, 8.5e4], - [9.9e4, 8.5e4], - [9.9e4, 8e4], - [1.05e5, 8e4], - [1.05e5, 7e4], - [9.5e4, 7e4], - ]], - hmean, - Δh, - ) - cfloe.v = -0.1 + east_large_coords = [[[9e4, -7e4], [9e4, -5e4], [1.4e5, -5e4], [1.4e5, -7e4], [9e4, -7e4]]] # Floe overlapping with western open boundary - wfloe = Floe( - [[ - [-9.75e4, 7e4], - [-9.75e4, 5e4], - [-10.05e4, 5e4], - [-10.05e4, 7e4], - [-9.75e4, 7e4], - ]], - hmean, - Δh, - ) - # Floe overlapping with northern periodic boundary - nfloe = Floe( - [[ - [5e4, 9.75e4], - [5e4, 10.05e4], - [7e4, 10.05e4], - [7e4, 9.75e4], - [5e4, 9.75e4], - ]], - hmean, - Δh, - ) + west_coords = [[[-9.75e4, 7e4], [-9.75e4, 5e4], [-10.05e4, 5e4], [-10.05e4, 7e4], [-9.75e4, 7e4]]] + # Floe overlapping with boundary with more than one region + cshape_coord = [[[9.5e4, 7e4], [9.5e4, 9e4], [1.05e5, 9e4], [1.05e5, 8.5e4], [9.9e4, 8.5e4], [9.9e4, 8e4], [1.05e5, 8e4], [1.05e5, 7e4], [9.5e4, 7e4]]] # Floe overlapping with topography element - tfloe = Floe( - [[ - [-0.5e4, 0.0], - [-0.5e4, 0.75e4], - [0.5e4, 0.75e4], - [0.5e4, 0.0], - [-0.5e4, 0.0], - ]], - hmean, - Δh, - ) - efloe_large.u = -0.4 - efloe_large.v = 0.2 - - consts = Constants() - - # Test floe overlapping slightly with collision boundary - Subzero.floe_domain_interaction!(efloe_small, domain, consts, Δt, max_overlap) - @test efloe_small.interactions[1, floeidx] == -3 - @test isapprox(efloe_small.interactions[1, xforce], -311304795.629, atol = 1e-3) - @test isapprox(efloe_small.interactions[1, yforce], -23618874.648, atol = 1e-3) - @test isapprox(efloe_small.interactions[1, overlap], 1704545.454, atol = 1e-3) - @test isapprox(efloe_small.interactions[1, xpoint], 100166.666, atol = 1e-3) - @test isapprox(efloe_small.interactions[1, ypoint], 21060.606, atol = 1e-3) - - Subzero.floe_domain_interaction!(cfloe, domain, consts, Δt, max_overlap) - @test cfloe.interactions[1, floeidx] == -3 - @test cfloe.interactions[2, floeidx] == -3 - @test isapprox(cfloe.interactions[1, xforce], -2876118708.17, atol = 1e-2) - @test isapprox(cfloe.interactions[2, xforce], -5752237416.35, atol = 1e-2) - @test isapprox(cfloe.interactions[1, yforce], 575223741.63, atol = 1e-2) - @test isapprox(cfloe.interactions[2, yforce], 1150447483.27, atol = 1e-2) - @test isapprox(cfloe.interactions[1, xpoint], 102500, atol = 1e-2) - @test isapprox(cfloe.interactions[2, xpoint], 102500, atol = 1e-2) - @test isapprox(cfloe.interactions[1, ypoint], 87500, atol = 1e-2) - @test isapprox(cfloe.interactions[2, ypoint], 75000, atol = 1e-2) - @test isapprox(cfloe.interactions[1, overlap], 25000000, atol = 1e-2) - @test isapprox(cfloe.interactions[2, overlap], 50000000, atol = 1e-2) + topo_overlap_coords = [[[-0.5e4, 0.0], [-0.5e4, 0.75e4], [0.5e4, 0.75e4], [0.5e4, 0.0], [-0.5e4, 0.0]]] + # Floe in a corner hitting more than one wall at once + corner_coords = [[[9.5e4, 7e4], [9e4, 7.5e4], [10e4, 1.05e5], [10.05e4, 9.5e4], [9.5e4, 7e4]]] + + # Test floe overlapping slightly with collision boundary (one region) + east_small_floe = Floe(east_small_coords, hmean, Δh) + east_small_floe.u, east_small_floe.v = 0.5, 0.25 + Subzero.floe_domain_interaction!(east_small_floe, topo_domain, consts, Δt, max_overlap) + @test east_small_floe.interactions[1, floeidx] == -3 + @test isapprox(east_small_floe.interactions[1, xforce], -311304795.629, atol = 1e-3) + @test isapprox(east_small_floe.interactions[1, yforce], -23618874.648, atol = 1e-3) + @test isapprox(east_small_floe.interactions[1, overlap], 1704545.454, atol = 1e-3) + @test isapprox(east_small_floe.interactions[1, xpoint], 100166.666, atol = 1e-3) + @test isapprox(east_small_floe.interactions[1, ypoint], 21060.606, atol = 1e-3) + + # Test floe overlapping slightly with collision boundary (two regions) + cshape_floe = Floe(cshape_coord, hmean, Δh) + cshape_floe.v = -0.1 + Subzero.floe_domain_interaction!(cshape_floe, topo_domain, consts, Δt, max_overlap) + @test cshape_floe.interactions[1, floeidx] == -3 + @test cshape_floe.interactions[2, floeidx] == -3 + @test isapprox(cshape_floe.interactions[1, xforce], -2876118708.17, atol = 1e-2) + @test isapprox(cshape_floe.interactions[2, xforce], -5752237416.35, atol = 1e-2) + @test isapprox(cshape_floe.interactions[1, yforce], 575223741.63, atol = 1e-2) + @test isapprox(cshape_floe.interactions[2, yforce], 1150447483.27, atol = 1e-2) + @test isapprox(cshape_floe.interactions[1, xpoint], 102500, atol = 1e-2) + @test isapprox(cshape_floe.interactions[2, xpoint], 102500, atol = 1e-2) + @test isapprox(cshape_floe.interactions[1, ypoint], 87500, atol = 1e-2) + @test isapprox(cshape_floe.interactions[2, ypoint], 75000, atol = 1e-2) + @test isapprox(cshape_floe.interactions[1, overlap], 25000000, atol = 1e-2) + @test isapprox(cshape_floe.interactions[2, overlap], 50000000, atol = 1e-2) # Test floe overlapping >75% with collision boundary - Subzero.floe_domain_interaction!(efloe_large, domain, consts, Δt, max_overlap) - @test isempty(efloe_large.interactions) - @test efloe_large.status.tag == Subzero.remove + east_large_floe = Floe(east_large_coords, hmean, Δh) + east_large_floe.u, east_large_floe.v = -0.4, 0.2 + Subzero.floe_domain_interaction!(east_large_floe, topo_domain, consts, Δt, max_overlap) + @test isempty(east_large_floe.interactions) + @test east_large_floe.status.tag == Subzero.remove + # Test floe overlapping >75% with collision boundary -> but max overlap is now 1 - Subzero.floe_domain_interaction!(efloe_large, domain, consts, Δt, 1.0) - @test !isempty(efloe_large.interactions) - @test efloe_large.num_inters > 0 + east_large_floe = Floe(east_large_coords, hmean, Δh) + east_large_floe.u, east_large_floe.v = -0.4, 0.2 + Subzero.floe_domain_interaction!(east_large_floe, topo_domain, consts, Δt, 1.0) + @test !isempty(east_large_floe.interactions) + @test east_large_floe.num_inters > 0 + # Test floe passing through open boundary is killed - Subzero.floe_domain_interaction!(wfloe, domain, consts, Δt, max_overlap) - @test wfloe.status.tag == Subzero.remove + west_floe = Floe(west_coords, hmean, Δh) + Subzero.floe_domain_interaction!(west_floe, topo_domain, consts, Δt, max_overlap) + @test west_floe.status.tag == Subzero.remove + # Test floes not not interact with periodic boundary - nfloe_copy = deepcopy(nfloe) - Subzero.floe_domain_interaction!(nfloe, domain, consts, Δt, max_overlap) - @test nfloe_copy.status.tag == nfloe.status.tag && nfloe_copy.interactions == nfloe.interactions + north_floe = Floe(north_coords, hmean, Δh) + Subzero.floe_domain_interaction!(north_floe, topo_domain, consts, Δt, max_overlap) + @test north_floe.status.tag == north_floe.status.tag && north_floe.interactions == north_floe.interactions + # Test floe overlapping with topography -> different from Subzero since topography are now boundaries - Subzero.floe_domain_interaction!(tfloe, domain, consts, Δt, max_overlap) - @test tfloe.interactions[1, xforce] < 0 - @test tfloe.interactions[1, yforce] < 0 - @test tfloe.interactions[1, floeidx] == -5 + topo_overlap_floe = Floe(topo_overlap_coords, hmean, Δh) + Subzero.floe_domain_interaction!(topo_overlap_floe, topo_domain, consts, Δt, max_overlap) + @test topo_overlap_floe.interactions[1, xforce] < 0 + @test topo_overlap_floe.interactions[1, yforce] < 0 + @test topo_overlap_floe.interactions[1, floeidx] == -5 - # Test floe hitting more than one wall at once -> different from Subzero - collision_domain = Domain( - CollisionBoundary(North, grid), - CollisionBoundary(South, grid), - CollisionBoundary(East, grid), - CollisionBoundary(West, grid), - ) - corner_floe = Floe( - [[ - [9.5e4, 7e4], - [9e4, 7.5e4], - [10e4, 1.05e5], - [10.05e4, 9.5e4], - [9.5e4, 7e4], - ]], - hmean, - Δh, - ) - Subzero.floe_domain_interaction!( - corner_floe, - collision_domain, - consts, - Δt, - max_overlap, - ) + # Test floe in corner hitting more than one wall at once + corner_floe = Floe(corner_coords, hmean, Δh) + Subzero.floe_domain_interaction!(corner_floe, collision_domain, consts, Δt, max_overlap) @test all(corner_floe.interactions[:, xforce] .<= 0) @test all(corner_floe.interactions[:, yforce] .<= 0) - # Test compression boundaries movement + + # Test compression boundaries movement - TODO: Move these to different file nc_boundary = MovingBoundary(North, grid, 0.0, -0.1) nc_coords = deepcopy(nc_boundary.coords) sc_boundary = MovingBoundary(South, grid, 0.0, 0.1) @@ -320,19 +213,8 @@ end @testset "Add Ghosts" begin - FT = Float64 - Lx = 1e5 - grid = RegRectilinearGrid( - (-Lx, Lx), - (-Lx, Lx), - 1e4, - 1e4, - ) - nboundary = PeriodicBoundary(North, grid) - sboundary = PeriodicBoundary(South, grid) - eboundary = PeriodicBoundary(East, grid) - wboundary = PeriodicBoundary(West, grid) - + hmean = 0.5 + ## Floe coordinates # corner floe - overlaps with north and east boundary coords1 = [[[9.9e4, 9.9e4], [9.9e4, 1.02e5], [1.02e5, 1.02e5], [1.02e5, 9.9e4], [9.9e4, 9.9e4]]] # overlaps with western boundary @@ -341,147 +223,100 @@ coords3 = [[[-2e4, 9.5e4], [-2e4, 1.1e5], [-1e4, 1.1e5], [-1e4, 9.5e4], [-2e4, 9.5e4]]] # doesn't overlap with any boundary coords4 = [[[0.0, 0.0], [0.0, 2e4], [2e4, 2e4], [2e4, 0.0], [0.0, 0.0]]] - - floe_arr = StructArray( - [Floe(c, 0.5, 0.0) for c in [coords1, coords2, coords3, coords4]], - ) - for i in eachindex(floe_arr) - floe_arr.id[i] = Float64(i) - end - - nonperiodic_domain = Domain( - OpenBoundary(North, grid), - OpenBoundary(South, grid), - OpenBoundary(East, grid), - OpenBoundary(West, grid), - ) - ew_periodic_domain = Domain( - OpenBoundary(North, grid), - OpenBoundary(South, grid), - PeriodicBoundary(East, grid), - PeriodicBoundary(West, grid), - ) - ns_periodic_domain = Domain( - PeriodicBoundary(North, grid), - PeriodicBoundary(South, grid), - OpenBoundary(East, grid), - OpenBoundary(West, grid), - ) - double_periodic_domain = Domain( - PeriodicBoundary(North, grid), - PeriodicBoundary(South, grid), - PeriodicBoundary(East, grid), - PeriodicBoundary(West, grid), - ) + # List of above coordinates + coord_list = [coords1, coords2, coords3, coords4] # Make sure nothing is added with non-periodic domain - new_floe_arr = deepcopy(floe_arr) - add_ghosts!(new_floe_arr, nonperiodic_domain) - @test new_floe_arr.coords == floe_arr.coords + floe_arr = initialize_floe_field(FT, coord_list, open_domain, hmean, Δh; supress_warnings) + add_ghosts!(floe_arr, open_domain) + @test floe_arr.coords == coord_list # Add ghost floes in east-west direction - new_floe_arr = deepcopy(floe_arr) - add_ghosts!(new_floe_arr, ew_periodic_domain) - @test -1e5 < new_floe_arr[1].centroid[1] < 1e5 - @test -1e5 < new_floe_arr[2].centroid[2] < 1e5 - @test new_floe_arr.coords[1] == Subzero.translate(floe_arr.coords[1], -2e5, 0.0) - @test new_floe_arr.coords[2:4] == floe_arr.coords[2:4] - @test new_floe_arr.coords[5] == floe_arr.coords[1] - @test new_floe_arr.coords[6] == Subzero.translate(floe_arr.coords[2], 2e5, 0.0) - @test new_floe_arr.id == [1, 2, 3, 4, 1, 2] - @test new_floe_arr.ghost_id == [0, 0, 0, 0, 1, 1] - @test new_floe_arr.ghosts[1] == [5] - @test new_floe_arr.ghosts[2] == [6] - @test new_floe_arr.ghosts[3:6] == [[], [], [], []] + floe_arr = initialize_floe_field(FT, coord_list, ew_periodic_domain, hmean, Δh; supress_warnings) + add_ghosts!(floe_arr, ew_periodic_domain) + @test -1e5 < floe_arr[1].centroid[1] < 1e5 + @test -1e5 < floe_arr[2].centroid[2] < 1e5 + @test floe_arr.coords[1] == Subzero.translate(coords1, -2e5, 0.0) + @test floe_arr.coords[2:4] == coord_list[2:4] + @test floe_arr.coords[5] == coords1 + @test floe_arr.coords[6] == Subzero.translate(coords2, 2e5, 0.0) + @test floe_arr.id == [1, 2, 3, 4, 1, 2] + @test floe_arr.ghost_id == [0, 0, 0, 0, 1, 1] + @test floe_arr.ghosts[1] == [5] + @test floe_arr.ghosts[2] == [6] + @test floe_arr.ghosts[3:6] == [[], [], [], []] # Add ghost floes in the north-south direction - new_floe_arr = deepcopy(floe_arr) - add_ghosts!(new_floe_arr, ns_periodic_domain) - @test -1e5 < new_floe_arr[1].centroid[2] < 1e5 - @test -1e5 < new_floe_arr[3].centroid[2] < 1e5 - @test new_floe_arr.coords[1] == Subzero.translate(floe_arr.coords[1], 0.0, -2e5) - @test new_floe_arr.coords[3] == Subzero.translate(floe_arr.coords[3], 0.0, -2e5) - @test new_floe_arr.coords[[2, 4]] == floe_arr.coords[[2, 4]] - @test new_floe_arr.coords[5] == floe_arr.coords[1] - @test new_floe_arr.coords[6] == floe_arr.coords[3] - @test new_floe_arr.id == [1, 2, 3, 4, 1, 3] - @test new_floe_arr.ghost_id == [0, 0, 0, 0, 1, 1] - @test new_floe_arr.ghosts[1] == [5] - @test new_floe_arr.ghosts[3] == [6] - @test new_floe_arr.ghosts[[2; 4:6]] == [[], [], [], []] + floe_arr = initialize_floe_field(FT, coord_list, ns_periodic_domain, hmean, Δh; supress_warnings) + add_ghosts!(floe_arr, ns_periodic_domain) + @test -1e5 < floe_arr[1].centroid[2] < 1e5 + @test -1e5 < floe_arr[3].centroid[2] < 1e5 + @test floe_arr.coords[1] == Subzero.translate(coords1, 0.0, -2e5) + @test floe_arr.coords[3] == Subzero.translate(coords3, 0.0, -2e5) + @test floe_arr.coords[[2, 4]] == coord_list[[2, 4]] + @test floe_arr.coords[5] == coords1 + @test floe_arr.coords[6] == coords3 + @test floe_arr.id == [1, 2, 3, 4, 1, 3] + @test floe_arr.ghost_id == [0, 0, 0, 0, 1, 1] + @test floe_arr.ghosts[1] == [5] + @test floe_arr.ghosts[3] == [6] + @test floe_arr.ghosts[[2; 4:6]] == [[], [], [], []] # Add ghosts in both east-west and north-south directions - new_floe_arr = deepcopy(floe_arr) - add_ghosts!(new_floe_arr, double_periodic_domain) - @test -1e5 < new_floe_arr.centroid[1][1] < 1e5 - @test -1e5 < new_floe_arr.centroid[1][2] < 1e5 - @test new_floe_arr.coords[1] == Subzero.translate(floe_arr.coords[1], -2e5, -2e5) - @test new_floe_arr.coords[3] == Subzero.translate(floe_arr.coords[3], 0.0, -2e5) - @test new_floe_arr.coords[[2, 4]] == floe_arr.coords[[2, 4]] - @test new_floe_arr.coords[5] == floe_arr.coords[1] - @test new_floe_arr.coords[6] == Subzero.translate(floe_arr.coords[2], 2e5, 0.0) - @test new_floe_arr.coords[7] == Subzero.translate(floe_arr.coords[1], 0.0, -2e5) - @test new_floe_arr.coords[8] == Subzero.translate(floe_arr.coords[1], -2e5, 0.0) - @test new_floe_arr.coords[9] == floe_arr.coords[3] - @test new_floe_arr.id == [1, 2, 3, 4, 1, 2, 1, 1, 3] - @test new_floe_arr.ghost_id == [0, 0, 0, 0, 1, 1, 2, 3, 1] - @test new_floe_arr.ghosts[1] == [5, 7, 8] - @test new_floe_arr.ghosts[2] == [6] - @test new_floe_arr.ghosts[3] == [9] - @test new_floe_arr.ghosts[4:9] == [[], [], [], [], [], []] - + floe_arr = initialize_floe_field(FT, coord_list, double_periodic_domain, hmean, Δh; supress_warnings) + add_ghosts!(floe_arr, double_periodic_domain) + @test -1e5 < floe_arr.centroid[1][1] < 1e5 + @test -1e5 < floe_arr.centroid[1][2] < 1e5 + @test floe_arr.coords[1] == Subzero.translate(coords1, -2e5, -2e5) + @test floe_arr.coords[3] == Subzero.translate(coords3, 0.0, -2e5) + @test floe_arr.coords[[2, 4]] == coord_list[[2, 4]] + @test floe_arr.coords[5] == coords1 + @test floe_arr.coords[6] == Subzero.translate(coords2, 2e5, 0.0) + @test floe_arr.coords[7] == Subzero.translate(coords1, 0.0, -2e5) + @test floe_arr.coords[8] == Subzero.translate(coords1, -2e5, 0.0) + @test floe_arr.coords[9] == coords3 + @test floe_arr.id == [1, 2, 3, 4, 1, 2, 1, 1, 3] + @test floe_arr.ghost_id == [0, 0, 0, 0, 1, 1, 2, 3, 1] + @test floe_arr.ghosts[1] == [5, 7, 8] + @test floe_arr.ghosts[2] == [6] + @test floe_arr.ghosts[3] == [9] + @test floe_arr.ghosts[4:9] == [[], [], [], [], [], []] end @testset "Ghost Collisions" begin - Δt = 10 - Lx = 1e5 - Ly = 1e5 - collision_settings = CollisionSettings() - grid = RegRectilinearGrid( - (-Lx, Lx), - (-Lx, Lx), - 1e4, - 1e4, - ) - double_periodic_domain = Domain( - PeriodicBoundary(North, grid), - PeriodicBoundary(South, grid), - PeriodicBoundary(East, grid), - PeriodicBoundary(West, grid), - ) - # Parent-parent collison (parents are touching) - coords1 = splitdims([Lx/2 Lx/2 3*Lx/4 3*Lx/4 Lx+10000 Lx+10000; Ly/2 Ly+10000 Ly+10000 3*Ly/4 3*Ly/4 Ly/2]) - th = 0:pi/50:2*pi - r = Ly/4+1000 - coords2 = invert([r * cos.(th) .+ (Lx-1), r * sin.(th) .+ (Ly-1)]) + hmean = 0.5 + ## Floe coordinates + # L-shaped polygon extending beyond top and right domain edges + lshape_coords = [splitdims([Lx/2 Lx/2 3*Lx/4 3*Lx/4 Lx+10000 Lx+10000; Ly/2 Ly+10000 Ly+10000 3*Ly/4 3*Ly/4 Ly/2])] + # Oval floe touching two edges of lshape_coords + th, r = 0:pi/50:2*pi, Ly/4+1000 + oval_coords = [invert([r * cos.(th) .+ (Lx-1), r * sin.(th) .+ (Ly-1)])] + # Tall rectangle in top right corner of domain + tall_rect_coords = [splitdims(vcat([5*Lx/8 5*Lx/8 3*Lx/4 3*Lx/4].+1000, [3*Ly/4 5*Ly/4 5*Ly/4 3*Ly/4]))] + # Long rectangle in bottom left corner of domain (ghost hits tall rectangle) + long_rect_coords = [splitdims(vcat(-[5*Lx/4 5*Lx/4 3*Lx/4-1000 3*Lx/4-1000], -[7*Lx/8 3*Lx/4-1000 3*Lx/4-1000 7*Lx/8]))] + # Shifted tall rectangle to ghost position in bottom right corner + shifted_down_tall_rect_coords = Subzero.translate(tall_rect_coords, 0.0, -2Ly) + # Shifted tall rectangle to ghost postion in top left corner + shifted_left_tall_rect_coords = Subzero.translate(tall_rect_coords, -2Lx, 0.0) + # Shifted long rectangle to ghost postion in bottom right corner + shifted_right_long_rect_coords = Subzero.translate(long_rect_coords, 2Lx, 0.0) + # Shifted long rectangle to ghost postion in top left corner + shifted_up_long_rect_coords = Subzero.translate(long_rect_coords, 0.0, 1.615Ly) + # Small rectangle in the corner that has 3 ghosts in all other corners + small_corner_rect_coords = [[[-1.1e5, -1.1e5], [-1.1e5, -9.5e4], [-9.5e4, -9.5e4], [-9.5e4, -1.1e5], [-1.1e5, -1.1e5]]] + # triangle in the middle of the domain with no ghosts - touches 3/4 corners + large_tri_coords = [[[-1e5, -1e5], [-1e5, 1e5], [1e5, -1e5], [-1e5, -1e5]]] + # rectangle along south boundary, ghost along north boundary + south_bound_rect_coords = [[[-9.8e4, -1.1e5], [-9.8e4, -9.5e4], [9.8e4, -9.5e4], [9.8e4, -1.1e5], [-9.8e4, -1.1e5]]] + + # Parent-parent collison (parents are touching) + floe_arr = initialize_floe_field(FT, [lshape_coords, oval_coords], double_periodic_domain, hmean, Δh; supress_warnings) + Subzero.timestep_collisions!(floe_arr, 2, double_periodic_domain, consts, Δt, collision_settings, spinlock) + xforce_vals, yforce_vals = abs(floe_arr[1].collision_force[1]), abs(floe_arr[1].collision_force[2]) + f1_torque, f2_torque = floe_arr[1].collision_trq, floe_arr[2].collision_trq - floe_arr = StructArray(Floe([c], 0.5, 0.0) for c in [coords1, coords2]) - for i in eachindex(floe_arr) - floe_arr.id[i] = Float64(i) - end - spinlock = Threads.SpinLock() - Subzero.timestep_collisions!( - floe_arr, - 2, - double_periodic_domain, - Subzero.Constants(), - Δt, - collision_settings, - spinlock, - ) - xforce_vals = abs(floe_arr[1].collision_force[1]) - yforce_vals = abs(floe_arr[1].collision_force[2]) - f1_torque = floe_arr[1].collision_trq - f2_torque = floe_arr[2].collision_trq add_ghosts!(floe_arr, double_periodic_domain) - Subzero.timestep_collisions!( - floe_arr, - 2, - double_periodic_domain, - Subzero.Constants(), - Δt, - collision_settings, - spinlock, - ) + Subzero.timestep_collisions!(floe_arr, 2, double_periodic_domain, consts, Δt, collision_settings, spinlock) # 1 and 2 are the "parent" floes - floe 1 and floe 2 interact @test xforce_vals == abs(floe_arr[1].collision_force[1]) == abs(floe_arr[2].collision_force[1]) @test yforce_vals == abs(floe_arr[2].collision_force[2]) == abs(floe_arr[2].collision_force[2]) @@ -493,53 +328,19 @@ @test floe_arr[3:end].collision_trq == zeros(6) # Ghost-Ghost collision (parents aren't touching, only ghosts touch) - coords1 = splitdims(vcat([5*Lx/8 5*Lx/8 3*Lx/4 3*Lx/4].+1000, [3*Ly/4 5*Ly/4 5*Ly/4 3*Ly/4])) - coords2 = splitdims(vcat(-[5*Lx/4 5*Lx/4 3*Lx/4-1000 3*Lx/4-1000], -[7*Lx/8 3*Lx/4-1000 3*Lx/4-1000 7*Lx/8])) - floe_arr = StructArray(Floe([c], 0.5, 0.0) for c in [coords1, coords2]) - for i in eachindex(floe_arr) - floe_arr.id[i] = i - end - trans_arr = StructArray([ - Floe( - Subzero.translate([coords1], - 0.0, -2Ly), - 0.5, - 0.0, - ), - Floe( - Subzero.translate([coords2], 2Lx, 0.0), - 0.5, - 0.0, - ), - ]) - for i in eachindex(trans_arr) - trans_arr.id[i] = i - end - Subzero.timestep_collisions!( - trans_arr, - 2, - double_periodic_domain, - Subzero.Constants(), - Δt, - collision_settings, - spinlock, - ) + tall_rect_coords = [splitdims(vcat([5*Lx/8 5*Lx/8 3*Lx/4 3*Lx/4].+1000, [3*Ly/4 5*Ly/4 5*Ly/4 3*Ly/4]))] + long_rect_coords = [splitdims(vcat(-[5*Lx/4 5*Lx/4 3*Lx/4-1000 3*Lx/4-1000], -[7*Lx/8 3*Lx/4-1000 3*Lx/4-1000 7*Lx/8]))] + floe_arr = initialize_floe_field(FT, [tall_rect_coords, long_rect_coords], double_periodic_domain, hmean, Δh; supress_warnings) + trans_arr = initialize_floe_field(FT, [shifted_down_tall_rect_coords, shifted_right_long_rect_coords], double_periodic_domain, hmean, Δh; supress_warnings) + + Subzero.timestep_collisions!(trans_arr, 2, double_periodic_domain, consts, Δt, collision_settings, spinlock) xforce_vals = abs(trans_arr[1].collision_force[1]) yforce_vals = abs(trans_arr[1].collision_force[2]) f1_torque = trans_arr[1].collision_trq f2_torque = trans_arr[2].collision_trq add_ghosts!(floe_arr, double_periodic_domain) - Subzero.timestep_collisions!( - floe_arr, - 2, - double_periodic_domain, - Subzero.Constants(), - Δt, - collision_settings, - spinlock, - ) - # floes 1 and 2 are the parents - floe 4 is floe 1's ghost and floe 3 is - # floe 2's ghost - floe 3 and 4 collide + Subzero.timestep_collisions!( floe_arr, 2, double_periodic_domain, consts, Δt, collision_settings, spinlock) + # floes 1 and 2 are the parents - floe 4 is floe 1's ghost and floe 3 is floe 2's ghost - floe 3 and 4 collide @test repeat([xforce_vals], 2) == abs.(first.(floe_arr.collision_force[1:2])) @test repeat([yforce_vals], 2) == abs.(last.(floe_arr.collision_force[1:2])) @test f1_torque == floe_arr[1].collision_trq @@ -549,40 +350,15 @@ @test floe_arr[2].interactions[:, [1:5; 7]] == floe_arr[3].interactions[:, [1:5; 7]] # Parent-Ghost Collision - coords1 = splitdims(vcat([5*Lx/8 5*Lx/8 3*Lx/4 3*Lx/4].+1000, [3*Ly/4 5*Ly/4 5*Ly/4 3*Ly/4])) - coords2 = splitdims(vcat(-[5*Lx/4 5*Lx/4 3*Lx/4-1000 3*Lx/4-1000], [7*Lx/8 3*Lx/4-1000 3*Lx/4-1000 7*Lx/8])) - floe_arr = StructArray(Floe([c], 0.5, 0.0) for c in [coords1, coords2]) - for i in eachindex(floe_arr) - floe_arr.id[i] = Float64(i) - end - trans_arr = StructArray([Floe(Subzero.translate([coords1], -2Lx, 0.0), 0.5, 0.0), - Floe([coords2], 0.5, 0.0)]) - for i in eachindex(trans_arr) - trans_arr.id[i] = Float64(i) - end - Subzero.timestep_collisions!( - trans_arr, - 2, - double_periodic_domain, - Subzero.Constants(), - Δt, - collision_settings, - spinlock, - ) + floe_arr = initialize_floe_field(FT, [tall_rect_coords, shifted_up_long_rect_coords], double_periodic_domain, hmean, Δh; supress_warnings) + trans_arr = initialize_floe_field(FT, [shifted_left_tall_rect_coords, shifted_up_long_rect_coords], double_periodic_domain, hmean, Δh; supress_warnings) + Subzero.timestep_collisions!(trans_arr, 2, double_periodic_domain, consts, Δt, collision_settings, spinlock) xforce_vals = abs(trans_arr[1].collision_force[1]) yforce_vals = abs(trans_arr[1].collision_force[2]) f1_torque = trans_arr[1].collision_trq f2_torque = trans_arr[2].collision_trq add_ghosts!(floe_arr, double_periodic_domain) - Subzero.timestep_collisions!( - floe_arr, - 2, - double_periodic_domain, - Subzero.Constants(), - Δt, - collision_settings, - spinlock, - ) + Subzero.timestep_collisions!(floe_arr, 2, double_periodic_domain, consts, Δt, collision_settings, spinlock) # Floe 1's ghost if floe 4 and floe 2's ghost is floe 3 and floe 3 and floe 1 interact @test repeat([xforce_vals], 2) == abs.(first.(floe_arr.collision_force[1:2])) @test repeat([yforce_vals], 2) == abs.(last.(floe_arr.collision_force[1:2])) @@ -592,77 +368,19 @@ @test isempty(floe_arr[4].interactions) # Parent and ghosts hitting the same floe - # small rectangle in the corner that has 3 ghosts in all other corners - small_rect = Floe( - [[ - [-1.1e5, -1.1e5], - [-1.1e5, -9.5e4], - [-9.5e4, -9.5e4], - [-9.5e4, -1.1e5], - [-1.1e5, -1.1e5], - ]], - 0.5, - 0.0, - ) - # triangle in the middle of the domain with no ghosts - touches 3/4 corners - large_tri = Floe( - [[ - [-1e5, -1e5], - [-1e5, 1e5], - [1e5, -1e5], - [-1e5, -1e5], - ]], - 0.5, - 0.0, - ) - # rectangle along south boundary, ghost along north boundary - bound_rect = Floe( - [[ - [-9.8e4, -1.1e5], - [-9.8e4, -9.5e4], - [9.8e4, -9.5e4], - [9.8e4, -1.1e5], - [-9.8e4, -1.1e5], - ]], - 0.5, - 0.0, - ) - floe_arr = StructArray([small_rect, large_tri]) - for i in eachindex(floe_arr) - floe_arr.id[i] = Float64(i) - end + floe_arr = initialize_floe_field(FT, [small_corner_rect_coords, large_tri_coords], double_periodic_domain, hmean, Δh; supress_warnings) add_ghosts!(floe_arr, double_periodic_domain) @test length(floe_arr) == 5 - Subzero.timestep_collisions!( - floe_arr, - 2, - double_periodic_domain, - Subzero.Constants(), - Δt, - collision_settings, - spinlock, - ) + Subzero.timestep_collisions!(floe_arr, 2, double_periodic_domain, consts, Δt, collision_settings, spinlock) @test size(floe_arr[1].interactions)[1] == 3 @test size(floe_arr[2].interactions)[1] == 3 @test floe_arr[1].interactions[1, Subzero.xforce] != floe_arr[1].interactions[2, Subzero.xforce] && floe_arr[1].interactions[1, Subzero.xforce] != floe_arr[1].interactions[3,Subzero.xforce] @test floe_arr[1].interactions[1, Subzero.yforce] != floe_arr[1].interactions[2, Subzero.yforce] && floe_arr[1].interactions[1, Subzero.yforce] != floe_arr[1].interactions[3, Subzero.yforce] - empty!(small_rect.ghosts) - floe_arr = StructArray([small_rect, bound_rect]) - for i in eachindex(floe_arr) - floe_arr.id[i] = Float64(i) - end + floe_arr = initialize_floe_field(FT, [small_corner_rect_coords, south_bound_rect_coords], double_periodic_domain, hmean, Δh; supress_warnings) add_ghosts!(floe_arr, double_periodic_domain) @test length(floe_arr) == 6 - Subzero.timestep_collisions!( - floe_arr, - 2, - double_periodic_domain, - Subzero.Constants(), - Δt, - collision_settings, - spinlock, - ) + Subzero.timestep_collisions!(floe_arr, 2, double_periodic_domain, consts, Δt, collision_settings, spinlock) @test size(floe_arr[1].interactions)[1] == 2 @test size(floe_arr[2].interactions)[1] == 2 @test floe_arr[1].interactions[1, Subzero.xpoint] != floe_arr[1].interactions[2, Subzero.xpoint] diff --git a/test/test_physical_processes/test_coupling.jl b/test/test_physical_processes/test_coupling.jl index 4cac171..e0ee7c8 100644 --- a/test/test_physical_processes/test_coupling.jl +++ b/test/test_physical_processes/test_coupling.jl @@ -14,13 +14,14 @@ -centroid1[2], ) origin_poly = Subzero.make_polygon(origin_coords) + origin_centroid = GO.centroid(origin_poly) 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 = GO.area(poly1) mc_x, mc_y, status = Subzero.generate_subfloe_points( MonteCarloPointsGenerator(), - origin_coords, - rmax, + origin_poly, + origin_centroid, area, Subzero.Status(), Xoshiro(1) @@ -35,8 +36,8 @@ # Test that random number generator is working mc_x2, mc_y2, status2 = Subzero.generate_subfloe_points( MonteCarloPointsGenerator(), - origin_coords, - rmax, + origin_poly, + origin_centroid, area, Subzero.Status(), Xoshiro(1) @@ -47,8 +48,8 @@ mc_x3, mc_y3, status3 = Subzero.generate_subfloe_points( MonteCarloPointsGenerator{Float32}(), - origin_coords, - rmax, + GO.tuples(origin_poly, Float32), + origin_centroid, area, Subzero.Status(), Xoshiro(1) @@ -59,17 +60,18 @@ # test generating sub-grid points for grid with Δx = Δy = 10 point_generator = SubGridPointsGenerator{Float64}(10/sqrt(2)) # Floe is smaller than grid cells --> centroid and vertices added - square = [[ + square = Subzero.make_polygon([[ [-2.5, -2.5], [-2.5, 2.5], [2.5, 2.5], [2.5, -2.5], [-2.5, -2.5], - ]] + ]]) + square_centroid = GO.centroid(square) xpoints, ypoints = Subzero.generate_subfloe_points( point_generator, square, - 0.0, # Not used + square_centroid, 0.0, # Not used Subzero.Status(), Xoshiro(), # Not random @@ -77,17 +79,18 @@ @test xpoints == [-2.5, -2.5, 2.5, 2.5, 0.0] @test ypoints == [-2.5, 2.5, 2.5, -2.5, 0.0] # Floe is larger than grid cell - tall_rect = [[ + tall_rect = Subzero.make_polygon([[ [-2.0, -10.0], [-2.0, 10.0], [2.0, 10.0], [2.0, -10.0], [-2.0, -10.0], - ]] + ]]) + tall_rect_centroid = GO.centroid(tall_rect) xpoints, ypoints = Subzero.generate_subfloe_points( point_generator, tall_rect, - 0.0, # Not used + tall_rect_centroid, 0.0, # Not used Subzero.Status(), Xoshiro(), # Not random @@ -100,17 +103,18 @@ ], atol = 1e-5)) - wide_rect = [[ + wide_rect = Subzero.make_polygon([[ [-10.0, -2.0], [-10.0, 2.0], [10.0, 2.0], [10.0, -2.0], [-10.0, -2.0], - ]] + ]]) + wide_rect_centroid = GO.centroid(wide_rect) xpoints, ypoints = Subzero.generate_subfloe_points( point_generator, wide_rect, - 0.0, # Not used + wide_rect_centroid, 0.0, # Not used Subzero.Status(), Xoshiro(), # Not random @@ -123,21 +127,23 @@ atol = 1e-5 )) @test ypoints == [-2; repeat([2], 5); repeat([-2], 4); repeat([0], 3)] - trapeziod = [[ + trapeziod = Subzero.make_polygon([[ [-8.0, -8.0], [-4.0, 8.0], [4.0, 8.0], [8.0, -8.0], [-8.0, -8.0], - ]] + ]]) + trapeziod_centroid = GO.centroid(trapeziod) xpoints, ypoints = Subzero.generate_subfloe_points( point_generator, trapeziod, - 0.0, # Not used + trapeziod_centroid, 0.0, # Not used Subzero.Status(), Xoshiro(), # Not random ) + ypoints .+= trapeziod_centroid[2] # y-point of centroid is not centered on the origin @test all(isapprox.( xpoints, [-8, -7.14251, -6.0, -4.85749, -4.0, 0.0, 4.0, 4.85749, 6.0, @@ -154,7 +160,6 @@ ], atol = 1e-5 )) - end @testset "Coupling Helper Functions" begin @@ -274,52 +279,18 @@ ) == (0:10:80, 1:9) #Test cell_coords - cell = Subzero.center_cell_coords( - 2, - 3, - grid, - periodic_bound, - periodic_bound, - ) - 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, - 1, - grid, - open_bound, - open_bound, - ) == [[[-10, -8], [-10, -6], [-9, -6], [-9, -8], [-10, -8]]] - @test Subzero.center_cell_coords( - 11, - 6, - grid, - periodic_bound, - periodic_bound, - ) == [[[9, 10], [9, 14], [11, 14], [11, 10], [9, 10]]] - @test Subzero.center_cell_coords( - 11, - 6, - grid, - open_bound, - open_bound, - ) == [[[9, 8], [9, 8], [10, 8], [10, 8], [9, 8]]] - @test Subzero.center_cell_coords( - 11, - 6, - grid, - open_bound, - periodic_bound, - ) == [[[9, 8], [9, 8], [11, 8], [11, 8], [9, 8]]] - @test Subzero.center_cell_coords( - 11, - 6, - grid, - periodic_bound, - open_bound, - ) == [[[9, 10], [9, 14], [10, 14], [10, 10], [9, 10]]] + @test GO.equals(Subzero.center_cell_coords(Float64, 2, 3, grid, periodic_bound, periodic_bound), + GI.Polygon([[[-9, -2], [-9, 2], [-7, 2], [-7, -2], [-9, -2]]])) + @test GO.equals(Subzero.center_cell_coords(Float64, 1, 1, grid, open_bound, open_bound), + GI.Polygon([[[-10, -8], [-10, -6], [-9, -6], [-9, -8], [-10, -8]]])) + @test GO.equals(Subzero.center_cell_coords(Float64, 11, 6, grid, periodic_bound, periodic_bound), + GI.Polygon([[[9, 10], [9, 14], [11, 14], [11, 10], [9, 10]]])) + @test GO.equals(Subzero.center_cell_coords(Float64, 11, 6, grid, open_bound, open_bound), + GI.Polygon([[[9, 8], [9, 8], [10, 8], [10, 8], [9, 8]]])) + @test GO.equals(Subzero.center_cell_coords(Float64, 11, 6, grid, open_bound, periodic_bound), + GI.Polygon([[[9, 8], [9, 8], [11, 8], [11, 8], [9, 8]]])) + @test GO.equals(Subzero.center_cell_coords(Float64, 11, 6, grid, periodic_bound, open_bound), + GI.Polygon([[[9, 10], [9, 14], [10, 14], [10, 10], [9, 10]]])) # Test aggregate_grid_stress! function test_floe_to_grid( diff --git a/test/test_physical_processes/test_fractures.jl b/test/test_physical_processes/test_fractures.jl index 52b6171..c8e1e00 100644 --- a/test/test_physical_processes/test_fractures.jl +++ b/test/test_physical_processes/test_fractures.jl @@ -1,23 +1,23 @@ @testset "Fractures" begin @testset "Fracture Criteria" begin + FT = Float64 # Test NoFracturee criteria @test NoFracture() isa NoFracture # Test HiblerYieldCurve criteria - @test_throws ArgumentError HiblerYieldCurve(2.25, 20.0, [[[0.0, 0.0]]]) @test HiblerYieldCurve( 2.25, 20.0, - [[[0.0, 0.0], [0, 1], [1 ,1], [1, 0]]] + Subzero.make_polygon([[[0.0, 0.0], [0, 1], [1 ,1], [1, 0]]]), ) isa HiblerYieldCurve - # Test calculate_hibler - hibler_verts = Subzero.calculate_hibler(0.5, 5e5, -1) - hibler_poly = Subzero.make_polygon(hibler_verts) + # Test _calculate_hibler + hibler_poly = Subzero._calculate_hibler(FT, 0.5, 5e5, -1) @test isapprox(GO.area(hibler_poly), 49054437859.374, atol = -1e3) @test all(isapprox.( GO.centroid(hibler_poly), (-1.25e5, -1.25e5), atol = 1e-3 )) + hibler_verts = Subzero.find_poly_coords(hibler_poly) x_verts, y_verts = first.(hibler_verts[1]), last.(hibler_verts[1]) @test all(isapprox.( extrema(x_verts), @@ -29,8 +29,8 @@ [-264743.588, 14727.999], atol = 1e-3 )) - hibler_verts = Subzero.calculate_hibler(0.25, 2.25e5, 20.0) - hibler_poly = Subzero.make_polygon(hibler_verts) + hibler_poly = Subzero._calculate_hibler(FT, 0.25, 2.25e5, 20.0) + hibler_verts = Subzero.find_poly_coords(hibler_poly) @test isapprox(GO.area(hibler_poly), 2483380916.630, atol = -1e3) @test all(isapprox.( GO.centroid(hibler_poly), @@ -53,7 +53,7 @@ @test typeof(Subzero.MohrsCone(Float64)) <: MohrsCone{Float64} # Float64 Mohr's Cone with q, σc, σ11 - mohrs_verts_64 = Subzero.calculate_mohrs(5.2, 2.5e5, -3.375e4) + mohrs_verts_64 = Subzero.find_poly_coords(Subzero._calculate_mohrs(FT, 5.2, 2.5e5, -3.375e4)) @test all(isapprox.( mohrs_verts_64[1], [ @@ -65,7 +65,7 @@ atol = 1e-3 )) # Float32 Mohr's Cone with q, σc, σ11 - mohrs_verts_32 = Subzero.calculate_mohrs(5.2, 2.5e5, 1.5e5) + mohrs_verts_32 = Subzero.find_poly_coords(Subzero._calculate_mohrs(FT, 5.2, 2.5e5, 1.5e5)) @test all(isapprox.( mohrs_verts_32[1], [ @@ -77,7 +77,8 @@ atol = 1e-3, )) # Float64 Mohr's Cone with σ1, σ2, σ11, σ22 - mohrs_verts_coords = Subzero.calculate_mohrs( + mohrs_verts_coords = Subzero._calculate_mohrs( + FT, 5.95e4, 5.95e4, -1.5e5, @@ -90,17 +91,17 @@ 0.0, )]) yield_curve = HiblerYieldCurve(floes) - verts = deepcopy(yield_curve.vertices) + old_poly = yield_curve.poly @test yield_curve isa HiblerYieldCurve @test yield_curve.pstar == 2.25e5 && yield_curve.c == 20 floes.height .= 0.5 Subzero.update_criteria!(yield_curve, floes) - @test verts != yield_curve.vertices + @test !GO.equals(old_poly, yield_curve.poly) # Test update criteria for Mohr's cone cone_curve = MohrsCone() - verts = cone_curve.vertices + old_poly = cone_curve.poly Subzero.update_criteria!(cone_curve, floes) - @test verts == cone_curve.vertices + @test GO.equals(old_poly, cone_curve.poly) end @testset "Fractures Floes" begin # Fracture tests depend on these floes and settings @@ -188,7 +189,7 @@ 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, + Subzero.make_polygon(colliding_coords), deforming_forces, FloeSettings(), 10, diff --git a/test/test_physical_processes/test_process_settings.jl b/test/test_physical_processes/test_process_settings.jl index e89ba3a..5302e85 100644 --- a/test/test_physical_processes/test_process_settings.jl +++ b/test/test_physical_processes/test_process_settings.jl @@ -103,7 +103,7 @@ criteria = HiblerYieldCurve( 0.0, 0.0, - [[[0.0, 0.0], [0, 1], [1 ,1], [1, 0]]], + Subzero.make_polygon([[[0.0, 0.0], [0, 1], [1 ,1], [1, 0]]]), ) custom_info = FractureSettings( fractures_on = true, diff --git a/test/test_physical_processes/test_ridge_raft.jl b/test/test_physical_processes/test_ridge_raft.jl index 42e0da9..480619a 100644 --- a/test/test_physical_processes/test_ridge_raft.jl +++ b/test/test_physical_processes/test_ridge_raft.jl @@ -1,4 +1,5 @@ @testset "Ridging and Rafting" begin + FT = Float64 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 @@ -24,7 +25,7 @@ return end function setup_floes_with_inters(coords, domain, consts, - collision_settings, lock, Δx = nothing, Δy = nothing, + collision_settings, lock, FT = Float64, Δx = nothing, Δy = nothing, ) floes = initialize_floe_field( Float64, @@ -38,6 +39,7 @@ Subzero.translate!(floes.coords[i], Δx[i], Δy[i]) floes.centroid[i][1] += Δx[i] floes.centroid[i][2] += Δy[i] + floes.poly[i] = Subzero._translate_poly(FT, floes.poly[i], Δx[i], Δy[i]) end end assign_random_velocities!(floes) @@ -410,7 +412,7 @@ [[[3e4, 3e4], [3e4, 5e4], [5e4, 5e4], [5e4, 3e4], [3e4, 3e4]]] ] floes_base = setup_floes_with_inters(coords, collision_domain, consts, - collision_settings, lock, [0.0, 0.5e4], [0.0, 0.5e4], + collision_settings, lock, FT, [0.0, 0.5e4], [0.0, 0.5e4], ) 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) diff --git a/test/test_physical_processes/test_simplification.jl b/test/test_physical_processes/test_simplification.jl index 033e828..8b799b3 100644 --- a/test/test_physical_processes/test_simplification.jl +++ b/test/test_physical_processes/test_simplification.jl @@ -74,7 +74,7 @@ @test f2.coords == coords2 # Test two floes intersecting -> will fuse into floe1 since same size - Subzero.Subzero.translate!(coords2, -13.0, 0.0) + Subzero.translate!(coords2, -13.0, 0.0) f2 = Floe(coords2, 0.75, 0.0) f1.id = 1 f2.id = 2 @@ -464,7 +464,7 @@ # 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 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 sum(GO.area, Subzero.intersect_polys(Subzero.make_polygon(floe_set2.coords[1]), open_domain_with_topo.topography.poly[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 diff --git a/test/test_physical_processes/test_welding.jl b/test/test_physical_processes/test_welding.jl index 0bc082a..bc67d4a 100644 --- a/test/test_physical_processes/test_welding.jl +++ b/test/test_physical_processes/test_welding.jl @@ -149,7 +149,7 @@ end coupling_settings = CouplingSettings() coords = [ [[[0.0, 0.0], [0.0, 5e4], [6e4, 5e4], [6e4, 0.0], [0.0, 0.0]]], - [[[4e4, 0.0], [4e4, 5e4], [1e5, 5e4], [1e5, 0.0], [0.0, 0.0]]], + [[[4e4, 0.0], [4e4, 5e4], [1e5, 5e4], [1e5, 0.0], [4e4, 0.0]]], [[[2e4, 4e4], [2e4, 8e4], [3e4, 8e4], [3e4, 4e4], [2e4, 4e4]]] ] floe_base = initialize_floe_field(