commit 2159025aedb6377215c2e589172c77c90cc1be36
parent 8d953de0b3264014f25228a388cbc7179cf26f1b
Author: Anders Damsgaard <andersd@riseup.net>
Date: Fri, 2 Feb 2018 15:06:29 -0500
Implement first simple porosity estimation
Diffstat:
2 files changed, 54 insertions(+), 5 deletions(-)
diff --git a/src/grid.jl b/src/grid.jl
@@ -1036,14 +1036,48 @@ function fitGridToGrains!(simulation::Simulation, grid::Any;
nothing
end
-function findPorosity(grid::Any)
+function findPorosity!(sim::Simulation, grid::Any; verbose::Bool=true)
+ if !isassigned(grid.grain_list)
+ info("Sorting grains in grid")
+ sortGrainsInGrid!(sim, grid, verbose=verbose)
+ end
+
+ sw = Vector{Float64}(2)
+ se = Vector{Float64}(2)
+ ne = Vector{Float64}(2)
+ nw = Vector{Float64}(2)
cell_area = 0.0
- for ix in size(grid.xh, 1)
- for iy in size(grid.xh, 2)
- cell_area =
+ for ix in 1:size(grid.xh, 1)
+ for iy in 1:size(grid.xh, 2)
+
+ @views sw .= grid.xq[ ix, iy], grid.yq[ ix, iy]
+ @views se .= grid.xq[ ix+1, iy], grid.yq[ ix+1, iy]
+ @views ne .= grid.xq[ ix+1, iy+1], grid.yq[ ix+1, iy+1]
+ @views nw .= grid.xq[ ix, iy+1], grid.yq[ ix, iy+1]
+ cell_area = areaOfQuadrilateral(sw, se, ne, nw)
+
+ # Subtract grain area from cell area
+ particle_area = 0.0
+ #=for ix_ = -1:1
+ for iy_ = -1:1
+
+ if ix + ix_ < 1 || ix + ix_ > size(grid.xh, 1) ||
+ iy + iy_ < 1 || iy + iy_ > size(grid.xh, 2)
+ continue
+ end
+
+ for i in grid.grain_list[ix + ix_, iy + iy_]=#
+ for i in grid.grain_list[ix, iy]
+ particle_area +=
+ grainHorizontalSurfaceArea(sim.grains[i])
+ end
+ #=end
+ end
+ end=#
+
+ grid.porosity[ix, iy] = (cell_area - particle_area)/cell_area
end
end
-
end
diff --git a/test/grid.jl b/test/grid.jl
@@ -1,4 +1,6 @@
#!/usr/bin/env julia
+using Compat.Test
+import Granular
# Check the grid interpolation and sorting functions
verbose = false
@@ -399,3 +401,16 @@ Granular.fitGridToGrains!(sim, sim.atmosphere, padding=.5, verbose=true)
@test sim.atmosphere.yq[1,1] ≈ -1.
@test sim.atmosphere.xq[end,end] ≈ 3.5
@test sim.atmosphere.yq[end,end] ≈ 5.5
+
+info("Testing porosity estimation")
+sim = Granular.createSimulation()
+dx = 1.0; dy = 1.0
+nx = 3; ny = 3
+sim.ocean = Granular.createRegularOceanGrid([nx, ny, 1], [nx*dx, ny*dy, 1.])
+Granular.addGrainCylindrical!(sim, [1.5, 1.5], 0.5*dx, 1.0)
+A_particle = π*(0.5*dx)^2
+A_cell = dx^2
+Granular.findPorosity!(sim, sim.ocean)
+@test sim.ocean.porosity ≈ [1. 1. 1.;
+ 1. (A_cell - A_particle)/A_cell 1.;
+ 1. 1. 1]
