More gravity experiments
This commit is contained in:
parent
aa8f2faf75
commit
b0bd606727
1 changed files with 675 additions and 93 deletions
|
|
@ -6,11 +6,13 @@
|
|||
#include <psemek/geom/camera.hpp>
|
||||
#include <psemek/geom/constants.hpp>
|
||||
#include <psemek/util/clock.hpp>
|
||||
#include <psemek/util/to_string.hpp>
|
||||
#include <psemek/audio/engine.hpp>
|
||||
#include <psemek/audio/constants.hpp>
|
||||
#include <psemek/audio/oscillator.hpp>
|
||||
#include <psemek/audio/effect/compressor.hpp>
|
||||
#include <psemek/prof/profiler.hpp>
|
||||
#include <psemek/util/moving_average.hpp>
|
||||
|
||||
#include <random>
|
||||
|
||||
|
|
@ -62,6 +64,8 @@ struct particle
|
|||
geom::vector<float, 2> delta_pos{0.f, 0.f};
|
||||
geom::vector<float, 2> delta_vel{0.f, 0.f};
|
||||
|
||||
geom::point<float, 2> old_pos{0.f, 0.f};
|
||||
|
||||
geom::vector<float, 2> acc{0.f, 0.f};
|
||||
|
||||
float T = 0.f;
|
||||
|
|
@ -70,88 +74,16 @@ struct particle
|
|||
float const G = 50.f;
|
||||
float const GG = 0*1000.f;
|
||||
float const GC = 0*100000.f;
|
||||
float const K = 10000.f;
|
||||
float const FR = 0.99f;
|
||||
|
||||
float const dt = 0.01f;
|
||||
|
||||
float const world_size = 10000.f;
|
||||
float const world_size = 10000000.f;
|
||||
geom::point world_center{0.f, 0.f};
|
||||
|
||||
struct sound_stream
|
||||
: audio::stream
|
||||
{
|
||||
sound_stream(std::vector<particle> const & particles)
|
||||
: particles_(particles)
|
||||
{}
|
||||
|
||||
std::optional<std::size_t> length() const override
|
||||
{
|
||||
return std::nullopt;
|
||||
}
|
||||
|
||||
std::size_t read(float * data, std::size_t sample_count) override
|
||||
{
|
||||
std::size_t time = played_.load() / 2;
|
||||
|
||||
if (velocity_.size() < particles_.size())
|
||||
{
|
||||
velocity_.resize(particles_.size(), 0.f);
|
||||
amplitude_.resize(particles_.size(), 0.f);
|
||||
oscillator_.resize(particles_.size(), 0.f);
|
||||
}
|
||||
|
||||
float mean_v = 0.f;
|
||||
for (std::size_t i = 0; i < particles_.size(); ++i)
|
||||
{
|
||||
float v = geom::length(particles_[i].vel);
|
||||
velocity_[i] = geom::lerp(velocity_[i], v, 0.1f);
|
||||
mean_v += velocity_[i];
|
||||
}
|
||||
mean_v /= particles_.size();
|
||||
|
||||
for (std::size_t i = 0; i < particles_.size(); ++i)
|
||||
{
|
||||
// oscillator_[i].frequency(500.f * velocity_[i] / mean_v);
|
||||
// oscillator_[i].frequency(geom::lerp(100.f, 10000.f, velocity_[i] / 1000.f));
|
||||
oscillator_[i].frequency(10.f * velocity_[i]);
|
||||
// oscillator_[i].frequency(0.1f * geom::sqr(velocity_[i]));
|
||||
}
|
||||
|
||||
|
||||
for (std::size_t i = 0; i < sample_count; i += 2)
|
||||
{
|
||||
float v = 0.f;
|
||||
for (std::size_t j = 0; j < particles_.size(); ++j)
|
||||
v += oscillator_[j].next().imag();
|
||||
|
||||
if (i >= 2)
|
||||
{
|
||||
v = (v + data[i - 1]) / 2.f;
|
||||
}
|
||||
|
||||
data[i + 0] = v;
|
||||
data[i + 1] = v;
|
||||
|
||||
time += 1;
|
||||
}
|
||||
|
||||
played_.fetch_add(sample_count);
|
||||
|
||||
return sample_count;
|
||||
}
|
||||
|
||||
// The number of samples already played from this stream
|
||||
std::size_t played() const override
|
||||
{
|
||||
return played_.load();
|
||||
}
|
||||
|
||||
private:
|
||||
std::vector<particle> const & particles_;
|
||||
std::vector<float> velocity_;
|
||||
std::vector<float> amplitude_;
|
||||
std::vector<audio::oscillator> oscillator_;
|
||||
std::atomic<std::size_t> played_;
|
||||
};
|
||||
int const SOLVE_ITERATIONS = 16;
|
||||
float const BIAS = 1.f / 8.f;
|
||||
|
||||
struct myapp : app::app
|
||||
{
|
||||
|
|
@ -163,15 +95,13 @@ struct myapp : app::app
|
|||
|
||||
vsync(true);
|
||||
|
||||
audio_.output()->stream(audio::compressor(std::make_shared<sound_stream>(particles_), audio::from_db(-4.f), 0.95f));
|
||||
|
||||
std::uniform_real_distribution<float> d{-50.f, 50.f};
|
||||
std::uniform_real_distribution<float> rr{0.5f, 2.f};
|
||||
std::uniform_real_distribution<float> rden{0.25f, 1.f};
|
||||
std::uniform_real_distribution<float> ra{0.f, 2.f * geom::pi};
|
||||
|
||||
float min_R = 0.f;
|
||||
float max_R = 50.f;
|
||||
float max_R = 100.f;
|
||||
std::uniform_real_distribution<float> rR{0.f, 1.f};
|
||||
|
||||
bool star = false;
|
||||
|
|
@ -179,8 +109,10 @@ struct myapp : app::app
|
|||
if (star)
|
||||
particles_.push_back({{0.f, 0.f}, {0.f, 0.f}, 0.f, 0.f, 10.f, geom::pi * 10000.f, 1.f});
|
||||
|
||||
// particles_.push_back({{-1.f, 0.f}, {0.f, -1.f}, 0.f, 0.f, 1.f, 1.f, 1.f});
|
||||
// particles_.push_back({{ 1.f, 0.f}, {0.f, 1.f}, 0.f, 0.f, 1.f, 1.f, 1.f});
|
||||
// particles_.push_back({{-2.f, 0.f}, {0.f, 0.f}, 0.f, 0.f, 1.f, 1.f, 1.f});
|
||||
// particles_.push_back({{ 2.f, 0.f}, {0.f, 0.f}, 0.f, 0.f, 1.f, 1.f, 1.f});
|
||||
// particles_.push_back({{ 0.f, 3.f}, {0.f, 0.f}, 0.f, 0.f, 1.f, 1.f, 1.f});
|
||||
// particles_.push_back({{ 0.f, -3.f}, {0.f, 0.f}, 0.f, 0.f, 1.f, 1.f, 1.f});
|
||||
|
||||
// float planet_R[] = {200.f, 300.f, 400.f};
|
||||
// float planet_a[] = {0.f, geom::rad(120.f), geom::rad(240.f)};
|
||||
|
|
@ -219,9 +151,9 @@ struct myapp : app::app
|
|||
R * std::cos(a),
|
||||
R * std::sin(a),
|
||||
};
|
||||
p += geom::vector{((i % 2) ? -1000.f : 1000.f), 0.f};
|
||||
// p += geom::vector{((i % 2) ? -1000.f : 1000.f), 0.f};
|
||||
|
||||
v[0] = {(i % 2) ? 100.f : -100.f};
|
||||
// v[0] = {(i % 2) ? 100.f : -100.f};
|
||||
|
||||
if (std::all_of(particles_.begin(), particles_.end(), [&](particle const & q){ return geom::distance(q.pos, p) > q.radius + r; }))
|
||||
break;
|
||||
|
|
@ -244,10 +176,15 @@ struct myapp : app::app
|
|||
particles_[i].vel = geom::ort(r / R) * V;
|
||||
(void)V;
|
||||
}
|
||||
|
||||
for (auto & p : particles_)
|
||||
p.old_pos = p.pos;
|
||||
}
|
||||
|
||||
void on_resize(int width, int height) override
|
||||
{
|
||||
app::on_resize(width, height);
|
||||
|
||||
gl::Viewport(0, 0, width, height);
|
||||
|
||||
window_size_ = {width, height};
|
||||
|
|
@ -331,12 +268,21 @@ struct myapp : app::app
|
|||
|
||||
if (key == SDLK_SPACE)
|
||||
{
|
||||
particles_.push_back({{200.f, 20.f}, {-1000.f, 0.f}, 0.f, 0.f, 1.f, 100.f, 1.f});
|
||||
paused_ = !paused_;
|
||||
}
|
||||
|
||||
if (key == SDLK_c)
|
||||
{
|
||||
particles_.push_back({{200.f, 0.f}, {-5000.f, 0.f}, 0.f, 0.f, 1.f, 100.f, 1.f});
|
||||
particles_.back().old_pos = particles_.back().pos;
|
||||
}
|
||||
}
|
||||
|
||||
void update() override
|
||||
{
|
||||
if (paused_)
|
||||
return;
|
||||
|
||||
for (std::size_t step = 0; step < 1; ++step)
|
||||
{
|
||||
{
|
||||
|
|
@ -378,6 +324,39 @@ struct myapp : app::app
|
|||
}
|
||||
}
|
||||
|
||||
// Force-based collisions
|
||||
// if (false)
|
||||
for (std::size_t i = 0; i < particles_.size(); ++i)
|
||||
{
|
||||
for (std::size_t j = i + 1; j < particles_.size(); ++j)
|
||||
{
|
||||
auto const r = particles_[i].pos - particles_[j].pos;
|
||||
float const l = length(r);
|
||||
|
||||
float R = particles_[i].radius + particles_[j].radius;
|
||||
auto n = r / l;
|
||||
|
||||
if (l < R)
|
||||
{
|
||||
auto f = K * n * (R - l);
|
||||
(void)n;
|
||||
|
||||
// auto const f = -2.f * G * particles_[i].mass * particles_[j].mass * r / std::pow(l, 3.f);
|
||||
// f -= 2.f * G * particles_[i].mass * particles_[j].mass * r / std::pow(l, 3.f);
|
||||
|
||||
particles_[i].acc += f / particles_[i].mass;
|
||||
particles_[j].acc -= f / particles_[j].mass;
|
||||
}
|
||||
|
||||
// auto vij = particles_[i].vel - particles_[j].vel;
|
||||
|
||||
// auto vn = geom::dot(vij, n) * n;
|
||||
|
||||
// particles_[i].acc -= vn * particles_[j].mass / (particles_[i].mass + particles_[j].mass);
|
||||
// particles_[j].acc += vn * particles_[i].mass / (particles_[i].mass + particles_[j].mass);
|
||||
}
|
||||
}
|
||||
|
||||
total_forces_ += clock.count();
|
||||
}
|
||||
|
||||
|
|
@ -386,6 +365,12 @@ struct myapp : app::app
|
|||
// log::info() << "Force: #" << i << " = " << std::setprecision(10) << particles_[i].acc;
|
||||
// }
|
||||
|
||||
if (is_key_down(SDLK_m))
|
||||
{
|
||||
for (auto & p : particles_)
|
||||
p.vel *= std::exp(-100.f*dt);
|
||||
}
|
||||
|
||||
if (force_target_ && false)
|
||||
{
|
||||
for (auto & p : particles_)
|
||||
|
|
@ -417,13 +402,466 @@ struct myapp : app::app
|
|||
if (force_target_)
|
||||
world_center = *force_target_;
|
||||
|
||||
// if (false)
|
||||
for (auto & p : particles_)
|
||||
{
|
||||
p.vel += p.acc * dt;
|
||||
p.vel *= FR;
|
||||
p.pos += p.vel * dt;
|
||||
p.angle += p.angle_vel * dt;
|
||||
}
|
||||
|
||||
// Verlet
|
||||
if (false)
|
||||
for (auto & p : particles_)
|
||||
{
|
||||
auto old = p.pos;
|
||||
p.pos += (p.pos - p.old_pos) + p.acc * dt * dt;
|
||||
p.old_pos = old;
|
||||
}
|
||||
|
||||
// New iterative algorithm
|
||||
if (false)
|
||||
{
|
||||
std::vector<geom::point<float, 2>> old_pos(particles_.size());
|
||||
for (std::size_t i = 0; i < particles_.size(); ++i)
|
||||
old_pos[i] = particles_[i].pos;
|
||||
|
||||
struct collision
|
||||
{
|
||||
std::size_t i, j;
|
||||
geom::vector<float, 2> n; // i -> j
|
||||
};
|
||||
|
||||
std::vector<collision> collisions;
|
||||
|
||||
for (std::size_t i = 0; i < particles_.size(); ++i)
|
||||
{
|
||||
for (std::size_t j = i + 1; j < particles_.size(); ++j)
|
||||
{
|
||||
auto const r = particles_[i].pos - particles_[j].pos;
|
||||
|
||||
float const R = particles_[i].radius + particles_[j].radius;
|
||||
float const l = length(r);
|
||||
|
||||
if (l < R)
|
||||
collisions.push_back({i, j, - r / l});
|
||||
}
|
||||
}
|
||||
|
||||
for (std::size_t iteration = 0; iteration < 16; ++iteration)
|
||||
{
|
||||
for (auto const & c : collisions)
|
||||
{
|
||||
auto const r = particles_[c.j].pos - particles_[c.i].pos;
|
||||
auto const R = particles_[c.j].radius + particles_[c.i].radius;
|
||||
auto const l = geom::dot(r, c.n);
|
||||
|
||||
if (l < R)
|
||||
{
|
||||
auto const M = particles_[c.j].mass + particles_[c.i].mass;
|
||||
|
||||
auto const d = (l - R) * c.n / M;
|
||||
|
||||
particles_[c.i].pos += d * particles_[c.j].mass;
|
||||
particles_[c.j].pos -= d * particles_[c.i].mass;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for (std::size_t i = 0; i < particles_.size(); ++i)
|
||||
particles_[i].vel += (particles_[i].pos - old_pos[i]) / dt;
|
||||
}
|
||||
|
||||
// First-event algorithm
|
||||
if (false)
|
||||
{
|
||||
// std::vector<float> particle_time(particles_.size(), 0.f);
|
||||
|
||||
for (std::size_t i = 0; i < particles_.size(); ++i)
|
||||
{
|
||||
bool collided = false;
|
||||
|
||||
for (std::size_t j = 0; j < particles_.size(); ++j)
|
||||
{
|
||||
if (i == j) continue;
|
||||
|
||||
// auto dp = (particles_[i].pos - particles_[i].vel * particle_time[i]) - (particles_[j].pos - particles_[j].vel * particle_time[j]);
|
||||
auto dp = particles_[i].pos - particles_[j].pos;
|
||||
auto dv = particles_[i].vel - particles_[j].vel;
|
||||
auto R = particles_[i].radius + particles_[j].radius;
|
||||
|
||||
float A = geom::dot(dv, dv);
|
||||
float B = geom::dot(dp, dv);
|
||||
float C = geom::dot(dp, dp) - geom::sqr(R);
|
||||
|
||||
if (B >= 0.f) return;
|
||||
|
||||
float collision_time;
|
||||
|
||||
if (C <= 0.f)
|
||||
{
|
||||
collision_time = 0.f;
|
||||
}
|
||||
else
|
||||
{
|
||||
float D = B * B - A * C;
|
||||
|
||||
if (D <= 0.f)
|
||||
continue;
|
||||
|
||||
collision_time = C / (- B + std::sqrt(D));
|
||||
}
|
||||
|
||||
if (collision_time >= dt)
|
||||
continue;
|
||||
|
||||
collided = true;
|
||||
particles_[i].pos += particles_[i].vel * collision_time;
|
||||
|
||||
// particles_[i].pos += particles_[i].vel * (e.t - particle_time[i]);
|
||||
// particles_[j].pos += particles_[j].vel * (e.t - particle_time[j]);
|
||||
|
||||
// particle_time[i] = e.t;
|
||||
// particle_time[j] = e.t;
|
||||
|
||||
auto n = geom::normalized(particles_[i].pos - particles_[j].pos);
|
||||
|
||||
float S = 0.5f * (1.f / particles_[i].mass + 1.f / particles_[j].mass);
|
||||
float T = dot(n, particles_[i].vel - particles_[j].vel);
|
||||
|
||||
auto np = n * (- T / S);
|
||||
|
||||
np *= std::exp(-10.f * dt);
|
||||
|
||||
particles_[i].vel += np / particles_[i].mass;
|
||||
particles_[j].vel -= np / particles_[j].mass;
|
||||
|
||||
break;
|
||||
}
|
||||
|
||||
if (!collided)
|
||||
particles_[i].pos += particles_[i].vel * dt;
|
||||
}
|
||||
}
|
||||
|
||||
// Event-based algorithm
|
||||
if (false)
|
||||
{
|
||||
// log::info() << "===== ITERATION =====";
|
||||
collisions_ = 0;
|
||||
struct collision_event
|
||||
{
|
||||
std::uint32_t i0, i1;
|
||||
float t;
|
||||
bool erased = false;
|
||||
};
|
||||
|
||||
struct comparator
|
||||
{
|
||||
bool operator()(collision_event const & e1, collision_event const & e2) const
|
||||
{
|
||||
return std::tie(e1.t, e1.i0, e1.i1) < std::tie(e2.t, e2.i0, e2.i1);
|
||||
}
|
||||
};
|
||||
|
||||
using events_container = std::set<collision_event, comparator>;
|
||||
events_container events;
|
||||
std::vector<events_container::node_type> erased_events;
|
||||
std::vector<std::vector<events_container::iterator>> event_list(particles_.size());
|
||||
|
||||
std::vector<float> particle_time(particles_.size(), 0.f);
|
||||
|
||||
std::unordered_map<geom::vector<int, 2>, std::vector<std::size_t>> cells;
|
||||
|
||||
float time = 0.f;
|
||||
|
||||
auto check_collision = [&](std::size_t i, std::size_t j)
|
||||
{
|
||||
if (i > j)
|
||||
std::swap(i, j);
|
||||
|
||||
// (pi + vi * (t - ti) - pj - vj * (t - tj))^2 = (ri+rj)^2
|
||||
|
||||
auto dp = (particles_[i].pos - particles_[i].vel * particle_time[i]) - (particles_[j].pos - particles_[j].vel * particle_time[j]);
|
||||
auto dv = particles_[i].vel - particles_[j].vel;
|
||||
auto R = particles_[i].radius + particles_[j].radius;
|
||||
|
||||
float A = geom::dot(dv, dv);
|
||||
float B = geom::dot(dp, dv);
|
||||
float C = geom::dot(dp, dp) - geom::sqr(R);
|
||||
|
||||
// log::info() << "DV = " << dv;
|
||||
// log::info() << "DP = " << dp;
|
||||
// log::info() << "B = " << B;
|
||||
|
||||
// if (C < 0.f)
|
||||
// log::info() << "DEFINITELY A COLLISION!";
|
||||
|
||||
// if (B > 0.f) return;
|
||||
|
||||
// auto result = geom::solve_quadratic(A, B, C);
|
||||
// if (!result) return;
|
||||
|
||||
// log::info() << result->first << " " << result->second;
|
||||
|
||||
// float min_time = 0.f;
|
||||
// min_time = std::max(particle_time[i], particle_time[j]);
|
||||
|
||||
// std::optional<float> collision_time;
|
||||
|
||||
// if (result->first >= min_time && result->first <= dt)
|
||||
// collision_time = result->first;
|
||||
// else if (result->second >= min_time && result->second <= dt)
|
||||
// collision_time = result->second;
|
||||
|
||||
// if (result->first >= 0.f && result->first <= dt)
|
||||
// if (result->first <= dt)
|
||||
// collision_time = result->first;
|
||||
// collision_time = std::max(min_time, result->first);
|
||||
// (void)min_time;
|
||||
|
||||
// if (!collision_time) return;
|
||||
|
||||
if (B >= -0.01f * geom::length(dv) * geom::length(dp)) return;
|
||||
|
||||
float collision_time;
|
||||
|
||||
if (C <= 0.f)
|
||||
{
|
||||
collision_time = time;
|
||||
}
|
||||
else
|
||||
{
|
||||
float D = B * B - A * C;
|
||||
|
||||
if (D <= 0.f)
|
||||
return;
|
||||
|
||||
collision_time = C / (- B + std::sqrt(D));
|
||||
}
|
||||
|
||||
if (collision_time >= dt) return;
|
||||
|
||||
auto dpn = (particles_[i].pos + particles_[i].vel * (collision_time - particle_time[i])) - (particles_[j].pos + particles_[j].vel * (collision_time - particle_time[j]));
|
||||
if (geom::dot(dv, dpn) >= -0.01f * geom::length(dv) * geom::length(dpn)) return;
|
||||
|
||||
auto insert_result = events.insert(collision_event{i, j, collision_time});
|
||||
if (insert_result.second)
|
||||
{
|
||||
event_list[i].push_back(insert_result.first);
|
||||
event_list[j].push_back(insert_result.first);
|
||||
}
|
||||
};
|
||||
|
||||
auto erase = [&](events_container::iterator it)
|
||||
{
|
||||
if (!it->erased)
|
||||
{
|
||||
auto n = events.extract(it);
|
||||
n.value().erased = true;
|
||||
erased_events.push_back(std::move(n));
|
||||
}
|
||||
};
|
||||
|
||||
auto clear_collisions = [&](std::size_t i)
|
||||
{
|
||||
for (auto it : event_list[i])
|
||||
erase(it);
|
||||
event_list[i].clear();
|
||||
};
|
||||
|
||||
auto foreach_cells = [&](std::size_t i, auto && callback)
|
||||
{
|
||||
geom::box<float, 2> bbox;
|
||||
|
||||
bbox |= particles_[i].pos - geom::vector{1.f, 1.f} * particles_[i].radius;
|
||||
bbox |= particles_[i].pos + geom::vector{1.f, 1.f} * particles_[i].radius;
|
||||
|
||||
auto npos = particles_[i].pos + particles_[i].vel * dt;//(dt - particle_time[i]);
|
||||
|
||||
bbox |= npos - geom::vector{1.f, 1.f} * particles_[i].radius;
|
||||
bbox |= npos + geom::vector{1.f, 1.f} * particles_[i].radius;
|
||||
|
||||
int xmin = std::floor(bbox[0].min);
|
||||
int xmax = std::floor(bbox[0].max);
|
||||
int ymin = std::floor(bbox[1].min);
|
||||
int ymax = std::floor(bbox[1].max);
|
||||
|
||||
for (int x = xmin; x <= xmax; ++x)
|
||||
{
|
||||
for (int y = ymin; y <= ymax; ++y)
|
||||
{
|
||||
callback(cells[{x, y}]);
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
auto clear_cells = [&](std::size_t i)
|
||||
{
|
||||
foreach_cells(i, [i](auto & cell){
|
||||
auto it = std::find(cell.begin(), cell.end(), i);
|
||||
if (it != cell.end())
|
||||
cell.erase(it);
|
||||
});
|
||||
};
|
||||
|
||||
auto fill_cells = [&](std::size_t i)
|
||||
{
|
||||
foreach_cells(i, [i](auto & cell){ cell.push_back(i); });
|
||||
};
|
||||
|
||||
auto find_collisions = [&](std::size_t i)
|
||||
{
|
||||
foreach_cells(i, [&, i](auto & cell){
|
||||
for (auto j : cell)
|
||||
if (i != j) check_collision(i, j);
|
||||
});
|
||||
};
|
||||
|
||||
auto handle_collision = [&](collision_event const & e)
|
||||
{
|
||||
clear_cells(e.i0);
|
||||
clear_cells(e.i1);
|
||||
|
||||
particles_[e.i0].pos += particles_[e.i0].vel * (e.t - particle_time[e.i0]);
|
||||
particles_[e.i1].pos += particles_[e.i1].vel * (e.t - particle_time[e.i1]);
|
||||
|
||||
// log::info() << "Collision " << e.i0 << " " << e.i1 << " " << geom::distance(particles_[e.i0].pos, particles_[e.i1].pos);
|
||||
|
||||
particle_time[e.i0] = e.t;
|
||||
particle_time[e.i1] = e.t;
|
||||
|
||||
auto n = geom::normalized(particles_[e.i0].pos - particles_[e.i1].pos);
|
||||
|
||||
float A = 0.5f * (1.f / particles_[e.i0].mass + 1.f / particles_[e.i1].mass);
|
||||
float B = dot(n, particles_[e.i0].vel - particles_[e.i1].vel);
|
||||
|
||||
auto BB = B / geom::length(particles_[e.i0].vel - particles_[e.i1].vel);
|
||||
(void)BB;
|
||||
auto np = n * (- B / A);
|
||||
|
||||
np *= std::exp(-10.f * dt);
|
||||
|
||||
particles_[e.i0].vel += np / particles_[e.i0].mass;
|
||||
particles_[e.i1].vel -= np / particles_[e.i1].mass;
|
||||
|
||||
fill_cells(e.i0);
|
||||
fill_cells(e.i1);
|
||||
};
|
||||
|
||||
// 0 - naive
|
||||
// 1 - optimized
|
||||
int algorithm = 1;
|
||||
|
||||
if (algorithm == 0)
|
||||
{
|
||||
collisions_ = 0;
|
||||
|
||||
while (time < dt)
|
||||
{
|
||||
particle_time.assign(particles_.size(), 0.f);
|
||||
|
||||
collision_event next{0, 0, std::numeric_limits<float>::infinity()};
|
||||
|
||||
for (std::size_t i = 0; i < particles_.size(); ++i)
|
||||
{
|
||||
for (std::size_t j = i + 1; j < particles_.size(); ++j)
|
||||
{
|
||||
auto rij = particles_[i].pos - particles_[j].pos;
|
||||
auto vij = particles_[i].vel - particles_[j].vel;
|
||||
auto R = particles_[i].radius + particles_[j].radius;
|
||||
|
||||
float B = geom::dot(rij, vij);
|
||||
|
||||
if (B >= 0.f) continue;
|
||||
|
||||
std::optional<collision_event> e;
|
||||
|
||||
float C = geom::dot(rij, rij) - geom::sqr(R);
|
||||
|
||||
if (C <= 0.f)
|
||||
{
|
||||
e = collision_event{i, j, 0.f};
|
||||
}
|
||||
else
|
||||
{
|
||||
float A = geom::dot(vij, vij);
|
||||
|
||||
float D = B * B - A * C;
|
||||
|
||||
if (D <= 0.f)
|
||||
continue;
|
||||
|
||||
float t = C / (- B + std::sqrt(D));
|
||||
e = collision_event{i, j, t};
|
||||
}
|
||||
|
||||
if (e && e->t < next.t)
|
||||
next = *e;
|
||||
}
|
||||
}
|
||||
|
||||
if (time + next.t < dt)
|
||||
{
|
||||
handle_collision(next);
|
||||
++collisions_;
|
||||
|
||||
for (std::size_t i = 0; i < particles_.size(); ++i)
|
||||
particles_[i].pos += particles_[i].vel * (next.t - particle_time[i]);
|
||||
|
||||
time += next.t;
|
||||
}
|
||||
else
|
||||
break;
|
||||
}
|
||||
|
||||
for (std::size_t i = 0; i < particles_.size(); ++i)
|
||||
particles_[i].pos += particles_[i].vel * (dt - time);
|
||||
}
|
||||
else if (algorithm == 1)
|
||||
{
|
||||
for (std::size_t i = 0; i < particles_.size(); ++i)
|
||||
fill_cells(i);
|
||||
|
||||
for (std::size_t i = 0; i < particles_.size(); ++i)
|
||||
find_collisions(i);
|
||||
|
||||
while (!events.empty())
|
||||
{
|
||||
collision_event e = *events.begin();
|
||||
time = e.t;
|
||||
|
||||
// if (e.t < time)
|
||||
// {
|
||||
// erase(events.begin());
|
||||
// continue;
|
||||
// }
|
||||
|
||||
// time = e.t;
|
||||
|
||||
// if (e.t < particle_time[e.i0])
|
||||
// log::info() << e.t << " < " << particle_time[e.i0];
|
||||
// if (e.t < particle_time[e.i1])
|
||||
// log::info() << e.t << " < " << particle_time[e.i1];
|
||||
|
||||
handle_collision(e);
|
||||
clear_collisions(e.i0);
|
||||
clear_collisions(e.i1);
|
||||
|
||||
find_collisions(e.i0);
|
||||
find_collisions(e.i1);
|
||||
|
||||
++collisions_;
|
||||
}
|
||||
|
||||
for (std::size_t i = 0; i < particles_.size(); ++i)
|
||||
{
|
||||
particles_[i].pos += particles_[i].vel * (dt - particle_time[i]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (false)
|
||||
{
|
||||
util::clock<> clock;
|
||||
|
||||
|
|
@ -555,8 +993,30 @@ struct myapp : app::app
|
|||
particles_[j].pos += kj * n;
|
||||
}
|
||||
|
||||
// inelastic collision without rotation
|
||||
if (l < R && dot(vij, n) < 0.f && false)
|
||||
// if (l < R)
|
||||
{
|
||||
float D = R - l;
|
||||
|
||||
float elasticity = 0.25f;
|
||||
|
||||
float A = (1.f / particles_[i].mass + 1.f / particles_[j].mass);
|
||||
float B = dot(n, vij);
|
||||
auto np = n * (- B / A) * (1.f + elasticity);
|
||||
|
||||
particles_[i].vel += np / particles_[i].mass;
|
||||
particles_[j].vel -= np / particles_[j].mass;
|
||||
|
||||
float ki = D * particles_[j].mass / (particles_[i].mass + particles_[j].mass);
|
||||
float kj = - D * particles_[i].mass / (particles_[i].mass + particles_[j].mass);
|
||||
|
||||
particles_[i].pos += ki * n;
|
||||
particles_[j].pos += kj * n;
|
||||
}
|
||||
|
||||
// collision replaced by force
|
||||
if (l < R)
|
||||
if (l < R && false)
|
||||
{
|
||||
auto f = 10000.f * n * (R - l);
|
||||
particles_[i].vel += dt * f / particles_[i].mass;
|
||||
|
|
@ -564,14 +1024,14 @@ struct myapp : app::app
|
|||
|
||||
auto vn = geom::dot(vij, n) * n;
|
||||
// vn *= 1.f - std::exp(- (1.f - l / R));
|
||||
vn *= 0.5f;
|
||||
// vn *= 0.5f;
|
||||
|
||||
particles_[i].vel -= vn * particles_[j].mass / (particles_[i].mass + particles_[j].mass);
|
||||
particles_[j].vel += vn * particles_[i].mass / (particles_[i].mass + particles_[j].mass);
|
||||
}
|
||||
|
||||
// Verlet collision
|
||||
if (l < R && false)
|
||||
if (l < R)
|
||||
{
|
||||
float D = (R - l) * 1.f;
|
||||
float ki = D * particles_[j].mass / (particles_[i].mass + particles_[j].mass);
|
||||
|
|
@ -583,6 +1043,19 @@ struct myapp : app::app
|
|||
particles_[i].pos += di;
|
||||
particles_[j].pos += dj;
|
||||
}
|
||||
|
||||
// impulse-based collision
|
||||
if (l < R && dot(n, particles_[i].vel - particles_[j].vel) < 0.f && false)
|
||||
{
|
||||
float bias = 0.5f;
|
||||
float constraint = l - R;
|
||||
float reduced_mass = 1.f / (1.f / particles_[i].mass + 1.f / particles_[j].mass);
|
||||
float lambda = (- geom::dot(n, particles_[i].vel - particles_[j].vel) - bias / dt * constraint) * reduced_mass;
|
||||
auto impulse = lambda * n;
|
||||
|
||||
particles_[i].vel += impulse / particles_[i].mass;
|
||||
particles_[j].vel -= impulse / particles_[j].mass;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -598,6 +1071,70 @@ struct myapp : app::app
|
|||
total_collisions_ += clock.count();
|
||||
}
|
||||
|
||||
// Verlet finalization
|
||||
if(false)
|
||||
for (auto & p : particles_)
|
||||
p.vel = (p.pos - p.old_pos) / dt;
|
||||
|
||||
// Iterative impulse-based collision
|
||||
if (false)
|
||||
{
|
||||
struct collision
|
||||
{
|
||||
std::size_t i, j;
|
||||
geom::vector<float, 2> normal;
|
||||
float value;
|
||||
float impulse = 0.f;
|
||||
};
|
||||
|
||||
std::vector<collision> collisions;
|
||||
|
||||
for (std::size_t i = 0; i < particles_.size(); ++i)
|
||||
{
|
||||
for (std::size_t j = i + 1; j < particles_.size(); ++j)
|
||||
{
|
||||
auto const r = particles_[i].pos - particles_[j].pos;
|
||||
|
||||
float const l = length(r);
|
||||
|
||||
float const R = particles_[i].radius + particles_[j].radius;
|
||||
|
||||
auto const n = r / l;
|
||||
|
||||
if (l < R)
|
||||
{
|
||||
collisions.push_back({i, j, n, l - R});
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for (std::size_t iteration = 0; iteration < SOLVE_ITERATIONS; ++iteration)
|
||||
{
|
||||
for (auto & c : collisions)
|
||||
{
|
||||
float bias = BIAS;
|
||||
|
||||
float reduced_mass = 1.f / (1.f / particles_[c.i].mass + 1.f / particles_[c.j].mass);
|
||||
|
||||
auto dv = geom::dot(c.normal, particles_[c.i].vel - particles_[c.j].vel);
|
||||
float elasticity = 0.75f;
|
||||
|
||||
float lambda = (- dv - bias / dt * c.value - dv * elasticity) * reduced_mass;
|
||||
float new_impulse = std::max(0.f, c.impulse + lambda);
|
||||
auto delta = new_impulse - c.impulse;
|
||||
c.impulse = new_impulse;
|
||||
|
||||
particles_[c.i].vel += delta * c.normal / particles_[c.i].mass;
|
||||
particles_[c.j].vel -= delta * c.normal / particles_[c.j].mass;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (false)
|
||||
for (auto & p : particles_)
|
||||
{
|
||||
p.pos += p.vel * dt;
|
||||
}
|
||||
|
||||
float Ep = 0.f;
|
||||
float Ek = 0.f;
|
||||
|
|
@ -618,6 +1155,8 @@ struct myapp : app::app
|
|||
omega += geom::det(particles_[i].mass * particles_[i].vel, particles_[i].pos - geom::point{0.f, 0.f});
|
||||
}
|
||||
|
||||
energy_ = Ek + Ep;
|
||||
|
||||
// log::info() << "Angular velocity: " << omega;
|
||||
|
||||
// log::info() << "Energy: " << Ek << " - " << (-Ep) << " = " << (Ek + Ep);
|
||||
|
|
@ -635,6 +1174,14 @@ struct myapp : app::app
|
|||
|
||||
painter_.circle(world_center, world_size, {255, 255, 255, 255}, 72);
|
||||
|
||||
// for (auto const & c : cells_)
|
||||
// {
|
||||
// float a = 0.25f;
|
||||
// if (!c.second.empty())
|
||||
// a = 0.5f;
|
||||
// painter_.rect({{{c.first[0], c.first[0] + 1.f}, {c.first[1], c.first[1] + 1.f}}}, gfx::to_coloru8(gfx::color_4f{1.f, 0.f, 0.f, a}));
|
||||
// }
|
||||
|
||||
for (auto & p : particles_)
|
||||
{
|
||||
// float c = 2.f / (std::exp(-p.T / 100000.f) + 1.f) - 1.f;
|
||||
|
|
@ -644,8 +1191,8 @@ struct myapp : app::app
|
|||
float s = window_size_[1] / camera_size_;
|
||||
|
||||
float r = std::max(p.radius * s, 1.5f) / s;
|
||||
painter_.circle(p.pos, r, {x, x, x, 255});
|
||||
painter_.line(p.pos, p.pos + r * geom::direction(p.angle), r / 16.f, {255, 0, 0, 255}, false);
|
||||
painter_.circle(p.pos, r, {x, x, x, 191});
|
||||
// painter_.line(p.pos, p.pos + r * geom::direction(p.angle), r / 16.f, {255, 0, 0, 255}, false);
|
||||
// painter_.circle(p.pos, r * 0.75f, {255, 255, 255, 255});
|
||||
}
|
||||
|
||||
|
|
@ -658,6 +1205,35 @@ struct myapp : app::app
|
|||
camera.box[2].max = 1.f;
|
||||
|
||||
painter_.render(camera.transform());
|
||||
|
||||
float rotation = 0.f;
|
||||
for (auto const & p : particles_)
|
||||
rotation += geom::det(p.pos - p.pos.zero(), p.vel) * p.mass;
|
||||
rotation_.push(rotation);
|
||||
|
||||
{
|
||||
gfx::painter::text_options opts;
|
||||
opts.scale = 3.f;
|
||||
opts.c = {0, 0, 0, 255};
|
||||
opts.x = gfx::painter::x_align::center;
|
||||
opts.y = gfx::painter::y_align::top;
|
||||
|
||||
auto put = [&](geom::point<float, 2> const & pos, std::string const & str)
|
||||
{
|
||||
return;
|
||||
painter_.text(pos, str, opts);
|
||||
painter_.text(pos + geom::vector{1.f, 0.f}, str, opts);
|
||||
};
|
||||
|
||||
put({width() / 2.f, 30.f}, util::to_string("ITERATIONS: ", SOLVE_ITERATIONS));
|
||||
put({width() / 2.f, 60.f}, BIAS == 1.f ? "BIAS: 1" : util::to_string("BIAS: 1/", std::round(1.f / BIAS)));
|
||||
|
||||
// painter_.text({10.f, 10.f}, util::to_string("Angular velocity: ", rotation_.average()), opts);
|
||||
// painter_.text({10.f, 30.f}, util::to_string("Collisions: ", collisions_, " = ", (collisions_ * 1.f / particles_.size() / particles_.size()), " N^2"), opts);
|
||||
// painter_.text({10.f, 50.f}, util::to_string("Energy: ", energy_), opts);
|
||||
}
|
||||
|
||||
painter_.render(geom::window_camera{width(), height()}.transform());
|
||||
}
|
||||
|
||||
~myapp()
|
||||
|
|
@ -676,7 +1252,7 @@ private:
|
|||
std::optional<geom::point<float, 2>> force_target_;
|
||||
|
||||
geom::point<float, 2> camera_center_ { 0.f, 0.f };
|
||||
float camera_size_ = 500.f;
|
||||
float camera_size_ = 150.f;
|
||||
float camera_ratio_ = 1.f;
|
||||
|
||||
std::optional<geom::point<int, 2>> mouse_;
|
||||
|
|
@ -691,6 +1267,12 @@ private:
|
|||
float total_collisions_ = 0.f;
|
||||
|
||||
audio::engine audio_;
|
||||
|
||||
util::moving_average<float> rotation_{300};
|
||||
int collisions_ = 0;
|
||||
float energy_ = 0.f;
|
||||
|
||||
bool paused_ = true;
|
||||
};
|
||||
|
||||
int main()
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue