#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace psemek; struct main_scene : app::ui_scene { main_scene(ui::controller & ui_controller); void update() override; void present() override; float simulation_radius = 1.f; float simulation_radius_tgt = 1.f; float potential = 1000.f; float step = 1e-4f; float multipole_threshold = 2.f; random::generator rng{random::device{}}; std::vector> points; std::vector> edges; std::vector degree; gfx::painter painter; util::clock<> clock; }; main_scene::main_scene(ui::controller & ui_controller) : app::ui_scene(ui_controller) { ui::default_element_factory element_factory; auto root = element_factory.make_screen(); auto wheel_event = std::make_shared(); wheel_event->on_mouse_wheel([this](ui::mouse_wheel const & event){ simulation_radius_tgt *= std::pow(0.8f, event.delta); return true; }); wheel_event->set_child(element_factory.make_screen()); root->add_child(wheel_event, ui::screen::x_policy::fill, ui::screen::y_policy::fill); auto panel = element_factory.make_window(""); panel->caption()->set_tagged_text("[bold][uline]Options[/uline][/bold]"); root->add_child(panel, ui::screen::x_policy::floating, ui::screen::y_policy::floating); auto layout = element_factory.make_grid_layout(); panel->set_child(layout); auto count_name_label = element_factory.make_label("Count:"); count_name_label->set_valign(ui::label::valignment::center); count_name_label->set_halign(ui::label::halignment::right); auto count_value_label = element_factory.make_label(""); count_value_label->set_valign(ui::label::valignment::center); count_value_label->set_halign(ui::label::halignment::center); auto count_slider = element_factory.make_slider(); count_slider->set_value_range({1, 5000}, false); count_slider->on_value_changed([this, count_value_label](int value){ count_value_label->set_text(util::to_string(value)); if (value < points.size()) points.resize(value); else { float radius = std::max(0.5f, std::sqrt((points.size() + value) / 1000.f)); random::uniform_sphere_point_distribution d({0.f, 0.f}, radius * 1.25f); while (points.size() < value) points.push_back(d(rng)); } }); count_slider->set_value(100); auto step_name_label = element_factory.make_label("Step:"); step_name_label->set_valign(ui::label::valignment::center); step_name_label->set_halign(ui::label::halignment::right); auto step_value_label = element_factory.make_label(""); step_value_label->set_valign(ui::label::valignment::center); step_value_label->set_halign(ui::label::halignment::center); auto step_slider = element_factory.make_slider(); step_slider->set_value_range({0, 12}, false); step_slider->on_value_changed([this, step_value_label](int value){ step = std::pow(2.f, value - 20); step_value_label->set_text(util::to_string(std::setprecision(1), std::scientific, step)); }); step_slider->set_value(11); auto precision_name_label = element_factory.make_label("Precision:"); precision_name_label->set_valign(ui::label::valignment::center); precision_name_label->set_halign(ui::label::halignment::right); auto precision_value_label = element_factory.make_label(""); precision_value_label->set_valign(ui::label::valignment::center); precision_value_label->set_halign(ui::label::halignment::center); auto precision_slider = element_factory.make_slider(); precision_slider->set_value_range({5, 50}, false); precision_slider->on_value_changed([this, precision_value_label](int value){ multipole_threshold = value * 0.1f; precision_value_label->set_text(util::to_string(std::setprecision(1), std::fixed, multipole_threshold)); }); precision_slider->set_value(15); layout->set_size(3, 3); layout->set_column_weight(0, 0.5f); layout->set_column_weight(1, 0.5f); layout->set(0, 0, count_name_label); layout->set(0, 1, count_value_label); layout->set(0, 2, count_slider); layout->set(1, 0, step_name_label); layout->set(1, 1, step_value_label); layout->set(1, 2, step_slider); layout->set(2, 0, precision_name_label); layout->set(2, 1, precision_value_label); layout->set(2, 2, precision_slider); ui::style style; style.font = ui::make_default_9x12_font(); style.bold_font = ui::make_default_10x12_bold_font(); style.text_scale = 1; style.bg_color = gfx::color_rgba{127, 127, 191, 255}; style.fg_color = gfx::color_rgba{127, 127, 255, 255}; style.action_color = gfx::color_rgba{0, 0, 255, 255}; style.highlight_color = gfx::color_rgba{0, 255, 255, 255}; style.border_width = 0; style.bevel_width = 0; root->set_style(std::make_shared(std::move(style))); set_ui(root); } struct node { static constexpr std::uint32_t null = -1; float size; float mass = 0.f; geom::point center = geom::point::zero(); std::uint32_t children[2][2] {{null, null}, {null, null}}; }; void main_scene::update() { ui_scene::update(); float const dt = clock.restart().count(); const int iterations = 1; auto const origin = geom::point::zero(); std::vector nodes; std::uint32_t root; auto add_node = util::recursive([&](auto && self, geom::box const & bbox, auto begin, auto end) -> std::uint32_t { if (begin == end) return node::null; std::uint32_t id = nodes.size(); nodes.emplace_back().size = bbox[0].length(); if (end - begin > 1) { auto middle_x = bbox[0].center(); auto middle_y = bbox[1].center(); auto mid_it = std::partition(begin, end, [middle_x](auto const & p){ return p[0] < middle_x; }); auto left_mid_it = std::partition(begin, mid_it, [middle_y](auto const & p){ return p[1] < middle_y; }); auto right_mid_it = std::partition(mid_it, end, [middle_y](auto const & p){ return p[1] < middle_y; }); nodes[id].children[0][0] = self(geom::box{{{bbox[0].min, middle_x}, {bbox[1].min, middle_y}}}, begin, left_mid_it); nodes[id].children[0][1] = self(geom::box{{{bbox[0].min, middle_x}, {middle_y, bbox[1].max}}}, left_mid_it, mid_it); nodes[id].children[1][0] = self(geom::box{{{middle_x, bbox[0].max}, {bbox[1].min, middle_y}}}, mid_it, right_mid_it); nodes[id].children[1][1] = self(geom::box{{{middle_x, bbox[0].max}, {middle_y, bbox[1].max}}}, right_mid_it, end); geom::vector sum{0.f, 0.f}; for (int x : {0, 1}) { for (int y : {0, 1}) { auto cid = nodes[id].children[x][y]; if (cid != node::null) { nodes[id].mass += nodes[cid].mass; sum += (nodes[cid].center - origin) * nodes[cid].mass; } } } nodes[id].center = origin + sum / nodes[id].mass; } else { nodes[id].mass = 1.f; nodes[id].center = *begin; } return id; }); { float max_size = 0.f; for (auto const & p : points) geom::make_max(max_size, geom::distance(p, origin)); prof::profiler prof("tree"); root = add_node(geom::box{{{-max_size, max_size}, {-max_size, max_size}}}, points.begin(), points.end()); } auto force_at = util::recursive([&](auto && self, auto const & p, auto id){ auto const & n = nodes[id]; if (n.center == p) return geom::vector{0.f, 0.f}; auto d = p - n.center; auto l = geom::length(d); if (n.mass == 1.f || l > n.size * multipole_threshold) return n.mass * d / std::pow(l, 3.f); geom::vector sum{0.f, 0.f}; for (int x : {0, 1}) { for (int y : {0, 1}) { auto cid = nodes[id].children[x][y]; if (cid != node::null) sum += self(p, cid); } } return sum; }); std::vector> delta; for (int iteration = 0; iteration < iterations; ++iteration) { prof::profiler prof("descent"); delta.assign(points.size(), geom::vector::zero()); for (std::size_t i = 0; i < points.size(); ++i) delta[i] += 2.f * potential * (origin - points[i]); // Exact O(n^2) method // for (std::size_t i = 0; i < points.size(); ++i) // { // for (std::size_t j = i + 1; j < points.size(); ++j) // { // auto d = points[i] - points[j]; // d /= std::pow(geom::length(d), 3.f); // delta[i] += d; // delta[j] -= d; // } // } // Barnes-Hut algorithm for (std::size_t i = 0; i < points.size(); ++i) delta[i] += force_at(points[i], root); for (std::size_t i = 0; i < points.size(); ++i) points[i] += delta[i] * step * dt; } simulation_radius += (simulation_radius_tgt - simulation_radius) * (1.f - std::exp(-10.f * dt)); prof::profiler prof("delaunay"); auto const dcel = cg::delaunay(geom::robust, points.begin(), points.end()); edges.clear(); degree.assign(points.size(), 0); for (std::size_t i = 0; i < dcel.edges.size(); ++i) { auto const e = dcel.edge(i); auto const v0 = e.origin().index(); auto const v1 = e.twin().origin().index(); if (v0 < v1) edges.push_back({v0, v1}); degree[v0]++; } } void main_scene::present() { gl::ClearColor(1.f, 1.f, 1.f, 0.f); gl::Clear(gl::COLOR_BUFFER_BIT | gl::DEPTH_BUFFER_BIT); geom::box view_area; { view_area[0] = {-simulation_radius, simulation_radius}; view_area[1] = {-simulation_radius, simulation_radius}; view_area[2] = {-1.f, 1.f}; float const aspect_ratio = width() * 1.f / height(); float extra_y = view_area[1].length() * 0.1f; view_area[1].min -= extra_y; view_area[1].max += extra_y; float extra_x = view_area[1].length() * aspect_ratio - view_area[0].length(); view_area[0].min -= extra_x / 2.f; view_area[0].max += extra_x / 2.f; } geom::matrix const transform = geom::orthographic{view_area}.homogeneous_matrix(); float const pixel_size = view_area[0].length() / width(); for (auto const & e : edges) painter.line(points[e[0]], points[e[1]], 2.f * pixel_size, {127, 127, 127, 255}, false); for (std::size_t i = 0; i < points.size(); ++i) { gfx::color_rgba color{0, 0, 0, 255}; int quality = 6; if (degree[i] == 5) { quality = 3; color = {255, 0, 0, 255}; } else if (degree[i] == 7) { quality = 4; color = {0, 0, 255, 255}; } painter.circle(points[i], 6.f * pixel_size, color, quality); } painter.render(transform); ui_scene::present(); } struct electron_crystal_app : app::app { electron_crystal_app(); ~electron_crystal_app(); async::event_loop event_loop; ui::controller ui_controller; }; electron_crystal_app::electron_crystal_app() : app("Electron crystal simulation", 4) , ui_controller(&event_loop) { push_scene(std::make_shared(ui_controller)); } electron_crystal_app::~electron_crystal_app() { prof::dump(); } int main() { return app::main(); }