From c056b8005058d425737784362cf4339ec8e1d2f8 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Mon, 28 Apr 2025 22:31:52 -0600
Subject: [PATCH 01/26] Added support for 3D conservative resampling to
 `regridding.weights()`.

---
 .../_weights/_weights_conservative_3d.py      | 568 ++++++++++++++++++
 .../_weights/_weights_conservative_3d_test.py |  79 +++
 2 files changed, 647 insertions(+)
 create mode 100644 regridding/_weights/_weights_conservative_3d.py
 create mode 100644 regridding/_weights/_weights_conservative_3d_test.py

diff --git a/regridding/_weights/_weights_conservative_3d.py b/regridding/_weights/_weights_conservative_3d.py
new file mode 100644
index 0000000..a053bda
--- /dev/null
+++ b/regridding/_weights/_weights_conservative_3d.py
@@ -0,0 +1,568 @@
+import sys
+import numpy as np
+import numpy.typing as npt
+import numba
+import regridding
+
+__all__ = [
+    "weights_conservative_3d",
+]
+
+axis_x = 0
+axis_y = 1
+axis_z = 2
+
+
+@numba.njit(cache=True)
+def weights_conservative_3d(
+    grid_input: tuple[np.ndarray, np.ndarray, np.ndarray],
+    grid_output: tuple[np.ndarray, np.ndarray, np.ndarray],
+) -> numba.typed.List[tuple[int, int, float]]:
+    """
+    For each cell of `grid_output`,
+    compute the fraction of volume shared with each cell of `grid_input`
+    and save as a list of weights.
+
+    Parameters
+    ----------
+    grid_input
+        The vertices of the old grid.
+        Every component must have the same 3D shape.
+    grid_output
+        The vertices of the new grid.
+        Every component must have the same 3D shape.
+    """
+
+    x_input, y_input, z_input = grid_input
+    x_output, y_output, z_output = grid_output
+
+    shape_input = x_input.shape
+    shape_output = x_output.shape
+
+    weights = numba.typed.List()
+    for x in range(0):
+        weights.append((0.0, 0.0, 0.0))
+
+    volume_input = _cell_volume(grid_input)
+
+    intercepts = _empty_intercepts(shape_input, shape_output)
+
+    for sweep_input in (False, True):
+        _sweep_grid(
+            grid_input=grid_input,
+            grid_output=grid_output,
+            volume_input=volume_input,
+            weights=weights,
+            intercepts=intercepts,
+            sweep_input=sweep_input,
+        )
+
+    # _sweep_intercepts(
+    #     grid_input=grid_input,
+    #     grid_output=grid_output,
+    #     volume_input=volume_input,
+    #     intercepts=intercepts,
+    #     weights=weights,
+    # )
+
+    return weights, intercepts
+
+
+@numba.njit(cache=True)
+def _sweep_grid(
+    grid_input: tuple[np.ndarray, np.ndarray, np.ndarray],
+    grid_output: tuple[np.ndarray, np.ndarray, np.ndarray],
+    volume_input: np.ndarray,
+    weights: numba.typed.List[tuple[int, int, float]],
+    intercepts: numba.typed.List,
+    sweep_input: bool,
+) -> None:
+    """
+    Sweep along each of the three axes and each of the three diagonals
+    of either `grid_input` or `grid_output` and compute the weights and
+    intercepts.
+
+    Parameters
+    ----------
+    grid_input
+        The vertices of the old grid.
+    grid_output
+        The vertices of the new grid.
+    volume_input
+        The area of the input grid which will used to compute the weight.
+    weights
+        The current list of weights.
+        New weights will be appended to this list.
+    intercepts
+        A sorted list of intercepts to be traversed later.
+        As new intercepts are found, they are inserted into the list.
+    sweep_input
+        A boolean flag indicating whether to iterate along `grid_input`
+        or `grid_output`.
+        If :obj:`True`, this function sweeps along `grid_input`.
+    """
+
+    grid_sweep, grid_static = _grid_sweep_static(grid_input, grid_output, sweep_input)
+
+    grid_static_x, grid_static_y, grid_static_z = grid_static
+
+    bbox_boundary = (
+        (grid_static_x.min(), grid_static_y.min(), grid_static_z.min()),
+        (grid_static_x.max(), grid_static_y.max(), grid_static_z.max()),
+    )
+
+    index_boundary, coord_boundary = _calc_boundary(grid_static)
+
+    axes_parallel = (
+        axis_x,
+        axis_y,
+        axis_z,
+    )
+
+    axes_diagonal = (
+        (axis_x, axis_y),
+        (axis_y, axis_z),
+        (axis_z, axis_x),
+    )
+
+    for axis in axes_parallel:
+        _sweep_along_axis(
+            grid_input=grid_input,
+            grid_output=grid_output,
+            volume_input=volume_input,
+            bbox_boundary=bbox_boundary,
+            index_boundary=index_boundary,
+            coord_boundary=coord_boundary,
+            weights=weights,
+            intercepts=intercepts,
+            sweep_input=sweep_input,
+            axes=axis,
+        )
+
+    for axis in axes_diagonal:
+        _sweep_along_diagonal(
+
+        )
+
+
+@numba.njit(cache=True, parallel=True)
+def _sweep_along_axis(
+    grid_input: tuple[np.ndarray, np.ndarray, np.ndarray],
+    grid_output: tuple[np.ndarray, np.ndarray, np.ndarray],
+    volume_input: np.ndarray,
+    bbox_boundary: tuple[tuple[float, float, float], tuple[float, float, float]],
+    index_boundary: tuple[np.ndarray, np.ndarray, np.ndarray],
+    coord_boundary: tuple[np.ndarray, np.ndarray, np.ndarray],
+    weights: numba.typed.List[tuple[int, int, float]],
+    intercepts: numba.typed.List,
+    sweep_input: bool,
+    axis: int,
+) -> None:
+    """
+    Sweep along one logical axis of either `grid_input` or `grid_output`.
+    At each vertex, compute and save the weights and intercepts.
+
+    Parameters
+    ----------
+    grid_input
+        The vertices of the old grid.
+    grid_output
+        The vertices of the new grid.
+    volume_input
+        The area of the input grid which will used to compute the weight.
+    bbox_boundary
+        Two points defining the bounding box of the static grid.
+    coord_boundary
+        A sequence of triangles defining the outer surface of the static grid.
+    index_boundary
+        The index of the static grid corresponding to each triangle in
+        `coord_boundary`.
+    weights
+        The current list of weights.
+        New weights will be appended to this list.
+    sweep_input
+        A boolean flag indicating whether to iterate along `grid_input`
+        or `grid_output`.
+        If :obj:`True`, this function sweeps along `grid_input`.
+    axis
+        The logical axis of the `grid_sweep` to iterate along.
+    """
+
+    grid_sweep, grid_static = _grid_sweep_static(grid_input, grid_output, sweep_input)
+
+    grid_sweep_x, grid_sweep_y, grid_sweep_z = grid_sweep
+    grid_static_x, grid_static_y, grid_static_z = grid_static
+
+    grid_sweep_x = _align_axis_right(grid_sweep_x, axis)
+    grid_sweep_y = _align_axis_right(grid_sweep_y, axis)
+    grid_sweep_z = _align_axis_right(grid_sweep_z, axis)
+
+    shape_sweep_x, shape_sweep_y, shape_sweep_z = grid_sweep.shape
+    shape_static_x, shape_static_y, shape_static_z = grid_static.shape
+
+    for i in numba.prange(shape_sweep_x):
+
+        i = numba.types.int64(i)
+
+        for j in range(shape_sweep_y):
+
+            k = 0
+
+            index_sweep = i, j, k
+
+            index_static = sys.maxsize, sys.maxsize, sys.maxsize
+
+            sweep_is_outside_static = True
+
+            point_1 = (
+                grid_sweep_x[index_sweep],
+                grid_sweep_y[index_sweep],
+                grid_sweep_z[index_sweep],
+            )
+
+            if regridding.geometry.point_is_inside_box_3d(
+                point=point_1,
+                box=bbox_boundary,
+            ):
+                if regridding.geometry.point_is_inside_polyhedron(
+                    point=point_1,
+                    polyhedron=coord_boundary,
+                ):
+                    index_static = _index_of_point_brute(
+                        point=point_1,
+                        grid=grid_static,
+                    )
+
+            while True:
+
+                point_2 = (
+                    grid_sweep_x[i, j, k + 1],
+                    grid_sweep_y[i, j, k + 1],
+                    grid_sweep_z[i, j, k + 1],
+                )
+
+                sweep_is_outside_static = index_static[0] == sys.maxsize
+
+                if sweep_is_outside_static:
+
+                    point_2, k, index_static = _step_outside_static(
+                        point_1=point_1,
+                        point_2=point_2,
+                        bbox_boundary=bbox_boundary,
+                        index_boundary=index_boundary,
+                        coord_boundary=coord_boundary,
+                        intercepts=intercepts,
+                    )
+
+                else:
+
+                    point_2, k, index_static = _step_inside_static(
+                        point_1=point_1,
+                        point_2=point_2,
+                        index_static=index_static,
+                        grid_static=grid_static,
+                        weights=weights,
+                        intercepts=intercepts,
+                    )
+
+                    if not 0 <= index_static[axis_x] < (shape_static_x - 1):
+                        break
+                    if not 0 <= index_static[axis_y] < (shape_static_y - 1):
+                        break
+                    if not 0 <= index_static[axis_z] < (shape_static_z - 1):
+                        break
+
+                if k >= (shape_sweep_z - 1):
+                    break
+
+
+
+
+
+
+
+
+@numba.njit(cache=True)
+def _step_outside_static(
+    point_1: tuple[float, float, float],
+    point_2: tuple[float, float, float],
+    index_static: tuple[int, int, int],
+    index_sweep: tuple[int, int, int],
+    grid_static: tuple[np.ndarray, np.ndarray, np.ndarray],
+    bbox_boundary: tuple[tuple[float, float, float], tuple[float, float, float]],
+    coord_boundary: tuple[np.ndarray, np.ndarray, np.ndarray],
+    index_boundary: tuple[np.ndarray, np.ndarray, np.ndarray],
+) -> tuple[
+    tuple[float, float, float],
+]:
+    """
+    Check if the current line segment crosses the boundary.
+    If it does, compute the intersection and return the intersection point
+    and the indices of the `grid_static` where the intersection occurs.
+
+    Parameters
+    ----------
+    point_1
+        The first point of the line segment.
+    point_2
+        The second point of the line segment.
+    index_static
+        The current 3D index in the static grid.
+    index_sweep
+        The current 3D index in the sweep grid.
+    grid_static
+        Coordinates of the static grid.
+    bbox_boundary
+        Two points defining the bounding box of the static grid.
+    coord_boundary
+        A sequence of triangles defining the outer surface of the static grid.
+    index_boundary
+        The index of the static grid corresponding to each triangle in
+        `coord_boundary`.
+    """
+
+
+@numba.njit(cache=True)
+def _step_inside_static(
+    point_1: tuple[float, float, float],
+    point_2: tuple[float, float, float],
+    index_static: tuple[int, int, int],
+    index_sweep: tuple[int, int, int],
+    grid_static: tuple[np.ndarray, np.ndarray, np.ndarray],
+    weights: numba.typed.List[tuple[int, int, float]],
+) -> None:
+    """
+    Check if the current line segment crosses any face of the current cell
+    in the static grid.
+    If it does, compute the intersection point and indices of the new
+    cell in the static grid.
+
+    Parameters
+    ----------
+    point_1
+        The first point of the line segment.
+    point_2
+        The second point of the line segment.
+    index_static
+        The current index in the static grid.
+    index_sweep
+        The current index in the sweep grid.
+    grid_static
+        Coordinates of the static grid.
+    weights
+        The current list of weights.
+        New weights will be appended to this list.
+    """
+
+
+@numba.njit(cache=True)
+def _grid_sweep_static(
+    grid_input: tuple[np.ndarray, np.ndarray, np.ndarray],
+    grid_output: tuple[np.ndarray, np.ndarray, np.ndarray],
+    sweep_input: bool,
+) -> tuple[
+    tuple[np.ndarray, np.ndarray, np.ndarray],
+    tuple[np.ndarray, np.ndarray, np.ndarray],
+]:
+    if sweep_input:
+        grid_sweep = grid_input
+        grid_static = grid_output
+    else:
+        grid_sweep = grid_output
+        grid_static = grid_input
+
+    return grid_sweep, grid_static
+
+
+@numba.njit(cache=True)
+def _cell_volume(
+    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
+) -> np.ndarray:
+    """
+    Compute the volume of each cell in a logically-rectangular grid.
+
+    Parameters
+    ----------
+    grid
+        A 3D grid of cell vertices.
+    """
+
+    x, y, z = grid
+
+    num_i, num_j, num_k = x.shape
+
+    result = np.zeros((num_i - 1, num_j - 1, num_k - 1))
+
+    for axis in (0, 1, 2):
+        _cell_volume_sweep(
+            grid=grid,
+            out=result,
+            axis=axis,
+        )
+
+    return result
+
+
+@numba.njit(cache=True)
+def _cell_volume_sweep(
+    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
+    out: np.ndarray,
+    axis: int,
+) -> None:
+    """
+    Compute the volume contribution of this axis.
+
+    Parameters
+    ----------
+    grid
+        A 3D grid of cell vertices.
+    out
+        An output to array to store the result.
+    axis
+        The axis along which to iterate.
+    """
+
+    x, y, z = grid
+
+    x = _align_axis_right(x, axis)
+    y = _align_axis_right(y, axis)
+    z = _align_axis_right(z, axis)
+
+    out = _align_axis_right(out, axis)
+
+    num_i, num_j, num_k = x.shape
+
+    for i in numba.prange(num_i):
+
+        i = numba.types.int64(i)
+
+        for j in range(num_j):
+
+            for k in range(num_k - 1):
+
+                i_left = i - 1
+                i_right = i
+
+                j_lower = j - 1
+                j_upper = j
+
+                k1 = k
+                k2 = k + 1
+
+                v1 = x[i, j, k1], y[i, j, k1], z[i, j, k1]
+                v2 = x[i, j, k2], y[i, j, k2], z[i, j, k2]
+
+                volume_left = 0
+                volume_lower = 0
+
+                if i_left >= 0:
+                    v0_left = x[i_left, j, k], y[i_left, j, k], z[i_left, j, k]
+                    volume_left = _volume(v0_left, v1, v2)
+
+                if j_lower >= 0:
+                    v0_lower = x[i, j_lower, k], y[i, j_lower, k], z[i, j_lower, k]
+                    volume_lower = -_volume(v0_lower, v1, v2)
+
+                if i_left >= 0:
+                    if j_lower >= 0:
+                        out[i_left, j_lower, k] -= volume_left
+                    if j_upper < (num_j - 1):
+                        out[i_left, j_upper, k] += volume_left
+
+                if j_lower >= 0:
+                    if i_left >= 0:
+                        out[i_left, j_lower, k] -= volume_lower
+                    if i_right < (num_i - 1):
+                        out[i_right, j_lower, k] += volume_lower
+
+
+@numba.njit(cache=True)
+def _volume(
+    v0: tuple[float, float, float],
+    v1: tuple[float, float, float],
+    v2: tuple[float, float, float],
+) -> float:
+    """
+    Compute the volume of a tetrahedron constructed from three
+    vertices and the origin.
+
+    Parameters
+    ----------
+    v0
+        First vertex of the tetrahedron.
+    v1
+        Second vertex of the tetrahedron.
+    v2
+        Third vertex of the tetrahedron.
+    """
+    triple_product = regridding.math.dot_3d(
+        a=v0,
+        b=regridding.math.cross_3d(v1, v2),
+    )
+
+    return triple_product / 6
+
+
+@numba.njit(cache=True)
+def _empty_intercepts(
+    shape_input: tuple[int, int, int],
+    shape_output: tuple[int, int, int],
+) -> numba.typed.List:
+    """
+    Create an empty list of intercepts between each plane in the
+    input grid and each plane in the output grid.
+
+    Parameters
+    ----------
+    shape_input
+        The shape of the input grid.
+    shape_output
+        The shape of the output grid.
+    """
+    intercepts = numba.typed.List()
+    for a in range(len(shape_input)):
+        intercepts_a = numba.typed.List()
+        for b in range(len(shape_output)):
+            intercepts_ab = numba.typed.List()
+            for i in range(shape_input[a]):
+                intercepts_abi = numba.typed.List()
+                for j in range(shape_output[b]):
+                    intercepts_abij = numba.typed.List()
+                    for x in range(0):
+                        intercepts_abij.append((0.0, 0.0, 0.0))
+                    intercepts_abi.append(intercepts_abij)
+                intercepts_ab.append(intercepts_abi)
+            intercepts_a.append(intercepts_ab)
+        intercepts.append(intercepts_a)
+    return intercepts
+
+
+@numba.njit(cache=True)
+def _align_axis_right(
+    a: np.ndarray,
+    axis: int,
+) -> np.ndarray:
+    """
+    Roll all the axes of a 3D array to the right until `axis` is the last axis..
+
+    This function is needed to permute the axes in such a way to retain
+    their right-handedness.
+
+    Parameters
+    ----------
+    a
+        The array to modify the axes of.
+    axis
+        The axis to set as the last axis.
+    """
+    axis = axis % a.ndim
+
+    if axis == axis_x:
+        result = a.transpose((axis_y, axis_z, axis_x))
+    elif axis == axis_y:
+        result = a.transpose((axis_z, axis_x, axis_y))
+    else:
+        result = a
+
+    return result
+
diff --git a/regridding/_weights/_weights_conservative_3d_test.py b/regridding/_weights/_weights_conservative_3d_test.py
new file mode 100644
index 0000000..7744159
--- /dev/null
+++ b/regridding/_weights/_weights_conservative_3d_test.py
@@ -0,0 +1,79 @@
+import pytest
+import numpy as np
+from . import _weights_conservative_3d
+
+_x = np.arange(3)
+_y = np.arange(4)
+_z = np.arange(5)
+
+_x, _y, _z = np.meshgrid(_x, _y, _z, indexing="ij")
+
+
+@pytest.mark.parametrize(
+    argnames="grid_input,grid_output",
+    argvalues=[
+        (
+            np.meshgrid(
+                np.arange(2),
+                np.arange(2),
+                np.arange(2),
+                indexing="ij",
+            ),
+            np.meshgrid(
+                np.arange(2),
+                np.arange(2),
+                np.arange(2),
+                indexing="ij",
+            ),
+        )
+    ]
+)
+def test_weights_conservative_3d(
+    grid_input: tuple[np.ndarray, np.ndarray, np.ndarray],
+    grid_output: tuple[np.ndarray, np.ndarray, np.ndarray],
+):
+    result = _weights_conservative_3d.weights_conservative_3d(
+        grid_input=grid_input,
+        grid_output=grid_output,
+    )
+
+
+@pytest.mark.parametrize(
+    argnames="grid,result_expected",
+    argvalues=[
+        (
+            np.meshgrid(
+                np.arange(2),
+                np.arange(2),
+                np.arange(2),
+                indexing="ij",
+            ),
+            np.ones((1, 1, 1))
+        ),
+        (
+            np.meshgrid(
+                np.arange(3),
+                np.arange(4),
+                np.arange(5),
+                indexing="ij",
+            ),
+            np.ones((2, 3, 4))
+        ),
+        (
+            np.meshgrid(
+                2 * np.arange(3),
+                2 * np.arange(4),
+                2 * np.arange(5),
+                indexing="ij",
+            ),
+            8 * np.ones((2, 3, 4))
+        ),
+    ]
+)
+def test_cell_volume(
+    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
+    result_expected: np.ndarray,
+):
+    result = _weights_conservative_3d._cell_volume(grid)
+    assert np.allclose(result, result_expected)
+    assert result.shape == result_expected.shape

From 90ea88ce2a481579e1e22b4d6c271a4dd4fcf575 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Tue, 29 Apr 2025 20:08:59 -0600
Subject: [PATCH 02/26] added `_index_of_point_brute()`

---
 .../_weights/_weights_conservative_3d.py      | 72 +++++++++++++++++++
 .../_weights/_weights_conservative_3d_test.py | 38 ++++++++++
 2 files changed, 110 insertions(+)

diff --git a/regridding/_weights/_weights_conservative_3d.py b/regridding/_weights/_weights_conservative_3d.py
index a053bda..8f35aac 100644
--- a/regridding/_weights/_weights_conservative_3d.py
+++ b/regridding/_weights/_weights_conservative_3d.py
@@ -503,6 +503,78 @@ def _volume(
     return triple_product / 6
 
 
+@numba.njit(cache=True)
+def _index_of_point_brute(
+    point: tuple[float, float, float],
+    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
+) -> tuple[int, int, int]:
+    """
+    Find the index of the cell in the grid which contains the given point.
+
+    This function uses brute force to search,
+    but this could be improved significantly by using the secant method
+    or possibly the bisection method.
+
+    Parameters
+    ----------
+    point
+        The query point.
+    grid
+        A logically-rectangular grid of cell vertices.
+    """
+
+    x, y, z = grid
+
+    shape_x, shape_y, shape_z = x.shape
+
+    for i in range(shape_x - 1):
+        for j in range(shape_y - 1):
+            for k in range(shape_z - 1):
+
+                i_000 = i + 0, j + 0, k + 0
+                i_001 = i + 0, j + 0, k + 1
+                i_010 = i + 0, j + 1, k + 0
+                i_011 = i + 0, j + 1, k + 1
+                i_100 = i + 1, j + 0, k + 0
+                i_101 = i + 1, j + 0, k + 1
+                i_110 = i + 1, j + 1, k + 0
+                i_111 = i + 1, j + 1, k + 1
+
+                indices = (
+                    (i_000, i_010, i_110),
+                    (i_110, i_100, i_000),
+                    (i_001, i_101, i_111),
+                    (i_111, i_011, i_001),
+                    (i_000, i_001, i_011),
+                    (i_011, i_010, i_000),
+                    (i_100, i_110, i_111),
+                    (i_111, i_101, i_100),
+                    (i_000, i_100, i_101),
+                    (i_101, i_001, i_000),
+                    (i_010, i_011, i_111),
+                    (i_111, i_110, i_010),
+                )
+
+                polyhedron = numba.typed.List()
+
+                for index in indices:
+
+                    t0, t1, t2 = index
+
+                    triangle = (
+                        (x[t0], y[t0], z[t0]),
+                        (x[t1], y[t1], z[t1]),
+                        (x[t2], y[t2], z[t2]),
+                    )
+
+                    polyhedron.append(triangle)
+
+                if regridding.geometry.point_is_inside_polyhedron(
+                    point=point, polyhedron=polyhedron
+                ):
+                    return i, j, k
+
+
 @numba.njit(cache=True)
 def _empty_intercepts(
     shape_input: tuple[int, int, int],
diff --git a/regridding/_weights/_weights_conservative_3d_test.py b/regridding/_weights/_weights_conservative_3d_test.py
index 7744159..ee2fdba 100644
--- a/regridding/_weights/_weights_conservative_3d_test.py
+++ b/regridding/_weights/_weights_conservative_3d_test.py
@@ -77,3 +77,41 @@ def test_cell_volume(
     result = _weights_conservative_3d._cell_volume(grid)
     assert np.allclose(result, result_expected)
     assert result.shape == result_expected.shape
+
+
+@pytest.mark.parametrize(
+    argnames="point,grid,result_expected",
+    argvalues=[
+        (
+            (0.5, 0.5, 0.5),
+            np.meshgrid(
+                np.arange(3),
+                np.arange(4),
+                np.arange(5),
+                indexing="ij",
+            ),
+            (0, 0, 0),
+        ),
+        (
+            (-0.5, -0.5, 3.5),
+            np.meshgrid(
+                -np.arange(3),
+                -np.arange(4),
+                np.arange(5),
+                indexing="ij",
+            ),
+            (0, 0, 3),
+        ),
+    ],
+)
+def test_index_of_point_brute(
+    point: tuple[float, float, float],
+    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
+    result_expected: np.ndarray,
+):
+    result = _weights_conservative_3d._index_of_point_brute(
+        point=point,
+        grid=grid,
+    )
+
+    assert result == result_expected

From 66887cb0b13cbf28769ed7928eb51bd470bbb811 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Tue, 29 Apr 2025 23:08:09 -0600
Subject: [PATCH 03/26] Added `_indices_boundary()` function.

---
 .../_weights/_weights_conservative_3d.py      | 123 ++++++++++++++++++
 1 file changed, 123 insertions(+)

diff --git a/regridding/_weights/_weights_conservative_3d.py b/regridding/_weights/_weights_conservative_3d.py
index 8f35aac..363e161 100644
--- a/regridding/_weights/_weights_conservative_3d.py
+++ b/regridding/_weights/_weights_conservative_3d.py
@@ -503,6 +503,129 @@ def _volume(
     return triple_product / 6
 
 
+@numba.njit(cache=True)
+def _indices_boundary(
+    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
+) -> numba.typed.List[
+    tuple[
+        tuple[int, int, int],
+        tuple[int, int, int],
+        tuple[int, int, int],
+    ],
+]:
+    """
+    For a given grid, find the set of indices that expresses the boundary
+    as a sequence of triangles.
+
+    Parameters
+    ----------
+    grid
+        A logically-rectangular grid of cell vertices.
+    """
+
+    x, y, z = grid
+
+    shape_x, shape_y, shape_z = x.shape
+
+    triangles = numba.typed.List()
+
+    for i0 in range(shape_x - 1):
+        for j0 in range(shape_y - 1):
+
+            k0 = 0
+
+            i1 = i0 + 1
+            j1 = j0 + 1
+
+            v_000 = i0, j0, k0
+            v_010 = i0, j1, k0
+            v_110 = i1, j1, k0
+            v_100 = i1, j0, k0
+
+            triangles.append((v_000, v_010, v_110))
+            triangles.append((v_110, v_100, v_000))
+
+    for i0 in range(shape_x - 1):
+        for j0 in range(shape_y - 1):
+
+            k1 = shape_z
+
+            i1 = i0 + 1
+            j1 = j0 + 1
+
+            v_001 = i0, j0, k1
+            v_101 = i1, j0, k1
+            v_111 = i1, j1, k1
+            v_011 = i0, j1, k1
+
+            triangles.append((v_001, v_101, v_111))
+            triangles.append((v_111, v_011, v_001))
+
+    for j0 in range(shape_y - 1):
+        for k0 in range(shape_z - 1):
+
+            i0 = 0
+
+            j1 = j0 + 1
+            k1 = k0 + 1
+
+            v_000 = i0, j0, k0
+            v_001 = i0, j0, k1
+            v_011 = i0, j1, k1
+            v_010 = i0, j1, k0
+
+            triangles.append((v_000, v_001, v_011))
+            triangles.append((v_011, v_010, v_000))
+
+    for j0 in range(shape_y - 1):
+        for k0 in range(shape_z - 1):
+
+            i1 = shape_x
+
+            j1 = j0 + 1
+            k1 = k0 + 1
+
+            v_100 = i1, j0, k0
+            v_110 = i1, j1, k0
+            v_111 = i1, j1, k1
+            v_101 = i1, j0, k1
+
+            triangles.append((v_100, v_110, v_111))
+            triangles.append((v_111, v_101, v_100))
+
+    for i0 in range(shape_x - 1):
+        for k0 in range(shape_z - 1):
+
+            j0 = 0
+
+            i1 = i0 + 1
+            k1 = k0 + 1
+
+            v_000 = i0, j0, k0
+            v_100 = i1, j0, k0
+            v_101 = i1, j0, k1
+            v_001 = i0, j0, k1
+
+            triangles.append((v_000, v_100, v_101))
+            triangles.append(v_101, v_001, v_000)
+
+    for i0 in range(shape_x - 1):
+        for k0 in range(shape_z - 1):
+
+            j1 = shape_y
+
+            i1 = i0 + 1
+            k1 = k0 + 1
+
+            v_010 = i0, j1, k0
+            v_011 = i0, j1, k1
+            v_111 = i1, j1, k1
+            v_110 = i1, j1, k0
+
+            triangles.append((v_010, v_011, v_111))
+            triangles.append((v_111, v_110, v_010))
+
+
 @numba.njit(cache=True)
 def _index_of_point_brute(
     point: tuple[float, float, float],

From 01348d705349d45ab577fc3da425f2d90c5d30d1 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Fri, 2 May 2025 11:31:17 -0600
Subject: [PATCH 04/26] Added `_arrays` module.

---
 .../_weights_conservative_3d/__init__.py      |   0
 .../_weights_conservative_3d/_arrays.py       | 139 ++++++++++++++++++
 .../_weights_conservative_3d/_arrays_test.py  |  89 +++++++++++
 3 files changed, 228 insertions(+)
 create mode 100644 regridding/_weights/_weights_conservative_3d/__init__.py
 create mode 100644 regridding/_weights/_weights_conservative_3d/_arrays.py
 create mode 100644 regridding/_weights/_weights_conservative_3d/_arrays_test.py

diff --git a/regridding/_weights/_weights_conservative_3d/__init__.py b/regridding/_weights/_weights_conservative_3d/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/regridding/_weights/_weights_conservative_3d/_arrays.py b/regridding/_weights/_weights_conservative_3d/_arrays.py
new file mode 100644
index 0000000..6962e6d
--- /dev/null
+++ b/regridding/_weights/_weights_conservative_3d/_arrays.py
@@ -0,0 +1,139 @@
+"""
+3-dimensional array manipulation routines.
+"""
+
+import numpy as np
+import numba
+
+__all__ = [
+    "axis_x",
+    "axis_y",
+    "axis_z",
+    "index_in_bounds",
+    "index_flat",
+    "index_3d",
+    "align_axis_right",
+
+]
+
+axis_x = 0
+axis_y = 1
+axis_z = 2
+
+
+@numba.njit(cache=True)
+def index_in_bounds(
+    index: tuple[int, int, int],
+    shape: tuple[int, int, int],
+):
+    """
+    Check if a 3D index is within array bounds specified by `shape`.
+
+    Parameters
+    ----------
+    index
+        The 3D index to check.
+    shape
+        The shape of the array to use as the bounds.
+    """
+    i, j, k = index
+
+    if i < 0:
+        return False
+    if j < 0:
+        return False
+    if k < 0:
+        return False
+
+    nx, ny, nz = shape
+
+    if i >= nx:
+        return False
+    if j >= ny:
+        return False
+    if k >= nz:
+        return False
+
+    return True
+
+
+@numba.njit(cache=True)
+def index_flat(
+    index: tuple[int, int, int],
+    shape: tuple[int, int, int],
+) -> int:
+    """
+    Convert a 3D index to a flat index.
+
+    Parameters
+    ----------
+    index
+        A 3D index to convert.
+    shape
+        The sizes of each axis of the array.
+    """
+
+    i, j, k = index
+
+    num_x, num_y, num_z = shape
+
+    result = i * num_y * num_z + j * num_z + k
+
+    return result
+
+
+@numba.njit(cache=True)
+def index_3d(
+    index: int,
+    shape: tuple[int, int, int],
+) -> tuple[int, int, int]:
+    """
+    Convert a flat index to a 3D index.
+
+    Parameters
+    ----------
+    index
+        A flat index to convert.
+    shape
+        The sizes of each axis of the array.
+    """
+
+    num_x, num_y, num_z = shape
+
+    k = index % num_z
+
+    j = (index // num_z) % num_y
+
+    i = index // (num_z * num_y)
+
+    return i, j, k
+
+
+@numba.njit(cache=True)
+def align_axis_right(
+    a: np.ndarray,
+    axis: int,
+) -> np.ndarray:
+    """
+    Roll all the axes of a 3D array to the right until `axis` is the last axis..
+
+    This function is needed to permute the axes in such a way to retain
+    their right-handedness.
+
+    Parameters
+    ----------
+    a
+        The array to modify the axes of.
+    axis
+        The axis to set as the last axis.
+    """
+    axis = axis % a.ndim
+
+    if axis == axis_x:
+        result = a.transpose((axis_y, axis_z, axis_x))
+    elif axis == axis_y:
+        result = a.transpose((axis_z, axis_x, axis_y))
+    else:
+        result = a
+
+    return result
diff --git a/regridding/_weights/_weights_conservative_3d/_arrays_test.py b/regridding/_weights/_weights_conservative_3d/_arrays_test.py
new file mode 100644
index 0000000..6f76f7a
--- /dev/null
+++ b/regridding/_weights/_weights_conservative_3d/_arrays_test.py
@@ -0,0 +1,89 @@
+import pytest
+import numpy as np
+from . import _arrays
+
+
+@pytest.mark.parametrize(
+    argnames="index,shape,result_expected",
+    argvalues=[
+        (
+            (0, 0, 0),
+            (11, 12, 13),
+            True,
+        ),
+        (
+            (5, 5, 5),
+            (11, 12, 13),
+            True,
+        ),
+        (
+            (-1, 0, 0),
+            (11, 12, 13),
+            False,
+        ),
+        (
+            (0, 12, 0),
+            (11, 12, 13),
+            False,
+        ),
+    ]
+)
+def test_index_in_bounds(
+    index: tuple[int, int, int],
+    shape: tuple[int, int, int],
+    result_expected: bool,
+):
+    result = _arrays.index_in_bounds(index=index, shape=shape)
+
+    assert result == result_expected
+
+
+@pytest.mark.parametrize(
+    argnames="index",
+    argvalues=[
+        (1, 1, 1),
+        (5, 6, 7),
+    ]
+)
+@pytest.mark.parametrize(
+    argnames="shape",
+    argvalues=[
+        (11, 12, 13),
+    ]
+)
+def test_index_flat(
+    index: tuple[int, int, int],
+    shape: tuple[int, int, int],
+):
+    result = _arrays.index_flat(index=index, shape=shape)
+    result_expected = np.ravel_multi_index(index, shape)
+
+    assert result == result_expected
+
+
+@pytest.mark.parametrize(
+    argnames="index",
+    argvalues=[
+        0,
+        11,
+        25
+    ]
+)
+@pytest.mark.parametrize(
+    argnames="shape",
+    argvalues=[
+        (11, 12, 13),
+    ]
+)
+def test_index_3d(
+    index: int,
+    shape: tuple[int, int, int],
+):
+    result = _arrays.index_3d(index=index, shape=shape)
+    result_expected = np.unravel_index(index, shape)
+
+    print(f"{result=}")
+    print(f"{result_expected=}")
+
+    assert result == result_expected
+    assert index == np.ravel_multi_index(result, shape)

From 42cd564f43e046ca75d089a8e6078b52b992695e Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Fri, 2 May 2025 11:32:04 -0600
Subject: [PATCH 05/26] Added `_volumes` module.

---
 .../_weights_conservative_3d/_volumes.py      | 138 ++++++++++++++++++
 .../_weights_conservative_3d/_volumes_test.py |  44 ++++++
 2 files changed, 182 insertions(+)
 create mode 100644 regridding/_weights/_weights_conservative_3d/_volumes.py
 create mode 100644 regridding/_weights/_weights_conservative_3d/_volumes_test.py

diff --git a/regridding/_weights/_weights_conservative_3d/_volumes.py b/regridding/_weights/_weights_conservative_3d/_volumes.py
new file mode 100644
index 0000000..1255c73
--- /dev/null
+++ b/regridding/_weights/_weights_conservative_3d/_volumes.py
@@ -0,0 +1,138 @@
+import numpy as np
+import numba
+import regridding as rg
+from . import _arrays
+
+__all__ = [
+    "volume_tetrahedron",
+    "volume_grid",
+]
+
+
+@numba.njit(cache=True)
+def volume_tetrahedron(
+    v0: tuple[float, float, float],
+    v1: tuple[float, float, float],
+    v2: tuple[float, float, float],
+) -> float:
+    """
+    Compute the volume of a tetrahedron constructed from three
+    vertices and the origin.
+
+    Parameters
+    ----------
+    v0
+        First vertex of the tetrahedron.
+    v1
+        Second vertex of the tetrahedron.
+    v2
+        Third vertex of the tetrahedron.
+    """
+    triple_product = rg.math.dot_3d(
+        a=v0,
+        b=rg.math.cross_3d(v1, v2),
+    )
+
+    return triple_product / 6
+
+
+@numba.njit(cache=True)
+def volume_grid(
+    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
+) -> np.ndarray:
+    """
+    Compute the volume of each cell in a logically-rectangular grid.
+
+    Parameters
+    ----------
+    grid
+        A 3D grid of cell vertices.
+    """
+
+    x, y, z = grid
+
+    num_i, num_j, num_k = x.shape
+
+    result = np.zeros((num_i - 1, num_j - 1, num_k - 1))
+
+    for axis in (0, 1, 2):
+        _volume_grid_sweep(
+            grid=grid,
+            out=result,
+            axis=axis,
+        )
+
+    return result
+
+
+@numba.njit(cache=True)
+def _volume_grid_sweep(
+    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
+    out: np.ndarray,
+    axis: int,
+) -> None:
+    """
+    Compute the volume contribution of this axis.
+
+    Parameters
+    ----------
+    grid
+        A 3D grid of cell vertices.
+    out
+        An output to array to store the result.
+    axis
+        The axis along which to iterate.
+    """
+
+    x, y, z = grid
+
+    x = _arrays.align_axis_right(x, axis)
+    y = _arrays.align_axis_right(y, axis)
+    z = _arrays.align_axis_right(z, axis)
+
+    out = _arrays.align_axis_right(out, axis)
+
+    num_i, num_j, num_k = x.shape
+
+    for i in numba.prange(num_i):
+
+        i = numba.types.int64(i)
+
+        for j in range(num_j):
+
+            for k in range(num_k - 1):
+
+                i_left = i - 1
+                i_right = i
+
+                j_lower = j - 1
+                j_upper = j
+
+                k1 = k
+                k2 = k + 1
+
+                v1 = x[i, j, k1], y[i, j, k1], z[i, j, k1]
+                v2 = x[i, j, k2], y[i, j, k2], z[i, j, k2]
+
+                volume_left = 0
+                volume_lower = 0
+
+                if i_left >= 0:
+                    v0_left = x[i_left, j, k], y[i_left, j, k], z[i_left, j, k]
+                    volume_left = volume_tetrahedron(v0_left, v1, v2)
+
+                if j_lower >= 0:
+                    v0_lower = x[i, j_lower, k], y[i, j_lower, k], z[i, j_lower, k]
+                    volume_lower = -volume_tetrahedron(v0_lower, v1, v2)
+
+                if i_left >= 0:
+                    if j_lower >= 0:
+                        out[i_left, j_lower, k] -= volume_left
+                    if j_upper < (num_j - 1):
+                        out[i_left, j_upper, k] += volume_left
+
+                if j_lower >= 0:
+                    if i_left >= 0:
+                        out[i_left, j_lower, k] -= volume_lower
+                    if i_right < (num_i - 1):
+                        out[i_right, j_lower, k] += volume_lower
diff --git a/regridding/_weights/_weights_conservative_3d/_volumes_test.py b/regridding/_weights/_weights_conservative_3d/_volumes_test.py
new file mode 100644
index 0000000..d72c309
--- /dev/null
+++ b/regridding/_weights/_weights_conservative_3d/_volumes_test.py
@@ -0,0 +1,44 @@
+import pytest
+import numpy as np
+from . import _volumes
+
+
+@pytest.mark.parametrize(
+    argnames="grid,result_expected",
+    argvalues=[
+        (
+            np.meshgrid(
+                np.arange(2),
+                np.arange(2),
+                np.arange(2),
+                indexing="ij",
+            ),
+            np.ones((1, 1, 1))
+        ),
+        (
+            np.meshgrid(
+                np.arange(3),
+                np.arange(4),
+                np.arange(5),
+                indexing="ij",
+            ),
+            np.ones((2, 3, 4))
+        ),
+        (
+            np.meshgrid(
+                2 * np.arange(3),
+                2 * np.arange(4),
+                2 * np.arange(5),
+                indexing="ij",
+            ),
+            8 * np.ones((2, 3, 4))
+        ),
+    ]
+)
+def test_volume_grid(
+    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
+    result_expected: np.ndarray,
+):
+    result = _volumes.volume_grid(grid)
+    assert np.allclose(result, result_expected)
+    assert result.shape == result_expected.shape

From 0991d1f1afa945cc81b9deab1f07c87560c8c2df Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Fri, 2 May 2025 13:54:43 -0600
Subject: [PATCH 06/26] Added `_grids` module.

---
 .../_weights_conservative_3d/_grids.py        | 29 +++++++++++++++++++
 1 file changed, 29 insertions(+)
 create mode 100644 regridding/_weights/_weights_conservative_3d/_grids.py

diff --git a/regridding/_weights/_weights_conservative_3d/_grids.py b/regridding/_weights/_weights_conservative_3d/_grids.py
new file mode 100644
index 0000000..0573bc4
--- /dev/null
+++ b/regridding/_weights/_weights_conservative_3d/_grids.py
@@ -0,0 +1,29 @@
+"""
+Utilities for inspecting and searching logically-rectangular grids of
+coordinates.
+"""
+
+import numba
+
+__all__ = [
+    "shape_centers",
+]
+
+
+@numba.njit(cache=True)
+def shape_centers(
+    shape: tuple[int, int, int],
+) -> tuple[int, int, int]:
+    """
+    Given the shape of the grid of cell vertices,
+    compute the shape of the grid of cell centers.
+
+    Parameters
+    ----------
+    shape
+        The shape of the grid of cell vertices.
+    """
+
+    nx, ny, nz = shape
+
+    return nx - 1, ny - 1, nz - 1

From b708151d909890b5b0969e707d59d32873e9aeac Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Fri, 2 May 2025 13:56:19 -0600
Subject: [PATCH 07/26] Added `_intercepts` module

---
 .../_weights_conservative_3d/_intercepts.py   | 216 ++++++++++++++++++
 .../_intercepts_test.py                       |  64 ++++++
 2 files changed, 280 insertions(+)
 create mode 100644 regridding/_weights/_weights_conservative_3d/_intercepts.py
 create mode 100644 regridding/_weights/_weights_conservative_3d/_intercepts_test.py

diff --git a/regridding/_weights/_weights_conservative_3d/_intercepts.py b/regridding/_weights/_weights_conservative_3d/_intercepts.py
new file mode 100644
index 0000000..6519f3b
--- /dev/null
+++ b/regridding/_weights/_weights_conservative_3d/_intercepts.py
@@ -0,0 +1,216 @@
+"""
+Create lists of intercepts generated during axis sweeps and determine their
+contribution to the final list of weights.
+"""
+
+import numpy as np
+import numba
+import regridding as rg
+from . import _arrays
+from . import _volumes
+from . import _grids
+
+
+@numba.njit(cache=True)
+def empty(
+    shape_input: tuple[int, int, int],
+    shape_output: tuple[int, int, int],
+) -> numba.typed.List:
+    """
+    Create an empty list of intercepts between each plane in the
+    input grid and each plane in the output grid.
+
+    There are three planes in the input grid and thre
+
+    Parameters
+    ----------
+    shape_input
+        The shape of the input grid.
+    shape_output
+        The shape of the output grid.
+    """
+    intercepts = numba.typed.List()
+    for a in range(len(shape_input)):
+        intercepts_a = numba.typed.List()
+        for b in range(len(shape_output)):
+            intercepts_ab = numba.typed.List()
+            for i in range(shape_input[a]):
+                intercepts_abi = numba.typed.List()
+                for j in range(shape_output[b]):
+                    intercepts_abij = numba.typed.List()
+                    for x in range(0):
+                        intercepts_abij.append((0, 0, (0.0, 0.0, 0.0)))
+                    intercepts_abi.append(intercepts_abij)
+                intercepts_ab.append(intercepts_abi)
+            intercepts_a.append(intercepts_ab)
+        intercepts.append(intercepts_a)
+    return intercepts
+
+
+@numba.njit(cache=True)
+def sweep(
+    intercepts: numba.typed.List,
+    weights: numba.typed.List,
+    grid_input: tuple[np.ndarray, np.ndarray, np.ndarray],
+    grid_output: tuple[np.ndarray, np.ndarray, np.ndarray],
+    volume_input: np.ndarray,
+) -> None:
+    """
+    Sweep through the list of intercepts and add the volume of each
+    segment to the list of weights.
+
+    Parameters
+    ----------
+    intercepts
+        A sequence of intercepts and the corresponding indices in the
+        input/output grids.
+    weights
+        The current list of weights.
+        New weights will be appended to this list.
+    grid_input
+        The vertices of the input grid.
+    grid_output
+        The vertices of the output grid.
+    volume_input
+        The volume of each cell in the input grid.
+    """
+
+    x_input, y_input, z_input = grid_input
+    x_output, y_output, z_output = grid_output
+
+    shape_input = x_input.shape
+    shape_output = x_output.shape
+
+    shape_centers_input = _grids.shape_centers(shape_input)
+    shape_centers_output = _grids.shape_centers(shape_output)
+
+    for a in range(len(intercepts)):
+
+        intercepts_a = intercepts[a]
+
+        normal_a = [1 if axis == a else 0 for axis in range(len(intercepts))]
+
+        for b in range(len(intercepts_a)):
+
+            intercepts_ab = intercepts_a[b]
+
+            normal_b = [1 if axis == b else 0 for axis in range(len(intercepts))]
+
+            for i in range(len(intercepts_ab)):
+
+                intercepts_abi = intercepts_ab[i]
+
+                for j in range(len(intercepts_abi)):
+
+                    intercepts_abij = intercepts_abi[j]
+
+                    if not intercepts_abij:
+                        continue
+
+                    i0_input, i0_output, p0 = intercepts_abij[0]
+
+                    i0_input = _arrays.index_3d(i0_input, shape_input)
+                    i0_output = _arrays.index_3d(i0_output, shape_output)
+
+                    for t in range(1, len(intercepts_abij)):
+
+                        i1_input, i1_output, p1 = intercepts_abij[t]
+
+                        i1_input = _arrays.index_3d(i1_input, shape_input)
+                        i1_output = _arrays.index_3d(i1_output, shape_output)
+
+                        i0_input_lower = rg.math.difference_3d(i0_input, normal_a)
+                        i0_input_upper = i0_input
+
+                        i0_output_lower = rg.math.difference_3d(i0_output, normal_b)
+                        i0_output_upper = i0_output
+
+                        apex_input = (
+                            x_input[i0_input],
+                            y_input[i0_input],
+                            z_input[i0_input],
+                        )
+                        apex_output = (
+                            x_output[i0_output],
+                            y_output[i0_output],
+                            z_output[i0_output],
+                        )
+
+                        vol_input = _volumes.volume_tetrahedron(apex_input, p0, p1)
+                        vol_output = _volumes.volume_tetrahedron(apex_output, p0, p1)
+
+                        volume_input_lower = volume_input[i0_input_lower]
+                        volume_input_upper = volume_input[i0_output_upper]
+
+                        input_lower_in_bounds = _arrays.index_in_bounds(
+                            index=i0_input_lower,
+                            shape=shape_centers_input,
+                        )
+                        input_upper_in_bounds = _arrays.index_in_bounds(
+                            index=i0_input_upper,
+                            shape=shape_centers_input,
+                        )
+                        output_lower_in_bounds = _arrays.index_in_bounds(
+                            index=i0_output_lower,
+                            shape=shape_centers_output,
+                        )
+                        output_upper_in_bounds = _arrays.index_in_bounds(
+                            index=i0_output_upper,
+                            shape=shape_centers_output,
+                        )
+
+                        i0_input_lower = _arrays.index_flat(
+                            index=i0_input_lower,
+                            shape=shape_input,
+                        )
+                        i0_input_upper = _arrays.index_flat(
+                            index=i0_input_upper,
+                            shape=shape_input,
+                        )
+
+                        i0_output_lower = _arrays.index_flat(
+                            index=i0_output_lower,
+                            shape=shape_output,
+                        )
+                        i0_output_upper = _arrays.index_flat(
+                            index=i0_output_upper,
+                            shape=shape_output,
+                        )
+
+                        if input_lower_in_bounds:
+                            if output_lower_in_bounds:
+                                weight_lower_lower = (
+                                    i0_input_lower,
+                                    i0_output_lower,
+                                    (-vol_input - vol_output) / volume_input_lower,
+                                )
+                                weights.append(weight_lower_lower)
+
+                            if output_upper_in_bounds:
+                                weight_lower_upper = (
+                                    i0_input_lower,
+                                    i0_output_upper,
+                                    (-vol_input + vol_output) / volume_input_lower,
+                                )
+                                weights.append(weight_lower_upper)
+
+                        if input_upper_in_bounds:
+                            if output_lower_in_bounds:
+                                weight_upper_lower = (
+                                    i0_input_upper,
+                                    i0_output_lower,
+                                    (vol_input - vol_output) / volume_input_upper,
+                                )
+                                weights.append(weight_upper_lower)
+
+                            if output_upper_in_bounds:
+                                weight_upper_upper = (
+                                    i0_input_upper,
+                                    i0_output_upper,
+                                    (vol_input + vol_output) / volume_input_upper,
+                                )
+                                weights.append(weight_upper_upper)
+
+                        p0 = p1
+                        i0_input = i1_input
+                        i0_output = i1_output
diff --git a/regridding/_weights/_weights_conservative_3d/_intercepts_test.py b/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
new file mode 100644
index 0000000..34d9b62
--- /dev/null
+++ b/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
@@ -0,0 +1,64 @@
+import pytest
+import numpy as np
+import numba
+from . import _intercepts
+
+
+@pytest.mark.parametrize(
+    argnames="grid_input,grid_output,volume_input",
+    argvalues=[
+        (
+            np.meshgrid(
+                np.arange(3),
+                np.arange(4),
+                np.arange(5),
+                indexing="ij",
+            ),
+            np.meshgrid(
+                np.arange(3) + 0.5,
+                np.arange(4) + 0.5,
+                np.arange(5) + 0.5,
+                indexing="ij",
+            ),
+            np.ones((2, 3, 4)),
+        )
+    ],
+)
+def test_sweep(
+    grid_input: tuple[np.ndarray, np.ndarray, np.ndarray],
+    grid_output: tuple[np.ndarray, np.ndarray, np.ndarray],
+    volume_input: np.ndarray,
+):
+
+    weights = _test_sweep_compiled(grid_input, grid_output, volume_input)
+
+    assert len(weights) > 0
+
+
+@numba.njit
+def _test_sweep_compiled(
+    grid_input: tuple[np.ndarray, np.ndarray, np.ndarray],
+    grid_output: tuple[np.ndarray, np.ndarray, np.ndarray],
+    volume_input: np.ndarray,
+):
+    intercepts = _intercepts.empty(
+        shape_input=grid_input[0].shape,
+        shape_output=grid_output[0].shape,
+    )
+
+    intercepts[0][0][0][0].append((0, 0, (0.05, 0.15, 0.35)))
+    intercepts[0][0][0][0].append((0, 1, (0.85, 0.25, 0.25)))
+
+    weights = numba.typed.List()
+    for _ in range(0):
+        weights.append((0, 0, 0.))
+
+    _intercepts.sweep(
+        intercepts=intercepts,
+        weights=weights,
+        grid_input=grid_input,
+        grid_output=grid_output,
+        volume_input=volume_input,
+    )
+
+    return weights

From 43c3bc3e470cd4747ad3b6e18d785815b5f516c2 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Fri, 2 May 2025 15:45:28 -0600
Subject: [PATCH 08/26] Moving stuff

---
 .../_weights_conservative_3d.py               | 612 ++++++++++--------
 .../_weights_conservative_3d_test.py          |  41 +-
 2 files changed, 345 insertions(+), 308 deletions(-)
 rename regridding/_weights/{ => _weights_conservative_3d}/_weights_conservative_3d.py (58%)
 rename regridding/_weights/{ => _weights_conservative_3d}/_weights_conservative_3d_test.py (63%)

diff --git a/regridding/_weights/_weights_conservative_3d.py b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
similarity index 58%
rename from regridding/_weights/_weights_conservative_3d.py
rename to regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
index 363e161..2a4f2c5 100644
--- a/regridding/_weights/_weights_conservative_3d.py
+++ b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
@@ -2,16 +2,19 @@
 import numpy as np
 import numpy.typing as npt
 import numba
-import regridding
+import regridding as rg
+from . import  (
+    _arrays,
+    _volumes,
+    _grids,
+    _intercepts,
+)
+from _arrays import axis_x, axis_y, axis_z
 
 __all__ = [
     "weights_conservative_3d",
 ]
 
-axis_x = 0
-axis_y = 1
-axis_z = 2
-
 
 @numba.njit(cache=True)
 def weights_conservative_3d(
@@ -43,9 +46,9 @@ def weights_conservative_3d(
     for x in range(0):
         weights.append((0.0, 0.0, 0.0))
 
-    volume_input = _cell_volume(grid_input)
+    volume_input = _volumes.volume_grid(grid_input)
 
-    intercepts = _empty_intercepts(shape_input, shape_output)
+    intercepts = _intercepts.empty(shape_input, shape_output)
 
     for sweep_input in (False, True):
         _sweep_grid(
@@ -111,7 +114,7 @@ def _sweep_grid(
         (grid_static_x.max(), grid_static_y.max(), grid_static_z.max()),
     )
 
-    index_boundary, coord_boundary = _calc_boundary(grid_static)
+    indices_boundary, triangles_boundary = _indices_and_triangles_of_boundary(grid_static)
 
     axes_parallel = (
         axis_x,
@@ -140,9 +143,7 @@ def _sweep_grid(
         )
 
     for axis in axes_diagonal:
-        _sweep_along_diagonal(
-
-        )
+        _sweep_along_diagonal()
 
 
 @numba.njit(cache=True, parallel=True)
@@ -151,7 +152,7 @@ def _sweep_along_axis(
     grid_output: tuple[np.ndarray, np.ndarray, np.ndarray],
     volume_input: np.ndarray,
     bbox_boundary: tuple[tuple[float, float, float], tuple[float, float, float]],
-    index_boundary: tuple[np.ndarray, np.ndarray, np.ndarray],
+    indices_boundary: tuple[np.ndarray, np.ndarray, np.ndarray],
     coord_boundary: tuple[np.ndarray, np.ndarray, np.ndarray],
     weights: numba.typed.List[tuple[int, int, float]],
     intercepts: numba.typed.List,
@@ -180,26 +181,37 @@ def _sweep_along_axis(
     weights
         The current list of weights.
         New weights will be appended to this list.
+    intercepts
+        A sorted list of intercepts to be traversed later.
+        As new intercepts are found, they are inserted into the list.
     sweep_input
         A boolean flag indicating whether to iterate along `grid_input`
         or `grid_output`.
         If :obj:`True`, this function sweeps along `grid_input`.
     axis
-        The logical axis of the `grid_sweep` to iterate along.
+        The logical axis of the sweep grid to iterate along.
     """
 
-    grid_sweep, grid_static = _grid_sweep_static(grid_input, grid_output, sweep_input)
-
-    grid_sweep_x, grid_sweep_y, grid_sweep_z = grid_sweep
-    grid_static_x, grid_static_y, grid_static_z = grid_static
+    grid_sweep, grid_static = _grid_sweep_static(
+        grid_input=grid_input,
+        grid_output=grid_output,
+        sweep_input=sweep_input,
+        axis_x=axis,
+    )
 
-    grid_sweep_x = _align_axis_right(grid_sweep_x, axis)
-    grid_sweep_y = _align_axis_right(grid_sweep_y, axis)
-    grid_sweep_z = _align_axis_right(grid_sweep_z, axis)
+    x_sweep, y_sweep, z_sweep = grid_sweep
 
     shape_sweep_x, shape_sweep_y, shape_sweep_z = grid_sweep.shape
     shape_static_x, shape_static_y, shape_static_z = grid_static.shape
 
+    weight = numba.typed.List()
+
+    for i in range(shape_sweep_x):
+        weight_i = numba.typed.List()
+        for _ in range(0):
+            weight_i.append((0., 0., 0.))
+        weight.append(weight_i)
+
     for i in numba.prange(shape_sweep_x):
 
         i = numba.types.int64(i)
@@ -215,16 +227,16 @@ def _sweep_along_axis(
             sweep_is_outside_static = True
 
             point_1 = (
-                grid_sweep_x[index_sweep],
-                grid_sweep_y[index_sweep],
-                grid_sweep_z[index_sweep],
+                x_sweep[index_sweep],
+                y_sweep[index_sweep],
+                z_sweep[index_sweep],
             )
 
-            if regridding.geometry.point_is_inside_box_3d(
+            if rg.geometry.point_is_inside_box_3d(
                 point=point_1,
                 box=bbox_boundary,
             ):
-                if regridding.geometry.point_is_inside_polyhedron(
+                if rg.geometry.point_is_inside_polyhedron(
                     point=point_1,
                     polyhedron=coord_boundary,
                 ):
@@ -232,17 +244,16 @@ def _sweep_along_axis(
                         point=point_1,
                         grid=grid_static,
                     )
+                    sweep_is_outside_static = False
 
-            while True:
+            while k < (shape_sweep_z - 1):
 
                 point_2 = (
-                    grid_sweep_x[i, j, k + 1],
-                    grid_sweep_y[i, j, k + 1],
-                    grid_sweep_z[i, j, k + 1],
+                    x_sweep[i, j, k + 1],
+                    y_sweep[i, j, k + 1],
+                    z_sweep[i, j, k + 1],
                 )
 
-                sweep_is_outside_static = index_static[0] == sys.maxsize
-
                 if sweep_is_outside_static:
 
                     point_2, k, index_static = _step_outside_static(
@@ -254,6 +265,8 @@ def _sweep_along_axis(
                         intercepts=intercepts,
                     )
 
+                    sweep_is_outside_static = index_static[0] == sys.maxsize
+
                 else:
 
                     point_2, k, index_static = _step_inside_static(
@@ -272,29 +285,33 @@ def _sweep_along_axis(
                     if not 0 <= index_static[axis_z] < (shape_static_z - 1):
                         break
 
-                if k >= (shape_sweep_z - 1):
-                    break
-
-
-
-
-
-
+    for i in range(shape_sweep_x):
+        weight_i = weight[i]
+        for w in range(len(weight_i)):
+            weights.append(weight_i[w])
 
 
 @numba.njit(cache=True)
 def _step_outside_static(
     point_1: tuple[float, float, float],
     point_2: tuple[float, float, float],
-    index_static: tuple[int, int, int],
     index_sweep: tuple[int, int, int],
     grid_static: tuple[np.ndarray, np.ndarray, np.ndarray],
     bbox_boundary: tuple[tuple[float, float, float], tuple[float, float, float]],
-    coord_boundary: tuple[np.ndarray, np.ndarray, np.ndarray],
-    index_boundary: tuple[np.ndarray, np.ndarray, np.ndarray],
-) -> tuple[
-    tuple[float, float, float],
-]:
+    indices_boundary: numba.typed.List[
+        tuple[int, int, int],
+    ],
+    triangles_boundary: numba.typed.List[
+        tuple[
+            tuple[float, float, float],
+            tuple[float, float, float],
+            tuple[float, float, float],
+        ],
+    ],
+    intercepts: numba.typed.List,
+    sweep_input: bool,
+    axis: int,
+) -> tuple[tuple[float, float, float], int, tuple[int, int, int]]:
     """
     Check if the current line segment crosses the boundary.
     If it does, compute the intersection and return the intersection point
@@ -306,26 +323,58 @@ def _step_outside_static(
         The first point of the line segment.
     point_2
         The second point of the line segment.
-    index_static
-        The current 3D index in the static grid.
-    index_sweep
-        The current 3D index in the sweep grid.
-    grid_static
-        Coordinates of the static grid.
     bbox_boundary
         Two points defining the bounding box of the static grid.
-    coord_boundary
-        A sequence of triangles defining the outer surface of the static grid.
-    index_boundary
-        The index of the static grid corresponding to each triangle in
-        `coord_boundary`.
+    indices_boundary
+        A sequence of triangle indices defining the outer boundary of the
+        static grid.
+    intercepts
+        A sorted list of intercepts to be traversed later.
+        As new intercepts are found, they are inserted into the list.
     """
 
+    line = (point_1, point_2)
+
+    x, y, z = grid_static
+
+    if rg.geometry.point_is_inside_box_3d(point_2, bbox_boundary):
+
+        for t in range(len(triangles_boundary)):
+
+            triangle = triangles_boundary[t]
+
+            tuv = rg.geometry.line_triangle_intersection_parameters(
+                line=line,
+                triangle=triangle
+            )
+
+            if rg.geometry.line_intersects_triangle(tuv):
+
+                index_static = indices_boundary[t]
+
+                intercept = rg.geometry.line_triangle_intersection(
+                    line=line,
+                    tuv=tuv,
+                )
+
+                index_input, index_output = _index_input_output(
+                    index_sweep=index_sweep,
+                    index_static=index_static,
+                    sweep_input=sweep_input,
+                    axis=axis,
+                )
+
+                if axis == axis_x:
+
+                    pass
+
 
 @numba.njit(cache=True)
 def _step_inside_static(
-    point_1: tuple[float, float, float],
-    point_2: tuple[float, float, float],
+    line: tuple[
+        tuple[float, float, float],
+        tuple[float, float, float],
+    ],
     index_static: tuple[int, int, int],
     index_sweep: tuple[int, int, int],
     grid_static: tuple[np.ndarray, np.ndarray, np.ndarray],
@@ -354,167 +403,194 @@ def _step_inside_static(
         New weights will be appended to this list.
     """
 
+    point_1, point_2 = line
+
+    i, j, k = index_sweep
+    l, m, n = index_static
+
+    x, y, z = grid_static
+
+    i_000 = l + 0, m + 0, n + 0
+    i_001 = l + 0, m + 0, n + 1
+    i_010 = l + 0, m + 1, n + 0
+    i_011 = l + 0, m + 1, n + 1
+    i_100 = l + 1, m + 0, n + 0
+    i_101 = l + 1, m + 0, n + 1
+    i_110 = l + 1, m + 1, n + 0
+    i_111 = l + 1, m + 1, n + 1
+
+    indices_triangles = (
+        (i_000, i_010, i_110),
+        (i_110, i_100, i_000),
+        (i_001, i_101, i_111),
+        (i_111, i_011, i_001),
+        (i_000, i_001, i_011),
+        (i_011, i_010, i_000),
+        (i_100, i_110, i_111),
+        (i_111, i_101, i_100),
+        (i_000, i_100, i_101),
+        (i_101, i_001, i_000),
+        (i_010, i_011, i_111),
+        (i_111, i_110, i_010),
+    )
 
-@numba.njit(cache=True)
-def _grid_sweep_static(
-    grid_input: tuple[np.ndarray, np.ndarray, np.ndarray],
-    grid_output: tuple[np.ndarray, np.ndarray, np.ndarray],
-    sweep_input: bool,
-) -> tuple[
-    tuple[np.ndarray, np.ndarray, np.ndarray],
-    tuple[np.ndarray, np.ndarray, np.ndarray],
-]:
-    if sweep_input:
-        grid_sweep = grid_input
-        grid_static = grid_output
-    else:
-        grid_sweep = grid_output
-        grid_static = grid_input
+    normals_index = (
+        (0, 0, -1),
+        (0, 0, -1),
+        (0, 0, +1),
+        (0, 0, +1),
+        (-1, 0, 0),
+        (-1, 0, 0),
+        (+1, 0, 0),
+        (+1, 0, 0),
+        (0, -1, 0),
+        (0, -1, 0),
+        (0, +1, 0),
+        (0, +1, 0),
+    )
 
-    return grid_sweep, grid_static
+    for t in range(len(indices_triangles)):
 
+        i0, i1, i2 = indices_triangles[t]
 
-@numba.njit(cache=True)
-def _cell_volume(
-    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
-) -> np.ndarray:
-    """
-    Compute the volume of each cell in a logically-rectangular grid.
+        v0 = (x[i0], y[i0], z[i0])
+        v1 = (x[i1], y[i1], z[i1])
+        v2 = (x[i2], y[i2], z[i2])
 
-    Parameters
-    ----------
-    grid
-        A 3D grid of cell vertices.
-    """
+        triangle = (v0, v1, v2)
 
-    x, y, z = grid
+        tuv = rg.geometry.line_triangle_intersection_parameters(
+            line=line,
+            triangle=triangle,
+        )
 
-    num_i, num_j, num_k = x.shape
+        if rg.geometry.line_intersects_triangle(tuv):
 
-    result = np.zeros((num_i - 1, num_j - 1, num_k - 1))
+            intercept = rg.geometry.line_triangle_intersection(
+                line=line,
+                tuv=tuv,
+            )
 
-    for axis in (0, 1, 2):
-        _cell_volume_sweep(
-            grid=grid,
-            out=result,
-            axis=axis,
-        )
+            line = (point_1, intercept)
 
-    return result
+            normal_index = normals_index[t]
 
+            index_static = rg.math.sum_3d(index_static, normal_index)
 
-@numba.njit(cache=True)
-def _cell_volume_sweep(
-    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
-    out: np.ndarray,
-    axis: int,
-) -> None:
-    """
-    Compute the volume contribution of this axis.
 
-    Parameters
-    ----------
-    grid
-        A 3D grid of cell vertices.
-    out
-        An output to array to store the result.
-    axis
-        The axis along which to iterate.
-    """
 
-    x, y, z = grid
+        else:
+            index_sweep = i, j, k + 1
 
-    x = _align_axis_right(x, axis)
-    y = _align_axis_right(y, axis)
-    z = _align_axis_right(z, axis)
 
-    out = _align_axis_right(out, axis)
 
-    num_i, num_j, num_k = x.shape
+        # index_left
+        #
+        # volume_sweep =
+        #
+        # i_normal, j_normal, k_normal = normals_index[t]
 
-    for i in numba.prange(num_i):
 
-        i = numba.types.int64(i)
 
-        for j in range(num_j):
 
-            for k in range(num_k - 1):
 
-                i_left = i - 1
-                i_right = i
 
-                j_lower = j - 1
-                j_upper = j
 
-                k1 = k
-                k2 = k + 1
 
-                v1 = x[i, j, k1], y[i, j, k1], z[i, j, k1]
-                v2 = x[i, j, k2], y[i, j, k2], z[i, j, k2]
+@numba.njit(cache=True)
+def _grid_sweep_static(
+    grid_input: tuple[np.ndarray, np.ndarray, np.ndarray],
+    grid_output: tuple[np.ndarray, np.ndarray, np.ndarray],
+    sweep_input: bool,
+    axis: int,
+) -> tuple[
+    tuple[np.ndarray, np.ndarray, np.ndarray],
+    tuple[np.ndarray, np.ndarray, np.ndarray],
+]:
+    """
+    Given the input and output grids,
+    prepare the sweep and static grids.
 
-                volume_left = 0
-                volume_lower = 0
+    Parameters
+    ----------
+    grid_input
+        The vertices of the old grid.
+    grid_output
+        The vertices of the new grid.
+    sweep_input
+        A boolean flag indicating whether to iterate along `grid_input`
+        or `grid_output`.
+        If :obj:`True`, this function sweeps along `grid_input`.
+    axis
+        The logical axis of the `grid_sweep` to iterate along.
+    """
+    if sweep_input:
+        grid_sweep = grid_input
+        grid_static = grid_output
+    else:
+        grid_sweep = grid_output
+        grid_static = grid_input
 
-                if i_left >= 0:
-                    v0_left = x[i_left, j, k], y[i_left, j, k], z[i_left, j, k]
-                    volume_left = _volume(v0_left, v1, v2)
+    grid_sweep_x, grid_sweep_y, grid_sweep_z = grid_sweep
 
-                if j_lower >= 0:
-                    v0_lower = x[i, j_lower, k], y[i, j_lower, k], z[i, j_lower, k]
-                    volume_lower = -_volume(v0_lower, v1, v2)
+    grid_sweep_x = _arrays.align_axis_right(grid_sweep_x, axis)
+    grid_sweep_y = _arrays.align_axis_right(grid_sweep_y, axis)
+    grid_sweep_z = _arrays.align_axis_right(grid_sweep_z, axis)
 
-                if i_left >= 0:
-                    if j_lower >= 0:
-                        out[i_left, j_lower, k] -= volume_left
-                    if j_upper < (num_j - 1):
-                        out[i_left, j_upper, k] += volume_left
+    grid_sweep = grid_sweep_x, grid_sweep_y, grid_sweep_z
 
-                if j_lower >= 0:
-                    if i_left >= 0:
-                        out[i_left, j_lower, k] -= volume_lower
-                    if i_right < (num_i - 1):
-                        out[i_right, j_lower, k] += volume_lower
+    return grid_sweep, grid_static
 
 
 @numba.njit(cache=True)
-def _volume(
-    v0: tuple[float, float, float],
-    v1: tuple[float, float, float],
-    v2: tuple[float, float, float],
-) -> float:
-    """
-    Compute the volume of a tetrahedron constructed from three
-    vertices and the origin.
+def _index_input_output(
+    index_sweep: tuple[int, int, int],
+    index_static: tuple[int, int, int],
+    sweep_input: bool,
+    axis: int,
+):
+
+    i_sweep, j_sweep, k_sweep = index_sweep
+
+    if axis == axis_x:
+        index_sweep = j_sweep, k_sweep, i_sweep
+    elif axis == axis_y:
+        index_sweep = k_sweep, i_sweep, j_sweep
+
+    if sweep_input:
+        index_input = index_sweep
+        index_output = index_static
+    else:
+        index_input = index_static
+        index_output = index_sweep
+
+    return index_input, index_output
+
+
 
-    Parameters
-    ----------
-    v0
-        First vertex of the tetrahedron.
-    v1
-        Second vertex of the tetrahedron.
-    v2
-        Third vertex of the tetrahedron.
-    """
-    triple_product = regridding.math.dot_3d(
-        a=v0,
-        b=regridding.math.cross_3d(v1, v2),
-    )
 
-    return triple_product / 6
 
 
 @numba.njit(cache=True)
-def _indices_boundary(
+def _indices_and_triangles_of_boundary(
     grid: tuple[np.ndarray, np.ndarray, np.ndarray],
-) -> numba.typed.List[
-    tuple[
-        tuple[int, int, int],
-        tuple[int, int, int],
+) -> tuple[
+    numba.typed.List[
         tuple[int, int, int],
     ],
+    numba.typed.List[
+        tuple[
+            tuple[float, float, float],
+            tuple[float, float, float],
+            tuple[float, float, float],
+        ],
+    ],
 ]:
     """
-    For a given grid, find the set of indices that expresses the boundary
+    For a given grid, find the indices of the boundary cells
+    and the vertices of all the triangles on the boundary.
+
+    find the set of indices that expresses the boundary
     as a sequence of triangles.
 
     Parameters
@@ -527,6 +603,8 @@ def _indices_boundary(
 
     shape_x, shape_y, shape_z = x.shape
 
+    indices = numba.typed.List()
+
     triangles = numba.typed.List()
 
     for i0 in range(shape_x - 1):
@@ -537,10 +615,20 @@ def _indices_boundary(
             i1 = i0 + 1
             j1 = j0 + 1
 
-            v_000 = i0, j0, k0
-            v_010 = i0, j1, k0
-            v_110 = i1, j1, k0
-            v_100 = i1, j0, k0
+            i_cell = i0, j0, k0
+
+            i_000 = i0, j0, k0
+            i_010 = i0, j1, k0
+            i_110 = i1, j1, k0
+            i_100 = i1, j0, k0
+
+            v_000 = x[i_000], y[i_000], z[i_000]
+            v_010 = x[i_010], y[i_010], z[i_010]
+            v_110 = x[i_110], y[i_110], z[i_110]
+            v_100 = x[i_100], y[i_100], z[i_100]
+
+            indices.append(i_cell)
+            indices.append(i_cell)
 
             triangles.append((v_000, v_010, v_110))
             triangles.append((v_110, v_100, v_000))
@@ -553,10 +641,20 @@ def _indices_boundary(
             i1 = i0 + 1
             j1 = j0 + 1
 
-            v_001 = i0, j0, k1
-            v_101 = i1, j0, k1
-            v_111 = i1, j1, k1
-            v_011 = i0, j1, k1
+            i_cell = i0, j0, k1 - 1
+
+            i_001 = i0, j0, k1
+            i_101 = i1, j0, k1
+            i_111 = i1, j1, k1
+            i_011 = i0, j1, k1
+
+            v_001 = x[i_001], y[i_001], z[i_001]
+            v_101 = x[i_101], y[i_101], z[i_101]
+            v_111 = x[i_111], y[i_111], z[i_111]
+            v_011 = x[i_011], y[i_011], z[i_011]
+
+            indices.append(i_cell)
+            indices.append(i_cell)
 
             triangles.append((v_001, v_101, v_111))
             triangles.append((v_111, v_011, v_001))
@@ -569,10 +667,20 @@ def _indices_boundary(
             j1 = j0 + 1
             k1 = k0 + 1
 
-            v_000 = i0, j0, k0
-            v_001 = i0, j0, k1
-            v_011 = i0, j1, k1
-            v_010 = i0, j1, k0
+            i_cell = i0, j0, k0
+
+            i_000 = i0, j0, k0
+            i_001 = i0, j0, k1
+            i_011 = i0, j1, k1
+            i_010 = i0, j1, k0
+
+            v_000 = x[i_000], y[i_000], z[i_000]
+            v_001 = x[i_001], y[i_001], z[i_001]
+            v_011 = x[i_011], y[i_011], z[i_011]
+            v_010 = x[i_010], y[i_010], z[i_010]
+
+            indices.append(i_cell)
+            indices.append(i_cell)
 
             triangles.append((v_000, v_001, v_011))
             triangles.append((v_011, v_010, v_000))
@@ -585,10 +693,20 @@ def _indices_boundary(
             j1 = j0 + 1
             k1 = k0 + 1
 
-            v_100 = i1, j0, k0
-            v_110 = i1, j1, k0
-            v_111 = i1, j1, k1
-            v_101 = i1, j0, k1
+            i_cell = i1 - 1, j0, k0
+
+            i_100 = i1, j0, k0
+            i_110 = i1, j1, k0
+            i_111 = i1, j1, k1
+            i_101 = i1, j0, k1
+
+            v_100 = x[i_100], y[i_100], z[i_100]
+            v_110 = x[i_110], y[i_110], z[i_110]
+            v_111 = x[i_111], y[i_111], z[i_111]
+            v_101 = x[i_101], y[i_101], z[i_101]
+
+            indices.append(i_cell)
+            indices.append(i_cell)
 
             triangles.append((v_100, v_110, v_111))
             triangles.append((v_111, v_101, v_100))
@@ -601,13 +719,23 @@ def _indices_boundary(
             i1 = i0 + 1
             k1 = k0 + 1
 
-            v_000 = i0, j0, k0
-            v_100 = i1, j0, k0
-            v_101 = i1, j0, k1
-            v_001 = i0, j0, k1
+            i_cell = i0, j0, k0
+
+            i_000 = i0, j0, k0
+            i_100 = i1, j0, k0
+            i_101 = i1, j0, k1
+            i_001 = i0, j0, k1
+
+            v_000 = x[i_000], y[i_000], z[i_000]
+            v_100 = x[i_100], y[i_100], z[i_100]
+            v_101 = x[i_101], y[i_101], z[i_101]
+            v_001 = x[i_001], y[i_001], z[i_001]
+
+            indices.append(i_cell)
+            indices.append(i_cell)
 
             triangles.append((v_000, v_100, v_101))
-            triangles.append(v_101, v_001, v_000)
+            triangles.append((v_101, v_001, v_000))
 
     for i0 in range(shape_x - 1):
         for k0 in range(shape_z - 1):
@@ -617,14 +745,26 @@ def _indices_boundary(
             i1 = i0 + 1
             k1 = k0 + 1
 
-            v_010 = i0, j1, k0
-            v_011 = i0, j1, k1
-            v_111 = i1, j1, k1
-            v_110 = i1, j1, k0
+            i_cell = i0, j1 - 1, k0
+
+            i_010 = i0, j1, k0
+            i_011 = i0, j1, k1
+            i_111 = i1, j1, k1
+            i_110 = i1, j1, k0
+
+            v_010 = x[i_010], y[i_010], z[i_010]
+            v_011 = x[i_011], y[i_011], z[i_011]
+            v_111 = x[i_111], y[i_111], z[i_111]
+            v_110 = x[i_110], y[i_110], z[i_110]
+
+            indices.append(i_cell)
+            indices.append(i_cell)
 
             triangles.append((v_010, v_011, v_111))
             triangles.append((v_111, v_110, v_010))
 
+    return indices, triangles
+
 
 @numba.njit(cache=True)
 def _index_of_point_brute(
@@ -692,72 +832,8 @@ def _index_of_point_brute(
 
                     polyhedron.append(triangle)
 
-                if regridding.geometry.point_is_inside_polyhedron(
-                    point=point, polyhedron=polyhedron
+                if rg.geometry.point_is_inside_polyhedron(
+                    point=point,
+                    polyhedron=polyhedron,
                 ):
                     return i, j, k
-
-
-@numba.njit(cache=True)
-def _empty_intercepts(
-    shape_input: tuple[int, int, int],
-    shape_output: tuple[int, int, int],
-) -> numba.typed.List:
-    """
-    Create an empty list of intercepts between each plane in the
-    input grid and each plane in the output grid.
-
-    Parameters
-    ----------
-    shape_input
-        The shape of the input grid.
-    shape_output
-        The shape of the output grid.
-    """
-    intercepts = numba.typed.List()
-    for a in range(len(shape_input)):
-        intercepts_a = numba.typed.List()
-        for b in range(len(shape_output)):
-            intercepts_ab = numba.typed.List()
-            for i in range(shape_input[a]):
-                intercepts_abi = numba.typed.List()
-                for j in range(shape_output[b]):
-                    intercepts_abij = numba.typed.List()
-                    for x in range(0):
-                        intercepts_abij.append((0.0, 0.0, 0.0))
-                    intercepts_abi.append(intercepts_abij)
-                intercepts_ab.append(intercepts_abi)
-            intercepts_a.append(intercepts_ab)
-        intercepts.append(intercepts_a)
-    return intercepts
-
-
-@numba.njit(cache=True)
-def _align_axis_right(
-    a: np.ndarray,
-    axis: int,
-) -> np.ndarray:
-    """
-    Roll all the axes of a 3D array to the right until `axis` is the last axis..
-
-    This function is needed to permute the axes in such a way to retain
-    their right-handedness.
-
-    Parameters
-    ----------
-    a
-        The array to modify the axes of.
-    axis
-        The axis to set as the last axis.
-    """
-    axis = axis % a.ndim
-
-    if axis == axis_x:
-        result = a.transpose((axis_y, axis_z, axis_x))
-    elif axis == axis_y:
-        result = a.transpose((axis_z, axis_x, axis_y))
-    else:
-        result = a
-
-    return result
-
diff --git a/regridding/_weights/_weights_conservative_3d_test.py b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d_test.py
similarity index 63%
rename from regridding/_weights/_weights_conservative_3d_test.py
rename to regridding/_weights/_weights_conservative_3d/_weights_conservative_3d_test.py
index ee2fdba..d665855 100644
--- a/regridding/_weights/_weights_conservative_3d_test.py
+++ b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d_test.py
@@ -38,46 +38,6 @@ def test_weights_conservative_3d(
     )
 
 
-@pytest.mark.parametrize(
-    argnames="grid,result_expected",
-    argvalues=[
-        (
-            np.meshgrid(
-                np.arange(2),
-                np.arange(2),
-                np.arange(2),
-                indexing="ij",
-            ),
-            np.ones((1, 1, 1))
-        ),
-        (
-            np.meshgrid(
-                np.arange(3),
-                np.arange(4),
-                np.arange(5),
-                indexing="ij",
-            ),
-            np.ones((2, 3, 4))
-        ),
-        (
-            np.meshgrid(
-                2 * np.arange(3),
-                2 * np.arange(4),
-                2 * np.arange(5),
-                indexing="ij",
-            ),
-            8 * np.ones((2, 3, 4))
-        ),
-    ]
-)
-def test_cell_volume(
-    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
-    result_expected: np.ndarray,
-):
-    result = _weights_conservative_3d._cell_volume(grid)
-    assert np.allclose(result, result_expected)
-    assert result.shape == result_expected.shape
-
 
 @pytest.mark.parametrize(
     argnames="point,grid,result_expected",
@@ -115,3 +75,4 @@ def test_index_of_point_brute(
     )
 
     assert result == result_expected
+

From 107d83d4ac6a0b6e6475aba96f4239da241e3b6c Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sun, 4 May 2025 07:02:53 -0600
Subject: [PATCH 09/26] Moving stuff around

---
 .../_weights_conservative_3d/_grids.py        | 105 +++++++++++++
 .../{_volumes_test.py => _grids_test.py}      |   0
 .../_weights_conservative_3d/_volumes.py      | 138 ------------------
 3 files changed, 105 insertions(+), 138 deletions(-)
 rename regridding/_weights/_weights_conservative_3d/{_volumes_test.py => _grids_test.py} (100%)
 delete mode 100644 regridding/_weights/_weights_conservative_3d/_volumes.py

diff --git a/regridding/_weights/_weights_conservative_3d/_grids.py b/regridding/_weights/_weights_conservative_3d/_grids.py
index 0573bc4..596965c 100644
--- a/regridding/_weights/_weights_conservative_3d/_grids.py
+++ b/regridding/_weights/_weights_conservative_3d/_grids.py
@@ -3,7 +3,10 @@
 coordinates.
 """
 
+import numpy as np
 import numba
+import regridding as rg
+from . import _arrays
 
 __all__ = [
     "shape_centers",
@@ -27,3 +30,105 @@ def shape_centers(
     nx, ny, nz = shape
 
     return nx - 1, ny - 1, nz - 1
+
+
+@numba.njit(cache=True)
+def grid_volume(
+    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
+) -> np.ndarray:
+    """
+    Compute the volume of each cell in a logically-rectangular grid.
+
+    Parameters
+    ----------
+    grid
+        A 3D grid of cell vertices.
+    """
+
+    x, y, z = grid
+
+    num_i, num_j, num_k = x.shape
+
+    result = np.zeros((num_i - 1, num_j - 1, num_k - 1))
+
+    for axis in (0, 1, 2):
+        _grid_volume_sweep(
+            grid=grid,
+            out=result,
+            axis=axis,
+        )
+
+    return result
+
+
+@numba.njit(cache=True)
+def _grid_volume_sweep(
+    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
+    out: np.ndarray,
+    axis: int,
+) -> None:
+    """
+    Compute the volume contribution of this axis.
+
+    Parameters
+    ----------
+    grid
+        A 3D grid of cell vertices.
+    out
+        An output to array to store the result.
+    axis
+        The axis along which to iterate.
+    """
+
+    x, y, z = grid
+
+    x = _arrays.align_axis_right(x, axis)
+    y = _arrays.align_axis_right(y, axis)
+    z = _arrays.align_axis_right(z, axis)
+
+    out = _arrays.align_axis_right(out, axis)
+
+    num_i, num_j, num_k = x.shape
+
+    for i in numba.prange(num_i):
+
+        i = numba.types.int64(i)
+
+        for j in range(num_j):
+
+            for k in range(num_k - 1):
+
+                i_left = i - 1
+                i_right = i
+
+                j_lower = j - 1
+                j_upper = j
+
+                k1 = k
+                k2 = k + 1
+
+                v1 = x[i, j, k1], y[i, j, k1], z[i, j, k1]
+                v2 = x[i, j, k2], y[i, j, k2], z[i, j, k2]
+
+                volume_left = 0
+                volume_lower = 0
+
+                if i_left >= 0:
+                    v0_left = x[i_left, j, k], y[i_left, j, k], z[i_left, j, k]
+                    volume_left = rg.geometry.volume_tetrahedron(v0_left, v1, v2)
+
+                if j_lower >= 0:
+                    v0_lower = x[i, j_lower, k], y[i, j_lower, k], z[i, j_lower, k]
+                    volume_lower = -rg.geometry.volume_tetrahedron(v0_lower, v1, v2)
+
+                if i_left >= 0:
+                    if j_lower >= 0:
+                        out[i_left, j_lower, k] -= volume_left
+                    if j_upper < (num_j - 1):
+                        out[i_left, j_upper, k] += volume_left
+
+                if j_lower >= 0:
+                    if i_left >= 0:
+                        out[i_left, j_lower, k] -= volume_lower
+                    if i_right < (num_i - 1):
+                        out[i_right, j_lower, k] += volume_lower
diff --git a/regridding/_weights/_weights_conservative_3d/_volumes_test.py b/regridding/_weights/_weights_conservative_3d/_grids_test.py
similarity index 100%
rename from regridding/_weights/_weights_conservative_3d/_volumes_test.py
rename to regridding/_weights/_weights_conservative_3d/_grids_test.py
diff --git a/regridding/_weights/_weights_conservative_3d/_volumes.py b/regridding/_weights/_weights_conservative_3d/_volumes.py
deleted file mode 100644
index 1255c73..0000000
--- a/regridding/_weights/_weights_conservative_3d/_volumes.py
+++ /dev/null
@@ -1,138 +0,0 @@
-import numpy as np
-import numba
-import regridding as rg
-from . import _arrays
-
-__all__ = [
-    "volume_tetrahedron",
-    "volume_grid",
-]
-
-
-@numba.njit(cache=True)
-def volume_tetrahedron(
-    v0: tuple[float, float, float],
-    v1: tuple[float, float, float],
-    v2: tuple[float, float, float],
-) -> float:
-    """
-    Compute the volume of a tetrahedron constructed from three
-    vertices and the origin.
-
-    Parameters
-    ----------
-    v0
-        First vertex of the tetrahedron.
-    v1
-        Second vertex of the tetrahedron.
-    v2
-        Third vertex of the tetrahedron.
-    """
-    triple_product = rg.math.dot_3d(
-        a=v0,
-        b=rg.math.cross_3d(v1, v2),
-    )
-
-    return triple_product / 6
-
-
-@numba.njit(cache=True)
-def volume_grid(
-    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
-) -> np.ndarray:
-    """
-    Compute the volume of each cell in a logically-rectangular grid.
-
-    Parameters
-    ----------
-    grid
-        A 3D grid of cell vertices.
-    """
-
-    x, y, z = grid
-
-    num_i, num_j, num_k = x.shape
-
-    result = np.zeros((num_i - 1, num_j - 1, num_k - 1))
-
-    for axis in (0, 1, 2):
-        _volume_grid_sweep(
-            grid=grid,
-            out=result,
-            axis=axis,
-        )
-
-    return result
-
-
-@numba.njit(cache=True)
-def _volume_grid_sweep(
-    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
-    out: np.ndarray,
-    axis: int,
-) -> None:
-    """
-    Compute the volume contribution of this axis.
-
-    Parameters
-    ----------
-    grid
-        A 3D grid of cell vertices.
-    out
-        An output to array to store the result.
-    axis
-        The axis along which to iterate.
-    """
-
-    x, y, z = grid
-
-    x = _arrays.align_axis_right(x, axis)
-    y = _arrays.align_axis_right(y, axis)
-    z = _arrays.align_axis_right(z, axis)
-
-    out = _arrays.align_axis_right(out, axis)
-
-    num_i, num_j, num_k = x.shape
-
-    for i in numba.prange(num_i):
-
-        i = numba.types.int64(i)
-
-        for j in range(num_j):
-
-            for k in range(num_k - 1):
-
-                i_left = i - 1
-                i_right = i
-
-                j_lower = j - 1
-                j_upper = j
-
-                k1 = k
-                k2 = k + 1
-
-                v1 = x[i, j, k1], y[i, j, k1], z[i, j, k1]
-                v2 = x[i, j, k2], y[i, j, k2], z[i, j, k2]
-
-                volume_left = 0
-                volume_lower = 0
-
-                if i_left >= 0:
-                    v0_left = x[i_left, j, k], y[i_left, j, k], z[i_left, j, k]
-                    volume_left = volume_tetrahedron(v0_left, v1, v2)
-
-                if j_lower >= 0:
-                    v0_lower = x[i, j_lower, k], y[i, j_lower, k], z[i, j_lower, k]
-                    volume_lower = -volume_tetrahedron(v0_lower, v1, v2)
-
-                if i_left >= 0:
-                    if j_lower >= 0:
-                        out[i_left, j_lower, k] -= volume_left
-                    if j_upper < (num_j - 1):
-                        out[i_left, j_upper, k] += volume_left
-
-                if j_lower >= 0:
-                    if i_left >= 0:
-                        out[i_left, j_lower, k] -= volume_lower
-                    if i_right < (num_i - 1):
-                        out[i_right, j_lower, k] += volume_lower

From 551006390cc76ad7246e7240729c0ae17fb09291 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sun, 4 May 2025 07:36:46 -0600
Subject: [PATCH 10/26] renaming bug

---
 regridding/_weights/_weights_conservative_3d/_grids_test.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/regridding/_weights/_weights_conservative_3d/_grids_test.py b/regridding/_weights/_weights_conservative_3d/_grids_test.py
index d72c309..0366604 100644
--- a/regridding/_weights/_weights_conservative_3d/_grids_test.py
+++ b/regridding/_weights/_weights_conservative_3d/_grids_test.py
@@ -1,6 +1,6 @@
 import pytest
 import numpy as np
-from . import _volumes
+from . import _grids
 
 
 @pytest.mark.parametrize(
@@ -39,6 +39,6 @@ def test_volume_grid(
     grid: tuple[np.ndarray, np.ndarray, np.ndarray],
     result_expected: np.ndarray,
 ):
-    result = _volumes.volume_grid(grid)
+    result = _grids.grid_volume(grid)
     assert np.allclose(result, result_expected)
     assert result.shape == result_expected.shape

From 1de6fd2afd439270ec8fe7151277a3bd39f4a943 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sun, 4 May 2025 07:37:28 -0600
Subject: [PATCH 11/26] added `_line_point_closest_approach()`

---
 .../_weights_conservative_3d/_intercepts.py   | 69 +++++++++++++++++--
 1 file changed, 65 insertions(+), 4 deletions(-)

diff --git a/regridding/_weights/_weights_conservative_3d/_intercepts.py b/regridding/_weights/_weights_conservative_3d/_intercepts.py
index 6519f3b..acf1421 100644
--- a/regridding/_weights/_weights_conservative_3d/_intercepts.py
+++ b/regridding/_weights/_weights_conservative_3d/_intercepts.py
@@ -7,7 +7,6 @@
 import numba
 import regridding as rg
 from . import _arrays
-from . import _volumes
 from . import _grids
 
 
@@ -47,6 +46,69 @@ def empty(
     return intercepts
 
 
+def _line_point_closest_approach_parameter(
+    line: tuple[
+        tuple[float, float, float],
+        tuple[float, float, float],
+    ],
+    point: tuple[float, float, float],
+) -> float:
+    """
+    Compute the parameter describing the point of closest approach of line
+    segment to a point.
+
+    Based on a `StackExchange answer by John Hughes <https://math.stackexchange.com/a/2193733/502546>`_.
+
+    Parameters
+    ----------
+    line
+        A line segment in 3D.
+    point
+        A point in 3D
+    """
+
+    a, b = line
+    p = point
+
+    v = rg.math.difference_3d(b, a)
+    u = rg.math.difference_3d(a, p)
+
+    numerator = -rg.math.dot_3d(v, u)
+    denominator = rg.math.dot_3d(v, v)
+
+    return numerator / denominator
+
+
+def _line_point_closest_approach(
+    line: tuple[
+        tuple[float, float, float],
+        tuple[float, float, float],
+    ],
+    t: float,
+) -> tuple[float, float, float]:
+    """
+    Compute the point of closest approach of line segment to a point using
+    the parameter computed using :func:`_line_point_closest_approach_parameter`.
+
+    Based on a `StackExchange answer by John Hughes <https://math.stackexchange.com/a/2193733/502546>`_.
+
+    Parameters
+    ----------
+    line
+        A line segment in 3D.
+    t
+        The parameter of closest approach.
+    """
+    a, b = line
+
+    s = 1 - t
+
+    a2 = rg.math.multiply_3d(s, a)
+    b2 = rg.math.multiply_3d(t, b)
+
+    return rg.math.sum_3d(a2, b2)
+
+
 @numba.njit(cache=True)
 def sweep(
     intercepts: numba.typed.List,
@@ -136,8 +198,8 @@ def sweep(
                             z_output[i0_output],
                         )
 
-                        vol_input = _volumes.volume_tetrahedron(apex_input, p0, p1)
-                        vol_output = _volumes.volume_tetrahedron(apex_output, p0, p1)
+                        vol_input = rg.geometry.volume_tetrahedron(apex_input, p0, p1)
+                        vol_output = rg.geometry.volume_tetrahedron(apex_output, p0, p1)
 
                         volume_input_lower = volume_input[i0_input_lower]
                         volume_input_upper = volume_input[i0_output_upper]
@@ -167,7 +229,6 @@ def sweep(
                             index=i0_input_upper,
                             shape=shape_input,
                         )
-
                         i0_output_lower = _arrays.index_flat(
                             index=i0_output_lower,
                             shape=shape_output,

From 664e6ff087a9459e144d979fde1bfebad1f3149e Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sun, 4 May 2025 07:37:40 -0600
Subject: [PATCH 12/26] added `_line_point_closest_approach()`

---
 .../_intercepts_test.py                       | 25 +++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/regridding/_weights/_weights_conservative_3d/_intercepts_test.py b/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
index 34d9b62..cba30a9 100644
--- a/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
+++ b/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
@@ -4,6 +4,31 @@
 from . import _intercepts
 
 
+@pytest.mark.parametrize(
+    argnames="line, point, result_expected",
+    argvalues=[
+        (
+            ((-1, -1, 0), (1, 1, 0)),
+            (-1, 1, 0),
+            (0, 0, 0),
+        )
+    ],
+)
+def test_line_point_closest_approach(
+    line: tuple[
+        tuple[float, float, float],
+        tuple[float, float, float],
+    ],
+    point: tuple[float, float, float],
+    result_expected: tuple[float, float, float],
+):
+    t = _intercepts._line_point_closest_approach_parameter(line, point)
+
+    result = _intercepts._line_point_closest_approach(line, t)
+
+    assert result == result_expected
+
+
 @pytest.mark.parametrize(
     argnames="grid_input,grid_output,volume_input",
     argvalues=[

From 20a0b2dfa3ff7299c35a23c9093800ec3435ac56 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Mon, 5 May 2025 16:07:32 -0600
Subject: [PATCH 13/26] Added `_arrays.axes`

---
 regridding/_weights/_weights_conservative_3d/_arrays.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/regridding/_weights/_weights_conservative_3d/_arrays.py b/regridding/_weights/_weights_conservative_3d/_arrays.py
index 6962e6d..03c973a 100644
--- a/regridding/_weights/_weights_conservative_3d/_arrays.py
+++ b/regridding/_weights/_weights_conservative_3d/_arrays.py
@@ -9,6 +9,7 @@
     "axis_x",
     "axis_y",
     "axis_z",
+    "axes",
     "index_in_bounds",
     "index_flat",
     "index_3d",
@@ -20,6 +21,8 @@
 axis_y = 1
 axis_z = 2
 
+axes = (axis_x, axis_y, axis_z)
+
 
 @numba.njit(cache=True)
 def index_in_bounds(

From e418099103ca8f2bee43964f316dd5e6990c3180 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Mon, 5 May 2025 16:30:40 -0600
Subject: [PATCH 14/26] slight improvements

---
 .../_weights/_weights_conservative_3d/_grids.py     | 13 +++++--------
 1 file changed, 5 insertions(+), 8 deletions(-)

diff --git a/regridding/_weights/_weights_conservative_3d/_grids.py b/regridding/_weights/_weights_conservative_3d/_grids.py
index 596965c..b8ed8ac 100644
--- a/regridding/_weights/_weights_conservative_3d/_grids.py
+++ b/regridding/_weights/_weights_conservative_3d/_grids.py
@@ -110,24 +110,21 @@ def _grid_volume_sweep(
                 v1 = x[i, j, k1], y[i, j, k1], z[i, j, k1]
                 v2 = x[i, j, k2], y[i, j, k2], z[i, j, k2]
 
-                volume_left = 0
-                volume_lower = 0
-
                 if i_left >= 0:
+
                     v0_left = x[i_left, j, k], y[i_left, j, k], z[i_left, j, k]
                     volume_left = rg.geometry.volume_tetrahedron(v0_left, v1, v2)
 
-                if j_lower >= 0:
-                    v0_lower = x[i, j_lower, k], y[i, j_lower, k], z[i, j_lower, k]
-                    volume_lower = -rg.geometry.volume_tetrahedron(v0_lower, v1, v2)
-
-                if i_left >= 0:
                     if j_lower >= 0:
                         out[i_left, j_lower, k] -= volume_left
                     if j_upper < (num_j - 1):
                         out[i_left, j_upper, k] += volume_left
 
                 if j_lower >= 0:
+
+                    v0_lower = x[i, j_lower, k], y[i, j_lower, k], z[i, j_lower, k]
+                    volume_lower = -rg.geometry.volume_tetrahedron(v0_lower, v1, v2)
+
                     if i_left >= 0:
                         out[i_left, j_lower, k] -= volume_lower
                     if i_right < (num_i - 1):

From ba123029997c27c65a731684ff6f481d2effc56f Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Mon, 5 May 2025 16:49:31 -0600
Subject: [PATCH 15/26] Moving stuff and first draft of
 `_step_inside_static()`.

---
 .../_weights_conservative_3d/_grids.py        | 270 ++++++++++
 .../_weights_conservative_3d/_grids_test.py   |  38 ++
 .../_weights_conservative_3d.py               | 483 +++++++-----------
 .../_weights_conservative_3d_test.py          |  35 --
 4 files changed, 492 insertions(+), 334 deletions(-)

diff --git a/regridding/_weights/_weights_conservative_3d/_grids.py b/regridding/_weights/_weights_conservative_3d/_grids.py
index b8ed8ac..5c7d849 100644
--- a/regridding/_weights/_weights_conservative_3d/_grids.py
+++ b/regridding/_weights/_weights_conservative_3d/_grids.py
@@ -10,6 +10,8 @@
 
 __all__ = [
     "shape_centers",
+    "boundary",
+    "index_of_point_brute",
 ]
 
 
@@ -129,3 +131,271 @@ def _grid_volume_sweep(
                         out[i_left, j_lower, k] -= volume_lower
                     if i_right < (num_i - 1):
                         out[i_right, j_lower, k] += volume_lower
+
+
+@numba.njit(cache=True)
+def boundary(
+    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
+) -> tuple[
+    numba.typed.List[
+        tuple[int, int, int],
+    ],
+    numba.typed.List[
+        tuple[
+            tuple[float, float, float],
+            tuple[float, float, float],
+            tuple[float, float, float],
+        ],
+    ],
+]:
+    """
+    For a given grid, find the indices of the boundary cells
+    and the vertices of all the triangles on the boundary.
+
+    find the set of indices that expresses the boundary
+    as a sequence of triangles.
+
+    Parameters
+    ----------
+    grid
+        A logically-rectangular grid of cell vertices.
+    """
+
+    x, y, z = grid
+
+    shape_x, shape_y, shape_z = x.shape
+
+    indices = numba.typed.List()
+
+    triangles = numba.typed.List()
+
+    for i0 in range(shape_x - 1):
+        for j0 in range(shape_y - 1):
+
+            k0 = 0
+
+            i1 = i0 + 1
+            j1 = j0 + 1
+
+            i_cell = i0, j0, k0
+
+            i_000 = i0, j0, k0
+            i_010 = i0, j1, k0
+            i_110 = i1, j1, k0
+            i_100 = i1, j0, k0
+
+            v_000 = x[i_000], y[i_000], z[i_000]
+            v_010 = x[i_010], y[i_010], z[i_010]
+            v_110 = x[i_110], y[i_110], z[i_110]
+            v_100 = x[i_100], y[i_100], z[i_100]
+
+            indices.append(i_cell)
+            indices.append(i_cell)
+
+            triangles.append((v_000, v_010, v_110))
+            triangles.append((v_110, v_100, v_000))
+
+    for i0 in range(shape_x - 1):
+        for j0 in range(shape_y - 1):
+
+            k1 = shape_z
+
+            i1 = i0 + 1
+            j1 = j0 + 1
+
+            i_cell = i0, j0, k1 - 1
+
+            i_001 = i0, j0, k1
+            i_101 = i1, j0, k1
+            i_111 = i1, j1, k1
+            i_011 = i0, j1, k1
+
+            v_001 = x[i_001], y[i_001], z[i_001]
+            v_101 = x[i_101], y[i_101], z[i_101]
+            v_111 = x[i_111], y[i_111], z[i_111]
+            v_011 = x[i_011], y[i_011], z[i_011]
+
+            indices.append(i_cell)
+            indices.append(i_cell)
+
+            triangles.append((v_001, v_101, v_111))
+            triangles.append((v_111, v_011, v_001))
+
+    for j0 in range(shape_y - 1):
+        for k0 in range(shape_z - 1):
+
+            i0 = 0
+
+            j1 = j0 + 1
+            k1 = k0 + 1
+
+            i_cell = i0, j0, k0
+
+            i_000 = i0, j0, k0
+            i_001 = i0, j0, k1
+            i_011 = i0, j1, k1
+            i_010 = i0, j1, k0
+
+            v_000 = x[i_000], y[i_000], z[i_000]
+            v_001 = x[i_001], y[i_001], z[i_001]
+            v_011 = x[i_011], y[i_011], z[i_011]
+            v_010 = x[i_010], y[i_010], z[i_010]
+
+            indices.append(i_cell)
+            indices.append(i_cell)
+
+            triangles.append((v_000, v_001, v_011))
+            triangles.append((v_011, v_010, v_000))
+
+    for j0 in range(shape_y - 1):
+        for k0 in range(shape_z - 1):
+
+            i1 = shape_x
+
+            j1 = j0 + 1
+            k1 = k0 + 1
+
+            i_cell = i1 - 1, j0, k0
+
+            i_100 = i1, j0, k0
+            i_110 = i1, j1, k0
+            i_111 = i1, j1, k1
+            i_101 = i1, j0, k1
+
+            v_100 = x[i_100], y[i_100], z[i_100]
+            v_110 = x[i_110], y[i_110], z[i_110]
+            v_111 = x[i_111], y[i_111], z[i_111]
+            v_101 = x[i_101], y[i_101], z[i_101]
+
+            indices.append(i_cell)
+            indices.append(i_cell)
+
+            triangles.append((v_100, v_110, v_111))
+            triangles.append((v_111, v_101, v_100))
+
+    for i0 in range(shape_x - 1):
+        for k0 in range(shape_z - 1):
+
+            j0 = 0
+
+            i1 = i0 + 1
+            k1 = k0 + 1
+
+            i_cell = i0, j0, k0
+
+            i_000 = i0, j0, k0
+            i_100 = i1, j0, k0
+            i_101 = i1, j0, k1
+            i_001 = i0, j0, k1
+
+            v_000 = x[i_000], y[i_000], z[i_000]
+            v_100 = x[i_100], y[i_100], z[i_100]
+            v_101 = x[i_101], y[i_101], z[i_101]
+            v_001 = x[i_001], y[i_001], z[i_001]
+
+            indices.append(i_cell)
+            indices.append(i_cell)
+
+            triangles.append((v_000, v_100, v_101))
+            triangles.append((v_101, v_001, v_000))
+
+    for i0 in range(shape_x - 1):
+        for k0 in range(shape_z - 1):
+
+            j1 = shape_y
+
+            i1 = i0 + 1
+            k1 = k0 + 1
+
+            i_cell = i0, j1 - 1, k0
+
+            i_010 = i0, j1, k0
+            i_011 = i0, j1, k1
+            i_111 = i1, j1, k1
+            i_110 = i1, j1, k0
+
+            v_010 = x[i_010], y[i_010], z[i_010]
+            v_011 = x[i_011], y[i_011], z[i_011]
+            v_111 = x[i_111], y[i_111], z[i_111]
+            v_110 = x[i_110], y[i_110], z[i_110]
+
+            indices.append(i_cell)
+            indices.append(i_cell)
+
+            triangles.append((v_010, v_011, v_111))
+            triangles.append((v_111, v_110, v_010))
+
+    return indices, triangles
+
+
+@numba.njit(cache=True)
+def index_of_point_brute(
+    point: tuple[float, float, float],
+    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
+) -> tuple[int, int, int]:
+    """
+    Find the index of the cell in the grid which contains the given point.
+
+    This function uses brute force to search,
+    but this could be improved significantly by using the secant method
+    or possibly the bisection method.
+
+    Parameters
+    ----------
+    point
+        The query point.
+    grid
+        A logically-rectangular grid of cell vertices.
+    """
+
+    x, y, z = grid
+
+    shape_x, shape_y, shape_z = x.shape
+
+    for i in range(shape_x - 1):
+        for j in range(shape_y - 1):
+            for k in range(shape_z - 1):
+
+                i_000 = i + 0, j + 0, k + 0
+                i_001 = i + 0, j + 0, k + 1
+                i_010 = i + 0, j + 1, k + 0
+                i_011 = i + 0, j + 1, k + 1
+                i_100 = i + 1, j + 0, k + 0
+                i_101 = i + 1, j + 0, k + 1
+                i_110 = i + 1, j + 1, k + 0
+                i_111 = i + 1, j + 1, k + 1
+
+                indices = (
+                    (i_000, i_010, i_110),
+                    (i_110, i_100, i_000),
+                    (i_001, i_101, i_111),
+                    (i_111, i_011, i_001),
+                    (i_000, i_001, i_011),
+                    (i_011, i_010, i_000),
+                    (i_100, i_110, i_111),
+                    (i_111, i_101, i_100),
+                    (i_000, i_100, i_101),
+                    (i_101, i_001, i_000),
+                    (i_010, i_011, i_111),
+                    (i_111, i_110, i_010),
+                )
+
+                polyhedron = numba.typed.List()
+
+                for index in indices:
+
+                    t0, t1, t2 = index
+
+                    triangle = (
+                        (x[t0], y[t0], z[t0]),
+                        (x[t1], y[t1], z[t1]),
+                        (x[t2], y[t2], z[t2]),
+                    )
+
+                    polyhedron.append(triangle)
+
+                if rg.geometry.point_is_inside_polyhedron(
+                    point=point,
+                    polyhedron=polyhedron,
+                ):
+                    return i, j, k
diff --git a/regridding/_weights/_weights_conservative_3d/_grids_test.py b/regridding/_weights/_weights_conservative_3d/_grids_test.py
index 0366604..cede370 100644
--- a/regridding/_weights/_weights_conservative_3d/_grids_test.py
+++ b/regridding/_weights/_weights_conservative_3d/_grids_test.py
@@ -42,3 +42,41 @@ def test_volume_grid(
     result = _grids.grid_volume(grid)
     assert np.allclose(result, result_expected)
     assert result.shape == result_expected.shape
+
+
+@pytest.mark.parametrize(
+    argnames="point,grid,result_expected",
+    argvalues=[
+        (
+            (0.5, 0.5, 0.5),
+            np.meshgrid(
+                np.arange(3),
+                np.arange(4),
+                np.arange(5),
+                indexing="ij",
+            ),
+            (0, 0, 0),
+        ),
+        (
+            (-0.5, -0.5, 3.5),
+            np.meshgrid(
+                -np.arange(3),
+                -np.arange(4),
+                np.arange(5),
+                indexing="ij",
+            ),
+            (0, 0, 3),
+        ),
+    ],
+)
+def test_index_of_point_brute(
+    point: tuple[float, float, float],
+    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
+    result_expected: np.ndarray,
+):
+    result = _grids._index_of_point_brute(
+        point=point,
+        grid=grid,
+    )
+
+    assert result == result_expected
diff --git a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
index 2a4f2c5..80bec54 100644
--- a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
+++ b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
@@ -44,7 +44,7 @@ def weights_conservative_3d(
 
     weights = numba.typed.List()
     for x in range(0):
-        weights.append((0.0, 0.0, 0.0))
+        weights.append((0, 0, 0.0))
 
     volume_input = _volumes.volume_grid(grid_input)
 
@@ -114,7 +114,7 @@ def _sweep_grid(
         (grid_static_x.max(), grid_static_y.max(), grid_static_z.max()),
     )
 
-    indices_boundary, triangles_boundary = _indices_and_triangles_of_boundary(grid_static)
+    indices_boundary, triangles_boundary = _grids.boundary(grid_static)
 
     axes_parallel = (
         axis_x,
@@ -134,8 +134,8 @@ def _sweep_grid(
             grid_output=grid_output,
             volume_input=volume_input,
             bbox_boundary=bbox_boundary,
-            index_boundary=index_boundary,
-            coord_boundary=coord_boundary,
+            indices_boundary=indices_boundary,
+            triangles_boundary=triangles_boundary,
             weights=weights,
             intercepts=intercepts,
             sweep_input=sweep_input,
@@ -209,7 +209,7 @@ def _sweep_along_axis(
     for i in range(shape_sweep_x):
         weight_i = numba.typed.List()
         for _ in range(0):
-            weight_i.append((0., 0., 0.))
+            weight_i.append((0, 0, 0.))
         weight.append(weight_i)
 
     for i in numba.prange(shape_sweep_x):
@@ -240,7 +240,7 @@ def _sweep_along_axis(
                     point=point_1,
                     polyhedron=coord_boundary,
                 ):
-                    index_static = _index_of_point_brute(
+                    index_static = _grids.index_of_point_brute(
                         point=point_1,
                         grid=grid_static,
                     )
@@ -375,40 +375,78 @@ def _step_inside_static(
         tuple[float, float, float],
         tuple[float, float, float],
     ],
-    index_static: tuple[int, int, int],
     index_sweep: tuple[int, int, int],
+    index_static: tuple[int, int, int],
+    grid_sweep: tuple[np.ndarray, np.ndarray, np.ndarray],
     grid_static: tuple[np.ndarray, np.ndarray, np.ndarray],
+    shape_cells_input: tuple[int, int, int],
+    shape_cells_output: tuple[int, int, int],
     weights: numba.typed.List[tuple[int, int, float]],
-) -> None:
+    intercepts: numba.typed.List,
+    sweep_input: bool,
+    axis: int,
+) -> tuple[
+    tuple[
+        tuple[float, float, float],
+        tuple[float, float, float],
+    ],
+    tuple[int, int, int],
+    tuple[int, int, int],
+]:
     """
     Check if the current line segment crosses any face of the current cell
     in the static grid.
     If it does, compute the intersection point and indices of the new
     cell in the static grid.
+    Otherwise, increment the index of the sweep grid.
 
     Parameters
     ----------
-    point_1
-        The first point of the line segment.
-    point_2
-        The second point of the line segment.
-    index_static
-        The current index in the static grid.
+    line
+        The current line segment of the sweep grid.
     index_sweep
         The current index in the sweep grid.
+        This is assumed to be a valid, positive index.
+    index_static
+        The current index in the static grid.
+        This is assumed to be a valid, positive index.
+    grid_sweep
+        The vertices of the sweep grid.
+        The last axis of this grid must be the sweep axis.
     grid_static
         Coordinates of the static grid.
+    shape_cells_input
+        The number of cells along each axis of the input grid.
+        This could be calculated from the other arguments,
+        but we don't to save computation time.
+    shape_cells_output
+        The number of cells along each axis of the output grid.
+        Like `shape_cells_input`,
+        this could be calculated from the other arguments,
+        but we don't to save computation time.
     weights
         The current list of weights.
         New weights will be appended to this list.
+    intercepts
+        A sorted list of intercepts to be traversed later.
+        As new intercepts are found, they are inserted into the list.
+    sweep_input
+        A boolean flag indicating whether to iterate along `grid_input`
+        or `grid_output`.
+        If :obj:`True`, this function sweeps along `grid_input`.
+    axis
+        The logical axis of the sweep grid to iterate along.
     """
 
-    point_1, point_2 = line
+    p1, p2 = line
 
     i, j, k = index_sweep
     l, m, n = index_static
 
-    x, y, z = grid_static
+    x_sweep, y_sweep, z_sweep = grid_sweep
+    x_static, y_static, z_static = grid_static
+
+    shape_sweep = x_sweep.shape
 
     i_000 = l + 0, m + 0, n + 0
     i_001 = l + 0, m + 0, n + 1
@@ -434,6 +472,21 @@ def _step_inside_static(
         (i_111, i_110, i_010),
     )
 
+    axis_index = (
+        2,
+        2,
+        2,
+        2,
+        0,
+        0,
+        0,
+        0,
+        1,
+        1,
+        1,
+        1,
+    )
+
     normals_index = (
         (0, 0, -1),
         (0, 0, -1),
@@ -453,9 +506,9 @@ def _step_inside_static(
 
         i0, i1, i2 = indices_triangles[t]
 
-        v0 = (x[i0], y[i0], z[i0])
-        v1 = (x[i1], y[i1], z[i1])
-        v2 = (x[i2], y[i2], z[i2])
+        v0 = (x_static[i0], y_static[i0], z_static[i0])
+        v1 = (x_static[i1], y_static[i1], z_static[i1])
+        v2 = (x_static[i2], y_static[i2], z_static[i2])
 
         triangle = (v0, v1, v2)
 
@@ -471,30 +524,133 @@ def _step_inside_static(
                 tuv=tuv,
             )
 
-            line = (point_1, intercept)
+            p2 = intercept
 
             normal_index = normals_index[t]
 
-            index_static = rg.math.sum_3d(index_static, normal_index)
+            index_static_new = rg.math.sum_3d(index_static, normal_index)
 
+            axis_static = axis_index[t]
 
+            if index_static_new[axis_static] > index_static[axis_static]:
+                index_static_intercept = index_static_new
+            else:
+                index_static_intercept = index_static
 
-        else:
-            index_sweep = i, j, k + 1
+            axes_orthogonal = [ax for ax in _arrays.axes if ax != axis]
 
+            for axis_orthogonal in axes_orthogonal:
 
+                if sweep_input:
+                    axis_input = axis_orthogonal
+                    axis_output = axis_static
+                else:
+                    axis_input = axis_static
+                    axis_output = axis_orthogonal
 
-        # index_left
-        #
-        # volume_sweep =
-        #
-        # i_normal, j_normal, k_normal = normals_index[t]
+                index_input, index_output = _index_input_output(
+                    index_sweep=index_sweep,
+                    index_static=index_static_intercept,
+                    sweep_input=sweep_input,
+                    axis=axis,
+                )
 
+                i_input = index_input[axis_input]
+                i_output = index_output[axis_output]
 
+                _intercepts.insert_intercept(
+                    intercepts=intercepts[axis_input][axis_output][i_input][i_output],
+                    intercept_new=(
+                        _arrays.index_flat(index_input, shape_cells_input),
+                        _arrays.index_flat(index_output, shape_cells_output),
+                        intercept,
+                    ),
+                )
 
+                break
 
+    else:
+        index_sweep = i, j, k + 1
 
+    i_left = i - 1
+    i_right = i
 
+    j_lower = j - 1
+    j_upper = j
+
+    index_sweep_left = i_left, j, k
+    index_sweep_lower = i, j_lower, k
+
+    p0_left = (
+        x_sweep[index_sweep_left],
+        y_sweep[index_sweep_left],
+        z_sweep[index_sweep_left],
+    )
+    p0_lower = (
+        x_sweep[index_sweep_lower],
+        y_sweep[index_sweep_lower],
+        z_sweep[index_sweep_lower],
+    )
+
+    volume_left = rg.geometry.volume_tetrahedron(p0_left, p1, p2)
+    volume_lower = -rg.geometry.volume_tetrahedron(p0_lower, p1, p1)
+
+    if i_left >= 0:
+        if j_lower >= 0:
+
+            index_sweep_00 = i_left, j_lower, k
+
+            index_input_00, index_output_00 = _index_input_output(
+                index_sweep=index_sweep_00,
+                index_static=index_static,
+                sweep_input=sweep_input,
+                axis=axis,
+            )
+
+            index_input_00 = _arrays.index_flat(index_input_00, shape_cells_input)
+            index_output_00 = _arrays.index_flat(index_output_00, shape_cells_output)
+
+            volume_lower_left = -volume_lower - volume_left
+            weights.append((index_input_00, index_output_00, volume_lower_left))
+
+        if j_upper < (shape_sweep[1] - 1):
+
+            index_sweep_01 = i_left, j_upper, k
+
+            index_input_01, index_output_01 = _index_input_output(
+                index_sweep=index_sweep_01,
+                index_static=index_static,
+                sweep_input=sweep_input,
+                axis=axis,
+            )
+
+            index_input_01 = _arrays.index_flat(index_input_01, shape_cells_input)
+            index_output_01 = _arrays.index_flat(index_output_01, shape_cells_output)
+
+            volume_upper_left = volume_left
+            weights.append((index_input_01, index_output_01, volume_upper_left))
+
+    if i_right < (shape_sweep[0] - 1):
+        if j_lower >= 0:
+
+            index_sweep_10 = i_right, j_lower, k
+
+            index_input_10, index_output_10 = _index_input_output(
+                index_sweep=index_sweep_10,
+                index_static=index_static,
+                sweep_input=sweep_input,
+                axis=axis,
+            )
+
+            index_input_10 = _arrays.index_flat(index_input_10, shape_cells_input)
+            index_output_10 = _arrays.index_flat(index_output_10, shape_cells_output)
+
+            volume_lower_right = volume_lower
+            weights.append((index_input_10, index_output_10, volume_lower_right))
+
+    line = p1, p2
+
+    return line, index_sweep, index_static
 
 
 @numba.njit(cache=True)
@@ -566,274 +722,3 @@ def _index_input_output(
 
     return index_input, index_output
 
-
-
-
-
-
-@numba.njit(cache=True)
-def _indices_and_triangles_of_boundary(
-    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
-) -> tuple[
-    numba.typed.List[
-        tuple[int, int, int],
-    ],
-    numba.typed.List[
-        tuple[
-            tuple[float, float, float],
-            tuple[float, float, float],
-            tuple[float, float, float],
-        ],
-    ],
-]:
-    """
-    For a given grid, find the indices of the boundary cells
-    and the vertices of all the triangles on the boundary.
-
-    find the set of indices that expresses the boundary
-    as a sequence of triangles.
-
-    Parameters
-    ----------
-    grid
-        A logically-rectangular grid of cell vertices.
-    """
-
-    x, y, z = grid
-
-    shape_x, shape_y, shape_z = x.shape
-
-    indices = numba.typed.List()
-
-    triangles = numba.typed.List()
-
-    for i0 in range(shape_x - 1):
-        for j0 in range(shape_y - 1):
-
-            k0 = 0
-
-            i1 = i0 + 1
-            j1 = j0 + 1
-
-            i_cell = i0, j0, k0
-
-            i_000 = i0, j0, k0
-            i_010 = i0, j1, k0
-            i_110 = i1, j1, k0
-            i_100 = i1, j0, k0
-
-            v_000 = x[i_000], y[i_000], z[i_000]
-            v_010 = x[i_010], y[i_010], z[i_010]
-            v_110 = x[i_110], y[i_110], z[i_110]
-            v_100 = x[i_100], y[i_100], z[i_100]
-
-            indices.append(i_cell)
-            indices.append(i_cell)
-
-            triangles.append((v_000, v_010, v_110))
-            triangles.append((v_110, v_100, v_000))
-
-    for i0 in range(shape_x - 1):
-        for j0 in range(shape_y - 1):
-
-            k1 = shape_z
-
-            i1 = i0 + 1
-            j1 = j0 + 1
-
-            i_cell = i0, j0, k1 - 1
-
-            i_001 = i0, j0, k1
-            i_101 = i1, j0, k1
-            i_111 = i1, j1, k1
-            i_011 = i0, j1, k1
-
-            v_001 = x[i_001], y[i_001], z[i_001]
-            v_101 = x[i_101], y[i_101], z[i_101]
-            v_111 = x[i_111], y[i_111], z[i_111]
-            v_011 = x[i_011], y[i_011], z[i_011]
-
-            indices.append(i_cell)
-            indices.append(i_cell)
-
-            triangles.append((v_001, v_101, v_111))
-            triangles.append((v_111, v_011, v_001))
-
-    for j0 in range(shape_y - 1):
-        for k0 in range(shape_z - 1):
-
-            i0 = 0
-
-            j1 = j0 + 1
-            k1 = k0 + 1
-
-            i_cell = i0, j0, k0
-
-            i_000 = i0, j0, k0
-            i_001 = i0, j0, k1
-            i_011 = i0, j1, k1
-            i_010 = i0, j1, k0
-
-            v_000 = x[i_000], y[i_000], z[i_000]
-            v_001 = x[i_001], y[i_001], z[i_001]
-            v_011 = x[i_011], y[i_011], z[i_011]
-            v_010 = x[i_010], y[i_010], z[i_010]
-
-            indices.append(i_cell)
-            indices.append(i_cell)
-
-            triangles.append((v_000, v_001, v_011))
-            triangles.append((v_011, v_010, v_000))
-
-    for j0 in range(shape_y - 1):
-        for k0 in range(shape_z - 1):
-
-            i1 = shape_x
-
-            j1 = j0 + 1
-            k1 = k0 + 1
-
-            i_cell = i1 - 1, j0, k0
-
-            i_100 = i1, j0, k0
-            i_110 = i1, j1, k0
-            i_111 = i1, j1, k1
-            i_101 = i1, j0, k1
-
-            v_100 = x[i_100], y[i_100], z[i_100]
-            v_110 = x[i_110], y[i_110], z[i_110]
-            v_111 = x[i_111], y[i_111], z[i_111]
-            v_101 = x[i_101], y[i_101], z[i_101]
-
-            indices.append(i_cell)
-            indices.append(i_cell)
-
-            triangles.append((v_100, v_110, v_111))
-            triangles.append((v_111, v_101, v_100))
-
-    for i0 in range(shape_x - 1):
-        for k0 in range(shape_z - 1):
-
-            j0 = 0
-
-            i1 = i0 + 1
-            k1 = k0 + 1
-
-            i_cell = i0, j0, k0
-
-            i_000 = i0, j0, k0
-            i_100 = i1, j0, k0
-            i_101 = i1, j0, k1
-            i_001 = i0, j0, k1
-
-            v_000 = x[i_000], y[i_000], z[i_000]
-            v_100 = x[i_100], y[i_100], z[i_100]
-            v_101 = x[i_101], y[i_101], z[i_101]
-            v_001 = x[i_001], y[i_001], z[i_001]
-
-            indices.append(i_cell)
-            indices.append(i_cell)
-
-            triangles.append((v_000, v_100, v_101))
-            triangles.append((v_101, v_001, v_000))
-
-    for i0 in range(shape_x - 1):
-        for k0 in range(shape_z - 1):
-
-            j1 = shape_y
-
-            i1 = i0 + 1
-            k1 = k0 + 1
-
-            i_cell = i0, j1 - 1, k0
-
-            i_010 = i0, j1, k0
-            i_011 = i0, j1, k1
-            i_111 = i1, j1, k1
-            i_110 = i1, j1, k0
-
-            v_010 = x[i_010], y[i_010], z[i_010]
-            v_011 = x[i_011], y[i_011], z[i_011]
-            v_111 = x[i_111], y[i_111], z[i_111]
-            v_110 = x[i_110], y[i_110], z[i_110]
-
-            indices.append(i_cell)
-            indices.append(i_cell)
-
-            triangles.append((v_010, v_011, v_111))
-            triangles.append((v_111, v_110, v_010))
-
-    return indices, triangles
-
-
-@numba.njit(cache=True)
-def _index_of_point_brute(
-    point: tuple[float, float, float],
-    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
-) -> tuple[int, int, int]:
-    """
-    Find the index of the cell in the grid which contains the given point.
-
-    This function uses brute force to search,
-    but this could be improved significantly by using the secant method
-    or possibly the bisection method.
-
-    Parameters
-    ----------
-    point
-        The query point.
-    grid
-        A logically-rectangular grid of cell vertices.
-    """
-
-    x, y, z = grid
-
-    shape_x, shape_y, shape_z = x.shape
-
-    for i in range(shape_x - 1):
-        for j in range(shape_y - 1):
-            for k in range(shape_z - 1):
-
-                i_000 = i + 0, j + 0, k + 0
-                i_001 = i + 0, j + 0, k + 1
-                i_010 = i + 0, j + 1, k + 0
-                i_011 = i + 0, j + 1, k + 1
-                i_100 = i + 1, j + 0, k + 0
-                i_101 = i + 1, j + 0, k + 1
-                i_110 = i + 1, j + 1, k + 0
-                i_111 = i + 1, j + 1, k + 1
-
-                indices = (
-                    (i_000, i_010, i_110),
-                    (i_110, i_100, i_000),
-                    (i_001, i_101, i_111),
-                    (i_111, i_011, i_001),
-                    (i_000, i_001, i_011),
-                    (i_011, i_010, i_000),
-                    (i_100, i_110, i_111),
-                    (i_111, i_101, i_100),
-                    (i_000, i_100, i_101),
-                    (i_101, i_001, i_000),
-                    (i_010, i_011, i_111),
-                    (i_111, i_110, i_010),
-                )
-
-                polyhedron = numba.typed.List()
-
-                for index in indices:
-
-                    t0, t1, t2 = index
-
-                    triangle = (
-                        (x[t0], y[t0], z[t0]),
-                        (x[t1], y[t1], z[t1]),
-                        (x[t2], y[t2], z[t2]),
-                    )
-
-                    polyhedron.append(triangle)
-
-                if rg.geometry.point_is_inside_polyhedron(
-                    point=point,
-                    polyhedron=polyhedron,
-                ):
-                    return i, j, k
diff --git a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d_test.py b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d_test.py
index d665855..551fd72 100644
--- a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d_test.py
+++ b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d_test.py
@@ -39,40 +39,5 @@ def test_weights_conservative_3d(
 
 
 
-@pytest.mark.parametrize(
-    argnames="point,grid,result_expected",
-    argvalues=[
-        (
-            (0.5, 0.5, 0.5),
-            np.meshgrid(
-                np.arange(3),
-                np.arange(4),
-                np.arange(5),
-                indexing="ij",
-            ),
-            (0, 0, 0),
-        ),
-        (
-            (-0.5, -0.5, 3.5),
-            np.meshgrid(
-                -np.arange(3),
-                -np.arange(4),
-                np.arange(5),
-                indexing="ij",
-            ),
-            (0, 0, 3),
-        ),
-    ],
-)
-def test_index_of_point_brute(
-    point: tuple[float, float, float],
-    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
-    result_expected: np.ndarray,
-):
-    result = _weights_conservative_3d._index_of_point_brute(
-        point=point,
-        grid=grid,
-    )
 
-    assert result == result_expected
 

From f60a5343ec1088e2b1741b4c12ef36b1c857db7a Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Wed, 7 May 2025 06:26:28 -0600
Subject: [PATCH 16/26] Lots of updates and refactoring

---
 .../_weights_conservative_3d/_grids.py        | 153 +++--
 .../_weights_conservative_3d/_intercepts.py   |  17 +
 .../_weights_conservative_3d.py               | 535 +++++++++++-------
 3 files changed, 453 insertions(+), 252 deletions(-)

diff --git a/regridding/_weights/_weights_conservative_3d/_grids.py b/regridding/_weights/_weights_conservative_3d/_grids.py
index 5c7d849..44e6802 100644
--- a/regridding/_weights/_weights_conservative_3d/_grids.py
+++ b/regridding/_weights/_weights_conservative_3d/_grids.py
@@ -10,7 +10,10 @@
 
 __all__ = [
     "shape_centers",
-    "boundary",
+    "cell_axes",
+    "cell_normals",
+    "cell_boundary",
+    "grid_boundary",
     "index_of_point_brute",
 ]
 
@@ -133,8 +136,111 @@ def _grid_volume_sweep(
                         out[i_right, j_lower, k] += volume_lower
 
 
+cell_axes = (
+    2,
+    2,
+    2,
+    2,
+    0,
+    0,
+    0,
+    0,
+    1,
+    1,
+    1,
+    1,
+)
+"""
+The index of the axis normal to each face in :func:`cell_boundary`.
+"""
+
+cell_normals = (
+    (0, 0, -1),
+    (0, 0, -1),
+    (0, 0, +1),
+    (0, 0, +1),
+    (-1, 0, 0),
+    (-1, 0, 0),
+    (+1, 0, 0),
+    (+1, 0, 0),
+    (0, -1, 0),
+    (0, -1, 0),
+    (0, +1, 0),
+    (0, +1, 0),
+)
+"""
+Vectors normal to each face in :func:`cell_boundary`.
+"""
+
+
+@numba.njit(cache=True)
+def cell_boundary(
+    index: tuple[int, int, int],
+    grid: tuple[np.ndarray, np.ndarray, np.ndarray],
+) -> numba.typed.List[
+    tuple[
+        tuple[float, float, float],
+        tuple[float, float, float],
+        tuple[float, float, float],
+    ],
+]:
+    """
+    For a given cell in a grid of cell vertices,
+    compute the 12 triangles composing the boundary.
+
+    Parameters
+    ----------
+    index
+        The index of the cell to compute the boundary of.
+    grid
+        A grid of cell vertices.
+
+    """
+    i, j, k = index
+
+    x, y, z = grid
+
+    i_000 = i + 0, j + 0, k + 0
+    i_001 = i + 0, j + 0, k + 1
+    i_010 = i + 0, j + 1, k + 0
+    i_011 = i + 0, j + 1, k + 1
+    i_100 = i + 1, j + 0, k + 0
+    i_101 = i + 1, j + 0, k + 1
+    i_110 = i + 1, j + 1, k + 0
+    i_111 = i + 1, j + 1, k + 1
+
+    indices_triangles = (
+        (i_000, i_010, i_110),
+        (i_110, i_100, i_000),
+        (i_001, i_101, i_111),
+        (i_111, i_011, i_001),
+        (i_000, i_001, i_011),
+        (i_011, i_010, i_000),
+        (i_100, i_110, i_111),
+        (i_111, i_101, i_100),
+        (i_000, i_100, i_101),
+        (i_101, i_001, i_000),
+        (i_010, i_011, i_111),
+        (i_111, i_110, i_010),
+    )
+
+    polyhedron = numba.typed.List()
+
+    for t in range(len(indices_triangles)):
+
+        i0, i1, i2 = indices_triangles[t]
+
+        v0 = x[i0], y[i0], z[i0]
+        v1 = x[i1], y[i1], z[i1]
+        v2 = x[i2], y[i2], z[i2]
+
+        polyhedron.append((v0, v1, v2))
+
+    return polyhedron
+
+
 @numba.njit(cache=True)
-def boundary(
+def grid_boundary(
     grid: tuple[np.ndarray, np.ndarray, np.ndarray],
 ) -> tuple[
     numba.typed.List[
@@ -356,46 +462,15 @@ def index_of_point_brute(
         for j in range(shape_y - 1):
             for k in range(shape_z - 1):
 
-                i_000 = i + 0, j + 0, k + 0
-                i_001 = i + 0, j + 0, k + 1
-                i_010 = i + 0, j + 1, k + 0
-                i_011 = i + 0, j + 1, k + 1
-                i_100 = i + 1, j + 0, k + 0
-                i_101 = i + 1, j + 0, k + 1
-                i_110 = i + 1, j + 1, k + 0
-                i_111 = i + 1, j + 1, k + 1
-
-                indices = (
-                    (i_000, i_010, i_110),
-                    (i_110, i_100, i_000),
-                    (i_001, i_101, i_111),
-                    (i_111, i_011, i_001),
-                    (i_000, i_001, i_011),
-                    (i_011, i_010, i_000),
-                    (i_100, i_110, i_111),
-                    (i_111, i_101, i_100),
-                    (i_000, i_100, i_101),
-                    (i_101, i_001, i_000),
-                    (i_010, i_011, i_111),
-                    (i_111, i_110, i_010),
-                )
-
-                polyhedron = numba.typed.List()
+                index = i, j, k
 
-                for index in indices:
-
-                    t0, t1, t2 = index
-
-                    triangle = (
-                        (x[t0], y[t0], z[t0]),
-                        (x[t1], y[t1], z[t1]),
-                        (x[t2], y[t2], z[t2]),
-                    )
-
-                    polyhedron.append(triangle)
+                polyhedron = cell_boundary(
+                    index=index,
+                    grid=grid,
+                )
 
                 if rg.geometry.point_is_inside_polyhedron(
                     point=point,
                     polyhedron=polyhedron,
                 ):
-                    return i, j, k
+                    return index
diff --git a/regridding/_weights/_weights_conservative_3d/_intercepts.py b/regridding/_weights/_weights_conservative_3d/_intercepts.py
index acf1421..8bc6ac3 100644
--- a/regridding/_weights/_weights_conservative_3d/_intercepts.py
+++ b/regridding/_weights/_weights_conservative_3d/_intercepts.py
@@ -46,6 +46,23 @@ def empty(
     return intercepts
 
 
+def insert_intercept(
+    intercepts: numba.typed.List[tuple[int, int, tuple[float, float, float]]],
+    intercept_new: tuple[int, int, tuple[float, float, float]],
+) -> None:
+    """
+    Insert a new intercept into the current list of intercepts at the correct
+    point so that the list maintains sorted order.
+
+    Parameters
+    ----------
+    intercepts
+        The current list of intercepts.
+    intercept_new
+        A new intercept to insert into the list.
+    """
+
+
 def _line_point_closest_approach_parameter(
     line: tuple[
         tuple[float, float, float],
diff --git a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
index 80bec54..c20360d 100644
--- a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
+++ b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
@@ -60,15 +60,15 @@ def weights_conservative_3d(
             sweep_input=sweep_input,
         )
 
-    # _sweep_intercepts(
-    #     grid_input=grid_input,
-    #     grid_output=grid_output,
-    #     volume_input=volume_input,
-    #     intercepts=intercepts,
-    #     weights=weights,
-    # )
+    _intercepts.sweep(
+        intercepts=intercepts,
+        weights=weights,
+        grid_input=grid_input,
+        grid_output=grid_output,
+        volume_input=volume_input,
+    )
 
-    return weights, intercepts
+    return weights
 
 
 @numba.njit(cache=True)
@@ -114,21 +114,18 @@ def _sweep_grid(
         (grid_static_x.max(), grid_static_y.max(), grid_static_z.max()),
     )
 
-    indices_boundary, triangles_boundary = _grids.boundary(grid_static)
-
-    axes_parallel = (
-        axis_x,
-        axis_y,
-        axis_z,
-    )
+    indices_boundary, triangles_boundary = _grids.grid_boundary(grid_static)
 
-    axes_diagonal = (
+    axes = (
+        (axis_x, ),
+        (axis_y, ),
+        (axis_z, ),
         (axis_x, axis_y),
         (axis_y, axis_z),
         (axis_z, axis_x),
     )
 
-    for axis in axes_parallel:
+    for axis in axes:
         _sweep_along_axis(
             grid_input=grid_input,
             grid_output=grid_output,
@@ -142,9 +139,6 @@ def _sweep_grid(
             axes=axis,
         )
 
-    for axis in axes_diagonal:
-        _sweep_along_diagonal()
-
 
 @numba.njit(cache=True, parallel=True)
 def _sweep_along_axis(
@@ -153,11 +147,11 @@ def _sweep_along_axis(
     volume_input: np.ndarray,
     bbox_boundary: tuple[tuple[float, float, float], tuple[float, float, float]],
     indices_boundary: tuple[np.ndarray, np.ndarray, np.ndarray],
-    coord_boundary: tuple[np.ndarray, np.ndarray, np.ndarray],
+    triangles_boundary: tuple[np.ndarray, np.ndarray, np.ndarray],
     weights: numba.typed.List[tuple[int, int, float]],
     intercepts: numba.typed.List,
     sweep_input: bool,
-    axis: int,
+    axis_sweep: tuple[int] | tuple[int, int],
 ) -> None:
     """
     Sweep along one logical axis of either `grid_input` or `grid_output`.
@@ -173,9 +167,9 @@ def _sweep_along_axis(
         The area of the input grid which will used to compute the weight.
     bbox_boundary
         Two points defining the bounding box of the static grid.
-    coord_boundary
+    triangles_boundary
         A sequence of triangles defining the outer surface of the static grid.
-    index_boundary
+    indices_boundary
         The index of the static grid corresponding to each triangle in
         `coord_boundary`.
     weights
@@ -188,7 +182,7 @@ def _sweep_along_axis(
         A boolean flag indicating whether to iterate along `grid_input`
         or `grid_output`.
         If :obj:`True`, this function sweeps along `grid_input`.
-    axis
+    axis_sweep
         The logical axis of the sweep grid to iterate along.
     """
 
@@ -196,13 +190,20 @@ def _sweep_along_axis(
         grid_input=grid_input,
         grid_output=grid_output,
         sweep_input=sweep_input,
-        axis_x=axis,
+        axis_sweep=axis_sweep,
     )
 
     x_sweep, y_sweep, z_sweep = grid_sweep
+    x_static, y_static, z_static = grid_static
+
+    shape_sweep = x_sweep.shape
+    shape_static = x_static.shape
 
-    shape_sweep_x, shape_sweep_y, shape_sweep_z = grid_sweep.shape
-    shape_static_x, shape_static_y, shape_static_z = grid_static.shape
+    shape_cells_sweep = _grids.shape_centers(shape_sweep)
+    shape_cells_static = _grids.shape_centers(shape_static)
+
+    shape_sweep_x, shape_sweep_y, shape_sweep_z = shape_sweep
+    shape_static_x, shape_static_y, shape_static_z = shape_static
 
     weight = numba.typed.List()
 
@@ -212,13 +213,30 @@ def _sweep_along_axis(
             weight_i.append((0, 0, 0.))
         weight.append(weight_i)
 
+    if len(axis_sweep) == 2:
+        d_start = -(shape_sweep_z - 2)
+        d_end = shape_sweep_y - 1
+        step = (0, 1, 1)
+    else:
+        d_start = 0
+        d_end = shape_sweep_y
+        step = (0, 0, 1)
+
+    j_end = d_end
+    k_end = shape_cells_sweep[2] - 1
+
     for i in numba.prange(shape_sweep_x):
 
         i = numba.types.int64(i)
 
-        for j in range(shape_sweep_y):
+        for d in range(d_start, d_end):
 
-            k = 0
+            if d < 0:
+                j = 0
+                k = -d
+            else:
+                j = d
+                k = 0
 
             index_sweep = i, j, k
 
@@ -238,7 +256,7 @@ def _sweep_along_axis(
             ):
                 if rg.geometry.point_is_inside_polyhedron(
                     point=point_1,
-                    polyhedron=coord_boundary,
+                    polyhedron=triangles_boundary,
                 ):
                     index_static = _grids.index_of_point_brute(
                         point=point_1,
@@ -246,12 +264,14 @@ def _sweep_along_axis(
                     )
                     sweep_is_outside_static = False
 
-            while k < (shape_sweep_z - 1):
+            while (j < j_end) and (k < k_end):
+
+                index_sweep_new = rg.math.sum_3d(index_sweep, step)
 
                 point_2 = (
-                    x_sweep[i, j, k + 1],
-                    y_sweep[i, j, k + 1],
-                    z_sweep[i, j, k + 1],
+                    x_sweep[index_sweep_new],
+                    y_sweep[index_sweep_new],
+                    z_sweep[index_sweep_new],
                 )
 
                 if sweep_is_outside_static:
@@ -260,8 +280,8 @@ def _sweep_along_axis(
                         point_1=point_1,
                         point_2=point_2,
                         bbox_boundary=bbox_boundary,
-                        index_boundary=index_boundary,
-                        coord_boundary=coord_boundary,
+                        indices_boundary=indices_boundary,
+                        triangles_boundary=triangles_boundary,
                         intercepts=intercepts,
                     )
 
@@ -269,20 +289,19 @@ def _sweep_along_axis(
 
                 else:
 
-                    point_2, k, index_static = _step_inside_static(
-                        point_1=point_1,
-                        point_2=point_2,
+                    line, index_sweep, index_static = _step_inside_static(
+                        line=line,
+                        index_sweep=index_sweep,
                         index_static=index_static,
                         grid_static=grid_static,
                         weights=weights,
                         intercepts=intercepts,
                     )
 
-                    if not 0 <= index_static[axis_x] < (shape_static_x - 1):
-                        break
-                    if not 0 <= index_static[axis_y] < (shape_static_y - 1):
-                        break
-                    if not 0 <= index_static[axis_z] < (shape_static_z - 1):
+                    if not _arrays.index_in_bounds(
+                        index=index_static,
+                        shape=shape_cells_static,
+                    ):
                         break
 
     for i in range(shape_sweep_x):
@@ -293,80 +312,137 @@ def _sweep_along_axis(
 
 @numba.njit(cache=True)
 def _step_outside_static(
-    point_1: tuple[float, float, float],
-    point_2: tuple[float, float, float],
+    line: tuple[
+        tuple[float, float, float],
+        tuple[float, float, float],
+    ],
     index_sweep: tuple[int, int, int],
+    index_static: tuple[int, int, int],
+    step_index: tuple[int, int, int],
+    grid_sweep: tuple[np.ndarray, np.ndarray, np.ndarray],
     grid_static: tuple[np.ndarray, np.ndarray, np.ndarray],
-    bbox_boundary: tuple[tuple[float, float, float], tuple[float, float, float]],
-    indices_boundary: numba.typed.List[
-        tuple[int, int, int],
-    ],
-    triangles_boundary: numba.typed.List[
-        tuple[
-            tuple[float, float, float],
-            tuple[float, float, float],
-            tuple[float, float, float],
-        ],
-    ],
+    bbox_static: tuple[tuple[float, float, float], tuple[float, float, float]],
+    shape_cells_input: tuple[int, int, int],
+    shape_cells_output: tuple[int, int, int],
+    weights: numba.typed.List[tuple[int, int, float]],
     intercepts: numba.typed.List,
     sweep_input: bool,
-    axis: int,
-) -> tuple[tuple[float, float, float], int, tuple[int, int, int]]:
+    axis_sweep: tuple[int] | tuple[int, int],
+) -> tuple[
+    tuple[
+        tuple[float, float, float],
+        tuple[float, float, float],
+    ],
+    tuple[int, int, int],
+    tuple[int, int, int],
+]:
     """
     Check if the current line segment crosses the boundary.
     If it does, compute the intersection and return the intersection point
     and the indices of the `grid_static` where the intersection occurs.
 
-    Parameters
-    ----------
-    point_1
-        The first point of the line segment.
-    point_2
-        The second point of the line segment.
-    bbox_boundary
-        Two points defining the bounding box of the static grid.
-    indices_boundary
-        A sequence of triangle indices defining the outer boundary of the
-        static grid.
-    intercepts
-        A sorted list of intercepts to be traversed later.
-        As new intercepts are found, they are inserted into the list.
     """
 
-    line = (point_1, point_2)
+    p1, p2 = line
 
     x, y, z = grid_static
 
-    if rg.geometry.point_is_inside_box_3d(point_2, bbox_boundary):
+    shape_cells_static = _grids.shape_centers(x.shape)
 
-        for t in range(len(triangles_boundary)):
+    if rg.geometry.point_is_inside_box_3d(p2, bbox_static):
 
-            triangle = triangles_boundary[t]
+        for axis in _arrays.axes:
 
-            tuv = rg.geometry.line_triangle_intersection_parameters(
-                line=line,
-                triangle=triangle
-            )
+            x_axis = _arrays.align_axis_right(x, axis)
+            y_axis = _arrays.align_axis_right(y, axis)
+            z_axis = _arrays.align_axis_right(z, axis)
 
-            if rg.geometry.line_intersects_triangle(tuv):
+            for direction_face in (-1, 1):
 
-                index_static = indices_boundary[t]
+                if direction_face > 0:
+                    k_face = ~0
+                else:
+                    k_face = 0
 
-                intercept = rg.geometry.line_triangle_intersection(
-                    line=line,
-                    tuv=tuv,
-                )
+                x_face = x_axis[..., k_face]
+                y_face = y_axis[..., k_face]
+                z_face = z_axis[..., k_face]
 
-                index_input, index_output = _index_input_output(
-                    index_sweep=index_sweep,
-                    index_static=index_static,
-                    sweep_input=sweep_input,
-                    axis=axis,
-                )
+                shape_face = x_face.shape
+
+                num_i, num_j = shape_face
+
+                for i0 in range(num_i - 1):
+                    for j0 in range(num_j - 1):
+
+                        i1 = i0 + 1
+                        j1 = j0 + 1
+
+                        for t in (-1, 1):
+
+                            if t > 0:
+                                ind = (
+                                    (i0, j0),
+                                    (i1, j0),
+                                    (i1, j1),
+                                )
+                            else:
+                                ind = (
+                                    (i1, j1),
+                                    (i0, j1),
+                                    (i0, j0),
+                                )
+
+                            if direction_face < 0:
+                                ind = ind[::-1]
 
-                if axis == axis_x:
+                            v0, v1, v2 = ind
 
-                    pass
+                            triangle = (
+                                (x_face[v0], y_face[v0], z_face[v0]),
+                                (x_face[v1], y_face[v1], z_face[v1]),
+                                (x_face[v2], y_face[v2], z_face[v2]),
+                            )
+
+                            tuv = rg.geometry.line_triangle_intersection_parameters(
+                                line=line,
+                                triangle=triangle
+                            )
+
+                            if rg.geometry.line_intersects_triangle(tuv):
+
+                                if direction_face > 0:
+                                    index_static = i0, j0, shape_cells_static[~0]
+                                else:
+                                    index_static = i0, j0, 0
+
+                                i, j, k = index_static
+
+                                if axis == axis_x:
+                                    index_static = j, k, i
+                                elif axis == axis_y:
+                                    index_static = k, i, j
+
+                                normal = [1 if ax == axis else 0 for ax in _arrays.axes]
+                                normal = rg.math.multiply_3d(normal, direction_face)
+
+                                p2 = _calc_and_save_intercept(
+                                    line=line,
+                                    tuv=tuv,
+                                    index_sweep=index_sweep,
+                                    index_static=index_static,
+                                    shape_cells_input=shape_cells_input,
+                                    shape_cells_output=shape_cells_output,
+                                    intercepts=intercepts,
+                                    sweep_input=sweep_input,
+                                    axis_sweep=axis_sweep,
+                                    axis_static=axis,
+                                    normal_static=normal,
+                                )
+
+                                return (p1, p2), index_sweep, index_static
+
+    return (p1, p2), index_sweep, index_static
 
 
 @numba.njit(cache=True)
@@ -377,6 +453,7 @@ def _step_inside_static(
     ],
     index_sweep: tuple[int, int, int],
     index_static: tuple[int, int, int],
+    step_index: tuple[int, int, int],
     grid_sweep: tuple[np.ndarray, np.ndarray, np.ndarray],
     grid_static: tuple[np.ndarray, np.ndarray, np.ndarray],
     shape_cells_input: tuple[int, int, int],
@@ -384,7 +461,7 @@ def _step_inside_static(
     weights: numba.typed.List[tuple[int, int, float]],
     intercepts: numba.typed.List,
     sweep_input: bool,
-    axis: int,
+    axis_sweep: tuple[int] | tuple[int, int],
 ) -> tuple[
     tuple[
         tuple[float, float, float],
@@ -410,6 +487,8 @@ def _step_inside_static(
     index_static
         The current index in the static grid.
         This is assumed to be a valid, positive index.
+    step_index
+        A step size in index space on the sweep grid.
     grid_sweep
         The vertices of the sweep grid.
         The last axis of this grid must be the sweep axis.
@@ -434,83 +513,20 @@ def _step_inside_static(
         A boolean flag indicating whether to iterate along `grid_input`
         or `grid_output`.
         If :obj:`True`, this function sweeps along `grid_input`.
-    axis
+    axis_sweep
         The logical axis of the sweep grid to iterate along.
     """
 
     p1, p2 = line
 
-    i, j, k = index_sweep
-    l, m, n = index_static
-
-    x_sweep, y_sweep, z_sweep = grid_sweep
-    x_static, y_static, z_static = grid_static
-
-    shape_sweep = x_sweep.shape
-
-    i_000 = l + 0, m + 0, n + 0
-    i_001 = l + 0, m + 0, n + 1
-    i_010 = l + 0, m + 1, n + 0
-    i_011 = l + 0, m + 1, n + 1
-    i_100 = l + 1, m + 0, n + 0
-    i_101 = l + 1, m + 0, n + 1
-    i_110 = l + 1, m + 1, n + 0
-    i_111 = l + 1, m + 1, n + 1
-
-    indices_triangles = (
-        (i_000, i_010, i_110),
-        (i_110, i_100, i_000),
-        (i_001, i_101, i_111),
-        (i_111, i_011, i_001),
-        (i_000, i_001, i_011),
-        (i_011, i_010, i_000),
-        (i_100, i_110, i_111),
-        (i_111, i_101, i_100),
-        (i_000, i_100, i_101),
-        (i_101, i_001, i_000),
-        (i_010, i_011, i_111),
-        (i_111, i_110, i_010),
-    )
-
-    axis_index = (
-        2,
-        2,
-        2,
-        2,
-        0,
-        0,
-        0,
-        0,
-        1,
-        1,
-        1,
-        1,
-    )
-
-    normals_index = (
-        (0, 0, -1),
-        (0, 0, -1),
-        (0, 0, +1),
-        (0, 0, +1),
-        (-1, 0, 0),
-        (-1, 0, 0),
-        (+1, 0, 0),
-        (+1, 0, 0),
-        (0, -1, 0),
-        (0, -1, 0),
-        (0, +1, 0),
-        (0, +1, 0),
+    polyhedron_cell = _grids.cell_boundary(
+        index=index_static,
+        grid=grid_static,
     )
 
-    for t in range(len(indices_triangles)):
-
-        i0, i1, i2 = indices_triangles[t]
-
-        v0 = (x_static[i0], y_static[i0], z_static[i0])
-        v1 = (x_static[i1], y_static[i1], z_static[i1])
-        v2 = (x_static[i2], y_static[i2], z_static[i2])
+    for t in range(len(polyhedron_cell)):
 
-        triangle = (v0, v1, v2)
+        triangle = polyhedron_cell[t]
 
         tuv = rg.geometry.line_triangle_intersection_parameters(
             line=line,
@@ -519,58 +535,133 @@ def _step_inside_static(
 
         if rg.geometry.line_intersects_triangle(tuv):
 
-            intercept = rg.geometry.line_triangle_intersection(
+            index_sweep_new = index_sweep
+
+            index_static_new, p2 = _calc_and_save_intercept(
                 line=line,
                 tuv=tuv,
+                index_sweep=index_sweep,
+                index_static=index_static,
+                shape_cells_input=shape_cells_input,
+                shape_cells_output=shape_cells_output,
+                intercepts=intercepts,
+                sweep_input=sweep_input,
+                axis_sweep=axis_sweep,
+                axis_static=_grids.cell_axes[t],
+                normal_static=_grids.cell_normals[t],
             )
 
-            p2 = intercept
+            break
+
+    else:
+        index_sweep_new = rg.math.sum_3d(index_sweep, step_index)
+        index_static_new = index_static
+
+    line = p1, p2
 
-            normal_index = normals_index[t]
+    if len(axis_sweep) == 1:
+        _calc_and_save_weights(
+            line=line,
+            index_sweep=index_sweep,
+            index_static=index_static,
+            grid_sweep=grid_sweep,
+            shape_cells_input=shape_cells_input,
+            shape_cells_output=shape_cells_output,
+            weights=weights,
+            sweep_input=sweep_input,
+            axis_sweep=axis_sweep,
+        )
 
-            index_static_new = rg.math.sum_3d(index_static, normal_index)
+    return line, index_sweep_new, index_static_new
 
-            axis_static = axis_index[t]
 
-            if index_static_new[axis_static] > index_static[axis_static]:
-                index_static_intercept = index_static_new
-            else:
-                index_static_intercept = index_static
+@numba.njit(cache=True)
+def _calc_and_save_intercept(
+    line: tuple[
+        tuple[float, float, float],
+        tuple[float, float, float],
+    ],
+    tuv: tuple[float, float, float],
+    index_sweep: tuple[int, int, int],
+    index_static: tuple[int, int, int],
+    shape_cells_input: tuple[int, int, int],
+    shape_cells_output: tuple[int, int, int],
+    intercepts: numba.typed.List,
+    sweep_input: bool,
+    axis_sweep: tuple[int] | tuple[int, int],
+    axis_static: int,
+    normal_static: tuple[int, int, int],
+) -> tuple[
+    tuple[int, int, int],
+    tuple[float, float, float],
+]:
+
+    intercept = rg.geometry.line_triangle_intersection(
+        line=line,
+        tuv=tuv,
+    )
 
-            axes_orthogonal = [ax for ax in _arrays.axes if ax != axis]
+    index_static_new = rg.math.sum_3d(index_static, normal_static)
 
-            for axis_orthogonal in axes_orthogonal:
+    if index_static_new[axis_static] > index_static[axis_static]:
+        index_static_intercept = index_static_new
+    else:
+        index_static_intercept = index_static
 
-                if sweep_input:
-                    axis_input = axis_orthogonal
-                    axis_output = axis_static
-                else:
-                    axis_input = axis_static
-                    axis_output = axis_orthogonal
-
-                index_input, index_output = _index_input_output(
-                    index_sweep=index_sweep,
-                    index_static=index_static_intercept,
-                    sweep_input=sweep_input,
-                    axis=axis,
-                )
+    axes_orthogonal = [ax for ax in _arrays.axes if ax not in axis_sweep]
 
-                i_input = index_input[axis_input]
-                i_output = index_output[axis_output]
+    for axis_orthogonal in axes_orthogonal:
 
-                _intercepts.insert_intercept(
-                    intercepts=intercepts[axis_input][axis_output][i_input][i_output],
-                    intercept_new=(
-                        _arrays.index_flat(index_input, shape_cells_input),
-                        _arrays.index_flat(index_output, shape_cells_output),
-                        intercept,
-                    ),
-                )
+        if sweep_input:
+            axis_input = axis_orthogonal
+            axis_output = axis_static
+        else:
+            axis_input = axis_static
+            axis_output = axis_orthogonal
 
-                break
+        index_input, index_output = _index_input_output(
+            index_sweep=index_sweep,
+            index_static=index_static_intercept,
+            sweep_input=sweep_input,
+            axis_sweep=axis_sweep,
+        )
 
-    else:
-        index_sweep = i, j, k + 1
+        i_input = index_input[axis_input]
+        i_output = index_output[axis_output]
+
+        _intercepts.insert_intercept(
+            intercepts=intercepts[axis_input][axis_output][i_input][i_output],
+            intercept_new=(
+                _arrays.index_flat(index_input, shape_cells_input),
+                _arrays.index_flat(index_output, shape_cells_output),
+                intercept,
+            ),
+        )
+
+    return index_static_new, intercept
+
+
+@numba.njit(cache=True)
+def _calc_and_save_weights(
+    line: tuple[
+        tuple[float, float, float],
+        tuple[float, float, float],
+    ],
+    index_sweep: tuple[int, int, int],
+    index_static: tuple[int, int, int],
+    grid_sweep: tuple[np.ndarray, np.ndarray, np.ndarray],
+    shape_cells_input: tuple[int, int, int],
+    shape_cells_output: tuple[int, int, int],
+    weights: numba.typed.List[tuple[int, int, float]],
+    sweep_input: bool,
+    axis_sweep: tuple[int] | tuple[int, int],
+):
+
+    x_sweep, y_sweep, z_sweep = grid_sweep
+
+    shape_sweep = x_sweep.shape
+
+    i, j, k = index_sweep
 
     i_left = i - 1
     i_right = i
@@ -592,6 +683,8 @@ def _step_inside_static(
         z_sweep[index_sweep_lower],
     )
 
+    p1, p2 = line
+
     volume_left = rg.geometry.volume_tetrahedron(p0_left, p1, p2)
     volume_lower = -rg.geometry.volume_tetrahedron(p0_lower, p1, p1)
 
@@ -604,7 +697,7 @@ def _step_inside_static(
                 index_sweep=index_sweep_00,
                 index_static=index_static,
                 sweep_input=sweep_input,
-                axis=axis,
+                axis_sweep=axis_sweep,
             )
 
             index_input_00 = _arrays.index_flat(index_input_00, shape_cells_input)
@@ -621,7 +714,7 @@ def _step_inside_static(
                 index_sweep=index_sweep_01,
                 index_static=index_static,
                 sweep_input=sweep_input,
-                axis=axis,
+                axis_sweep=axis_sweep,
             )
 
             index_input_01 = _arrays.index_flat(index_input_01, shape_cells_input)
@@ -639,7 +732,7 @@ def _step_inside_static(
                 index_sweep=index_sweep_10,
                 index_static=index_static,
                 sweep_input=sweep_input,
-                axis=axis,
+                axis_sweep=axis_sweep,
             )
 
             index_input_10 = _arrays.index_flat(index_input_10, shape_cells_input)
@@ -648,17 +741,13 @@ def _step_inside_static(
             volume_lower_right = volume_lower
             weights.append((index_input_10, index_output_10, volume_lower_right))
 
-    line = p1, p2
-
-    return line, index_sweep, index_static
-
 
 @numba.njit(cache=True)
 def _grid_sweep_static(
     grid_input: tuple[np.ndarray, np.ndarray, np.ndarray],
     grid_output: tuple[np.ndarray, np.ndarray, np.ndarray],
     sweep_input: bool,
-    axis: int,
+    axis_sweep: tuple[int] | tuple[int, int],
 ) -> tuple[
     tuple[np.ndarray, np.ndarray, np.ndarray],
     tuple[np.ndarray, np.ndarray, np.ndarray],
@@ -677,7 +766,7 @@ def _grid_sweep_static(
         A boolean flag indicating whether to iterate along `grid_input`
         or `grid_output`.
         If :obj:`True`, this function sweeps along `grid_input`.
-    axis
+    axis_sweep
         The logical axis of the `grid_sweep` to iterate along.
     """
     if sweep_input:
@@ -689,6 +778,7 @@ def _grid_sweep_static(
 
     grid_sweep_x, grid_sweep_y, grid_sweep_z = grid_sweep
 
+    axis = axis_sweep[~0]
     grid_sweep_x = _arrays.align_axis_right(grid_sweep_x, axis)
     grid_sweep_y = _arrays.align_axis_right(grid_sweep_y, axis)
     grid_sweep_z = _arrays.align_axis_right(grid_sweep_z, axis)
@@ -703,8 +793,27 @@ def _index_input_output(
     index_sweep: tuple[int, int, int],
     index_static: tuple[int, int, int],
     sweep_input: bool,
-    axis: int,
+    axis_sweep: tuple[int] | tuple[int, int],
 ):
+    """
+    Convert indices on the sweep and static grids into indices on the
+    input/output grids.
+
+    Parameters
+    ----------
+    index_sweep
+        The index on the sweep grid to convert.
+    index_static
+        The index on the static grid to convert.
+    sweep_input
+        A boolean flag indicating whether to iterate along `grid_input`
+        or `grid_output`.
+        If :obj:`True`, this function sweeps along `grid_input`.
+    axis_sweep
+        The logical axis of the `grid_sweep` we're iterating along.
+    """
+
+    axis = axis_sweep[~0]
 
     i_sweep, j_sweep, k_sweep = index_sweep
 

From 6c22d03922c96c2c54e583e45b8a81bfb7d12a99 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Thu, 8 May 2025 11:15:08 -0600
Subject: [PATCH 17/26] Added `_intercepts.insert_intercept()`.

---
 .../_weights_conservative_3d/_intercepts.py   | 66 +++++++++++++++++
 .../_intercepts_test.py                       | 71 +++++++++++++++++++
 2 files changed, 137 insertions(+)

diff --git a/regridding/_weights/_weights_conservative_3d/_intercepts.py b/regridding/_weights/_weights_conservative_3d/_intercepts.py
index 8bc6ac3..8710b1b 100644
--- a/regridding/_weights/_weights_conservative_3d/_intercepts.py
+++ b/regridding/_weights/_weights_conservative_3d/_intercepts.py
@@ -46,6 +46,7 @@ def empty(
     return intercepts
 
 
+@numba.njit(cache=True)
 def insert_intercept(
     intercepts: numba.typed.List[tuple[int, int, tuple[float, float, float]]],
     intercept_new: tuple[int, int, tuple[float, float, float]],
@@ -54,6 +55,28 @@ def insert_intercept(
     Insert a new intercept into the current list of intercepts at the correct
     point so that the list maintains sorted order.
 
+    Parameters
+    ----------
+    intercepts
+        The current list of intercepts.
+    intercept_new
+        A new intercept to insert into the list.
+    """
+    index = _bisect_intercepts(intercepts, intercept_new)
+
+    intercepts.insert(index, intercept_new)
+
+
+@numba.njit(cache=True)
+def _bisect_intercepts(
+    intercepts: numba.typed.List[tuple[int, int, tuple[float, float, float]]],
+    intercept_new: tuple[int, int, tuple[float, float, float]],
+) -> int:
+    """
+    Given an ordered sequence of intercepts,
+    find the index for which a new intercept should be inserted to maintain
+    the ordering
+
     Parameters
     ----------
     intercepts
@@ -62,7 +85,50 @@ def insert_intercept(
         A new intercept to insert into the list.
     """
 
+    num_intercepts = len(intercepts)
 
+    _, _, intercept_new = intercept_new
+
+    index_left = 0
+    index_right = num_intercepts- 1
+
+    _, _, intercept_left = intercepts[index_left]
+    _, _, intercept_right = intercepts[index_right]
+
+    t_start = _line_point_closest_approach_parameter(
+        line=(intercept_left, intercept_right),
+        point=intercept_new,
+    )
+
+    if t_start < 0:
+        return 0
+
+    if t_start >= 1:
+        return num_intercepts
+
+    while (index_right - index_left) > 1:
+
+        index_middle = (index_left + index_right) // 2
+
+        _, _, intercept_middle = intercepts[index_middle]
+
+        t_left = _line_point_closest_approach_parameter(
+            line=(intercept_left, intercept_middle),
+            point=intercept_new,
+        )
+
+        if 0 <= t_left < 1:
+            index_right = index_middle
+            intercept_right = intercept_middle
+
+        else:
+            index_left = index_middle
+            intercept_left = intercept_middle
+
+    return index_right
+
+
+@numba.njit(cache=True)
 def _line_point_closest_approach_parameter(
     line: tuple[
         tuple[float, float, float],
diff --git a/regridding/_weights/_weights_conservative_3d/_intercepts_test.py b/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
index cba30a9..fccbfcf 100644
--- a/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
+++ b/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
@@ -4,6 +4,77 @@
 from . import _intercepts
 
 
+intercepts1 = numba.typed.List([
+    (0, 0, (0, 0, 0)),
+    (0, 1, (1, 1, 1)),
+    (0, 2, (2, 2, 2)),
+])
+
+intercepts2 = numba.typed.List([
+    (0, 0, (0, 0, 0)),
+    (0, 1, (-1, 1, -1)),
+    (0, 2, (-2, 2, -2)),
+])
+
+
+@pytest.mark.parametrize(
+    argnames="intercepts,intercept_new,result_expected",
+    argvalues=[
+        (
+            intercepts1,
+            (0, 0, (-1, -1, -1)),
+            0,
+        ),
+        (
+            intercepts1,
+            (0, 0, (0.5, 0.25, 0.5)),
+            1,
+        ),
+        (
+            intercepts1,
+            (0, 0, (1.5, 1.25, 1.75)),
+            2,
+        ),
+        (
+            intercepts1,
+            (0, 0, (2.5, 2.25, 2.75)),
+            3,
+        ),
+        (
+            intercepts2,
+            (0, 0, (1, -1, 1)),
+            0,
+        ),
+        (
+            intercepts2,
+            (0, 0, (-0.5, 0.25, -0.5)),
+            1,
+        ),
+        (
+            intercepts2,
+            (0, 0, (-1.5, 1.25, -1.75)),
+            2,
+        ),
+        (
+            intercepts2,
+            (0, 0, (-2.5, 2.25, -2.75)),
+            3,
+        ),
+    ]
+)
+def test_bisect_intercepts(
+    intercepts: numba.typed.List[tuple[int, int, tuple[float, float, float]]],
+    intercept_new: tuple[int, int, tuple[float, float, float]],
+    result_expected: int,
+):
+    result = _intercepts._bisect_intercepts(
+        intercepts=intercepts,
+        intercept_new=intercept_new,
+    )
+
+    assert result == result_expected
+
+
 @pytest.mark.parametrize(
     argnames="line, point, result_expected",
     argvalues=[

From 0c2c08bdd0f959ff0680ebee850659acb3851117 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Wed, 14 May 2025 06:47:07 -0600
Subject: [PATCH 18/26] Fixed a ton of bugs so that the program runs now. Still
 doesn't produce correct results.

---
 regridding/_weights/_weights_conservative.py  |  21 +-
 .../_weights_conservative_3d/__init__.py      |   8 +
 .../_weights_conservative_3d/_grids.py        |   6 +-
 .../_weights_conservative_3d/_intercepts.py   |  43 ++-
 .../_intercepts_test.py                       |  28 ++
 .../_weights_conservative_3d.py               | 294 ++++++++++++------
 6 files changed, 278 insertions(+), 122 deletions(-)

diff --git a/regridding/_weights/_weights_conservative.py b/regridding/_weights/_weights_conservative.py
index 67fc9c1..04f183d 100644
--- a/regridding/_weights/_weights_conservative.py
+++ b/regridding/_weights/_weights_conservative.py
@@ -6,6 +6,7 @@
 from regridding import _util
 from ._weights_conservative_1d import weights_conservative_1d
 from regridding._conservative_ramshaw import _conservative_ramshaw
+from ._weights_conservative_3d import weights_conservative_3d
 
 
 def _weights_conservative(
@@ -127,9 +128,25 @@ def _weights_conservative(
                     weights_input=weights_input_index,
                 )
 
-            else:  # pragma: nocover
+            elif len(axis_input) == 3:
+            x_input, y_input, z_input = coordinates_input
+            x_output, y_output, z_output = coordinates_output
+            weights[index] = weights_conservative_3d(
+                grid_input=(
+                    x_input[index_vertices_input],
+                    y_input[index_vertices_input],
+                    z_input[index_vertices_input],
+                ),
+                grid_output=(
+                    x_output[index_vertices_output],
+                    y_output[index_vertices_output],
+                    z_output[index_vertices_output],
+                ),
+            )
+
+        else:  # pragma: nocover
                 raise NotImplementedError(
-                    "Regridding operations greater than 2D are not supported"
+                    "Regridding operations greater than 3D are not supported"
                 )
 
     return weights, shape_values_input, shape_values_output
diff --git a/regridding/_weights/_weights_conservative_3d/__init__.py b/regridding/_weights/_weights_conservative_3d/__init__.py
index e69de29..8a91a8b 100644
--- a/regridding/_weights/_weights_conservative_3d/__init__.py
+++ b/regridding/_weights/_weights_conservative_3d/__init__.py
@@ -0,0 +1,8 @@
+"""
+A subpackage for performing 3D conservative resampling
+"""
+from ._weights_conservative_3d import weights_conservative_3d
+
+__all__ = [
+    "weights_conservative_3d",
+]
diff --git a/regridding/_weights/_weights_conservative_3d/_grids.py b/regridding/_weights/_weights_conservative_3d/_grids.py
index 44e6802..ca16d46 100644
--- a/regridding/_weights/_weights_conservative_3d/_grids.py
+++ b/regridding/_weights/_weights_conservative_3d/_grids.py
@@ -304,7 +304,7 @@ def grid_boundary(
     for i0 in range(shape_x - 1):
         for j0 in range(shape_y - 1):
 
-            k1 = shape_z
+            k1 = shape_z - 1
 
             i1 = i0 + 1
             j1 = j0 + 1
@@ -356,7 +356,7 @@ def grid_boundary(
     for j0 in range(shape_y - 1):
         for k0 in range(shape_z - 1):
 
-            i1 = shape_x
+            i1 = shape_x - 1
 
             j1 = j0 + 1
             k1 = k0 + 1
@@ -408,7 +408,7 @@ def grid_boundary(
     for i0 in range(shape_x - 1):
         for k0 in range(shape_z - 1):
 
-            j1 = shape_y
+            j1 = shape_y - 1
 
             i1 = i0 + 1
             k1 = k0 + 1
diff --git a/regridding/_weights/_weights_conservative_3d/_intercepts.py b/regridding/_weights/_weights_conservative_3d/_intercepts.py
index 8710b1b..3642dbf 100644
--- a/regridding/_weights/_weights_conservative_3d/_intercepts.py
+++ b/regridding/_weights/_weights_conservative_3d/_intercepts.py
@@ -87,10 +87,13 @@ def _bisect_intercepts(
 
     num_intercepts = len(intercepts)
 
+    if num_intercepts < 2:
+        return num_intercepts
+
     _, _, intercept_new = intercept_new
 
     index_left = 0
-    index_right = num_intercepts- 1
+    index_right = num_intercepts - 1
 
     _, _, intercept_left = intercepts[index_left]
     _, _, intercept_right = intercepts[index_right]
@@ -304,53 +307,59 @@ def sweep(
                             shape=shape_centers_output,
                         )
 
-                        i0_input_lower = _arrays.index_flat(
+                        n0_input_lower = _arrays.index_flat(
                             index=i0_input_lower,
-                            shape=shape_input,
+                            shape=shape_centers_input,
                         )
-                        i0_input_upper = _arrays.index_flat(
+                        n0_input_upper = _arrays.index_flat(
                             index=i0_input_upper,
-                            shape=shape_input,
+                            shape=shape_centers_input,
                         )
-                        i0_output_lower = _arrays.index_flat(
+                        n0_output_lower = _arrays.index_flat(
                             index=i0_output_lower,
-                            shape=shape_output,
+                            shape=shape_centers_output,
                         )
-                        i0_output_upper = _arrays.index_flat(
+                        n0_output_upper = _arrays.index_flat(
                             index=i0_output_upper,
-                            shape=shape_output,
+                            shape=shape_centers_output,
                         )
 
                         if input_lower_in_bounds:
+
+                            volume_input_lower = volume_input[i0_input_lower]
+
                             if output_lower_in_bounds:
                                 weight_lower_lower = (
-                                    i0_input_lower,
-                                    i0_output_lower,
+                                    n0_input_lower,
+                                    n0_output_lower,
                                     (-vol_input - vol_output) / volume_input_lower,
                                 )
                                 weights.append(weight_lower_lower)
 
                             if output_upper_in_bounds:
                                 weight_lower_upper = (
-                                    i0_input_lower,
-                                    i0_output_upper,
+                                    n0_input_lower,
+                                    n0_output_upper,
                                     (-vol_input + vol_output) / volume_input_lower,
                                 )
                                 weights.append(weight_lower_upper)
 
                         if input_upper_in_bounds:
+
+                            volume_input_upper = volume_input[i0_input_upper]
+
                             if output_lower_in_bounds:
                                 weight_upper_lower = (
-                                    i0_input_upper,
-                                    i0_output_lower,
+                                    n0_input_upper,
+                                    n0_output_lower,
                                     (vol_input - vol_output) / volume_input_upper,
                                 )
                                 weights.append(weight_upper_lower)
 
                             if output_upper_in_bounds:
                                 weight_upper_upper = (
-                                    i0_input_upper,
-                                    i0_output_upper,
+                                    n0_input_upper,
+                                    n0_output_upper,
                                     (vol_input + vol_output) / volume_input_upper,
                                 )
                                 weights.append(weight_upper_upper)
diff --git a/regridding/_weights/_weights_conservative_3d/_intercepts_test.py b/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
index fccbfcf..a28932c 100644
--- a/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
+++ b/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
@@ -17,6 +17,34 @@
 ])
 
 
+@pytest.mark.parametrize(
+    argnames="intercepts,intercept_new",
+    argvalues=[
+        (
+            numba.typed.List([
+                (0, 0, (0, 0, 0)),
+                (0, 1, (-1, 1, -1)),
+                (0, 2, (-2, 2, -2)),
+            ]),
+            (0, 0, (-1, -1, -1)),
+        ),
+        (
+            _intercepts.empty((11,11,11), (11,11,11))[0][0][0][0],
+            (0, 0, (-1, -1, -1)),
+        ),
+    ],
+)
+def test_insert_intercept(
+    intercepts: numba.typed.List[tuple[int, int, tuple[float, float, float]]],
+    intercept_new: tuple[int, int, tuple[float, float, float]],
+):
+    length_old = len(intercepts)
+
+    _intercepts.insert_intercept(intercepts, intercept_new)
+
+    assert len(intercepts) == length_old + 1
+
+
 @pytest.mark.parametrize(
     argnames="intercepts,intercept_new,result_expected",
     argvalues=[
diff --git a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
index c20360d..dc6b5e2 100644
--- a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
+++ b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
@@ -1,15 +1,23 @@
 import sys
 import numpy as np
-import numpy.typing as npt
 import numba
+from numba import literal_unroll
 import regridding as rg
-from . import  (
+from . import (
     _arrays,
-    _volumes,
     _grids,
     _intercepts,
 )
-from _arrays import axis_x, axis_y, axis_z
+from ._arrays import axis_x, axis_y, axis_z
+
+axes = (
+    (axis_x, ),
+    (axis_y, ),
+    (axis_z, ),
+    (axis_x, axis_y),
+    (axis_y, axis_z),
+    (axis_z, axis_x),
+)
 
 __all__ = [
     "weights_conservative_3d",
@@ -46,7 +54,8 @@ def weights_conservative_3d(
     for x in range(0):
         weights.append((0, 0, 0.0))
 
-    volume_input = _volumes.volume_grid(grid_input)
+    volume_input = _grids.grid_volume(grid_input)
+
 
     intercepts = _intercepts.empty(shape_input, shape_output)
 
@@ -105,49 +114,67 @@ def _sweep_grid(
         If :obj:`True`, this function sweeps along `grid_input`.
     """
 
-    grid_sweep, grid_static = _grid_sweep_static(grid_input, grid_output, sweep_input)
+    shape_input = grid_input[0].shape
+    shape_output = grid_output[0].shape
 
-    grid_static_x, grid_static_y, grid_static_z = grid_static
+    shape_cells_input = _grids.shape_centers(shape_input)
+    shape_cells_output = _grids.shape_centers(shape_output)
 
-    bbox_boundary = (
-        (grid_static_x.min(), grid_static_y.min(), grid_static_z.min()),
-        (grid_static_x.max(), grid_static_y.max(), grid_static_z.max()),
+    grid_sweep, grid_static = _grid_sweep_static(
+        grid_input=grid_input,
+        grid_output=grid_output,
+        sweep_input=sweep_input,
     )
 
-    indices_boundary, triangles_boundary = _grids.grid_boundary(grid_static)
+    x_sweep, y_sweep, z_sweep = grid_sweep
+    x_static, y_static, z_static = grid_static
 
-    axes = (
-        (axis_x, ),
-        (axis_y, ),
-        (axis_z, ),
-        (axis_x, axis_y),
-        (axis_y, axis_z),
-        (axis_z, axis_x),
+    bbox_static = (
+        (x_static.min(), y_static.min(), z_static.min()),
+        (x_static.max(), y_static.max(), z_static.max()),
     )
 
+    _, boundary_static = _grids.grid_boundary(grid_static)
+
     for axis in axes:
+    # for axis in literal_unroll(axes):
+
+        axis_last = axis[~0]
+
         _sweep_along_axis(
-            grid_input=grid_input,
-            grid_output=grid_output,
+            grid_sweep=(
+                _arrays.align_axis_right(x_sweep, axis_last),
+                _arrays.align_axis_right(y_sweep, axis_last),
+                _arrays.align_axis_right(z_sweep, axis_last),
+            ),
+            grid_static=grid_static,
+            bbox_static=bbox_static,
+            boundary_static=boundary_static,
             volume_input=volume_input,
-            bbox_boundary=bbox_boundary,
-            indices_boundary=indices_boundary,
-            triangles_boundary=triangles_boundary,
+            shape_cells_input=shape_cells_input,
+            shape_cells_output=shape_cells_output,
             weights=weights,
             intercepts=intercepts,
             sweep_input=sweep_input,
-            axes=axis,
+            axis_sweep=axis,
         )
 
 
 @numba.njit(cache=True, parallel=True)
 def _sweep_along_axis(
-    grid_input: tuple[np.ndarray, np.ndarray, np.ndarray],
-    grid_output: tuple[np.ndarray, np.ndarray, np.ndarray],
+    grid_sweep: tuple[np.ndarray, np.ndarray, np.ndarray],
+    grid_static: tuple[np.ndarray, np.ndarray, np.ndarray],
+    bbox_static: tuple[tuple[float, float, float], tuple[float, float, float]],
+    boundary_static: numba.typed.List[
+        tuple[
+            tuple[float, float, float],
+            tuple[float, float, float],
+            tuple[float, float, float],
+        ],
+    ],
     volume_input: np.ndarray,
-    bbox_boundary: tuple[tuple[float, float, float], tuple[float, float, float]],
-    indices_boundary: tuple[np.ndarray, np.ndarray, np.ndarray],
-    triangles_boundary: tuple[np.ndarray, np.ndarray, np.ndarray],
+    shape_cells_input: tuple[int, int, int],
+    shape_cells_output: tuple[int, int, int],
     weights: numba.typed.List[tuple[int, int, float]],
     intercepts: numba.typed.List,
     sweep_input: bool,
@@ -159,19 +186,26 @@ def _sweep_along_axis(
 
     Parameters
     ----------
-    grid_input
-        The vertices of the old grid.
-    grid_output
-        The vertices of the new grid.
-    volume_input
-        The area of the input grid which will used to compute the weight.
-    bbox_boundary
+    grid_sweep
+        The vertices of the sweep grid.
+        The last axis of this grid must be the sweep axis.
+    grid_static
+        Coordinates of the static grid.
+    bbox_static
         Two points defining the bounding box of the static grid.
-    triangles_boundary
+    boundary_static
         A sequence of triangles defining the outer surface of the static grid.
-    indices_boundary
-        The index of the static grid corresponding to each triangle in
-        `coord_boundary`.
+    volume_input
+        The volume of each cell in the input grid.
+    shape_cells_input
+        The number of cells along each axis of the input grid.
+        This could be calculated from the other arguments,
+        but we don't to save computation time.
+    shape_cells_output
+        The number of cells along each axis of the output grid.
+        Like `shape_cells_input`,
+        this could be calculated from the other arguments,
+        but we don't to save computation time.
     weights
         The current list of weights.
         New weights will be appended to this list.
@@ -216,16 +250,16 @@ def _sweep_along_axis(
     if len(axis_sweep) == 2:
         d_start = -(shape_sweep_z - 2)
         d_end = shape_sweep_y - 1
-        step = (0, 1, 1)
+        step_index = (0, 1, 1)
     else:
         d_start = 0
         d_end = shape_sweep_y
-        step = (0, 0, 1)
+        step_index = (0, 0, 1)
 
     j_end = d_end
-    k_end = shape_cells_sweep[2] - 1
+    k_end = shape_cells_sweep[2]
 
-    for i in numba.prange(shape_sweep_x):
+    for i in range(shape_sweep_x):
 
         i = numba.types.int64(i)
 
@@ -252,11 +286,11 @@ def _sweep_along_axis(
 
             if rg.geometry.point_is_inside_box_3d(
                 point=point_1,
-                box=bbox_boundary,
+                box=bbox_static,
             ):
                 if rg.geometry.point_is_inside_polyhedron(
                     point=point_1,
-                    polyhedron=triangles_boundary,
+                    polyhedron=boundary_static,
                 ):
                     index_static = _grids.index_of_point_brute(
                         point=point_1,
@@ -264,9 +298,9 @@ def _sweep_along_axis(
                     )
                     sweep_is_outside_static = False
 
-            while (j < j_end) and (k < k_end):
+            while (index_sweep[1] < j_end) and (index_sweep[2] < k_end):
 
-                index_sweep_new = rg.math.sum_3d(index_sweep, step)
+                index_sweep_new = rg.math.sum_3d(index_sweep, step_index)
 
                 point_2 = (
                     x_sweep[index_sweep_new],
@@ -274,15 +308,21 @@ def _sweep_along_axis(
                     z_sweep[index_sweep_new],
                 )
 
+                line = point_1, point_2
                 if sweep_is_outside_static:
 
-                    point_2, k, index_static = _step_outside_static(
-                        point_1=point_1,
-                        point_2=point_2,
-                        bbox_boundary=bbox_boundary,
-                        indices_boundary=indices_boundary,
-                        triangles_boundary=triangles_boundary,
+                    line, index_sweep, index_static = _step_outside_static(
+                        line=line,
+                        index_sweep=index_sweep,
+                        index_static=index_static,
+                        step_index=step_index,
+                        grid_static=grid_static,
+                        bbox_static=bbox_static,
+                        shape_cells_input=shape_cells_input,
+                        shape_cells_output=shape_cells_output,
                         intercepts=intercepts,
+                        sweep_input=sweep_input,
+                        axis_sweep=axis_sweep,
                     )
 
                     sweep_is_outside_static = index_static[0] == sys.maxsize
@@ -293,9 +333,16 @@ def _sweep_along_axis(
                         line=line,
                         index_sweep=index_sweep,
                         index_static=index_static,
+                        step_index=step_index,
+                        grid_sweep=grid_sweep,
                         grid_static=grid_static,
+                        volume_input=volume_input,
+                        shape_cells_input=shape_cells_input,
+                        shape_cells_output=shape_cells_output,
                         weights=weights,
                         intercepts=intercepts,
+                        sweep_input=sweep_input,
+                        axis_sweep=axis_sweep,
                     )
 
                     if not _arrays.index_in_bounds(
@@ -304,6 +351,8 @@ def _sweep_along_axis(
                     ):
                         break
 
+                point_1 = line[1]
+
     for i in range(shape_sweep_x):
         weight_i = weight[i]
         for w in range(len(weight_i)):
@@ -319,12 +368,10 @@ def _step_outside_static(
     index_sweep: tuple[int, int, int],
     index_static: tuple[int, int, int],
     step_index: tuple[int, int, int],
-    grid_sweep: tuple[np.ndarray, np.ndarray, np.ndarray],
     grid_static: tuple[np.ndarray, np.ndarray, np.ndarray],
     bbox_static: tuple[tuple[float, float, float], tuple[float, float, float]],
     shape_cells_input: tuple[int, int, int],
     shape_cells_output: tuple[int, int, int],
-    weights: numba.typed.List[tuple[int, int, float]],
     intercepts: numba.typed.List,
     sweep_input: bool,
     axis_sweep: tuple[int] | tuple[int, int],
@@ -341,6 +388,40 @@ def _step_outside_static(
     If it does, compute the intersection and return the intersection point
     and the indices of the `grid_static` where the intersection occurs.
 
+    Parameters
+    ----------
+    line
+        The current line segment of the sweep grid.
+    index_sweep
+        The index of the current vertex in the sweep grid
+        This is assumed to be a valid, positive index.
+    index_static
+        The index of the current cell in the static grid.
+        This is assumed to be a valid, positive index.
+    step_index
+        A step size in index space on the sweep grid.
+    grid_static
+        The vertices of the static grid.
+    bbox_static
+        Two points defining the bounding box of the static grid.
+    shape_cells_input
+        The number of cells along each axis of the input grid.
+        This could be calculated from the other arguments,
+        but we don't to save computation time.
+    shape_cells_output
+        The number of cells along each axis of the output grid.
+        Like `shape_cells_input`,
+        this could be calculated from the other arguments,
+        but we don't to save computation time.
+    intercepts
+        A sorted list of intercepts to be traversed later.
+        As new intercepts are found, they are inserted into the list.
+    sweep_input
+        A boolean flag indicating whether to iterate along `grid_input`
+        or `grid_output`.
+        If :obj:`True`, this function sweeps along `grid_input`.
+    axis_sweep
+        The logical axis of the sweep grid to iterate along.
     """
 
     p1, p2 = line
@@ -412,21 +493,21 @@ def _step_outside_static(
                             if rg.geometry.line_intersects_triangle(tuv):
 
                                 if direction_face > 0:
-                                    index_static = i0, j0, shape_cells_static[~0]
+                                    index_static = i0, j0, shape_cells_static[axis] - 1
                                 else:
                                     index_static = i0, j0, 0
 
                                 i, j, k = index_static
 
                                 if axis == axis_x:
-                                    index_static = j, k, i
-                                elif axis == axis_y:
                                     index_static = k, i, j
+                                elif axis == axis_y:
+                                    index_static = j, k, i
 
                                 normal = [1 if ax == axis else 0 for ax in _arrays.axes]
-                                normal = rg.math.multiply_3d(normal, direction_face)
+                                normal = rg.math.multiply_3d(direction_face, normal)
 
-                                p2 = _calc_and_save_intercept(
+                                _, p2 = _calc_and_save_intercept(
                                     line=line,
                                     tuv=tuv,
                                     index_sweep=index_sweep,
@@ -442,6 +523,8 @@ def _step_outside_static(
 
                                 return (p1, p2), index_sweep, index_static
 
+    index_sweep = rg.math.sum_3d(index_sweep, step_index)
+
     return (p1, p2), index_sweep, index_static
 
 
@@ -456,6 +539,7 @@ def _step_inside_static(
     step_index: tuple[int, int, int],
     grid_sweep: tuple[np.ndarray, np.ndarray, np.ndarray],
     grid_static: tuple[np.ndarray, np.ndarray, np.ndarray],
+    volume_input: np.ndarray,
     shape_cells_input: tuple[int, int, int],
     shape_cells_output: tuple[int, int, int],
     weights: numba.typed.List[tuple[int, int, float]],
@@ -482,7 +566,7 @@ def _step_inside_static(
     line
         The current line segment of the sweep grid.
     index_sweep
-        The current index in the sweep grid.
+        The index of the current vertex in the sweep grid
         This is assumed to be a valid, positive index.
     index_static
         The current index in the static grid.
@@ -493,7 +577,9 @@ def _step_inside_static(
         The vertices of the sweep grid.
         The last axis of this grid must be the sweep axis.
     grid_static
-        Coordinates of the static grid.
+        The vertices of the static grid.
+    volume_input
+        The volume of each cell in the input grid.
     shape_cells_input
         The number of cells along each axis of the input grid.
         This could be calculated from the other arguments,
@@ -535,23 +621,25 @@ def _step_inside_static(
 
         if rg.geometry.line_intersects_triangle(tuv):
 
-            index_sweep_new = index_sweep
-
-            index_static_new, p2 = _calc_and_save_intercept(
-                line=line,
-                tuv=tuv,
-                index_sweep=index_sweep,
-                index_static=index_static,
-                shape_cells_input=shape_cells_input,
-                shape_cells_output=shape_cells_output,
-                intercepts=intercepts,
-                sweep_input=sweep_input,
-                axis_sweep=axis_sweep,
-                axis_static=_grids.cell_axes[t],
-                normal_static=_grids.cell_normals[t],
-            )
+            if tuv[0] > 1e-8:
+
+                index_sweep_new = index_sweep
+
+                index_static_new, p2 = _calc_and_save_intercept(
+                    line=line,
+                    tuv=tuv,
+                    index_sweep=index_sweep,
+                    index_static=index_static,
+                    shape_cells_input=shape_cells_input,
+                    shape_cells_output=shape_cells_output,
+                    intercepts=intercepts,
+                    sweep_input=sweep_input,
+                    axis_sweep=axis_sweep,
+                    axis_static=_grids.cell_axes[t],
+                    normal_static=_grids.cell_normals[t],
+                )
 
-            break
+                break
 
     else:
         index_sweep_new = rg.math.sum_3d(index_sweep, step_index)
@@ -565,6 +653,7 @@ def _step_inside_static(
             index_sweep=index_sweep,
             index_static=index_static,
             grid_sweep=grid_sweep,
+            volume_input=volume_input,
             shape_cells_input=shape_cells_input,
             shape_cells_output=shape_cells_output,
             weights=weights,
@@ -601,12 +690,16 @@ def _calc_and_save_intercept(
         tuv=tuv,
     )
 
+    shape_input = [n + 1 for n in shape_cells_input]
+    shape_output = [n + 1 for n in shape_cells_output]
+
     index_static_new = rg.math.sum_3d(index_static, normal_static)
+    # print(f"{index_static=}")
+    # print(f"{index_static_new=}")
 
+    index_static_intercept = index_static
     if index_static_new[axis_static] > index_static[axis_static]:
         index_static_intercept = index_static_new
-    else:
-        index_static_intercept = index_static
 
     axes_orthogonal = [ax for ax in _arrays.axes if ax not in axis_sweep]
 
@@ -632,8 +725,8 @@ def _calc_and_save_intercept(
         _intercepts.insert_intercept(
             intercepts=intercepts[axis_input][axis_output][i_input][i_output],
             intercept_new=(
-                _arrays.index_flat(index_input, shape_cells_input),
-                _arrays.index_flat(index_output, shape_cells_output),
+                _arrays.index_flat(index_input, shape_input),
+                _arrays.index_flat(index_output, shape_output),
                 intercept,
             ),
         )
@@ -650,6 +743,7 @@ def _calc_and_save_weights(
     index_sweep: tuple[int, int, int],
     index_static: tuple[int, int, int],
     grid_sweep: tuple[np.ndarray, np.ndarray, np.ndarray],
+    volume_input: np.ndarray,
     shape_cells_input: tuple[int, int, int],
     shape_cells_output: tuple[int, int, int],
     weights: numba.typed.List[tuple[int, int, float]],
@@ -700,11 +794,13 @@ def _calc_and_save_weights(
                 axis_sweep=axis_sweep,
             )
 
+            volume_input_lower_left = volume_input[index_input_00]
+
             index_input_00 = _arrays.index_flat(index_input_00, shape_cells_input)
             index_output_00 = _arrays.index_flat(index_output_00, shape_cells_output)
 
-            volume_lower_left = -volume_lower - volume_left
             weights.append((index_input_00, index_output_00, volume_lower_left))
+            volume_lower_left = (-volume_lower - volume_left) / volume_input_lower_left
 
         if j_upper < (shape_sweep[1] - 1):
 
@@ -717,11 +813,13 @@ def _calc_and_save_weights(
                 axis_sweep=axis_sweep,
             )
 
+            volume_input_upper_left = volume_input[index_input_01]
+
             index_input_01 = _arrays.index_flat(index_input_01, shape_cells_input)
             index_output_01 = _arrays.index_flat(index_output_01, shape_cells_output)
 
-            volume_upper_left = volume_left
             weights.append((index_input_01, index_output_01, volume_upper_left))
+            volume_upper_left = volume_left / volume_input_upper_left
 
     if i_right < (shape_sweep[0] - 1):
         if j_lower >= 0:
@@ -735,11 +833,12 @@ def _calc_and_save_weights(
                 axis_sweep=axis_sweep,
             )
 
+            volume_input_lower_right = volume_input[index_input_10]
+
             index_input_10 = _arrays.index_flat(index_input_10, shape_cells_input)
             index_output_10 = _arrays.index_flat(index_output_10, shape_cells_output)
 
-            volume_lower_right = volume_lower
-            weights.append((index_input_10, index_output_10, volume_lower_right))
+            volume_lower_right = volume_lower / volume_input_lower_right
 
 
 @numba.njit(cache=True)
@@ -747,7 +846,6 @@ def _grid_sweep_static(
     grid_input: tuple[np.ndarray, np.ndarray, np.ndarray],
     grid_output: tuple[np.ndarray, np.ndarray, np.ndarray],
     sweep_input: bool,
-    axis_sweep: tuple[int] | tuple[int, int],
 ) -> tuple[
     tuple[np.ndarray, np.ndarray, np.ndarray],
     tuple[np.ndarray, np.ndarray, np.ndarray],
@@ -766,8 +864,6 @@ def _grid_sweep_static(
         A boolean flag indicating whether to iterate along `grid_input`
         or `grid_output`.
         If :obj:`True`, this function sweeps along `grid_input`.
-    axis_sweep
-        The logical axis of the `grid_sweep` to iterate along.
     """
     if sweep_input:
         grid_sweep = grid_input
@@ -776,14 +872,12 @@ def _grid_sweep_static(
         grid_sweep = grid_output
         grid_static = grid_input
 
-    grid_sweep_x, grid_sweep_y, grid_sweep_z = grid_sweep
-
-    axis = axis_sweep[~0]
-    grid_sweep_x = _arrays.align_axis_right(grid_sweep_x, axis)
-    grid_sweep_y = _arrays.align_axis_right(grid_sweep_y, axis)
-    grid_sweep_z = _arrays.align_axis_right(grid_sweep_z, axis)
-
-    grid_sweep = grid_sweep_x, grid_sweep_y, grid_sweep_z
+    # axis = axis_sweep[~0]
+    # grid_sweep_x = _arrays.align_axis_right(grid_sweep_x, axis)
+    # grid_sweep_y = _arrays.align_axis_right(grid_sweep_y, axis)
+    # grid_sweep_z = _arrays.align_axis_right(grid_sweep_z, axis)
+    #
+    # grid_sweep = grid_sweep_x, grid_sweep_y, grid_sweep_z
 
     return grid_sweep, grid_static
 
@@ -818,9 +912,9 @@ def _index_input_output(
     i_sweep, j_sweep, k_sweep = index_sweep
 
     if axis == axis_x:
-        index_sweep = j_sweep, k_sweep, i_sweep
-    elif axis == axis_y:
         index_sweep = k_sweep, i_sweep, j_sweep
+    elif axis == axis_y:
+        index_sweep = j_sweep, k_sweep, i_sweep
 
     if sweep_input:
         index_input = index_sweep

From 2ee6087ad5c73c19045a9f6b67ee2eb147758ff9 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Wed, 14 May 2025 06:53:16 -0600
Subject: [PATCH 19/26] fixe bugs introduced by previous commit

---
 .../_weights_conservative_3d.py                      | 12 +++---------
 1 file changed, 3 insertions(+), 9 deletions(-)

diff --git a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
index dc6b5e2..bcd3f46 100644
--- a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
+++ b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
@@ -799,8 +799,8 @@ def _calc_and_save_weights(
             index_input_00 = _arrays.index_flat(index_input_00, shape_cells_input)
             index_output_00 = _arrays.index_flat(index_output_00, shape_cells_output)
 
-            weights.append((index_input_00, index_output_00, volume_lower_left))
             volume_lower_left = (-volume_lower - volume_left) / volume_input_lower_left
+            weights.append((index_input_00, index_output_00, volume_lower_left))
 
         if j_upper < (shape_sweep[1] - 1):
 
@@ -818,8 +818,8 @@ def _calc_and_save_weights(
             index_input_01 = _arrays.index_flat(index_input_01, shape_cells_input)
             index_output_01 = _arrays.index_flat(index_output_01, shape_cells_output)
 
-            weights.append((index_input_01, index_output_01, volume_upper_left))
             volume_upper_left = volume_left / volume_input_upper_left
+            weights.append((index_input_01, index_output_01, volume_upper_left))
 
     if i_right < (shape_sweep[0] - 1):
         if j_lower >= 0:
@@ -839,6 +839,7 @@ def _calc_and_save_weights(
             index_output_10 = _arrays.index_flat(index_output_10, shape_cells_output)
 
             volume_lower_right = volume_lower / volume_input_lower_right
+            weights.append((index_input_10, index_output_10, volume_lower_right))
 
 
 @numba.njit(cache=True)
@@ -872,13 +873,6 @@ def _grid_sweep_static(
         grid_sweep = grid_output
         grid_static = grid_input
 
-    # axis = axis_sweep[~0]
-    # grid_sweep_x = _arrays.align_axis_right(grid_sweep_x, axis)
-    # grid_sweep_y = _arrays.align_axis_right(grid_sweep_y, axis)
-    # grid_sweep_z = _arrays.align_axis_right(grid_sweep_z, axis)
-    #
-    # grid_sweep = grid_sweep_x, grid_sweep_y, grid_sweep_z
-
     return grid_sweep, grid_static
 
 

From e0a82849121e5b97636509a7c1b6c055c8c88b06 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Mon, 2 Feb 2026 22:46:27 -0700
Subject: [PATCH 20/26] First working version

---
 .../_weights_conservative_3d/_arrays.py       |  8 +++++++-
 .../_weights_conservative_3d/_grids.py        |  4 +++-
 .../_weights_conservative_3d/_intercepts.py   | 19 ++++++++-----------
 .../_intercepts_test.py                       | 17 ++++++++++++++++-
 4 files changed, 34 insertions(+), 14 deletions(-)

diff --git a/regridding/_weights/_weights_conservative_3d/_arrays.py b/regridding/_weights/_weights_conservative_3d/_arrays.py
index 03c973a..5f19bf4 100644
--- a/regridding/_weights/_weights_conservative_3d/_arrays.py
+++ b/regridding/_weights/_weights_conservative_3d/_arrays.py
@@ -23,6 +23,12 @@
 
 axes = (axis_x, axis_y, axis_z)
 
+vector_unit = (
+    (1, 0, 0),
+    (0, 1, 0),
+    (0, 0, 1),
+)
+
 
 @numba.njit(cache=True)
 def index_in_bounds(
@@ -118,7 +124,7 @@ def align_axis_right(
     axis: int,
 ) -> np.ndarray:
     """
-    Roll all the axes of a 3D array to the right until `axis` is the last axis..
+    Roll all the axes of a 3D array to the right until `axis` is the last axis.
 
     This function is needed to permute the axes in such a way to retain
     their right-handedness.
diff --git a/regridding/_weights/_weights_conservative_3d/_grids.py b/regridding/_weights/_weights_conservative_3d/_grids.py
index ca16d46..49c673c 100644
--- a/regridding/_weights/_weights_conservative_3d/_grids.py
+++ b/regridding/_weights/_weights_conservative_3d/_grids.py
@@ -2,7 +2,7 @@
 Utilities for inspecting and searching logically-rectangular grids of
 coordinates.
 """
-
+import sys
 import numpy as np
 import numba
 import regridding as rg
@@ -474,3 +474,5 @@ def index_of_point_brute(
                     polyhedron=polyhedron,
                 ):
                     return index
+
+    return sys.maxsize, sys.maxsize, sys.maxsize
diff --git a/regridding/_weights/_weights_conservative_3d/_intercepts.py b/regridding/_weights/_weights_conservative_3d/_intercepts.py
index 3642dbf..cdd0ab6 100644
--- a/regridding/_weights/_weights_conservative_3d/_intercepts.py
+++ b/regridding/_weights/_weights_conservative_3d/_intercepts.py
@@ -236,13 +236,13 @@ def sweep(
 
         intercepts_a = intercepts[a]
 
-        normal_a = [1 if axis == a else 0 for axis in range(len(intercepts))]
+        offset_a = [1 if axis == a else 0 for axis in range(len(intercepts))]
 
         for b in range(len(intercepts_a)):
 
             intercepts_ab = intercepts_a[b]
 
-            normal_b = [1 if axis == b else 0 for axis in range(len(intercepts))]
+            offset_b = [1 if axis == b else 0 for axis in range(len(intercepts))]
 
             for i in range(len(intercepts_ab)):
 
@@ -267,10 +267,10 @@ def sweep(
                         i1_input = _arrays.index_3d(i1_input, shape_input)
                         i1_output = _arrays.index_3d(i1_output, shape_output)
 
-                        i0_input_lower = rg.math.difference_3d(i0_input, normal_a)
+                        i0_input_lower = rg.math.difference_3d(i0_input, offset_a)
                         i0_input_upper = i0_input
 
-                        i0_output_lower = rg.math.difference_3d(i0_output, normal_b)
+                        i0_output_lower = rg.math.difference_3d(i0_output, offset_b)
                         i0_output_upper = i0_output
 
                         apex_input = (
@@ -285,10 +285,7 @@ def sweep(
                         )
 
                         vol_input = rg.geometry.volume_tetrahedron(apex_input, p0, p1)
-                        vol_output = rg.geometry.volume_tetrahedron(apex_output, p0, p1)
-
-                        volume_input_lower = volume_input[i0_input_lower]
-                        volume_input_upper = volume_input[i0_output_upper]
+                        vol_output = -rg.geometry.volume_tetrahedron(apex_output, p0, p1)
 
                         input_lower_in_bounds = _arrays.index_in_bounds(
                             index=i0_input_lower,
@@ -332,7 +329,7 @@ def sweep(
                                 weight_lower_lower = (
                                     n0_input_lower,
                                     n0_output_lower,
-                                    (-vol_input - vol_output) / volume_input_lower,
+                                    (vol_input + vol_output) / volume_input_lower,
                                 )
                                 weights.append(weight_lower_lower)
 
@@ -340,7 +337,7 @@ def sweep(
                                 weight_lower_upper = (
                                     n0_input_lower,
                                     n0_output_upper,
-                                    (-vol_input + vol_output) / volume_input_lower,
+                                    (-vol_input - vol_output) / volume_input_lower,
                                 )
                                 weights.append(weight_lower_upper)
 
@@ -352,7 +349,7 @@ def sweep(
                                 weight_upper_lower = (
                                     n0_input_upper,
                                     n0_output_lower,
-                                    (vol_input - vol_output) / volume_input_upper,
+                                    (-vol_input - vol_output) / volume_input_upper,
                                 )
                                 weights.append(weight_upper_lower)
 
diff --git a/regridding/_weights/_weights_conservative_3d/_intercepts_test.py b/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
index a28932c..18268b8 100644
--- a/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
+++ b/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
@@ -48,6 +48,21 @@ def test_insert_intercept(
 @pytest.mark.parametrize(
     argnames="intercepts,intercept_new,result_expected",
     argvalues=[
+        (
+            numba.typed.List([(0, 0, (0.9, 0, 0)), (0, 1, (1.2, 0, 0))]),
+            (0, 2, (1.1, 0, 0)),
+            1,
+        ),
+        (
+            numba.typed.List([(0, 0, (0, 0.9, 0)), (0, 1, (0, 1.2, 0))]),
+            (0, 2, (0, 1.1, 0)),
+            1,
+        ),
+        (
+            numba.typed.List([(0, 0, (0, 0, 0.9)), (0, 1, (0, 0, 1.2))]),
+            (0, 2, (0, 0, 1.1)),
+            1,
+        ),
         (
             intercepts1,
             (0, 0, (-1, -1, -1)),
@@ -88,7 +103,7 @@ def test_insert_intercept(
             (0, 0, (-2.5, 2.25, -2.75)),
             3,
         ),
-    ]
+    ],
 )
 def test_bisect_intercepts(
     intercepts: numba.typed.List[tuple[int, int, tuple[float, float, float]]],

From 7cb58768474335606c9ec28e5e350a5e789a96ac Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Mon, 2 Feb 2026 22:56:32 -0700
Subject: [PATCH 21/26] actual first working version

---
 .../_weights_conservative_3d.py               | 333 ++++++++++++------
 1 file changed, 221 insertions(+), 112 deletions(-)

diff --git a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
index bcd3f46..cec2d33 100644
--- a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
+++ b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
@@ -11,9 +11,9 @@
 from ._arrays import axis_x, axis_y, axis_z
 
 axes = (
-    (axis_x, ),
-    (axis_y, ),
-    (axis_z, ),
+    (axis_x,),
+    (axis_y,),
+    (axis_z,),
     (axis_x, axis_y),
     (axis_y, axis_z),
     (axis_z, axis_x),
@@ -136,8 +136,7 @@ def _sweep_grid(
 
     _, boundary_static = _grids.grid_boundary(grid_static)
 
-    for axis in axes:
-    # for axis in literal_unroll(axes):
+    for axis in literal_unroll(axes):
 
         axis_last = axis[~0]
 
@@ -160,7 +159,10 @@ def _sweep_grid(
         )
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(
+    cache=True,
+    # parallel=True,
+)
 def _sweep_along_axis(
     grid_sweep: tuple[np.ndarray, np.ndarray, np.ndarray],
     grid_static: tuple[np.ndarray, np.ndarray, np.ndarray],
@@ -220,13 +222,6 @@ def _sweep_along_axis(
         The logical axis of the sweep grid to iterate along.
     """
 
-    grid_sweep, grid_static = _grid_sweep_static(
-        grid_input=grid_input,
-        grid_output=grid_output,
-        sweep_input=sweep_input,
-        axis_sweep=axis_sweep,
-    )
-
     x_sweep, y_sweep, z_sweep = grid_sweep
     x_static, y_static, z_static = grid_static
 
@@ -309,6 +304,7 @@ def _sweep_along_axis(
                 )
 
                 line = point_1, point_2
+
                 if sweep_is_outside_static:
 
                     line, index_sweep, index_static = _step_outside_static(
@@ -316,8 +312,8 @@ def _sweep_along_axis(
                         index_sweep=index_sweep,
                         index_static=index_static,
                         step_index=step_index,
+                        grid_sweep=grid_sweep,
                         grid_static=grid_static,
-                        bbox_static=bbox_static,
                         shape_cells_input=shape_cells_input,
                         shape_cells_output=shape_cells_output,
                         intercepts=intercepts,
@@ -359,7 +355,7 @@ def _sweep_along_axis(
             weights.append(weight_i[w])
 
 
-@numba.njit(cache=True)
+@numba.njit(cache=True, error_model="numpy")
 def _step_outside_static(
     line: tuple[
         tuple[float, float, float],
@@ -368,8 +364,8 @@ def _step_outside_static(
     index_sweep: tuple[int, int, int],
     index_static: tuple[int, int, int],
     step_index: tuple[int, int, int],
+    grid_sweep: tuple[np.ndarray, np.ndarray, np.ndarray],
     grid_static: tuple[np.ndarray, np.ndarray, np.ndarray],
-    bbox_static: tuple[tuple[float, float, float], tuple[float, float, float]],
     shape_cells_input: tuple[int, int, int],
     shape_cells_output: tuple[int, int, int],
     intercepts: numba.typed.List,
@@ -400,6 +396,9 @@ def _step_outside_static(
         This is assumed to be a valid, positive index.
     step_index
         A step size in index space on the sweep grid.
+    grid_sweep
+        The vertices of the sweep grid.
+        The last axis of this grid must be the sweep axis.
     grid_static
         The vertices of the static grid.
     bbox_static
@@ -430,105 +429,105 @@ def _step_outside_static(
 
     shape_cells_static = _grids.shape_centers(x.shape)
 
-    if rg.geometry.point_is_inside_box_3d(p2, bbox_static):
+    for axis in _arrays.axes:
 
-        for axis in _arrays.axes:
+        x_axis = _arrays.align_axis_right(x, axis)
+        y_axis = _arrays.align_axis_right(y, axis)
+        z_axis = _arrays.align_axis_right(z, axis)
 
-            x_axis = _arrays.align_axis_right(x, axis)
-            y_axis = _arrays.align_axis_right(y, axis)
-            z_axis = _arrays.align_axis_right(z, axis)
+        for direction_face in (-1, 1):
 
-            for direction_face in (-1, 1):
-
-                if direction_face > 0:
-                    k_face = ~0
-                else:
-                    k_face = 0
+            if direction_face > 0:
+                k_face = ~0
+            else:
+                k_face = 0
 
-                x_face = x_axis[..., k_face]
-                y_face = y_axis[..., k_face]
-                z_face = z_axis[..., k_face]
+            x_face = x_axis[..., k_face]
+            y_face = y_axis[..., k_face]
+            z_face = z_axis[..., k_face]
 
-                shape_face = x_face.shape
+            shape_face = x_face.shape
 
-                num_i, num_j = shape_face
+            num_i, num_j = shape_face
 
-                for i0 in range(num_i - 1):
-                    for j0 in range(num_j - 1):
+            for i0 in range(num_i - 1):
+                for j0 in range(num_j - 1):
 
-                        i1 = i0 + 1
-                        j1 = j0 + 1
+                    i1 = i0 + 1
+                    j1 = j0 + 1
 
-                        for t in (-1, 1):
+                    for t in (-1, 1):
 
-                            if t > 0:
-                                ind = (
-                                    (i0, j0),
-                                    (i1, j0),
-                                    (i1, j1),
-                                )
-                            else:
-                                ind = (
-                                    (i1, j1),
-                                    (i0, j1),
-                                    (i0, j0),
-                                )
+                        if t > 0:
+                            ind = (
+                                (i0, j0),
+                                (i1, j0),
+                                (i1, j1),
+                            )
+                        else:
+                            ind = (
+                                (i1, j1),
+                                (i0, j1),
+                                (i0, j0),
+                            )
 
-                            if direction_face < 0:
-                                ind = ind[::-1]
+                        if direction_face < 0:
+                            ind = ind[::-1]
 
-                            v0, v1, v2 = ind
+                        v0, v1, v2 = ind
 
-                            triangle = (
-                                (x_face[v0], y_face[v0], z_face[v0]),
-                                (x_face[v1], y_face[v1], z_face[v1]),
-                                (x_face[v2], y_face[v2], z_face[v2]),
-                            )
+                        triangle = (
+                            (x_face[v0], y_face[v0], z_face[v0]),
+                            (x_face[v1], y_face[v1], z_face[v1]),
+                            (x_face[v2], y_face[v2], z_face[v2]),
+                        )
 
-                            tuv = rg.geometry.line_triangle_intersection_parameters(
-                                line=line,
-                                triangle=triangle
-                            )
+                        tuv = rg.geometry.line_triangle_intersection_parameters(
+                            line=line,
+                            triangle=triangle
+                        )
 
-                            if rg.geometry.line_intersects_triangle(tuv):
+                        if rg.geometry.line_intersects_triangle(tuv):
 
-                                if direction_face > 0:
-                                    index_static = i0, j0, shape_cells_static[axis] - 1
-                                else:
-                                    index_static = i0, j0, 0
+                            if direction_face > 0:
+                                index_static = i0, j0, shape_cells_static[axis] - 1
+                            else:
+                                index_static = i0, j0, 0
 
-                                i, j, k = index_static
+                            i, j, k = index_static
 
-                                if axis == axis_x:
-                                    index_static = k, i, j
-                                elif axis == axis_y:
-                                    index_static = j, k, i
+                            if axis == axis_x:
+                                index_static = k, i, j
+                            elif axis == axis_y:
+                                index_static = j, k, i
 
-                                normal = [1 if ax == axis else 0 for ax in _arrays.axes]
-                                normal = rg.math.multiply_3d(direction_face, normal)
+                            normal = [1 if ax == axis else 0 for ax in _arrays.axes]
+                            normal = rg.math.multiply_3d(direction_face, normal)
 
-                                _, p2 = _calc_and_save_intercept(
-                                    line=line,
-                                    tuv=tuv,
-                                    index_sweep=index_sweep,
-                                    index_static=index_static,
-                                    shape_cells_input=shape_cells_input,
-                                    shape_cells_output=shape_cells_output,
-                                    intercepts=intercepts,
-                                    sweep_input=sweep_input,
-                                    axis_sweep=axis_sweep,
-                                    axis_static=axis,
-                                    normal_static=normal,
-                                )
+                            _, p2 = _calc_and_save_intercept(
+                                line=line,
+                                tuv=tuv,
+                                index_sweep=index_sweep,
+                                index_static=index_static,
+                                grid_sweep=grid_sweep,
+                                grid_static=grid_static,
+                                shape_cells_input=shape_cells_input,
+                                shape_cells_output=shape_cells_output,
+                                intercepts=intercepts,
+                                sweep_input=sweep_input,
+                                axis_sweep=axis_sweep,
+                                axis_static=axis,
+                                offset_static=normal,
+                            )
 
-                                return (p1, p2), index_sweep, index_static
+                            return (p1, p2), index_sweep, index_static
 
     index_sweep = rg.math.sum_3d(index_sweep, step_index)
 
     return (p1, p2), index_sweep, index_static
 
 
-@numba.njit(cache=True)
+@numba.njit(cache=True, error_model="numpy")
 def _step_inside_static(
     line: tuple[
         tuple[float, float, float],
@@ -569,7 +568,7 @@ def _step_inside_static(
         The index of the current vertex in the sweep grid
         This is assumed to be a valid, positive index.
     index_static
-        The current index in the static grid.
+        The index of the current cell in the static grid.
         This is assumed to be a valid, positive index.
     step_index
         A step size in index space on the sweep grid.
@@ -630,13 +629,15 @@ def _step_inside_static(
                     tuv=tuv,
                     index_sweep=index_sweep,
                     index_static=index_static,
+                    grid_sweep=grid_sweep,
+                    grid_static=grid_static,
                     shape_cells_input=shape_cells_input,
                     shape_cells_output=shape_cells_output,
                     intercepts=intercepts,
                     sweep_input=sweep_input,
                     axis_sweep=axis_sweep,
                     axis_static=_grids.cell_axes[t],
-                    normal_static=_grids.cell_normals[t],
+                    offset_static=_grids.cell_normals[t],
                 )
 
                 break
@@ -673,13 +674,15 @@ def _calc_and_save_intercept(
     tuv: tuple[float, float, float],
     index_sweep: tuple[int, int, int],
     index_static: tuple[int, int, int],
+    grid_sweep: tuple[np.ndarray, np.ndarray, np.ndarray],
+    grid_static: tuple[np.ndarray, np.ndarray, np.ndarray],
     shape_cells_input: tuple[int, int, int],
     shape_cells_output: tuple[int, int, int],
     intercepts: numba.typed.List,
     sweep_input: bool,
     axis_sweep: tuple[int] | tuple[int, int],
     axis_static: int,
-    normal_static: tuple[int, int, int],
+    offset_static: tuple[int, int, int],
 ) -> tuple[
     tuple[int, int, int],
     tuple[float, float, float],
@@ -690,21 +693,109 @@ def _calc_and_save_intercept(
         tuv=tuv,
     )
 
+    point_sweep, point_sweep_new = line
+
+    x_sweep, y_sweep, z_sweep = grid_sweep
+    x_static, y_static, z_static = grid_static
+
+    vertex_static = (
+        x_static[index_static],
+        y_static[index_static],
+        z_static[index_static],
+    )
+
     shape_input = [n + 1 for n in shape_cells_input]
     shape_output = [n + 1 for n in shape_cells_output]
 
-    index_static_new = rg.math.sum_3d(index_static, normal_static)
-    # print(f"{index_static=}")
-    # print(f"{index_static_new=}")
+    index_static_new = rg.math.sum_3d(index_static, offset_static)
 
-    index_static_intercept = index_static
+    index_static_normal = index_static_new
+    if 0 <= index_static_normal[axis_static] < shape_output[axis_static]:
+        vertex_static_normal = (
+            x_static[index_static_normal],
+            y_static[index_static_normal],
+            z_static[index_static_normal],
+        )
+        normal_static = rg.math.difference_3d(vertex_static_normal, vertex_static)
+    else:
+        index_static_normal = rg.math.difference_3d(index_static, offset_static)
+        vertex_static_normal = (
+            x_static[index_static_normal],
+            y_static[index_static_normal],
+            z_static[index_static_normal],
+        )
+        normal_static = rg.math.difference_3d(vertex_static, vertex_static_normal)
+
+    index_apex_static = index_static
     if index_static_new[axis_static] > index_static[axis_static]:
-        index_static_intercept = index_static_new
+        index_apex_static = index_static_new
+    else:
+        normal_static = rg.math.negate_3d(normal_static)
 
     axes_orthogonal = [ax for ax in _arrays.axes if ax not in axis_sweep]
 
     for axis_orthogonal in axes_orthogonal:
 
+        offset_axis = 2 - axis_sweep[~0]
+
+        axis_orthogonal_apparent = (axis_orthogonal + offset_axis) % 3
+
+        offset_orthogonal = _arrays.vector_unit[axis_orthogonal_apparent]
+
+        index_orthogonal = rg.math.difference_3d(index_sweep, offset_orthogonal)
+
+        if index_orthogonal[axis_orthogonal_apparent] >= 0:
+            point_orthogonal = (
+                x_sweep[index_orthogonal],
+                y_sweep[index_orthogonal],
+                z_sweep[index_orthogonal],
+            )
+            normal_sweep = rg.math.difference_3d(point_sweep, point_orthogonal)
+        else:
+            index_orthogonal = rg.math.sum_3d(index_sweep, offset_orthogonal)
+            point_orthogonal = (
+                x_sweep[index_orthogonal],
+                y_sweep[index_orthogonal],
+                z_sweep[index_orthogonal],
+            )
+            normal_sweep = rg.math.difference_3d(point_orthogonal, point_sweep)
+
+        if sweep_input:
+            normal_input = normal_sweep
+            normal_output = normal_static
+        else:
+            normal_input = normal_static
+            normal_output = normal_sweep
+
+        direction = rg.math.cross_3d(normal_input, normal_output)
+
+        index_apex_sweep = index_sweep
+        if len(axis_sweep) == 1:
+
+            axis_parallel_actual = [
+                ax for ax in _arrays.axes
+                if ax is not axis_orthogonal_apparent and ax is not _arrays.axis_z
+            ][0]
+
+            offset_parallel = _arrays.vector_unit[axis_parallel_actual]
+
+            index_parallel = rg.math.difference_3d(index_sweep, offset_parallel)
+
+            if index_parallel[axis_parallel_actual] >= 0:
+
+                point_parallel = (
+                    x_sweep[index_parallel],
+                    y_sweep[index_parallel],
+                    z_sweep[index_parallel],
+                )
+
+                parallel = rg.math.difference_3d(point_sweep, point_parallel)
+
+                sgn = rg.math.dot_3d(parallel, direction)
+
+                if sgn < 0:
+                    index_apex_sweep = index_parallel
+
         if sweep_input:
             axis_input = axis_orthogonal
             axis_output = axis_static
@@ -713,8 +804,8 @@ def _calc_and_save_intercept(
             axis_output = axis_orthogonal
 
         index_input, index_output = _index_input_output(
-            index_sweep=index_sweep,
-            index_static=index_static_intercept,
+            index_sweep=index_apex_sweep,
+            index_static=index_apex_static,
             sweep_input=sweep_input,
             axis_sweep=axis_sweep,
         )
@@ -722,15 +813,34 @@ def _calc_and_save_intercept(
         i_input = index_input[axis_input]
         i_output = index_output[axis_output]
 
-        _intercepts.insert_intercept(
-            intercepts=intercepts[axis_input][axis_output][i_input][i_output],
-            intercept_new=(
-                _arrays.index_flat(index_input, shape_input),
-                _arrays.index_flat(index_output, shape_output),
-                intercept,
-            ),
+        intercept_new = (
+            _arrays.index_flat(index_input, shape_input),
+            _arrays.index_flat(index_output, shape_output),
+            intercept,
         )
 
+        intercepts_input_output = intercepts[axis_input][axis_output][i_input][i_output]
+        num_intercept = len(intercepts_input_output)
+        if num_intercept == 0:
+            index = 0
+
+        elif num_intercept == 1:
+            _, _, intercept_0 = intercepts_input_output[0]
+            direction_intercepts = rg.math.difference_3d(intercept, intercept_0)
+            sgn_intercepts = rg.math.dot_3d(direction_intercepts, direction)
+            if sgn_intercepts > 0:
+                index = 1
+            else:
+                index = 0
+
+        else:
+            index = _intercepts._bisect_intercepts(
+                intercepts=intercepts_input_output,
+                intercept_new=intercept_new,
+            )
+
+        intercepts_input_output.insert(index, intercept_new)
+
     return index_static_new, intercept
 
 
@@ -779,8 +889,8 @@ def _calc_and_save_weights(
 
     p1, p2 = line
 
-    volume_left = rg.geometry.volume_tetrahedron(p0_left, p1, p2)
-    volume_lower = -rg.geometry.volume_tetrahedron(p0_lower, p1, p1)
+    volume_left = -rg.geometry.volume_tetrahedron(p0_left, p1, p2)
+    volume_lower = rg.geometry.volume_tetrahedron(p0_lower, p1, p2)
 
     if i_left >= 0:
         if j_lower >= 0:
@@ -799,7 +909,7 @@ def _calc_and_save_weights(
             index_input_00 = _arrays.index_flat(index_input_00, shape_cells_input)
             index_output_00 = _arrays.index_flat(index_output_00, shape_cells_output)
 
-            volume_lower_left = (-volume_lower - volume_left) / volume_input_lower_left
+            volume_lower_left = (volume_lower + volume_left) / volume_input_lower_left
             weights.append((index_input_00, index_output_00, volume_lower_left))
 
         if j_upper < (shape_sweep[1] - 1):
@@ -818,7 +928,7 @@ def _calc_and_save_weights(
             index_input_01 = _arrays.index_flat(index_input_01, shape_cells_input)
             index_output_01 = _arrays.index_flat(index_output_01, shape_cells_output)
 
-            volume_upper_left = volume_left / volume_input_upper_left
+            volume_upper_left = -volume_left / volume_input_upper_left
             weights.append((index_input_01, index_output_01, volume_upper_left))
 
     if i_right < (shape_sweep[0] - 1):
@@ -838,7 +948,7 @@ def _calc_and_save_weights(
             index_input_10 = _arrays.index_flat(index_input_10, shape_cells_input)
             index_output_10 = _arrays.index_flat(index_output_10, shape_cells_output)
 
-            volume_lower_right = volume_lower / volume_input_lower_right
+            volume_lower_right = -volume_lower / volume_input_lower_right
             weights.append((index_input_10, index_output_10, volume_lower_right))
 
 
@@ -882,7 +992,7 @@ def _index_input_output(
     index_static: tuple[int, int, int],
     sweep_input: bool,
     axis_sweep: tuple[int] | tuple[int, int],
-):
+) -> tuple[tuple[int, int, int], tuple[int, int, int]]:
     """
     Convert indices on the sweep and static grids into indices on the
     input/output grids.
@@ -918,4 +1028,3 @@ def _index_input_output(
         index_output = index_sweep
 
     return index_input, index_output
-

From b0209a7f4434ae058401ad69df7b59b36058c47b Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Tue, 3 Feb 2026 08:08:15 -0700
Subject: [PATCH 22/26] fix test

---
 regridding/_weights/_weights_conservative_3d/_grids_test.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/regridding/_weights/_weights_conservative_3d/_grids_test.py b/regridding/_weights/_weights_conservative_3d/_grids_test.py
index cede370..71ed33f 100644
--- a/regridding/_weights/_weights_conservative_3d/_grids_test.py
+++ b/regridding/_weights/_weights_conservative_3d/_grids_test.py
@@ -74,7 +74,7 @@ def test_index_of_point_brute(
     grid: tuple[np.ndarray, np.ndarray, np.ndarray],
     result_expected: np.ndarray,
 ):
-    result = _grids._index_of_point_brute(
+    result = _grids.index_of_point_brute(
         point=point,
         grid=grid,
     )

From 5e7fc16c8c294d4eb6f71a006786241146fafa2c Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Tue, 3 Feb 2026 10:39:06 -0700
Subject: [PATCH 23/26] black

---
 .../_weights_conservative_3d/__init__.py      |  1 +
 .../_weights_conservative_3d/_arrays.py       |  1 -
 .../_weights_conservative_3d/_arrays_test.py  | 14 +++----
 .../_weights_conservative_3d/_grids.py        |  5 +--
 .../_weights_conservative_3d/_grids_test.py   |  8 ++--
 .../_weights_conservative_3d/_intercepts.py   | 12 +++++-
 .../_intercepts_test.py                       | 40 +++++++++++--------
 .../_weights_conservative_3d.py               |  8 ++--
 .../_weights_conservative_3d_test.py          |  7 +---
 9 files changed, 50 insertions(+), 46 deletions(-)

diff --git a/regridding/_weights/_weights_conservative_3d/__init__.py b/regridding/_weights/_weights_conservative_3d/__init__.py
index 8a91a8b..6c5336f 100644
--- a/regridding/_weights/_weights_conservative_3d/__init__.py
+++ b/regridding/_weights/_weights_conservative_3d/__init__.py
@@ -1,6 +1,7 @@
 """
 A subpackage for performing 3D conservative resampling
 """
+
 from ._weights_conservative_3d import weights_conservative_3d
 
 __all__ = [
diff --git a/regridding/_weights/_weights_conservative_3d/_arrays.py b/regridding/_weights/_weights_conservative_3d/_arrays.py
index 5f19bf4..07be0fc 100644
--- a/regridding/_weights/_weights_conservative_3d/_arrays.py
+++ b/regridding/_weights/_weights_conservative_3d/_arrays.py
@@ -14,7 +14,6 @@
     "index_flat",
     "index_3d",
     "align_axis_right",
-
 ]
 
 axis_x = 0
diff --git a/regridding/_weights/_weights_conservative_3d/_arrays_test.py b/regridding/_weights/_weights_conservative_3d/_arrays_test.py
index 6f76f7a..ad5f761 100644
--- a/regridding/_weights/_weights_conservative_3d/_arrays_test.py
+++ b/regridding/_weights/_weights_conservative_3d/_arrays_test.py
@@ -26,7 +26,7 @@
             (11, 12, 13),
             False,
         ),
-    ]
+    ],
 )
 def test_index_in_bounds(
     index: tuple[int, int, int],
@@ -43,13 +43,13 @@ def test_index_in_bounds(
     argvalues=[
         (1, 1, 1),
         (5, 6, 7),
-    ]
+    ],
 )
 @pytest.mark.parametrize(
     argnames="shape",
     argvalues=[
         (11, 12, 13),
-    ]
+    ],
 )
 def test_index_flat(
     index: tuple[int, int, int],
@@ -63,17 +63,13 @@ def test_index_flat(
 
 @pytest.mark.parametrize(
     argnames="index",
-    argvalues=[
-        0,
-        11,
-        25
-    ]
+    argvalues=[0, 11, 25],
 )
 @pytest.mark.parametrize(
     argnames="shape",
     argvalues=[
         (11, 12, 13),
-    ]
+    ],
 )
 def test_index_3d(
     index: int,
diff --git a/regridding/_weights/_weights_conservative_3d/_grids.py b/regridding/_weights/_weights_conservative_3d/_grids.py
index 49c673c..a0297cb 100644
--- a/regridding/_weights/_weights_conservative_3d/_grids.py
+++ b/regridding/_weights/_weights_conservative_3d/_grids.py
@@ -2,6 +2,7 @@
 Utilities for inspecting and searching logically-rectangular grids of
 coordinates.
 """
+
 import sys
 import numpy as np
 import numba
@@ -243,9 +244,7 @@ def cell_boundary(
 def grid_boundary(
     grid: tuple[np.ndarray, np.ndarray, np.ndarray],
 ) -> tuple[
-    numba.typed.List[
-        tuple[int, int, int],
-    ],
+    numba.typed.List[tuple[int, int, int]],
     numba.typed.List[
         tuple[
             tuple[float, float, float],
diff --git a/regridding/_weights/_weights_conservative_3d/_grids_test.py b/regridding/_weights/_weights_conservative_3d/_grids_test.py
index 71ed33f..b09d81b 100644
--- a/regridding/_weights/_weights_conservative_3d/_grids_test.py
+++ b/regridding/_weights/_weights_conservative_3d/_grids_test.py
@@ -13,7 +13,7 @@
                 np.arange(2),
                 indexing="ij",
             ),
-            np.ones((1, 1, 1))
+            np.ones((1, 1, 1)),
         ),
         (
             np.meshgrid(
@@ -22,7 +22,7 @@
                 np.arange(5),
                 indexing="ij",
             ),
-            np.ones((2, 3, 4))
+            np.ones((2, 3, 4)),
         ),
         (
             np.meshgrid(
@@ -31,9 +31,9 @@
                 2 * np.arange(5),
                 indexing="ij",
             ),
-            8 * np.ones((2, 3, 4))
+            8 * np.ones((2, 3, 4)),
         ),
-    ]
+    ],
 )
 def test_volume_grid(
     grid: tuple[np.ndarray, np.ndarray, np.ndarray],
diff --git a/regridding/_weights/_weights_conservative_3d/_intercepts.py b/regridding/_weights/_weights_conservative_3d/_intercepts.py
index cdd0ab6..0e5e69f 100644
--- a/regridding/_weights/_weights_conservative_3d/_intercepts.py
+++ b/regridding/_weights/_weights_conservative_3d/_intercepts.py
@@ -284,8 +284,16 @@ def sweep(
                             z_output[i0_output],
                         )
 
-                        vol_input = rg.geometry.volume_tetrahedron(apex_input, p0, p1)
-                        vol_output = -rg.geometry.volume_tetrahedron(apex_output, p0, p1)
+                        vol_input = rg.geometry.volume_tetrahedron(
+                            apex_input,
+                            p0,
+                            p1,
+                        )
+                        vol_output = -rg.geometry.volume_tetrahedron(
+                            apex_output,
+                            p0,
+                            p1,
+                        )
 
                         input_lower_in_bounds = _arrays.index_in_bounds(
                             index=i0_input_lower,
diff --git a/regridding/_weights/_weights_conservative_3d/_intercepts_test.py b/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
index 18268b8..1f293ac 100644
--- a/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
+++ b/regridding/_weights/_weights_conservative_3d/_intercepts_test.py
@@ -4,32 +4,38 @@
 from . import _intercepts
 
 
-intercepts1 = numba.typed.List([
-    (0, 0, (0, 0, 0)),
-    (0, 1, (1, 1, 1)),
-    (0, 2, (2, 2, 2)),
-])
+intercepts1 = numba.typed.List(
+    [
+        (0, 0, (0, 0, 0)),
+        (0, 1, (1, 1, 1)),
+        (0, 2, (2, 2, 2)),
+    ]
+)
 
-intercepts2 = numba.typed.List([
-    (0, 0, (0, 0, 0)),
-    (0, 1, (-1, 1, -1)),
-    (0, 2, (-2, 2, -2)),
-])
+intercepts2 = numba.typed.List(
+    [
+        (0, 0, (0, 0, 0)),
+        (0, 1, (-1, 1, -1)),
+        (0, 2, (-2, 2, -2)),
+    ]
+)
 
 
 @pytest.mark.parametrize(
     argnames="intercepts,intercept_new",
     argvalues=[
         (
-            numba.typed.List([
-                (0, 0, (0, 0, 0)),
-                (0, 1, (-1, 1, -1)),
-                (0, 2, (-2, 2, -2)),
-            ]),
+            numba.typed.List(
+                [
+                    (0, 0, (0, 0, 0)),
+                    (0, 1, (-1, 1, -1)),
+                    (0, 2, (-2, 2, -2)),
+                ]
+            ),
             (0, 0, (-1, -1, -1)),
         ),
         (
-            _intercepts.empty((11,11,11), (11,11,11))[0][0][0][0],
+            _intercepts.empty((11, 11, 11), (11, 11, 11))[0][0][0][0],
             (0, 0, (-1, -1, -1)),
         ),
     ],
@@ -190,7 +196,7 @@ def _test_sweep_compiled(
 
     weights = numba.typed.List()
     for _ in range(0):
-        weights.append((0, 0, 0.))
+        weights.append((0, 0, 0.0))
 
     _intercepts.sweep(
         intercepts=intercepts,
diff --git a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
index cec2d33..1f436e3 100644
--- a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
+++ b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
@@ -56,7 +56,6 @@ def weights_conservative_3d(
 
     volume_input = _grids.grid_volume(grid_input)
 
-
     intercepts = _intercepts.empty(shape_input, shape_output)
 
     for sweep_input in (False, True):
@@ -239,7 +238,7 @@ def _sweep_along_axis(
     for i in range(shape_sweep_x):
         weight_i = numba.typed.List()
         for _ in range(0):
-            weight_i.append((0, 0, 0.))
+            weight_i.append((0, 0, 0.0))
         weight.append(weight_i)
 
     if len(axis_sweep) == 2:
@@ -484,7 +483,7 @@ def _step_outside_static(
 
                         tuv = rg.geometry.line_triangle_intersection_parameters(
                             line=line,
-                            triangle=triangle
+                            triangle=triangle,
                         )
 
                         if rg.geometry.line_intersects_triangle(tuv):
@@ -773,7 +772,8 @@ def _calc_and_save_intercept(
         if len(axis_sweep) == 1:
 
             axis_parallel_actual = [
-                ax for ax in _arrays.axes
+                ax
+                for ax in _arrays.axes
                 if ax is not axis_orthogonal_apparent and ax is not _arrays.axis_z
             ][0]
 
diff --git a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d_test.py b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d_test.py
index 551fd72..af73249 100644
--- a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d_test.py
+++ b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d_test.py
@@ -26,7 +26,7 @@
                 indexing="ij",
             ),
         )
-    ]
+    ],
 )
 def test_weights_conservative_3d(
     grid_input: tuple[np.ndarray, np.ndarray, np.ndarray],
@@ -36,8 +36,3 @@ def test_weights_conservative_3d(
         grid_input=grid_input,
         grid_output=grid_output,
     )
-
-
-
-
-

From 7298d2c70ac0f91cf80871670fbf560a1faaa793 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Wed, 4 Mar 2026 21:58:32 -0700
Subject: [PATCH 24/26] Fixed a bug where the first point of the "parallel"
 vector was incorrect.

---
 .../_weights_conservative_3d.py                      | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
index 1f436e3..236ba80 100644
--- a/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
+++ b/regridding/_weights/_weights_conservative_3d/_weights_conservative_3d.py
@@ -743,13 +743,19 @@ def _calc_and_save_intercept(
 
         index_orthogonal = rg.math.difference_3d(index_sweep, offset_orthogonal)
 
+        apex_sweep = (
+            x_sweep[index_sweep],
+            y_sweep[index_sweep],
+            z_sweep[index_sweep],
+        )
+
         if index_orthogonal[axis_orthogonal_apparent] >= 0:
             point_orthogonal = (
                 x_sweep[index_orthogonal],
                 y_sweep[index_orthogonal],
                 z_sweep[index_orthogonal],
             )
-            normal_sweep = rg.math.difference_3d(point_sweep, point_orthogonal)
+            normal_sweep = rg.math.difference_3d(apex_sweep, point_orthogonal)
         else:
             index_orthogonal = rg.math.sum_3d(index_sweep, offset_orthogonal)
             point_orthogonal = (
@@ -757,7 +763,7 @@ def _calc_and_save_intercept(
                 y_sweep[index_orthogonal],
                 z_sweep[index_orthogonal],
             )
-            normal_sweep = rg.math.difference_3d(point_orthogonal, point_sweep)
+            normal_sweep = rg.math.difference_3d(point_orthogonal, apex_sweep)
 
         if sweep_input:
             normal_input = normal_sweep
@@ -789,7 +795,7 @@ def _calc_and_save_intercept(
                     z_sweep[index_parallel],
                 )
 
-                parallel = rg.math.difference_3d(point_sweep, point_parallel)
+                parallel = rg.math.difference_3d(apex_sweep, point_parallel)
 
                 sgn = rg.math.dot_3d(parallel, direction)
 

From 2028cc97efe225c994d1d8cfe9aefb241a7e7caa Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Thu, 7 May 2026 13:44:57 -0600
Subject: [PATCH 25/26] indent error

---
 regridding/_weights/_weights_conservative.py | 28 ++++++++++----------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/regridding/_weights/_weights_conservative.py b/regridding/_weights/_weights_conservative.py
index 04f183d..2046653 100644
--- a/regridding/_weights/_weights_conservative.py
+++ b/regridding/_weights/_weights_conservative.py
@@ -129,20 +129,20 @@ def _weights_conservative(
                 )
 
             elif len(axis_input) == 3:
-            x_input, y_input, z_input = coordinates_input
-            x_output, y_output, z_output = coordinates_output
-            weights[index] = weights_conservative_3d(
-                grid_input=(
-                    x_input[index_vertices_input],
-                    y_input[index_vertices_input],
-                    z_input[index_vertices_input],
-                ),
-                grid_output=(
-                    x_output[index_vertices_output],
-                    y_output[index_vertices_output],
-                    z_output[index_vertices_output],
-                ),
-            )
+                x_input, y_input, z_input = coordinates_input
+                x_output, y_output, z_output = coordinates_output
+                weights[index] = weights_conservative_3d(
+                    grid_input=(
+                        x_input[index_vertices_input],
+                        y_input[index_vertices_input],
+                        z_input[index_vertices_input],
+                    ),
+                    grid_output=(
+                        x_output[index_vertices_output],
+                        y_output[index_vertices_output],
+                        z_output[index_vertices_output],
+                    ),
+                )
 
         else:  # pragma: nocover
                 raise NotImplementedError(

From 4c600e9c2d630d2c7dabe6461e65299e8d843e1b Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Thu, 7 May 2026 13:48:27 -0600
Subject: [PATCH 26/26] more indent issues

---
 regridding/_weights/_weights_conservative.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/regridding/_weights/_weights_conservative.py b/regridding/_weights/_weights_conservative.py
index 2046653..d5dcee1 100644
--- a/regridding/_weights/_weights_conservative.py
+++ b/regridding/_weights/_weights_conservative.py
@@ -144,7 +144,7 @@ def _weights_conservative(
                     ),
                 )
 
-        else:  # pragma: nocover
+            else:  # pragma: nocover
                 raise NotImplementedError(
                     "Regridding operations greater than 3D are not supported"
                 )