Granular.jl

Julia package for granular dynamics simulation
git clone git://src.adamsgaard.dk/Granular.jl
Log | Files | Refs | README | LICENSE

commit abbde74a6163d80e6cf65e3298cb8891cc0e7e4b
parent ddf1f4d1c0e72b369fcae11a7cd37e9267265f65
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Mon,  1 May 2017 10:05:42 -0400

add flag to enable or disable ice floe

Diffstat:
Msrc/contact_search.jl | 10+++++++++-
Msrc/datatypes.jl | 2++
Msrc/grid.jl | 21+++++++++++++++++----
Msrc/icefloe.jl | 4++++
Msrc/io.jl | 1+
Msrc/ocean.jl | 5+++++
Msrc/simulation.jl | 11++++++++++-
Msrc/temporal_integration.jl | 6++++++
Mtest/grid.jl | 7+++----
Mtest/vtk.jl | 2+-
10 files changed, 58 insertions(+), 11 deletions(-)

diff --git a/src/contact_search.jl b/src/contact_search.jl @@ -61,7 +61,9 @@ function findContactsAllToAll!(simulation::Simulation) if i < j if simulation.ice_floes[i].fixed && - simulation.ice_floes[j].fixed + simulation.ice_floes[j].fixed || + !simulation.ice_floes[i].enabled || + !simulation.ice_floes[j].enabled continue end @@ -105,6 +107,12 @@ function findContactsOceanGrid!(simulation::Simulation) for idx_j in simulation.ocean.ice_floe_list[i, j] if idx_i < idx_j + if simulation.ice_floes[i].fixed && + simulation.ice_floes[j].fixed || + !simulation.ice_floes[i].enabled || + !simulation.ice_floes[j].enabled + continue + end # Inter-grain position vector and grain overlap position_ij = interIceFloePositionVector(simulation, diff --git a/src/datatypes.jl b/src/datatypes.jl @@ -34,6 +34,7 @@ type IceFloeCylindrical # Kinematic constraint flags fixed::Bool rotating::Bool + enabled::Bool # Rheological parameters contact_stiffness_normal::float @@ -80,6 +81,7 @@ type IceFloeArrays # Kinematic constraint flags fixed rotating + enabled # Rheological parameters contact_stiffness_normal diff --git a/src/grid.jl b/src/grid.jl @@ -50,6 +50,12 @@ function sortIceFloesInOceanGrid!(simulation::Simulation; verbose=false) i, j = findCellContainingPoint(simulation.ocean, simulation.ice_floes[idx].lin_pos) + # remove ice floe if it is outside of the grid + if i == 0 && j == 0 +export removeIceFloe! + + end + # add cell to ice floe simulation.ice_floes[idx].ocean_grid_pos = [i, j] @@ -60,11 +66,19 @@ end export findCellContainingPoint """ + findCellContainingPoint(ocean, point[, method]) + Returns the `i`, `j` index of the ocean grid cell containing the `point`. The function uses either an area-based approach (`method = "Area"`), or a conformal mapping approach (`method = "Conformal"`). The area-based approach is -more robust. This function returns the coordinates of the cell or raises an -error. +more robust. This function returns the coordinates of the cell. If no match is +found the function returns `(0,0)`. + +# Arguments +* `ocean::Ocean`: object containing ocean data. +* `point::Array{float, 1}`: two-dimensional vector of point to check. +* `method::String`: approach to use for determining if point is inside cell or + not, can be "Conformal" (default) or "Areal". """ function findCellContainingPoint(ocean::Ocean, point::Array{float, 1}; method::String="Conformal") @@ -75,8 +89,7 @@ function findCellContainingPoint(ocean::Ocean, point::Array{float, 1}; end end end - # throw an error for now, maybe remove ice floe later on - error("point not found in grid") + return 0, 0 end export getNonDimensionalCellCoordinates diff --git a/src/icefloe.jl b/src/icefloe.jl @@ -28,6 +28,7 @@ function addIceFloeCylindrical(simulation::Simulation, pressure::float = 0., fixed::Bool = false, rotating::Bool = true, + enabled::Bool = true, verbose::Bool = true, ocean_grid_pos::Array{Int, 1} = [0, 0]) @@ -90,6 +91,7 @@ function addIceFloeCylindrical(simulation::Simulation, fixed, rotating, + enabled, contact_stiffness_normal, contact_stiffness_tangential, @@ -183,6 +185,7 @@ function convertIceFloeDataToArrays(simulation::Simulation) Array(Int, length(simulation.ice_floes)), Array(Int, length(simulation.ice_floes)), + Array(Int, length(simulation.ice_floes)), Array(Float64, length(simulation.ice_floes)), Array(Float64, length(simulation.ice_floes)), @@ -221,6 +224,7 @@ function convertIceFloeDataToArrays(simulation::Simulation) ifarr.fixed[i] = Int(simulation.ice_floes[i].fixed) ifarr.rotating[i] = Int(simulation.ice_floes[i].rotating) + ifarr.enabled[i] = Int(simulation.ice_floes[i].enabled) ifarr.contact_stiffness_normal[i] = simulation.ice_floes[i].contact_stiffness_normal diff --git a/src/io.jl b/src/io.jl @@ -78,6 +78,7 @@ function writeIceFloeVTK(simulation::Simulation, WriteVTK.vtk_point_data(vtkfile, ifarr.fixed, "Fixed in space [-]") WriteVTK.vtk_point_data(vtkfile, ifarr.rotating, "Free to rotate [-]") + WriteVTK.vtk_point_data(vtkfile, ifarr.enabled, "Enabled [-]") WriteVTK.vtk_point_data(vtkfile, ifarr.contact_stiffness_normal, "Contact stiffness (normal) [N m^-1]") diff --git a/src/ocean.jl b/src/ocean.jl @@ -255,6 +255,11 @@ function addOceanDrag!(simulation::Simulation) u, v, h, e = interpolateOceanState(simulation.ocean, simulation.time) for ice_floe in simulation.ice_floes + + if !ice_floe.enabled + continue + end + i, j = ice_floe.ocean_grid_pos k = 1 diff --git a/src/simulation.jl b/src/simulation.jl @@ -144,7 +144,6 @@ confirmation message will be printed to stdout`." function addIceFloe!(simulation::Simulation, icefloe::IceFloeCylindrical, verbose::Bool = False) - # Append icefloe to global icefloe array push!(simulation.ice_floes, icefloe) if verbose @@ -162,6 +161,16 @@ function removeIceFloe!(simulation::Simulation, i::Integer) delete!(simulation.ice_floes, i) end +export disableIceFloe! +"Disable ice floe with index `i` in the `simulation` object." +function disableIceFloe!(simulation::Simulation, i::Integer) + if i < 1 + error("Index must be greater than 0 (i = $i)") + end + + simulation.ice_floes +end + export zeroForcesAndTorques! "Sets the `force` and `torque` values of all ice floes to zero." function zeroForcesAndTorques!(simulation::Simulation) diff --git a/src/temporal_integration.jl b/src/temporal_integration.jl @@ -21,10 +21,16 @@ function updateIceFloeKinematics!(simulation::Simulation; if method == "Two-term Taylor" for ice_floe in simulation.ice_floes + if !ice_floe.enabled + continue + end updateIceFloeKinematicsTwoTermTaylor!(ice_floe, simulation) end elseif method == "Three-term Taylor" for ice_floe in simulation.ice_floes + if !ice_floe.enabled + continue + end updateIceFloeKinematicsThreeTermTaylor!(ice_floe, simulation) end else diff --git a/test/grid.jl b/test/grid.jl @@ -82,16 +82,15 @@ ocean.u[1, 2, 1, 1] = 0.0 info("Testing cell binning - Area-based approach") @test SeaIce.findCellContainingPoint(ocean, [6.2,53.4], method="Area") == (1, 1) @test SeaIce.findCellContainingPoint(ocean, [7.2,53.4], method="Area") == (2, 1) -@test_throws ErrorException SeaIce.findCellContainingPoint(ocean, [0.2, 53.4], - method="Area") +@test SeaIce.findCellContainingPoint(ocean, [0.2,53.4], method="Area") == (0, 0) info("Testing cell binning - Conformal mapping") @test SeaIce.findCellContainingPoint(ocean, [6.2,53.4], method="Conformal") == (1, 1) @test SeaIce.findCellContainingPoint(ocean, [7.2,53.4], method="Conformal") == (2, 1) -@test_throws ErrorException SeaIce.findCellContainingPoint(ocean, [0.2, 53.4], - method="Conformal") +@test SeaIce.findCellContainingPoint(ocean, [0.2, 53.4], method="Conformal") == + (0, 0) sim = SeaIce.createSimulation() sim.ocean = SeaIce.readOceanNetCDF("Baltic/00010101.ocean_month.nc", diff --git a/test/vtk.jl b/test/vtk.jl @@ -18,7 +18,7 @@ elseif Base.is_apple() else error("checksum verification of VTK file not supported on this platform") end -@test readstring(`$(cmd) test.icefloes.1.vtu`) == "1c0c2bdd265abdda22ef3727e7cac829e2321462d494be2e23364653f9529c87 test.icefloes.1.vtu\n" +@test readstring(`$(cmd) test.icefloes.1.vtu`) == "72f4e4b854d7e92afd8cde0b79a4af6a29e49714b751ffc30a4ff3867f44b50 test.icefloes.1.vtu\n" @test readstring(`$(cmd) test.ocean.1.vts`) == "f0117e414c4e71a0c55980f63865eb03b6c597fa2546983258b8a57eb4ff2a25 test.ocean.1.vts\n" rm("test.icefloes.1.vtu")