Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor Makie plotting #13

Merged
merged 8 commits into from
Nov 21, 2019
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 10 additions & 0 deletions examples/1d_gas.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
using HardSphereDynamics
include("makie_visualization.jl")

fluid = HardSphereFluid(1, 10, 0.03) # create hard spheres in 1D
positions, velocities, times = evolve!(fluid, δt, final_time)

positions2 = HardSphereDynamics.to_3D(positions)
velocities2 = HardSphereDynamics.to_3D(velocities)

visualize_3d(positions2, velocities2, [ball.r for ball in fluid.balls], sleep_step=0.01)
35 changes: 0 additions & 35 deletions examples/makie_2d.jl

This file was deleted.

88 changes: 0 additions & 88 deletions examples/makie_3d.jl

This file was deleted.

107 changes: 107 additions & 0 deletions examples/makie_visualization.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,107 @@
using HardSphereDynamics

using Makie
using GeometryTypes
using AbstractPlotting
using Colors
using StatsBase
using StaticArrays

using LinearAlgebra



function visualize_2d(fluid, positions, velocities, times, sleep_step=0.001)

data = Makie.Node(Point2f0.(positions[1]))
# limits = FRect2D(fluid.box.lower, fluid.box.upper)
limits = FRect(0.0f0, 0.0f0, 1.0f0, 1.0f0)

# color by speed:
cs = Makie.Node(norm.(velocities[1]))
mean_c = mean(cs[])
crange = (0.0, 2 * mean_c)


scene = Scene(resolution = (1000, 1000))
s = Makie.scatter!(scene, data, marker=Circle(Point2f0(0), 10.0f0),
color=cs, colorrange=crange, colormap=:viridis,
markersize=[ball.r for ball in fluid.balls],
limits=limits)

display(s)

# update positions and velocities to animate:
for t in 1:length(positions)
data[] = positions[t]
cs[] = norm.(velocities[t])
sleep(sleep_step)
end

end

visualize_3d(fluid::HardSphereFluid, positions, velocities; sleep_step=0.001) =
visualize_3d(positions, velocities, [ball.r for ball in fluid.balls], sleep_step=sleep_step)


function visualize_3d(positions, velocities, radii;
lower = -0.5*ones(SVector{3,Float32}),
upper = 0.5*ones(SVector{3,Float32}),
sleep_step=0.001)

data = Makie.Node(Point3f0.(positions[1]))
limits = FRect3D(lower, upper .- lower) # 2nd argument are widths in each direction

# color by speed:
cs = Makie.Node(norm.(velocities[1]))
mean_c = mean(cs[])
crange = (0.0, 2 * mean_c)

scene = Scene(resolution = (1000, 1000))
s = Makie.meshscatter!(scene, data,
color=cs, colorrange=crange, colormap=:viridis,
markersize=radii,
limits=limits)

display(s)

# update positions and velocities to animagte:
for t in 1:length(positions)
data[] = positions[t]
cs[] = norm.(velocities[t])
sleep(sleep_step)
end

return data, cs

end

to_2d(v::SVector{1,T}) where {T} = SVector(zero(T), v[1])
to_3d(v::SVector{1,T}) where {T} = SVector(zero(T), zero(T), v[1])
to_3d(v::SVector{2,T}) where {T} = SVector(zero(T), v[1], v[2])

to_2d(v::Vector) = to_2D.(v)
to_3d(v::Vector) = to_3D.(v)





## Run:

d = 3
n = 200 # number of spheres
r = 0.05 # radius

δt = 0.01
final_time = 100

fluid = HardSphereFluid(d, n, r) # create hard spheres in unit box in d dimensions
positions, velocities, times = evolve!(fluid, δt, final_time)

visualize_3d(fluid, positions, velocities)


## Maxwell--Boltzmann distribution:
# using Plots
# Plots.histogram(reduce(vcat, [norm.(v) for v in velocities]), normed=true)
5 changes: 3 additions & 2 deletions src/box.jl
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,9 @@ end
RectangularBox(lower::SVector{N,T}, upper::SVector{N,T}) where {N,T} = RectangularBox{N,T}(lower, upper)



function unit_hypercube(N, T)
return RectangularBox(-0.5 * ones(SVector{N,T}),
+0.5 * ones(SVector{N,T}))
+0.5 * ones(SVector{N,T})
+0.5 * ones(SVector{N,T})
)
end