tadd method to scale grid according to granular assemblage - Granular.jl - Juli… | |
git clone git://src.adamsgaard.dk/Granular.jl | |
Log | |
Files | |
Refs | |
README | |
LICENSE | |
--- | |
commit 6942faf25139645b324b47eea7b671e7dde7aba6 | |
parent b9871a74c1289f2d04517636331f75e4432bc580 | |
Author: Anders Damsgaard <[email protected]> | |
Date: Mon, 6 Nov 2017 13:43:03 -0500 | |
add method to scale grid according to granular assemblage | |
Diffstat: | |
M src/atmosphere.jl | 10 ++++++---- | |
M src/grain.jl | 2 +- | |
M src/grid.jl | 81 ++++++++++++++++++++++++++++++ | |
M src/ocean.jl | 10 ++++++---- | |
M src/packing.jl | 9 +++++++++ | |
M test/grid.jl | 67 ++++++++++++++++++++++++++++-… | |
6 files changed, 164 insertions(+), 15 deletions(-) | |
--- | |
diff --git a/src/atmosphere.jl b/src/atmosphere.jl | |
t@@ -110,12 +110,14 @@ function createRegularAtmosphereGrid(n::Vector{Int}, | |
time::Array{Float64, 1} = zeros(1), | |
name::String = "unnamed") | |
- xq = repmat(linspace(origo[1], L[1], n[1] + 1), 1, n[2] + 1) | |
- yq = repmat(linspace(origo[2], L[2], n[2] + 1)', n[1] + 1, 1) | |
+ xq = repmat(linspace(origo[1], origo[1] + L[1], n[1] + 1), 1, n[2] + 1) | |
+ yq = repmat(linspace(origo[2], origo[2] + L[2], n[2] + 1)', n[1] + 1, 1) | |
dx = L./n | |
- xh = repmat(linspace(origo[1] + .5*dx[1], L[1] - .5*dx[1], n[1]), 1, n[2]) | |
- yh = repmat(linspace(origo[2] + .5*dx[2], L[2] - .5*dx[2], n[2])', n[1], 1) | |
+ xh = repmat(linspace(origo[1] + .5*dx[1], origo[1] + L[1] - .5*dx[1], | |
+ n[1]), 1, n[2]) | |
+ yh = repmat(linspace(origo[2] + .5*dx[2], origo[1] + L[2] - .5*dx[2], | |
+ n[2])', n[1], 1) | |
zl = -linspace(.5*dx[3], L[3] - .5*dx[3], n[3]) | |
diff --git a/src/grain.jl b/src/grain.jl | |
t@@ -81,7 +81,7 @@ are optional, and come with default values. The only requir… | |
* `atmosphere_drag_coeff_horiz::Float64 = 2.5e-4`: horizontal drag coefficient | |
for atmosphere against grain bottom [-]. | |
* `pressure::Float64 = 0.`: current compressive stress on grain [Pa]. | |
-* `fixed::Bool = false`: grain is fixed in space. | |
+* `fixed::Bool = false`: grain is fixed to a constant velocity (e.g. zero). | |
* `rotating::Bool = true`: grain is allowed to rotate. | |
* `enabled::Bool = true`: grain interacts with other grains. | |
* `verbose::Bool = true`: display diagnostic information during the function | |
diff --git a/src/grid.jl b/src/grid.jl | |
t@@ -792,3 +792,84 @@ function moveGrainsAcrossPeriodicBoundaries!(sim::Simulat… | |
end | |
nothing | |
end | |
+ | |
+export fitGridToGrains! | |
+""" | |
+ fitGridToGrains!(simulation, grid[, padding]) | |
+ | |
+Fit the ocean or atmosphere grid for a simulation to the current grains and | |
+their positions. | |
+""" | |
+function fitGridToGrains!(simulation::Simulation, grid::Any; | |
+ padding::Float64 = 0., verbose::Bool = true) | |
+ | |
+ if typeof(grid) != Ocean && typeof(grid) != Atmosphere | |
+ error("grid must be of Ocean or Atmosphere type") | |
+ end | |
+ | |
+ min_x = Inf | |
+ min_y = Inf | |
+ max_x = -Inf | |
+ max_y = -Inf | |
+ max_radius = 0. | |
+ | |
+ if length(simulation.grains) < 1 | |
+ error("Grains need to be initialized before calling fitGridToGrains") | |
+ end | |
+ | |
+ for grain in simulation.grains | |
+ | |
+ if grain.lin_pos[1] - grain.contact_radius < min_x | |
+ min_x = grain.lin_pos[1] - grain.contact_radius | |
+ end | |
+ | |
+ if grain.lin_pos[1] + grain.contact_radius > max_x | |
+ max_x = grain.lin_pos[1] + grain.contact_radius | |
+ end | |
+ | |
+ if grain.lin_pos[2] - grain.contact_radius < min_y | |
+ min_y = grain.lin_pos[2] - grain.contact_radius | |
+ end | |
+ | |
+ if grain.lin_pos[2] + grain.contact_radius > max_y | |
+ max_y = grain.lin_pos[2] + grain.contact_radius | |
+ end | |
+ | |
+ if grain.contact_radius > max_radius | |
+ max_radius = grain.contact_radius | |
+ end | |
+ end | |
+ min_x -= padding | |
+ min_y -= padding | |
+ max_x += padding | |
+ max_y += padding | |
+ | |
+ const L::Vector{Float64} = [max_x - min_x, max_y - min_y] | |
+ const dx::Float64 = 2.*max_radius | |
+ const n = convert(Vector{Int}, floor.(L./dx)) | |
+ if 1 in n | |
+ error("Grid is too small compared to grain size (n = $n). " * | |
+ "Use all-to-all contact search instead.") | |
+ end | |
+ | |
+ if typeof(grid) == Ocean | |
+ simulation.ocean = createRegularOceanGrid(vcat(n, 1), vcat(L, 1.), | |
+ origo=[min_x, min_y], | |
+ time=[0.], name="fitted") | |
+ elseif typeof(grid) == Atmosphere | |
+ simulation.atmosphere = createRegularAtmosphereGrid(vcat(n, 1), | |
+ vcat(L, 1.), | |
+ origo=[min_x, | |
+ min_y], | |
+ time=[0.], | |
+ name="fitted") | |
+ end | |
+ | |
+ if verbose | |
+ info("Created regular $(typeof(grid)) grid from " * | |
+ "[$min_x, $min_y] to [$max_x, $max_y] " * | |
+ "with a cell size of $dx ($n).") | |
+ end | |
+ | |
+ nothing | |
+end | |
diff --git a/src/ocean.jl b/src/ocean.jl | |
t@@ -220,12 +220,14 @@ function createRegularOceanGrid(n::Array{Int, 1}, | |
time::Vector{Float64} = zeros(1), | |
name::String = "unnamed") | |
- xq = repmat(linspace(origo[1], L[1], n[1] + 1), 1, n[2] + 1) | |
- yq = repmat(linspace(origo[2], L[2], n[2] + 1)', n[1] + 1, 1) | |
+ xq = repmat(linspace(origo[1], origo[1] + L[1], n[1] + 1), 1, n[2] + 1) | |
+ yq = repmat(linspace(origo[2], origo[2] + L[2], n[2] + 1)', n[1] + 1, 1) | |
dx = L./n | |
- xh = repmat(linspace(origo[1] + .5*dx[1], L[1] - .5*dx[1], n[1]), 1, n[2]) | |
- yh = repmat(linspace(origo[2] + .5*dx[2], L[2] - .5*dx[2], n[2])', n[1], 1) | |
+ xh = repmat(linspace(origo[1] + .5*dx[1], origo[1] + L[1] - .5*dx[1], | |
+ n[1]), 1, n[2]) | |
+ yh = repmat(linspace(origo[2] + .5*dx[2], origo[2] + L[2] - .5*dx[2], | |
+ n[2])', n[1], 1) | |
zl = -linspace(.5*dx[3], L[3] - .5*dx[3], n[3]) | |
zi = -linspace(0., L[3], n[3] + 1) | |
diff --git a/src/packing.jl b/src/packing.jl | |
t@@ -1,4 +1,13 @@ | |
## Functions for creating grain packings | |
+ | |
+export regularPacking! | |
+""" | |
+""" | |
+function regularPacking!(simulation::Simulation, | |
+ ) | |
+ | |
+end | |
+ | |
""" | |
Return random point in spherical annulus between `(r_i + r_j)` and `2.*(r_i + | |
r_j)` around `x_i`. Note: there is slightly higher point density towards (r_i… | |
diff --git a/test/grid.jl b/test/grid.jl | |
t@@ -1,7 +1,7 @@ | |
#!/usr/bin/env julia | |
# Check the grid interpolation and sorting functions | |
-verbose = true | |
+verbose = false | |
info("#### $(basename(@__FILE__)) ####") | |
t@@ -324,7 +324,7 @@ Granular.addGrainCylindrical!(sim, [2.6, 2.5], .1, 1., ver… | |
Granular.sortGrainsInGrid!(sim, sim.ocean, verbose=verbose) | |
Granular.setTimeStep!(sim) | |
Granular.setTotalTime!(sim, 1.0) | |
-Granular.run!(sim, single_step=true, verbose=verbose) | |
+Granular.run!(sim, single_step=true, verbose=false) | |
Test.@test sim.atmosphere.collocated_with_ocean_grid == true | |
Test.@test sim.grains[1].ocean_grid_pos == [1, 1] | |
Test.@test sim.grains[2].ocean_grid_pos == [1, 1] | |
t@@ -339,8 +339,63 @@ Test.@test sim.atmosphere.grain_list[1, 1] == [1, 2] | |
Test.@test sim.atmosphere.grain_list[2, 2] == [] | |
Test.@test sim.atmosphere.grain_list[3, 3] == [3] | |
-info("Testing ocean drag") | |
+info("Testing automatic grid-size adjustment") | |
+# ocean grid | |
sim = Granular.createSimulation() | |
-sim.ocean.u[:,:,1,1] = 5. | |
-Granular.addGrainCylindrical!(sim, [2.5, 3.5], 1., 1., verbose=verbose) | |
-Granular.addGrainCylindrical!(sim, [2.6, 2.5], 1., 1., verbose=verbose) | |
+Test.@test_throws ErrorException Granular.fitGridToGrains!(sim, sim.ocean) | |
+sim = Granular.createSimulation() | |
+Granular.addGrainCylindrical!(sim, [0.0, 1.5], .5, 1., verbose=verbose) | |
+Granular.addGrainCylindrical!(sim, [2.5, 5.5], 1., 1., verbose=verbose) | |
+Granular.fitGridToGrains!(sim, sim.ocean, verbose=true) | |
+Test.@test sim.ocean.xq[1,1] ≈ -.5 | |
+Test.@test sim.ocean.yq[1,1] ≈ 1.0 | |
+Test.@test sim.ocean.xq[end,end] ≈ 3.5 | |
+Test.@test sim.ocean.yq[end,end] ≈ 6.5 | |
+ | |
+sim = Granular.createSimulation() | |
+Granular.addGrainCylindrical!(sim, [0.5, 1.5], .5, 1., verbose=verbose) | |
+Granular.addGrainCylindrical!(sim, [2.5, 4.5], .5, 1., verbose=verbose) | |
+Granular.fitGridToGrains!(sim, sim.ocean, verbose=true) | |
+Test.@test sim.ocean.xq[1,1] ≈ 0. | |
+Test.@test sim.ocean.yq[1,1] ≈ 1. | |
+Test.@test sim.ocean.xq[end,end] ≈ 3. | |
+Test.@test sim.ocean.yq[end,end] ≈ 5. | |
+ | |
+sim = Granular.createSimulation() | |
+Granular.addGrainCylindrical!(sim, [0.5, 0.0], .5, 1., verbose=verbose) | |
+Granular.addGrainCylindrical!(sim, [2.0, 4.0], 1., 1., verbose=verbose) | |
+Granular.fitGridToGrains!(sim, sim.ocean, padding=.5, verbose=true) | |
+Test.@test sim.ocean.xq[1,1] ≈ -.5 | |
+Test.@test sim.ocean.yq[1,1] ≈ -1. | |
+Test.@test sim.ocean.xq[end,end] ≈ 3.5 | |
+Test.@test sim.ocean.yq[end,end] ≈ 5.5 | |
+ | |
+# atmosphere grid | |
+sim = Granular.createSimulation() | |
+Test.@test_throws ErrorException Granular.fitGridToGrains!(sim, sim.atmosphere) | |
+sim = Granular.createSimulation() | |
+Granular.addGrainCylindrical!(sim, [0.0, 1.5], .5, 1., verbose=verbose) | |
+Granular.addGrainCylindrical!(sim, [2.5, 5.5], 1., 1., verbose=verbose) | |
+Granular.fitGridToGrains!(sim, sim.atmosphere, verbose=true) | |
+Test.@test sim.atmosphere.xq[1,1] ≈ -.5 | |
+Test.@test sim.atmosphere.yq[1,1] ≈ 1.0 | |
+Test.@test sim.atmosphere.xq[end,end] ≈ 3.5 | |
+Test.@test sim.atmosphere.yq[end,end] ≈ 6.5 | |
+ | |
+sim = Granular.createSimulation() | |
+Granular.addGrainCylindrical!(sim, [0.5, 1.5], .5, 1., verbose=verbose) | |
+Granular.addGrainCylindrical!(sim, [2.5, 4.5], .5, 1., verbose=verbose) | |
+Granular.fitGridToGrains!(sim, sim.atmosphere, verbose=true) | |
+Test.@test sim.atmosphere.xq[1,1] ≈ 0. | |
+Test.@test sim.atmosphere.yq[1,1] ≈ 1. | |
+Test.@test sim.atmosphere.xq[end,end] ≈ 3. | |
+Test.@test sim.atmosphere.yq[end,end] ≈ 5. | |
+ | |
+sim = Granular.createSimulation() | |
+Granular.addGrainCylindrical!(sim, [0.5, 0.0], .5, 1., verbose=verbose) | |
+Granular.addGrainCylindrical!(sim, [2.0, 4.0], 1., 1., verbose=verbose) | |
+Granular.fitGridToGrains!(sim, sim.atmosphere, padding=.5, verbose=true) | |
+Test.@test sim.atmosphere.xq[1,1] ≈ -.5 | |
+Test.@test sim.atmosphere.yq[1,1] ≈ -1. | |
+Test.@test sim.atmosphere.xq[end,end] ≈ 3.5 | |
+Test.@test sim.atmosphere.yq[end,end] ≈ 5.5 |