Add floe creation bounds (#89)

* New example file

* Add floe creation bounds within domain

* Add comment

examples/forcing_contained_floes.jl

@@ -0,0 +1,106 @@
+using JLD2, Random, Statistics, Subzero
+
+# User Inputs
+const FT = Float64
+const Lx = 1e5
+const Ly = 1e5
+const Δgrid = 2e3
+const hmean = 0.25
+const Δh = 0.0
+const Δt = 20
+
+# Model instantiation
+grid = RegRectilinearGrid(
+    (0.0, Lx),
+    (0.0, Ly),
+    Δgrid,
+    Δgrid,
+)
+# Set ocean u velocities
+ocean_uvels = zeros(FT, (grid.Nx + 1, grid.Ny + 1))
+for i in CartesianIndices(ocean_uvels)
+    r, c = Tuple(i)
+    if r ≤ 5
+        ocean_uvels[i] = 0.2
+    elseif r ≥ grid.Nx - 4
+        ocean_uvels[i] = -0.2
+    end
+end
+ocean_uvels[20:40, 20:30] .= 0.15
+# Set ocean v velocities
+ocean_vvels = zeros(FT, (grid.Nx + 1, grid.Ny + 1))
+for i in CartesianIndices(ocean_vvels)
+    r, c = Tuple(i)
+    if c ≤ 5
+        ocean_vvels[i] = 0.2
+    elseif c ≥ grid.Ny - 4
+        ocean_vvels[i] = -0.2
+    end
+end
+# Create ocean
+ocean = Ocean(
+    ocean_uvels,
+    ocean_vvels,
+    zeros(FT, (grid.Nx + 1, grid.Ny + 1)),
+)
+
+# Create atmosphere
+atmos = Atmos(grid, 0.0, 0.0, -1.0)
+
+# Domain creation
+nboundary = OpenBoundary(North, grid)
+sboundary = OpenBoundary(South, grid)
+eboundary = OpenBoundary(East, grid)
+wboundary = OpenBoundary(West, grid)
+
+domain = Domain(nboundary, sboundary, eboundary, wboundary)
+
+floe_settings = FloeSettings(
+    subfloe_point_generator = SubGridPointsGenerator(grid, 2),
+)
+# Floe creation - bound floes within smaller part of the domain
+floe_bounds = [[[0.1Lx, 0.1Ly], [0.1Lx, 0.9Ly], [0.9Lx, 0.9Ly], [0.9Lx, 0.1Ly], [0.1Lx, 0.1Ly]]]
+floe_arr = initialize_floe_field(
+    FT,
+    300,
+    [0.4],
+    domain,
+    hmean,
+    Δh;
+    floe_bounds = floe_bounds,
+    rng = Xoshiro(1),
+    floe_settings = floe_settings
+)
+
+# Model creation
+model = Model(grid, ocean, atmos, domain, floe_arr)
+
+# Simulation setup
+modulus = 1.5e3*(mean(sqrt.(floe_arr.area)) + minimum(sqrt.(floe_arr.area)))
+consts = Constants(E = modulus)
+
+# Run simulation
+run_time!(simulation) =  @time run!(simulation)
+dir = "output/contained"
+
+# Output setup
+initwriter = InitialStateOutputWriter(dir = dir, overwrite = true)
+floewriter = FloeOutputWriter(50, dir = dir, overwrite = true)
+writers = OutputWriters(initwriter, floewriter)
+simulation = Simulation(
+    model = model,
+    consts = consts,
+    Δt = Δt,
+    nΔt = 15000,
+    verbose = true,
+    writers = writers,
+    rng = Xoshiro(1),
+)
+run_time!(simulation)
+
+plot_sim(
+    joinpath(dir, "floes.jld2"),
+    joinpath(dir, "initial_state.jld2"),
+    Δt,
+    joinpath(dir, "contained.mp4"),
+)

src/simulation_components/floe.jl

@@ -606,9 +606,12 @@ Inputs:
     hmean           <Float> average floe height
     Δh              <Float> height range - floes will range in height from
                         hmean - Δh to hmean + Δh
-    floe_settings       <FloeSettings> settings needed to initialize floes
-    rng                 <RNG> random number generator to generate random floe
-                            attributes - default uses Xoshiro256++
+    floe_bounds     <PolyVec> coordinates of boundary within which to populate floes. This
+                        can be smaller that the domain, but will be limited to open space
+                        within the domain
+    floe_settings   <FloeSettings> settings needed to initialize floes
+    rng             <RNG> random number generator to generate random floe
+                        attributes - default uses Xoshiro256++
 Output:
     floe_arr <StructArray> list of floes created using Voronoi Tesselation
         of the domain with given concentrations.
@@ -620,40 +623,39 @@ 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),
     rng = Xoshiro(),
 ) where {FT <: AbstractFloat}
     floe_arr = StructArray{Floe{FT}}(undef, 0)
-    # Split domain into cells with given concentrations
-    nrows, ncols = size(concentrations[:, :])
-    Lx = domain.east.val - domain.west.val
-    Ly = domain.north.val - domain.south.val
-    rowlen = Ly / nrows
-    collen = Lx / ncols
-    # Availible space in whole domain
-    open_water = LG.Polygon(rect_coords(
-        domain.west.val,
-        domain.east.val,
-        domain.south.val,
-        domain.north.val
-    ))
+    # Availible space in domain
+    open_water = LG.Polygon(floe_bounds)
+    open_water = LG.intersection(open_water, LG.Polygon(rect_coords(domain.west.val, domain.east.val, domain.south.val, domain.north.val)))
     if !isempty(domain.topography)
         open_water = LG.difference(
             open_water, 
             LG.MultiPolygon(domain.topography.coords)
         )
     end
+    (bounds_xmin, bounds_xmax), (bounds_ymin, bounds_ymax) = GeometryBasics.GeoInterface.extent(open_water)
     open_water_area = LG.area(open_water)
 
+    # Split domain into cells with given concentrations
+    nrows, ncols = size(concentrations[:, :])
+    Lx = bounds_xmax - bounds_xmin
+    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)
             c = concentrations[i, j]
             if c > 0
                 c = c > 1 ? 1 : c
                 # Grid cell bounds
-                xmin = domain.west.val + collen * (j - 1)
-                ymin = domain.south.val + rowlen * (i - 1)
+                xmin = bounds_xmin + collen * (j - 1) # TODO: change this
+                ymin = bounds_ymin + rowlen * (i - 1)
                 cell_bounds = rect_coords(
                     xmin,
                     xmin + collen,