From 063c487ad24c85cc1922cca58294285a8c5f566e Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Sun, 16 Nov 2025 22:16:24 -0800
Subject: [PATCH 1/7] modify voxels in place to avoid mem copy

---
 app/geometry/sphere.py | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/app/geometry/sphere.py b/app/geometry/sphere.py
index efcf00a..e945971 100644
--- a/app/geometry/sphere.py
+++ b/app/geometry/sphere.py
@@ -116,25 +116,24 @@ def _maybe_expand_voxel_space(self, voxel_space, max_index, axis_number):
         return voxel_space
 
     def _deposit_sphere(self, radius, voxel_space, nozzle_height, target_volume):
-        copy_voxel_space = voxel_space.copy()
-
-        copy_voxel_space = self.deform_voxel_space_for_big_spheres(
-            copy_voxel_space, radius
+        # Modify voxel_space in-place to avoid redundant copying
+        voxel_space = self.deform_voxel_space_for_big_spheres(
+            voxel_space, radius
         )
 
         lower_indexes, upper_indexes = self.find_sphere_limits(radius, nozzle_height)
 
-        copy_voxel_space = self.fill_voxels(
-            copy_voxel_space, radius, lower_indexes, upper_indexes
+        voxel_space = self.fill_voxels(
+            voxel_space, radius, lower_indexes, upper_indexes
         )
 
         current_volume = GeometryMath.calculate_filled_volume(
-            copy_voxel_space, self.voxel_size
+            voxel_space, self.voxel_size
         )
 
         volume_overshoot = current_volume / target_volume - 1.0
 
-        return volume_overshoot, copy_voxel_space
+        return volume_overshoot, voxel_space
 
     def _increase_solver_tolerance(self, radius_a, radius_b):
         return radius_b - radius_a < self.voxel_size * 0.5

From 5707e139c89f2e26bbdcfcb7229e3459e8ec3fb7 Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Sun, 16 Nov 2025 22:18:48 -0800
Subject: [PATCH 2/7] vectorize voxel fill operations using numpy

---
 app/geometry/sphere.py | 21 ++++++++++++---------
 1 file changed, 12 insertions(+), 9 deletions(-)

diff --git a/app/geometry/sphere.py b/app/geometry/sphere.py
index e945971..3be99f0 100644
--- a/app/geometry/sphere.py
+++ b/app/geometry/sphere.py
@@ -41,15 +41,18 @@ def fill_voxels(self, voxel_space, radius, lower_indexes, upper_indexes):
             voxel_space, lower_indexes, upper_indexes
         )
 
-        for voxel in empty_voxels:
-            voxel_coordinate = GeometryMath.find_coordinates(voxel, self.voxel_size)
-
-            distance_to_centre = GeometryMath.distance(
-                voxel_coordinate, self.centre_coordinates
-            )
-
-            if distance_to_centre <= radius + self.voxel_size * 1e-8:
-                voxel_space[tuple(voxel)] = 1
+        # Convert to numpy array for vectorized operations
+        empty_voxels_np = np.array(empty_voxels)
+        # Calculate coordinates for all voxels at once
+        voxel_coords = self.voxel_size * (2 * (empty_voxels_np + 1) - 1) * 0.5
+        # Calculate squared distances to centre for all voxels
+        centre = np.array(self.centre_coordinates)
+        dists = np.linalg.norm(voxel_coords - centre, axis=1)
+        # Mask for voxels within radius
+        mask = dists <= radius + self.voxel_size * 1e-8
+        # Assign value 1 to all selected voxels
+        for idx in empty_voxels_np[mask]:
+            voxel_space[tuple(idx)] = 1
 
         return voxel_space
 

From 072267cff9dce0597e5bd67b70638c732414d47a Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Sun, 16 Nov 2025 22:43:42 -0800
Subject: [PATCH 3/7] check for no empty voxels

---
 app/geometry/sphere.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/app/geometry/sphere.py b/app/geometry/sphere.py
index 3be99f0..ba451d2 100644
--- a/app/geometry/sphere.py
+++ b/app/geometry/sphere.py
@@ -41,6 +41,9 @@ def fill_voxels(self, voxel_space, radius, lower_indexes, upper_indexes):
             voxel_space, lower_indexes, upper_indexes
         )
 
+        if not empty_voxels:
+            return voxel_space
+
         # Convert to numpy array for vectorized operations
         empty_voxels_np = np.array(empty_voxels)
         # Calculate coordinates for all voxels at once

From 4bf5918e6acca381c961c649e922776489bc9793 Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Sun, 16 Nov 2025 23:12:27 -0800
Subject: [PATCH 4/7] feat: enhance voxel space expansion logging and update
 bisection method to handle updated arguments

---
 app/geometry/sphere.py          |  9 +++++----
 app/solvers/bisection_method.py | 35 ++++++++++++++++++++-------------
 2 files changed, 26 insertions(+), 18 deletions(-)

diff --git a/app/geometry/sphere.py b/app/geometry/sphere.py
index ba451d2..6091038 100644
--- a/app/geometry/sphere.py
+++ b/app/geometry/sphere.py
@@ -111,11 +111,12 @@ def _maybe_expand_voxel_space(self, voxel_space, max_index, axis_number):
             return voxel_space
 
         number_to_be_added = max_index - index_size + 1
+        import logging
+        logging.info(f"Expanding voxel space on axis {axis_number}: current size {index_size}, required {max_index}, adding {number_to_be_added}")
 
-        new_size = list(size)
-        new_size[axis_number] = number_to_be_added
-
-        mat_add = np.zeros(new_size, dtype=np.int8)
+        mat_add_size = list(size)
+        mat_add_size[axis_number] = number_to_be_added
+        mat_add = np.zeros(mat_add_size, dtype=np.int8)
 
         voxel_space = np.concatenate((voxel_space, mat_add), axis=axis_number)
 
diff --git a/app/solvers/bisection_method.py b/app/solvers/bisection_method.py
index 953ce58..5bfbaa9 100644
--- a/app/solvers/bisection_method.py
+++ b/app/solvers/bisection_method.py
@@ -23,35 +23,42 @@ def execute(
         else:
             point_b = initial_point
 
-        _, _, point_a, point_b = self._loop_fb(fun, point_b, increment, args)
+        fb, out, point_a, point_b, updated_args = self._loop_fb(fun, point_b, increment, args)
 
-        return self._loop_fc(
-            fun, point_a, point_b, tolerance, fun_increase_tolerance, args
+        # _loop_fc returns (fc, out) where out may be the raw output from `fun`.
+        fc, out = self._loop_fc(
+            fun, point_a, point_b, tolerance, fun_increase_tolerance, updated_args
         )
 
+        # If `out` is a tuple like (volume_overshoot, voxel_space), unwrap it
+        if isinstance(out, tuple) and len(out) == 2:
+            return fc, out[1]
+
+        return fc, out
+
     def _loop_fb(self, fun, point_b, inc, args):
         fb = -1.0
-
         point_a = 0.0
-
+        updated_args = list(args)
         while fb < 0.0:
-
-            fb, out = fun(point_b, *args)
-
+            fb, out = fun(point_b, *updated_args)
+            # If the function returns an updated voxel_space as the second value,
+            # `out` will be that ndarray (it will have a `shape` attribute).
+            if hasattr(out, "shape"):
+                updated_args[0] = out
             if fb < 0:
                 point_a = point_b
                 point_b += inc
-
-        return fb, out, point_a, point_b
+        return fb, out, point_a, point_b, updated_args
 
     def _loop_fc(self, fun, point_a, point_b, tolerance, fun_increase_tolerance, args):
         fc = 2.0 * tolerance
-
+        updated_args = list(args)
         while abs(fc) > tolerance:
             point_c = (point_a + point_b) * 0.5
-
-            fc, out = fun(point_c, *args)
-
+            fc, out = fun(point_c, *updated_args)
+            if hasattr(out, "shape"):
+                updated_args[0] = out
             if fun_increase_tolerance(point_a, point_b):
                 logger.debug(
                     "[BisectionMethod]: increasing tolerance because point_b and point_a are too close"

From ffd29cb7842f750871e25aa3f363af4930f6dcdf Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Sun, 16 Nov 2025 23:24:24 -0800
Subject: [PATCH 5/7] preallocate voxel space

---
 app/geometry/sphere.py      |  9 +++++++--
 app/geometry/voxel_space.py | 40 ++++++++++++++++++++++++++++++++++++-
 2 files changed, 46 insertions(+), 3 deletions(-)

diff --git a/app/geometry/sphere.py b/app/geometry/sphere.py
index 6091038..94554e7 100644
--- a/app/geometry/sphere.py
+++ b/app/geometry/sphere.py
@@ -112,10 +112,15 @@ def _maybe_expand_voxel_space(self, voxel_space, max_index, axis_number):
 
         number_to_be_added = max_index - index_size + 1
         import logging
-        logging.info(f"Expanding voxel space on axis {axis_number}: current size {index_size}, required {max_index}, adding {number_to_be_added}")
 
+        # Add a buffer to reduce the number of reallocations. Buffer is 20%
+        # of current size (rounded), but at least the minimum required.
+        buffer_layers = max(int(index_size * 0.2), 1)
+        layers_to_add = max(number_to_be_added, buffer_layers)
+
+        # Create only the additional block to concatenate
         mat_add_size = list(size)
-        mat_add_size[axis_number] = number_to_be_added
+        mat_add_size[axis_number] = layers_to_add
         mat_add = np.zeros(mat_add_size, dtype=np.int8)
 
         voxel_space = np.concatenate((voxel_space, mat_add), axis=axis_number)
diff --git a/app/geometry/voxel_space.py b/app/geometry/voxel_space.py
index 819d576..2bb0f0a 100644
--- a/app/geometry/voxel_space.py
+++ b/app/geometry/voxel_space.py
@@ -47,8 +47,46 @@ def __init__(
         self._consider_acceleration = self._simulation.consider_acceleration
 
     def initialize_space(self):
+        # Try to intelligently preallocate Z dimension to avoid repeated expansions.
+        # Estimate the maximum per-step sphere radius from filaments and add a safety margin.
+        z_dim = self.dimensions["z"]
+
+        try:
+            max_radius = 0.0
+            step_size = self._simulation.step_size
+            for filament in self._instruction.filaments_coordinates:
+                coord_old = filament[0]
+                coord_new = filament[1]
+                volume = filament[2]
+
+                # compute filament length and estimated number of steps
+                length = GeometryMath.distance(coord_old, coord_new)
+                n_steps = int(round(length / step_size))
+                if n_steps < 1:
+                    n_steps = 1
+
+                per_step_volume = volume / n_steps
+                # estimate radius of sphere that would contain this volume
+                radius = (3.0 * per_step_volume / (4.0 * math.pi)) ** (1.0 / 3.0)
+                if radius > max_radius:
+                    max_radius = radius
+
+            # include configured sphere z offset
+            max_extra_z = int(math.ceil((max_radius + self._simulation.sphere_z_offset) / self._simulation.voxel_size))
+
+            # safety margin (20% of current z) to reduce chance of further expansion
+            safety_margin = int(max(1, z_dim * 0.2))
+
+            if max_extra_z > 0:
+                z_dim = z_dim + max_extra_z + safety_margin
+
+            logging.getLogger(__name__).info(f"Preallocating voxel space z-dimension: base={self.dimensions['z']}, extra={max_extra_z}, safety={safety_margin}, final={z_dim}")
+        except Exception:
+            # Fallback to original allocation if any issue occurs during estimation
+            z_dim = self.dimensions["z"]
+
         self.space = np.zeros(
-            (self.dimensions["x"], self.dimensions["y"], self.dimensions["z"]),
+            (self.dimensions["x"], self.dimensions["y"], z_dim),
             dtype=np.int8,
         )
 

From 426b2754dee74ec576afa122f8f43f62912c87a2 Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Sun, 16 Nov 2025 23:27:43 -0800
Subject: [PATCH 6/7] keep track of filled voxels

---
 app/geometry/geometry_math.py   | 15 +++++++
 app/geometry/sphere.py          | 72 +++++++++++++++++++++------------
 app/geometry/voxel_space.py     | 22 ++++++++--
 app/solvers/bisection_method.py |  9 +++--
 4 files changed, 84 insertions(+), 34 deletions(-)

diff --git a/app/geometry/geometry_math.py b/app/geometry/geometry_math.py
index b0dc014..4855fcf 100644
--- a/app/geometry/geometry_math.py
+++ b/app/geometry/geometry_math.py
@@ -49,6 +49,21 @@ def find_coordinates(indexes, voxel_size):
 
     @staticmethod
     def calculate_filled_volume(voxel_space, voxel_size):
+        # Accept either a raw ndarray or a VoxelSpace-like object with a
+        # `_filled_voxels_count` running counter. Prefer the running counter
+        # when available to avoid expensive full-array scans.
+        if hasattr(voxel_space, "_filled_voxels_count"):
+            return int(voxel_space._filled_voxels_count) * voxel_size ** 3
+
+        # If an object with `.space` is provided, try to use its counter first.
+        if hasattr(voxel_space, "space") and hasattr(voxel_space.space, "shape"):
+            # If the container itself tracks a counter, use it.
+            if hasattr(voxel_space, "_filled_voxels_count"):
+                return int(voxel_space._filled_voxels_count) * voxel_size ** 3
+            # Fallback to scanning the ndarray
+            return np.count_nonzero(voxel_space.space) * voxel_size ** 3
+
+        # Otherwise assume voxel_space is a numpy array
         return np.count_nonzero(voxel_space) * voxel_size**3
 
     def find_index(coordinate, voxel_size):
diff --git a/app/geometry/sphere.py b/app/geometry/sphere.py
index 94554e7..37b4828 100644
--- a/app/geometry/sphere.py
+++ b/app/geometry/sphere.py
@@ -21,7 +21,10 @@ def deposit_sphere(
     ):
         initial_radius = self.estimate_initial_radius(sphere_volume)
 
-        _, voxel_space = BisectionMethod().execute(
+        # `voxel_space` is expected to be a VoxelSpace object here. The
+        # bisection solver and inner functions will mutate and ultimately
+        # return that object so callers can read `.space` and counters.
+        _, voxel_space_out = BisectionMethod().execute(
             self._deposit_sphere,
             initial_point=initial_radius,
             tolerance=solver_tolerance,
@@ -34,30 +37,51 @@ def deposit_sphere(
             ),
         )
 
-        return voxel_space
+        return voxel_space_out
 
-    def fill_voxels(self, voxel_space, radius, lower_indexes, upper_indexes):
+    def fill_voxels(self, voxel_space_obj, radius, lower_indexes, upper_indexes):
+        # `voxel_space_obj` is a VoxelSpace instance. Operate on its `.space` ndarray
         empty_voxels = GeometryMath.find_empty_voxels_in_space(
-            voxel_space, lower_indexes, upper_indexes
+            voxel_space_obj.space, lower_indexes, upper_indexes
         )
 
         if not empty_voxels:
-            return voxel_space
+            return voxel_space_obj
 
         # Convert to numpy array for vectorized operations
         empty_voxels_np = np.array(empty_voxels)
         # Calculate coordinates for all voxels at once
         voxel_coords = self.voxel_size * (2 * (empty_voxels_np + 1) - 1) * 0.5
-        # Calculate squared distances to centre for all voxels
+        # Calculate distances to centre for all voxels
         centre = np.array(self.centre_coordinates)
         dists = np.linalg.norm(voxel_coords - centre, axis=1)
         # Mask for voxels within radius
         mask = dists <= radius + self.voxel_size * 1e-8
-        # Assign value 1 to all selected voxels
-        for idx in empty_voxels_np[mask]:
-            voxel_space[tuple(idx)] = 1
 
-        return voxel_space
+        if not np.any(mask):
+            return voxel_space_obj
+
+        # Filter indices that should be filled
+        fill_indices = empty_voxels_np[mask]
+
+        # Advanced index assignment (vectorized)
+        xi = fill_indices[:, 0].astype(int)
+        yj = fill_indices[:, 1].astype(int)
+        zk = fill_indices[:, 2].astype(int)
+
+        # Before setting, count how many of these are actually zero (defensive)
+        # They should be zero because `find_empty_voxels_in_space` returned empties,
+        # but reconfirm to be robust in case of race or prior modifications.
+        current_vals = voxel_space_obj.space[xi, yj, zk]
+        new_mask = current_vals == 0
+        n_new = int(np.count_nonzero(new_mask))
+        if n_new > 0:
+            voxel_space_obj.space[xi[new_mask], yj[new_mask], zk[new_mask]] = 1
+            # Update running counter on the VoxelSpace object
+            if hasattr(voxel_space_obj, "_filled_voxels_count"):
+                voxel_space_obj._filled_voxels_count += n_new
+
+        return voxel_space_obj
 
     def estimate_initial_radius(self, volume):
         return (3.0 * volume / (4.0 * math.pi)) ** (1.0 / 3.0)
@@ -91,6 +115,7 @@ def find_sphere_limits(self, radius, nozzle_height):
     """
 
     def deform_voxel_space_for_big_spheres(self, voxel_space, radius):
+        # voxel_space is expected to be a VoxelSpace instance; check and operate on its `.space`
         max_indexes = [
             self._find_index(coord + radius) for coord in self.centre_coordinates
         ]
@@ -102,16 +127,16 @@ def deform_voxel_space_for_big_spheres(self, voxel_space, radius):
 
         return voxel_space
 
-    def _maybe_expand_voxel_space(self, voxel_space, max_index, axis_number):
-        size = voxel_space.shape
+    def _maybe_expand_voxel_space(self, voxel_space_obj, max_index, axis_number):
+        # Operate on `voxel_space_obj.space` and update in-place by reassigning.
+        size = voxel_space_obj.space.shape
 
         index_size = size[axis_number]
 
         if max_index < index_size:
-            return voxel_space
+            return voxel_space_obj
 
         number_to_be_added = max_index - index_size + 1
-        import logging
 
         # Add a buffer to reduce the number of reallocations. Buffer is 20%
         # of current size (rounded), but at least the minimum required.
@@ -123,28 +148,23 @@ def _maybe_expand_voxel_space(self, voxel_space, max_index, axis_number):
         mat_add_size[axis_number] = layers_to_add
         mat_add = np.zeros(mat_add_size, dtype=np.int8)
 
-        voxel_space = np.concatenate((voxel_space, mat_add), axis=axis_number)
+        voxel_space_obj.space = np.concatenate((voxel_space_obj.space, mat_add), axis=axis_number)
 
-        return voxel_space
+        return voxel_space_obj
 
     def _deposit_sphere(self, radius, voxel_space, nozzle_height, target_volume):
-        # Modify voxel_space in-place to avoid redundant copying
-        voxel_space = self.deform_voxel_space_for_big_spheres(
-            voxel_space, radius
-        )
+        # `voxel_space` is a VoxelSpace instance; deform and fill it in-place.
+        voxel_space = self.deform_voxel_space_for_big_spheres(voxel_space, radius)
 
         lower_indexes, upper_indexes = self.find_sphere_limits(radius, nozzle_height)
 
-        voxel_space = self.fill_voxels(
-            voxel_space, radius, lower_indexes, upper_indexes
-        )
+        voxel_space = self.fill_voxels(voxel_space, radius, lower_indexes, upper_indexes)
 
-        current_volume = GeometryMath.calculate_filled_volume(
-            voxel_space, self.voxel_size
-        )
+        current_volume = GeometryMath.calculate_filled_volume(voxel_space, self.voxel_size)
 
         volume_overshoot = current_volume / target_volume - 1.0
 
+        # Return the computed overshoot and the (possibly mutated) VoxelSpace object
         return volume_overshoot, voxel_space
 
     def _increase_solver_tolerance(self, radius_a, radius_b):
diff --git a/app/geometry/voxel_space.py b/app/geometry/voxel_space.py
index 2bb0f0a..d4c75cf 100644
--- a/app/geometry/voxel_space.py
+++ b/app/geometry/voxel_space.py
@@ -45,6 +45,8 @@ def __init__(
         self._simulation = simulation_config
         self._printer = printer
         self._consider_acceleration = self._simulation.consider_acceleration
+        # Running count of filled voxels (to avoid repeated `np.count_nonzero` calls)
+        self._filled_voxels_count = 0
 
     def initialize_space(self):
         # Try to intelligently preallocate Z dimension to avoid repeated expansions.
@@ -89,6 +91,8 @@ def initialize_space(self):
             (self.dimensions["x"], self.dimensions["y"], z_dim),
             dtype=np.int8,
         )
+        # reset the running counter when allocating space
+        self._filled_voxels_count = 0
 
     def print(self):
         number_printed_filaments = 0
@@ -173,8 +177,9 @@ def _deposit_filament(
         filament_initial_coordinates,
         volumes,
     ):
-        total_deposited_volume = GeometryMath.calculate_filled_volume(
-            self.space, self._simulation.voxel_size
+        # Use the running counter to compute total deposited volume quickly
+        total_deposited_volume = (
+            self._filled_voxels_count * self._simulation.voxel_size ** 3
         )
 
         for step_n in range(0, number_simulation_steps):
@@ -221,11 +226,20 @@ def _deposit_filament(
                 f"Depositing filament: step = {step_n + 1}/{number_simulation_steps}"
             )
 
-            self.space = sphere.deposit_sphere(
-                voxel_space=self.space,
+            # Pass the VoxelSpace instance so sphere code can update the
+            # running counter and expand the `.space` ndarray in-place.
+            voxel_space_out = sphere.deposit_sphere(
+                voxel_space=self,
                 nozzle_height=nozzle_height,
                 sphere_volume=sphere_volume,
                 voxel_space_target_volume=volume_target,
                 solver_tolerance=self._simulation.solver_tolerance,
                 radius_increment=self._simulation.radius_increment,
             )
+
+            # Ensure our internal `.space` and counter reflect any mutations
+            # returned by the sphere logic (voxel_space_out is a VoxelSpace)
+            if hasattr(voxel_space_out, "space"):
+                self.space = voxel_space_out.space
+            if hasattr(voxel_space_out, "_filled_voxels_count"):
+                self._filled_voxels_count = voxel_space_out._filled_voxels_count
diff --git a/app/solvers/bisection_method.py b/app/solvers/bisection_method.py
index 5bfbaa9..80b5e09 100644
--- a/app/solvers/bisection_method.py
+++ b/app/solvers/bisection_method.py
@@ -42,9 +42,10 @@ def _loop_fb(self, fun, point_b, inc, args):
         updated_args = list(args)
         while fb < 0.0:
             fb, out = fun(point_b, *updated_args)
-            # If the function returns an updated voxel_space as the second value,
-            # `out` will be that ndarray (it will have a `shape` attribute).
-            if hasattr(out, "shape"):
+            # If the function returns an updated voxel_space, propagate it.
+            # The returned `out` may be an ndarray (has `shape`) or a
+            # VoxelSpace-like object with a `space` attribute. Accept both.
+            if hasattr(out, "shape") or hasattr(out, "space"):
                 updated_args[0] = out
             if fb < 0:
                 point_a = point_b
@@ -57,7 +58,7 @@ def _loop_fc(self, fun, point_a, point_b, tolerance, fun_increase_tolerance, arg
         while abs(fc) > tolerance:
             point_c = (point_a + point_b) * 0.5
             fc, out = fun(point_c, *updated_args)
-            if hasattr(out, "shape"):
+            if hasattr(out, "shape") or hasattr(out, "space"):
                 updated_args[0] = out
             if fun_increase_tolerance(point_a, point_b):
                 logger.debug(

From 47d7438a28ab2edfe6357701b9cd2b5873393e7b Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Mon, 17 Nov 2025 07:40:13 -0800
Subject: [PATCH 7/7] fix py tests

---
 app/geometry/sphere.py             | 58 +++++++++++++++++-------------
 app/physics/volume.py              |  2 +-
 tests/geometry/voxel_space_test.py | 11 ++++--
 3 files changed, 44 insertions(+), 27 deletions(-)

diff --git a/app/geometry/sphere.py b/app/geometry/sphere.py
index 37b4828..a3e181f 100644
--- a/app/geometry/sphere.py
+++ b/app/geometry/sphere.py
@@ -40,9 +40,12 @@ def deposit_sphere(
         return voxel_space_out
 
     def fill_voxels(self, voxel_space_obj, radius, lower_indexes, upper_indexes):
-        # `voxel_space_obj` is a VoxelSpace instance. Operate on its `.space` ndarray
+        # Accept either a VoxelSpace-like object (has `.space`) or a raw ndarray.
+        is_voxel_space = hasattr(voxel_space_obj, "space")
+        space = voxel_space_obj.space if is_voxel_space else voxel_space_obj
+
         empty_voxels = GeometryMath.find_empty_voxels_in_space(
-            voxel_space_obj.space, lower_indexes, upper_indexes
+            space, lower_indexes, upper_indexes
         )
 
         if not empty_voxels:
@@ -70,18 +73,20 @@ def fill_voxels(self, voxel_space_obj, radius, lower_indexes, upper_indexes):
         zk = fill_indices[:, 2].astype(int)
 
         # Before setting, count how many of these are actually zero (defensive)
-        # They should be zero because `find_empty_voxels_in_space` returned empties,
-        # but reconfirm to be robust in case of race or prior modifications.
-        current_vals = voxel_space_obj.space[xi, yj, zk]
+        current_vals = space[xi, yj, zk]
         new_mask = current_vals == 0
         n_new = int(np.count_nonzero(new_mask))
         if n_new > 0:
-            voxel_space_obj.space[xi[new_mask], yj[new_mask], zk[new_mask]] = 1
-            # Update running counter on the VoxelSpace object
-            if hasattr(voxel_space_obj, "_filled_voxels_count"):
+            # Assign into the ndarray `space` (in-place)
+            space[xi[new_mask], yj[new_mask], zk[new_mask]] = 1
+            # Update running counter only when we were given a VoxelSpace object
+            if is_voxel_space and hasattr(voxel_space_obj, "_filled_voxels_count"):
                 voxel_space_obj._filled_voxels_count += n_new
 
-        return voxel_space_obj
+        # Return the same type we were given
+        if is_voxel_space:
+            return voxel_space_obj
+        return space
 
     def estimate_initial_radius(self, volume):
         return (3.0 * volume / (4.0 * math.pi)) ** (1.0 / 3.0)
@@ -115,22 +120,23 @@ def find_sphere_limits(self, radius, nozzle_height):
     """
 
     def deform_voxel_space_for_big_spheres(self, voxel_space, radius):
-        # voxel_space is expected to be a VoxelSpace instance; check and operate on its `.space`
+        # Accept either a VoxelSpace object or a raw ndarray and return the same type.
         max_indexes = [
             self._find_index(coord + radius) for coord in self.centre_coordinates
         ]
 
+        vs = voxel_space
         for axis_number in range(0, 3):
-            voxel_space = self._maybe_expand_voxel_space(
-                voxel_space, max_indexes[axis_number], axis_number
-            )
+            vs = self._maybe_expand_voxel_space(vs, max_indexes[axis_number], axis_number)
 
-        return voxel_space
+        return vs
 
     def _maybe_expand_voxel_space(self, voxel_space_obj, max_index, axis_number):
-        # Operate on `voxel_space_obj.space` and update in-place by reassigning.
-        size = voxel_space_obj.space.shape
+        # Accept either a VoxelSpace-like object (has `.space`) or a raw ndarray.
+        is_voxel_space = hasattr(voxel_space_obj, "space")
+        space = voxel_space_obj.space if is_voxel_space else voxel_space_obj
 
+        size = space.shape
         index_size = size[axis_number]
 
         if max_index < index_size:
@@ -148,24 +154,28 @@ def _maybe_expand_voxel_space(self, voxel_space_obj, max_index, axis_number):
         mat_add_size[axis_number] = layers_to_add
         mat_add = np.zeros(mat_add_size, dtype=np.int8)
 
-        voxel_space_obj.space = np.concatenate((voxel_space_obj.space, mat_add), axis=axis_number)
+        new_space = np.concatenate((space, mat_add), axis=axis_number)
+
+        if is_voxel_space:
+            voxel_space_obj.space = new_space
+            return voxel_space_obj
 
-        return voxel_space_obj
+        return new_space
 
     def _deposit_sphere(self, radius, voxel_space, nozzle_height, target_volume):
-        # `voxel_space` is a VoxelSpace instance; deform and fill it in-place.
-        voxel_space = self.deform_voxel_space_for_big_spheres(voxel_space, radius)
+        # Accept and return either a VoxelSpace object or a raw ndarray.
+        vs = self.deform_voxel_space_for_big_spheres(voxel_space, radius)
 
         lower_indexes, upper_indexes = self.find_sphere_limits(radius, nozzle_height)
 
-        voxel_space = self.fill_voxels(voxel_space, radius, lower_indexes, upper_indexes)
+        vs = self.fill_voxels(vs, radius, lower_indexes, upper_indexes)
 
-        current_volume = GeometryMath.calculate_filled_volume(voxel_space, self.voxel_size)
+        current_volume = GeometryMath.calculate_filled_volume(vs, self.voxel_size)
 
         volume_overshoot = current_volume / target_volume - 1.0
 
-        # Return the computed overshoot and the (possibly mutated) VoxelSpace object
-        return volume_overshoot, voxel_space
+        # Return the computed overshoot and the (possibly mutated) voxel container
+        return volume_overshoot, vs
 
     def _increase_solver_tolerance(self, radius_a, radius_b):
         return radius_b - radius_a < self.voxel_size * 0.5
diff --git a/app/physics/volume.py b/app/physics/volume.py
index 977a788..4ca2005 100644
--- a/app/physics/volume.py
+++ b/app/physics/volume.py
@@ -27,7 +27,7 @@ def get_volumes_for_filament(
         if not consider_acceleration:
             # Simple case: evenly distribute the volume across steps
             volume_per_step = total_volume / number_simulation_steps
-            return [volume_per_step * (i + 1) for i in range(number_simulation_steps)]
+            return [volume_per_step for i in range(number_simulation_steps)]
         
         # Complex case: use acceleration-specific implementation
         from app.physics.acceleration.volume import AccelerationVolume
diff --git a/tests/geometry/voxel_space_test.py b/tests/geometry/voxel_space_test.py
index d6c3895..7bfa546 100644
--- a/tests/geometry/voxel_space_test.py
+++ b/tests/geometry/voxel_space_test.py
@@ -13,6 +13,10 @@ def __init__(self):
         self._number_printed_filaments = 0
         self._filaments_coordinates = list()
         self._default_nozzle_speed = 10.0
+        # Patch: provide a _printer attribute for compatibility
+        class _Printer:
+            feedstock_filament_diameter = 1.75
+        self._printer = _Printer()
 
     def read(self):
         pass
@@ -47,6 +51,7 @@ def __init__(self):
         self.step_size = 0.2
         self.radius_increment = 0.1
         self.solver_tolerance = 0.001
+        self.consider_acceleration = False
 
 
 class FakePrinter:
@@ -114,17 +119,19 @@ def __init__(self):
         self.step_size = 0.2
         self.radius_increment = 0.1
         self.solver_tolerance = 0.001
+        self.consider_acceleration = False
+        self.sphere_z_offset = 0.0
 
 
 class TestPrint:
     def test_should_print(self):
+        printer = FakePrinter()
         test_instruction = Gcode(
-            gcode_path="tests/fixtures/gcode_example.gcode", default_nozzle_speed=40.0
+            gcode_path="tests/fixtures/gcode_example.gcode", default_nozzle_speed=40.0, printer=printer
         )
         test_instruction.read()
 
         simulation_config = FakeSimulationPrint()
-        printer = FakePrinter()
 
         voxel_space = VoxelSpace(
             instruction=test_instruction,