#include #include #include #include #include #include #include #include #include /* No optimizations, 125: [2020 Aug 29 09:54:27.348][ main][ info] Avg forces time: 0.000641577 [2020 Aug 29 09:54:27.348][ main][ info] Avg collision time: 0.000789326 No optimizations, 250: [2020 Aug 29 09:49:48.476][ main][ info] Avg forces time: 0.00257929 [2020 Aug 29 09:49:48.476][ main][ info] Avg collision time: 0.00314664 No optimizations, 500: [2020 Aug 29 09:44:02.200][ main][ info] Avg forces time: 0.00299429 [2020 Aug 29 09:44:02.200][ main][ info] Avg collision time: 0.00362124 No optimizations, 1000: [2020 Aug 29 09:42:34.143][ main][ info] Avg forces time: 0.0072538 [2020 Aug 29 09:42:34.143][ main][ info] Avg collision time: 0.00873274 No optimizations, 2000: [2020 Aug 29 09:44:35.334][ main][ info] Avg forces time: 0.029108 [2020 Aug 29 09:44:35.334][ main][ info] Avg collision time: 0.0345713 No optimizations, 4000: [2020 Aug 29 09:50:15.647][ main][ info] Avg forces time: 0.118145 [2020 Aug 29 09:50:15.647][ main][ info] Avg collision time: 0.1402 */ using namespace psemek; struct particle { geom::point pos; geom::vector vel; float radius; float mass; float T = 0.f; }; struct myapp : app::app { myapp() : app("Test app") , rng_{std::random_device{}()} { gl::ClearColor(1.f, 1.f, 1.f, 1.f); std::uniform_real_distribution d{-200.f, 200.f}; std::uniform_real_distribution rr{0.5f, 2.f}; for (int i = 0; i < 1000; ++i) { float r = rr(rng_); float m = r * r * r; float v = 0.1f; geom::point p; while (true) { p = {((i % 2) ? -300.f : 300.f) + d(rng_), d(rng_)}; if (std::all_of(particles_.begin(), particles_.end(), [&](particle const & q){ return geom::distance(q.pos, p) > q.radius + r; })) break; } particles_.push_back({p, {v * d(rng_), v * d(rng_)}, r, m}); // particles_.push_back({{ 400.f + d(rng), d(rng)}, {0.f, 0.f}, 1.f, 1.f}); // particles_.push_back({{d(rng), d(rng)}, {0.f, 0.f}, 1.f, 1.f}); // particles_.push_back({{d(rng), d(rng)}, {0.f, 0.f}, 1.f, 1.f}); } } void on_resize(int width, int height) override { gl::Viewport(0, 0, width, height); window_size_ = {width, height}; camera_ratio_ = static_cast(width) / height; } void on_mouse_wheel(int delta) override { camera_size_ *= std::pow(0.8f, delta); } void on_left_button_down() override { if (mouse_) { geom::scale const flip_y({1.f, -1.f}); float const scale = camera_size_ / window_size_[1]; geom::point const screen_center { window_size_[0] / 2, window_size_[1] / 2 }; auto target = camera_center_ + flip_y(geom::cast(*mouse_ - screen_center) * scale); geom::point pos{100.f, 0.f}; geom::vector vel = geom::normalized(target - pos) * 500.f; particles_.push_back({pos, vel, 1.f, 1.f}); } } void on_right_button_down() override { camera_drag_ = mouse_; } void on_right_button_up() override { camera_drag_ = std::nullopt; } void on_mouse_move(int x, int y, int, int) override { mouse_ = {x, y}; if (camera_drag_) { geom::scale const flip_y({1.f, -1.f}); float const scale = camera_size_ / window_size_[1]; camera_center_ += flip_y(geom::cast(*camera_drag_ - *mouse_) * scale); camera_drag_ = mouse_; } } void update() override { // std::shuffle(particles_.begin(), particles_.end(), rng_); float const dt = 0.01f; float const G = 100.f; for (std::size_t step = 0; step < 1; ++step) { { util::clock<> clock; 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); auto const f = G * particles_[i].mass * particles_[j].mass * r / std::pow(l, 3.f); particles_[i].vel -= f * dt / particles_[i].mass; particles_[j].vel += f * dt / particles_[j].mass; } } total_forces_ += clock.count(); } for (auto & p : particles_) { p.pos += p.vel * dt; } { util::clock<> clock; for (std::size_t iteration = 0; iteration < 10; ++iteration) { 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) { auto f = E * n * std::pow(l / R, -2.f); particles_[i].vel += f * dt / particles_[i].mass; particles_[j].vel -= f * dt / particles_[j].mass; auto const vij = particles_[i].vel - particles_[j].vel; auto dp = - K * dt * n * dot(n, vij) * 2.f / (1.f / particles_[i].mass + 1.f / particles_[j].mass); float Ei0 = geom::length_sqr(particles_[i].vel) * particles_[i].mass * 0.5f; particles_[i].vel += dp / particles_[i].mass; float Ei1 = geom::length_sqr(particles_[i].vel) * particles_[i].mass * 0.5f; particles_[i].T += Ei0 - Ei1; float Ej0 = geom::length_sqr(particles_[j].vel) * particles_[j].mass * 0.5f; particles_[j].vel -= dp / particles_[j].mass; float Ej1 = geom::length_sqr(particles_[j].vel) * particles_[j].mass * 0.5f; particles_[j].T += Ej0 - Ej1; float dT = particles_[i].T - particles_[j].T; particles_[i].T -= C * dT * dt; particles_[j].T += C * dT * dt; } */ auto const vij = particles_[i].vel - particles_[j].vel; if (l < R && dot(vij, n) < 0.f) { auto const d = n * (R - l) / 2.f * 0.5f;// * (particles_[i].mass + particles_[j].mass); particles_[i].pos += d; particles_[j].pos -= d; auto dp = - n * dot(n, vij) * 2.f / (1.f / particles_[i].mass + 1.f / particles_[j].mass); particles_[i].vel += dp / particles_[i].mass; particles_[j].vel -= dp / particles_[j].mass; } } } } total_collisions_ += clock.count(); } float Ep = 0.f; float Ek = 0.f; for (std::size_t i = 0; i < particles_.size(); ++i) { Ek += geom::length_sqr(particles_[i].vel) * particles_[i].mass / 2.f; for (std::size_t j = i + 1; j < particles_.size(); ++j) { Ep -= G * particles_[i].mass * particles_[j].mass / distance(particles_[i].pos, particles_[j].pos); } } log::info() << "Energy: " << Ek << " - " << (-Ep) << " = " << (Ek + Ep); ++frame_count_; } } void present() override { gl::Clear(gl::COLOR_BUFFER_BIT); for (auto & p : particles_) { float c = 2.f / (std::exp(-p.T / 100000.f) + 1.f) - 1.f; painter_.circle(p.pos, p.radius, {static_cast(c * 255.f), 0, 0, 255}); } geom::orthographic_camera camera; camera.box[0].min = camera_center_[0] - camera_size_ * camera_ratio_ / 2.f; camera.box[0].max = camera_center_[0] + camera_size_ * camera_ratio_ / 2.f; camera.box[1].min = camera_center_[1] - camera_size_ / 2.f; camera.box[1].max = camera_center_[1] + camera_size_ / 2.f; camera.box[2].min = -1.f; camera.box[2].max = 1.f; painter_.render(camera.transform()); } ~myapp() { log::info() << "Avg forces time: " << (total_forces_ / frame_count_); log::info() << "Avg collision time: " << (total_collisions_ / frame_count_); } private: std::default_random_engine rng_; geom::vector window_size_; geom::point camera_center_ { 0.f, 0.f }; float camera_size_ = 500.f; float camera_ratio_ = 1.f; std::optional> mouse_; std::optional> camera_drag_; gfx::painter painter_; std::vector particles_; int frame_count_ = 0; float total_forces_ = 0.f; float total_collisions_ = 0.f; }; int main() { return app::main(); }