reformat

benchmark/state_space.cpp

@@ -177,7 +177,7 @@ static auto explore_rush_hour_puzzle_space(benchmark::State& state) -> void
     constexpr std::tuple<uint8_t, uint8_t, uint8_t, uint8_t> target_block_pos_range = {
         0,
         goal_y,
-        board_width - 1,
+        goal_x,
         goal_y
     };
 

flake.nix

@@ -113,7 +113,7 @@ rec {
             abbr -a run "${buildRelease} && ./cmake-build-release/masssprings"
             abbr -a run-clusters "${buildRelease} && ./cmake-build-release/masssprings --output=clusters.puzzle --space=rh --w=6 --h=6 --gx=4 --gy=2 --blocks=4"
             abbr -a runtests "${buildDebug} && ./cmake-build-debug/tests"
-            abbr -a runbenchs "mv benchs.json benchs.old.json; ${buildRelease} && sudo cpupower frequency-set --governor performance && ./cmake-build-release/benchmarks --benchmark_out=benchs.json --benchmark_out_format=console; sudo cpupower frequency-set --governor powersave"
+            abbr -a runbenchs "mv -f benchs.json benchs.old.json; ${buildRelease} && sudo cpupower frequency-set --governor performance && ./cmake-build-release/benchmarks --benchmark_out=benchs.json --benchmark_out_format=console; sudo cpupower frequency-set --governor powersave"
             abbr -a rungdb "${buildDebug} && gdb --tui ./cmake-build-debug/masssprings"
             abbr -a runvalgrind "${buildDebug} && valgrind --leak-check=full --show-reachable=no --show-leak-kinds=definite,indirect,possible --track-origins=no --suppressions=valgrind.supp --log-file=valgrind.log ./cmake-build-debug/masssprings && cat valgrind.log"
             abbr -a runperf "${buildRelease} && perf record -g ./cmake-build-release/masssprings && hotspot ./perf.data"

include/octree.hpp

@@ -1,6 +1,8 @@
 #ifndef OCTREE_HPP_
 #define OCTREE_HPP_
 
+#include "util.hpp"
+
 #include <array>
 #include <vector>
 
@@ -21,7 +23,7 @@ class octree
         bool leaf = true;
 
     public:
-        [[nodiscard]] auto child_count() const -> int;
+        node(const Vector3& _box_min, const Vector3& _box_max) : box_min(_box_min), box_max(_box_max) {}
     };
 
 public:
@@ -29,23 +31,88 @@ public:
 
     std::vector<node> nodes;
 
+    // This approach is actually slower than the array of nodes
+    // beacuse we access all the attributes in the same function
+
+    // std::vector<Vector3> mass_centers;
+    // std::vector<float> mass_totals;
+    // std::vector<Vector3> box_mins;
+    // std::vector<Vector3> box_maxs;
+    // std::vector<std::array<int, 8>> childrens;
+    // std::vector<int> mass_ids;
+    // std::vector<uint8_t> leafs; // bitpacked std::vector<bool> is a lot slower
+
 public:
     octree() = default;
 
-    // octree(const octree& copy) = delete;
+    octree(const octree& copy) = delete;
-    // auto operator=(const octree& copy) -> octree& = delete;
+    auto operator=(const octree& copy) -> octree& = delete;
-    // octree(octree&& move) = delete;
+    octree(octree&& move) = delete;
-    // auto operator=(octree&& move) -> octree& = delete;
+    auto operator=(octree&& move) -> octree& = delete;
 
 public:
-    [[nodiscard]] auto get_octant(int node_idx, const Vector3& pos) const -> int;
+    auto clear() -> void;
-    [[nodiscard]] auto get_child_bounds(int node_idx, int octant) const
+    auto reserve(size_t count) -> void;
-        -> std::pair<Vector3, Vector3>;
+    [[nodiscard]] auto empty() const -> bool;
-    auto create_empty_leaf(const Vector3& box_min, const Vector3& box_max) -> int;
+    [[nodiscard]] INLINE static inline auto get_octant(const Vector3& box_min,
+                                                       const Vector3& box_max,
+                                                       const Vector3& pos) -> int;
+    [[nodiscard]] INLINE static inline auto get_child_bounds(const Vector3& box_min,
+                                                             const Vector3& box_max,
+                                                             int octant) -> std::pair<Vector3, Vector3>;
+    INLINE inline auto create_empty_leaf(const Vector3& box_min, const Vector3& box_max) -> int;
+    [[nodiscard]] auto get_child_count(int node_idx) const -> int;
     auto insert(int node_idx, int mass_id, const Vector3& pos, float mass, int depth) -> void;
     static auto build_octree(octree& t, const std::vector<Vector3>& positions) -> void;
 
     [[nodiscard]] auto calculate_force(int node_idx, const Vector3& pos) const -> Vector3;
 };
 
+INLINE inline auto octree::get_octant(const Vector3& box_min, const Vector3& box_max, const Vector3& pos) -> int
+{
+    auto [cx, cy, cz] = (box_min + box_max) / 2.0f;
+
+    // The octant is encoded as a 3-bit integer "zyx". The node area is split
+    // along all 3 axes, if a position is right of an axis, this bit is set to 1.
+    // If a position is right of the x-axis and y-axis and left of the z-axis, the
+    // encoded octant is "011".
+    return (pos.x >= cx) | ((pos.y >= cy) << 1) | ((pos.z >= cz) << 2);
+}
+
+INLINE inline auto octree::get_child_bounds(const Vector3& box_min,
+                                            const Vector3& box_max,
+                                            const int octant) -> std::pair<Vector3, Vector3>
+{
+    auto [cx, cy, cz] = (box_min + box_max) / 2.0f;
+
+    Vector3 min = Vector3Zero();
+    Vector3 max = Vector3Zero();
+
+    // If (octant & 1), the octant is to the right of the node region's x-axis
+    // (see GetOctant). This means the left bound is the x-axis and the right
+    // bound the node's region max.
+    min.x = octant & 1 ? cx : box_min.x;
+    max.x = octant & 1 ? box_max.x : cx;
+    min.y = octant & 2 ? cy : box_min.y;
+    max.y = octant & 2 ? box_max.y : cy;
+    min.z = octant & 4 ? cz : box_min.z;
+    max.z = octant & 4 ? box_max.z : cz;
+
+    return std::make_pair(min, max);
+}
+
+INLINE inline auto octree::create_empty_leaf(const Vector3& box_min, const Vector3& box_max) -> int
+{
+    nodes.emplace_back(box_min, box_max);
+    return static_cast<int>(nodes.size() - 1);
+
+    // mass_centers.emplace_back(Vector3Zero());
+    // mass_totals.emplace_back(0);
+    // box_mins.emplace_back(box_min);
+    // box_maxs.emplace_back(box_max);
+    // childrens.push_back({-1, -1, -1, -1, -1, -1, -1, -1});
+    // mass_ids.emplace_back(-1);
+    // leafs.emplace_back(true);
+}
+
 #endif

src/main.cpp

@@ -169,7 +169,12 @@ auto rush_hour_puzzle_space() -> int
         puzzle::block(0, 0, 1, 3, false, false)
     };
     const puzzle::block target_block = puzzle::block(0, 0, 2, 1, true, false);
-    const std::tuple<uint8_t, uint8_t, uint8_t, uint8_t> target_block_pos_range = {0, goal_y, board_width - 1, goal_y};
+    const std::tuple<uint8_t, uint8_t, uint8_t, uint8_t> target_block_pos_range = {
+        0,
+        goal_y,
+        board_width - target_block.get_width() - 1,
+        goal_y
+    };
 
     infoln("Exploring Rush-Hour puzzle space:");
     infoln("- Size:       {}x{}", board_width, board_height);

src/mass_spring_system.cpp

@@ -4,55 +4,6 @@
 #include <cfloat>
 #include <cstring>
 
-auto mass_spring_system::calculate_spring_force(const size_t s) -> void
-{
-    const spring _s = springs[s];
-    const Vector3 a_pos = positions[_s.a];
-    const Vector3 b_pos = positions[_s.b];
-    const Vector3 a_vel = velocities[_s.a];
-    const Vector3 b_vel = velocities[_s.b];
-
-    const Vector3 delta_pos = a_pos - b_pos;
-    const Vector3 delta_vel = a_vel - b_vel;
-
-    const float sq_len = Vector3DotProduct(delta_pos, delta_pos);
-    const float inv_len = 1.0f / sqrt(sq_len);
-    const float len = sq_len * inv_len;
-
-    const float hooke = SPRING_CONSTANT * (len - REST_LENGTH);
-    const float dampening = DAMPENING_CONSTANT * Vector3DotProduct(delta_vel, delta_pos) * inv_len;
-
-    const Vector3 a_force = Vector3Scale(delta_pos, -(hooke + dampening) * inv_len);
-    const Vector3 b_force = a_force * -1.0f;
-
-    forces[_s.a] += a_force;
-    forces[_s.b] += b_force;
-}
-
-auto mass_spring_system::integrate_velocity(const size_t m, const float dt) -> void
-{
-    const Vector3 acc = forces[m] / MASS;
-    velocities[m] += acc * dt;
-}
-
-auto mass_spring_system::integrate_position(const size_t m, const float dt) -> void
-{
-    previous_positions[m] = positions[m];
-    positions[m] += velocities[m] * dt;
-}
-
-auto mass_spring_system::verlet_update(const size_t m, const float dt) -> void
-{
-    const Vector3 acc = (forces[m] / MASS) * dt * dt;
-    const Vector3 pos = positions[m];
-
-    Vector3 delta_pos = pos - previous_positions[m];
-    delta_pos *= 1.0 - VERLET_DAMPENING; // Minimal dampening
-
-    positions[m] += delta_pos + acc;
-    previous_positions[m] = pos;
-}
-
 auto mass_spring_system::clear() -> void
 {
     positions.clear();
@@ -60,7 +11,7 @@ auto mass_spring_system::clear() -> void
     velocities.clear();
     forces.clear();
     springs.clear();
-    tree.nodes.clear();
+    tree.clear();
 }
 
 auto mass_spring_system::add_mass() -> void
@@ -87,12 +38,12 @@ auto mass_spring_system::add_spring(size_t a, size_t b) -> void
     const Vector3& mass_b = positions[b];
 
     Vector3 offset{static_cast<float>(GetRandomValue(-100, 100)),
-                   static_cast<float>(GetRandomValue(-100, 100)),
+        static_cast<float>(GetRandomValue(-100, 100)),
-                   static_cast<float>(GetRandomValue(-100, 100))};
+        static_cast<float>(GetRandomValue(-100, 100))};
     offset = Vector3Normalize(offset) * REST_LENGTH;
 
     // If the offset moves the mass closer to the current center of mass, flip it
-    if (!tree.nodes.empty()) {
+    if (!tree.empty()) {
         const Vector3 mass_center_direction =
             Vector3Subtract(positions[a], tree.nodes[0].mass_center);
         const float mass_center_distance = Vector3Length(mass_center_direction);
@@ -114,19 +65,44 @@ auto mass_spring_system::add_spring(size_t a, size_t b) -> void
 
 auto mass_spring_system::clear_forces() -> void
 {
-#ifdef TRACY
+    #ifdef TRACY
     ZoneScoped;
-#endif
+    #endif
 
     memset(forces.data(), 0, forces.size() * sizeof(Vector3));
 }
 
+auto mass_spring_system::calculate_spring_force(const size_t s) -> void
+{
+    const spring _s = springs[s];
+    const Vector3 a_pos = positions[_s.a];
+    const Vector3 b_pos = positions[_s.b];
+    const Vector3 a_vel = velocities[_s.a];
+    const Vector3 b_vel = velocities[_s.b];
+
+    const Vector3 delta_pos = a_pos - b_pos;
+    const Vector3 delta_vel = a_vel - b_vel;
+
+    const float sq_len = Vector3DotProduct(delta_pos, delta_pos);
+    const float inv_len = 1.0f / sqrt(sq_len);
+    const float len = sq_len * inv_len;
+
+    const float hooke = SPRING_CONSTANT * (len - REST_LENGTH);
+    const float dampening = DAMPENING_CONSTANT * Vector3DotProduct(delta_vel, delta_pos) * inv_len;
+
+    const Vector3 a_force = Vector3Scale(delta_pos, -(hooke + dampening) * inv_len);
+    const Vector3 b_force = a_force * -1.0f;
+
+    forces[_s.a] += a_force;
+    forces[_s.b] += b_force;
+}
+
 auto mass_spring_system::calculate_spring_forces(
     const std::optional<BS::thread_pool<>* const> thread_pool) -> void
 {
-#ifdef TRACY
+    #ifdef TRACY
     ZoneScoped;
-#endif
+    #endif
 
     const auto solve_spring_force = [&](const int i) { calculate_spring_force(i); };
 
@@ -144,9 +120,9 @@ auto mass_spring_system::calculate_spring_forces(
 auto mass_spring_system::calculate_repulsion_forces(
     const std::optional<BS::thread_pool<>* const> thread_pool) -> void
 {
-#ifdef TRACY
+    #ifdef TRACY
     ZoneScoped;
-#endif
+    #endif
 
     const auto solve_octree = [&](const int i)
     {
@@ -166,12 +142,36 @@ auto mass_spring_system::calculate_repulsion_forces(
     }
 }
 
+auto mass_spring_system::integrate_velocity(const size_t m, const float dt) -> void
+{
+    const Vector3 acc = forces[m] / MASS;
+    velocities[m] += acc * dt;
+}
+
+auto mass_spring_system::integrate_position(const size_t m, const float dt) -> void
+{
+    previous_positions[m] = positions[m];
+    positions[m] += velocities[m] * dt;
+}
+
+auto mass_spring_system::verlet_update(const size_t m, const float dt) -> void
+{
+    const Vector3 acc = (forces[m] / MASS) * dt * dt;
+    const Vector3 pos = positions[m];
+
+    Vector3 delta_pos = pos - previous_positions[m];
+    delta_pos *= 1.0 - VERLET_DAMPENING; // Minimal dampening
+
+    positions[m] += delta_pos + acc;
+    previous_positions[m] = pos;
+}
+
 auto mass_spring_system::update(const float dt,
                                 const std::optional<BS::thread_pool<>* const> thread_pool) -> void
 {
-#ifdef TRACY
+    #ifdef TRACY
     ZoneScoped;
-#endif
+    #endif
 
     const auto update = [&](const int i) { verlet_update(i, dt); };
 

src/octree.cpp

@@ -4,10 +4,43 @@
 #include <cfloat>
 #include <raymath.h>
 
-auto octree::node::child_count() const -> int
+auto octree::clear() -> void
+{
+    nodes.clear();
+
+    // mass_centers.clear();
+    // mass_totals.clear();
+    // box_mins.clear();
+    // box_maxs.clear();
+    // childrens.clear();
+    // mass_ids.clear();
+    // leafs.clear();
+}
+
+auto octree::reserve(const size_t count) -> void
+{
+    nodes.reserve(count);
+
+    // mass_centers.reserve(count);
+    // mass_totals.reserve(count);
+    // box_mins.reserve(count);
+    // box_maxs.reserve(count);
+    // childrens.reserve(count);
+    // mass_ids.reserve(count);
+    // leafs.reserve(count);
+}
+
+auto octree::empty() const -> bool
+{
+    return nodes.empty();
+
+    // return mass_centers.empty();
+}
+
+auto octree::get_child_count(const int node_idx) const -> int
 {
     int child_count = 0;
-    for (const int child : children) {
+    for (const int child : nodes[node_idx].children) {
         if (child != -1) {
             ++child_count;
         }
@@ -15,88 +48,33 @@ auto octree::node::child_count() const -> int
     return child_count;
 }
 
-auto octree::get_octant(const int node_idx, const Vector3& pos) const -> int
+auto octree::insert(const int node_idx,
-{
+                    const int mass_id,
-    const node& n = nodes[node_idx];
+                    const Vector3& pos,
-    auto [cx, cy, cz] = Vector3((n.box_min.x + n.box_max.x) / 2.0f, (n.box_min.y + n.box_max.y) / 2.0f,
+                    const float mass,
-                                (n.box_min.z + n.box_max.z) / 2.0f);
-
-    // The octant is encoded as a 3-bit integer "zyx". The node area is split
-    // along all 3 axes, if a position is right of an axis, this bit is set to 1.
-    // If a position is right of the x-axis and y-axis and left of the z-axis, the
-    // encoded octant is "011".
-    int octant = 0;
-    if (pos.x >= cx) {
-        octant |= 1;
-    }
-    if (pos.y >= cy) {
-        octant |= 2;
-    }
-    if (pos.z >= cz) {
-        octant |= 4;
-    }
-
-    return octant;
-}
-
-auto octree::get_child_bounds(const int node_idx, const int octant) const -> std::pair<Vector3, Vector3>
-{
-    const node& n = nodes[node_idx];
-    auto [cx, cy, cz] = Vector3((n.box_min.x + n.box_max.x) / 2.0f, (n.box_min.y + n.box_max.y) / 2.0f,
-                                (n.box_min.z + n.box_max.z) / 2.0f);
-
-    Vector3 min = Vector3Zero();
-    Vector3 max = Vector3Zero();
-
-    // If (octant & 1), the octant is to the right of the node region's x-axis
-    // (see GetOctant). This means the left bound is the x-axis and the right
-    // bound the node's region max.
-    min.x = octant & 1 ? cx : n.box_min.x;
-    max.x = octant & 1 ? n.box_max.x : cx;
-    min.y = octant & 2 ? cy : n.box_min.y;
-    max.y = octant & 2 ? n.box_max.y : cy;
-    min.z = octant & 4 ? cz : n.box_min.z;
-    max.z = octant & 4 ? n.box_max.z : cz;
-
-    return std::make_pair(min, max);
-}
-
-auto octree::create_empty_leaf(const Vector3& box_min, const Vector3& box_max) -> int
-{
-    node n;
-    n.box_min = box_min;
-    n.box_max = box_max;
-    nodes.emplace_back(n);
-
-    return static_cast<int>(nodes.size() - 1);
-}
-
-auto octree::insert(const int node_idx, const int mass_id, const Vector3& pos, const float mass,
                     const int depth) -> void
 {
-    // infoln("Inserting position ({}, {}, {}) into octree at node {} (depth {})", pos.x, pos.y,
-    // pos.z, node_idx, depth);
     if (depth > MAX_DEPTH) {
-        throw std::runtime_error(std::format("MAX_DEPTH! node={} box_min=({},{},{}) box_max=({},{},{}) pos=({},{},{})",
+        throw std::runtime_error("Octree insertion reachead max depth");
-                                             node_idx, nodes[node_idx].box_min.x, nodes[node_idx].box_min.y,
-                                             nodes[node_idx].box_min.z, nodes[node_idx].box_max.x,
-                                             nodes[node_idx].box_max.y, nodes[node_idx].box_max.z, pos.x, pos.y,
-                                             pos.z));
     }
 
     // NOTE: Do not store a nodes[node_idx] reference as the nodes vector might reallocate during
     // this function
 
     // We can place the particle in the empty leaf
-    if (nodes[node_idx].leaf && nodes[node_idx].mass_id == -1) {
+    const bool leaf = nodes[node_idx].leaf;
+    if (leaf && nodes[node_idx].mass_id == -1) {
         nodes[node_idx].mass_id = mass_id;
         nodes[node_idx].mass_center = pos;
         nodes[node_idx].mass_total = mass;
         return;
     }
 
+    const Vector3& box_min = nodes[node_idx].box_min;
+    const Vector3& box_max = nodes[node_idx].box_max;
+
     // The leaf is occupied, we need to subdivide
-    if (nodes[node_idx].leaf) {
+    if (leaf) {
         const int existing_id = nodes[node_idx].mass_id;
         const Vector3 existing_pos = nodes[node_idx].mass_center;
         const float existing_mass = nodes[node_idx].mass_total;
@@ -119,13 +97,13 @@ auto octree::insert(const int node_idx, const int mass_id, const Vector3& pos, c
         // Convert the leaf to an internal node
         nodes[node_idx].mass_id = -1;
         nodes[node_idx].leaf = false;
-        nodes[node_idx].mass_total = 0.0;
         nodes[node_idx].mass_center = Vector3Zero();
+        nodes[node_idx].mass_total = 0.0;
 
         // Re-insert the existing mass into a new empty leaf (see above)
-        const int oct = get_octant(node_idx, existing_pos);
+        const int oct = get_octant(box_min, box_max, existing_pos);
         if (nodes[node_idx].children[oct] == -1) {
-            const auto& [min, max] = get_child_bounds(node_idx, oct);
+            const auto& [min, max] = get_child_bounds(box_min, box_max, oct);
             const int child_idx = create_empty_leaf(min, max);
             nodes[node_idx].children[oct] = child_idx;
         }
@@ -133,9 +111,9 @@ auto octree::insert(const int node_idx, const int mass_id, const Vector3& pos, c
     }
 
     // Insert the new mass
-    const int oct = get_octant(node_idx, pos);
+    const int oct = get_octant(box_min, box_max, pos);
     if (nodes[node_idx].children[oct] == -1) {
-        const auto& [min, max] = get_child_bounds(node_idx, oct);
+        const auto& [min, max] = get_child_bounds(box_min, box_max, oct);
         const int child_idx = create_empty_leaf(min, max);
         nodes[node_idx].children[oct] = child_idx;
     }
@@ -143,9 +121,7 @@ auto octree::insert(const int node_idx, const int mass_id, const Vector3& pos, c
 
     // Update the center of mass
     const float new_mass = nodes[node_idx].mass_total + mass;
-    nodes[node_idx].mass_center.x = (nodes[node_idx].mass_center.x * nodes[node_idx].mass_total + pos.x) / new_mass;
+    nodes[node_idx].mass_center = (nodes[node_idx].mass_center * nodes[node_idx].mass_total + pos) / new_mass;
-    nodes[node_idx].mass_center.y = (nodes[node_idx].mass_center.y * nodes[node_idx].mass_total + pos.y) / new_mass;
-    nodes[node_idx].mass_center.z = (nodes[node_idx].mass_center.z * nodes[node_idx].mass_total + pos.z) / new_mass;
     nodes[node_idx].mass_total = new_mass;
 }
 
@@ -155,8 +131,8 @@ auto octree::build_octree(octree& t, const std::vector<Vector3>& positions) -> v
     ZoneScoped;
     #endif
 
-    t.nodes.clear();
+    t.clear();
-    t.nodes.reserve(positions.size() * 2);
+    t.reserve(positions.size() * 2);
 
     // Compute bounding box around all masses
     Vector3 min{FLT_MAX, FLT_MAX, FLT_MAX};
@@ -194,28 +170,36 @@ auto octree::calculate_force(const int node_idx, const Vector3& pos) const -> Ve
     }
 
     const node& n = nodes[node_idx];
-    if (std::abs(n.mass_total) <= 0.001f) {
+
+    const float mass_total = n.mass_total;
+    const Vector3& mass_center = n.mass_center;
+    const Vector3& box_min = n.box_min;
+    const Vector3& box_max = n.box_max;
+    const std::array<int, 8>& children = n.children;
+    const bool leaf = n.leaf;
+
+    if (std::abs(mass_total) <= 0.001f) {
         return Vector3Zero();
     }
 
-    const Vector3 diff = Vector3Subtract(pos, n.mass_center);
+    const Vector3 diff = Vector3Subtract(pos, mass_center);
     float dist_sq = diff.x * diff.x + diff.y * diff.y + diff.z * diff.z;
 
     // Softening
     dist_sq += SOFTENING;
 
     // Barnes-Hut
-    const float size = n.box_max.x - n.box_min.x;
+    const float size = box_max.x - box_min.x;
-    if (n.leaf || size * size / dist_sq < THETA * THETA) {
+    if (leaf || size * size / dist_sq < THETA * THETA) {
         const float dist = std::sqrt(dist_sq);
-        const float force_mag = BH_FORCE * n.mass_total / dist_sq;
+        const float force_mag = BH_FORCE * mass_total / dist_sq;
 
         return Vector3Scale(diff, force_mag / dist);
     }
 
     // Collect child forces
     Vector3 force = Vector3Zero();
-    for (const int child : n.children) {
+    for (const int child : children) {
         if (child >= 0) {
             const Vector3 child_force = calculate_force(child, pos);
 

src/puzzle.cpp

@@ -1127,7 +1127,7 @@ auto puzzle::explore_puzzle_space(const boost::unordered_flat_set<block, block_h
                                                       });
                              });
 
-    infoln("Found {} of {} clusters with a solution", visited_clusters.size(), total);
+    // infoln("Found {} of {} clusters with a solution", visited_clusters.size(), total);
 
     return visited_clusters;
 }

src/threaded_physics.cpp

@@ -154,7 +154,7 @@ auto threaded_physics::physics_thread(physics_state& state, const std::optional<
                 loop_iterations = 0;
                 ups_accumulator = std::chrono::duration<double>(0);
             }
-            if (mass_springs.tree.nodes.empty()) {
+            if (mass_springs.tree.empty()) {
                 state.mass_center = Vector3Zero();
             } else {
                 state.mass_center = mass_springs.tree.nodes[0].mass_center;