diff --git a/src/bonsai/bonsai/core/model.py b/src/bonsai/bonsai/core/model.py
index 7f4237fb45..fe289cbda1 100644
--- a/src/bonsai/bonsai/core/model.py
+++ b/src/bonsai/bonsai/core/model.py
@@ -21,7 +21,7 @@
 from __future__ import annotations
 
 import math
-from typing import TYPE_CHECKING, Literal, Optional
+from typing import TYPE_CHECKING, Any, Literal, Optional
 
 if TYPE_CHECKING:
     import bpy
@@ -34,6 +34,24 @@ if TYPE_CHECKING:
     OffsetType = Literal["CENTER", "EXTERIOR", "INTERIOR"]
 
 
+# Arc sample count for fillet preview polylines. 24 samples produces a visually
+# smooth arc at common viewport scales without bloating the GPU batch.
+FILLET_DEFAULT_ARC_RESOLUTION = 24
+# Dot-product floor for treating two wall-axis segments as parallel — below
+# this the projected intersection is too sensitive to floating-point noise
+# to be useful as a junction apex. Calibrated to ~2° from parallel.
+PARALLEL_DOT_THRESHOLD = 0.9994
+# Perpendicular distance (SI metres) under which two parallel wall axes are
+# considered to share the same infinite line. Calibrated to absorb sub-50mm
+# placement drift between authored-joined walls without merging genuinely
+# offset parallel walls.
+COLLINEAR_LINE_TOLERANCE = 0.05
+# Default proximity (SI metres) for classifying a layer offset against the
+# canonical EXTERIOR / CENTER / INTERIOR baselines. Tight enough that ordinary
+# millimetre-scale modelling intent always falls into the nearest baseline.
+BASELINE_OFFSET_TOLERANCE = 0.001
+
+
 def unjoin_walls(
     ifc: type[tool.Ifc],
     blender: type[tool.Blender],
@@ -179,16 +197,16 @@ class RequireLayeredElement(Exception):
 
 
 # --- Wall geometry math (pure) ------------------------------------------------
-# Tuple in / tuple out so these helpers run under ``pytest test/core/`` without
-# ``bpy`` or ``mathutils``. Callers convert ``mathutils.Vector`` at the boundary.
+# Tuple in / tuple out so these helpers run without ``bpy`` or ``mathutils``.
+# Callers convert ``mathutils.Vector`` at the boundary.
 
 
-def baseline_from_offset(offset: float, thickness: float, tolerance: float = 0.001) -> str:
+def baseline_from_offset(offset: float, thickness: float, tolerance: float = BASELINE_OFFSET_TOLERANCE) -> str:
     """Classify a numeric layer offset as EXTERIOR / CENTER / INTERIOR.
 
-    Mirrors the math in ``tool.Model.offset_wall`` for both POSITIVE and NEGATIVE
-    direction_sense walls. Returns the closest canonical baseline; falls back to
-    ``"CENTER"`` when nothing is within ``tolerance``."""
+    Handles both POSITIVE and NEGATIVE direction_sense walls. Returns the
+    closest canonical baseline; falls back to ``"CENTER"`` when nothing is
+    within ``tolerance``."""
     candidates = (
         ("EXTERIOR", 0.0),
         ("CENTER", -thickness / 2),
@@ -211,7 +229,7 @@ def project_axis_intersection(
     Each segment is a pair of 3-tuples. Returns the intersection as a 3-tuple
     (Z is the average of the four input Zs, for visual placement) or ``None`` if
     the segments are parallel within ``parallel_threshold`` (a dot-product magnitude
-    threshold — e.g. ``cos(2°) ≈ 0.9994`` treats walls within 2° of parallel as parallel)."""
+    threshold — see ``PARALLEL_DOT_THRESHOLD`` for the calibrated value)."""
     p1, p2 = seg_a
     p3, p4 = seg_b
     d1x, d1y = p2[0] - p1[0], p2[1] - p1[1]
@@ -233,21 +251,100 @@ def project_axis_intersection(
     return (ix, iy, iz)
 
 
-def displacement_from_x_angle(height: float, x_angle: float) -> float:
-    """Top-edge horizontal displacement for a wall of given vertical ``height`` and
-    slope ``x_angle`` (radians). Drives the slope dimension gizmo's display value.
+def opening_is_past_cut(min_t: float, cut_percentage: float) -> bool:
+    """True when the opening's near edge sits past the cut on the t axis.
 
-    Inverse of :func:`x_angle_from_displacement`."""
+    Strict inequality is load-bearing: a boundary touch or NaN keeps the
+    opening on both walls — the safe default when extent resolution fails."""
+    return min_t > cut_percentage
+
+
+def opening_is_before_cut(max_t: float, cut_percentage: float) -> bool:
+    """True when the opening's far edge sits before the cut on the t axis."""
+    return max_t < cut_percentage
+
+
+def opening_straddles_cut(min_t: float, max_t: float, cut_percentage: float) -> bool:
+    """True when the opening's extent crosses the cut on the t axis."""
+    return min_t < cut_percentage < max_t
+
+
+WallJoinState = Literal["joined", "collinear", "intersect", "none"]
+
+
+def classify_wall_join_state(
+    seg_a: tuple[tuple[float, float, float], tuple[float, float, float]],
+    seg_b: tuple[tuple[float, float, float], tuple[float, float, float]],
+    are_joined: bool,
+    parallel_threshold: float,
+    collinear_tolerance: float,
+) -> tuple[WallJoinState, Optional[tuple[float, float, float]]]:
+    """Classify a wall pair's geometric state — ``(state, intersection)``.
+
+    Priority: ``"joined"`` (caller-supplied flag) → ``"collinear"`` →
+    ``"intersect"`` (projected point returned) → ``"none"`` (parallel,
+    non-collinear)."""
+    if are_joined:
+        return "joined", None
+    if are_axes_collinear(seg_a, seg_b, parallel_threshold, collinear_tolerance):
+        return "collinear", None
+    intersection = project_axis_intersection(seg_a, seg_b, parallel_threshold)
+    if intersection is None:
+        return "none", None
+    return "intersect", intersection
+
+
+def wall_join_preview_lines(
+    seg_a: tuple[tuple[float, float, float], tuple[float, float, float]],
+    seg_b: tuple[tuple[float, float, float], tuple[float, float, float]],
+    intersection: tuple[float, float, float],
+) -> list[tuple[tuple[float, float, float], tuple[float, float, float]]]:
+    """Two segments showing each wall axis extending to ``intersection``.
+
+    Each segment runs from the input axis's nearest endpoint to the
+    intersection, held at that wall's own Z. Returned in input order
+    ``[floor_a, floor_b]``."""
+    ix, iy, _ = intersection
+
+    def _nearest(seg: tuple[tuple[float, float, float], tuple[float, float, float]]) -> tuple[float, float, float]:
+        return min(seg, key=lambda p: (p[0] - ix) ** 2 + (p[1] - iy) ** 2)
+
+    near_a = _nearest(seg_a)
+    near_b = _nearest(seg_b)
+    return [
+        (near_a, (ix, iy, near_a[2])),
+        (near_b, (ix, iy, near_b[2])),
+    ]
+
+
+def resolve_extend_walls_target(
+    target_obj: Any,
+    objs: list[Any],
+    reverse: bool,
+) -> tuple[Any, list[Any]]:
+    """Pick which object is the extend-target and which are extended.
+
+    Default direction: ``objs`` are extended to meet ``target_obj``.
+    Reversed direction (``reverse=True``) swaps the pair — equivalent to
+    having passed them in the opposite order. The swap is well-defined only
+    for the 1+1 case (one target + one other); for ``n>1`` it would be
+    ambiguous, so the default direction is preserved instead."""
+    if reverse and target_obj is not None and len(objs) == 1:
+        return objs[0], [target_obj]
+    return target_obj, objs
+
+
+def displacement_from_x_angle(height: float, x_angle: float) -> float:
+    """Top-edge horizontal displacement for a wall of given vertical ``height``
+    and slope ``x_angle`` (radians). Inverse of ``x_angle_from_displacement``."""
     return height * math.tan(x_angle)
 
 
 def x_angle_from_displacement(height: float, displacement: float) -> float:
     """Recover slope ``x_angle`` (radians) from a top-edge horizontal displacement.
 
-    ``height`` is clamped to ``max(height, 1e-6)`` so vertical walls of effectively
-    zero height map cleanly to ``±π/2`` via ``atan2`` rather than dividing by zero.
-
-    Inverse of :func:`displacement_from_x_angle`."""
+    ``height`` is clamped to ``max(height, 1e-6)`` so zero-height walls map
+    cleanly to ``±π/2`` instead of dividing by zero."""
     return math.atan2(displacement, max(height, 1e-6))
 
 
@@ -260,22 +357,38 @@ def vertical_height_from_extrusion_depth(extrusion_depth: float, x_angle: float)
     return extrusion_depth * abs(math.cos(x_angle))
 
 
+def extrusion_depth_from_vertical_height(vertical_height: float, x_angle: float) -> float:
+    """``vertical_height / cos(x_angle)`` with ``cos`` clamped at ``1e-6`` to
+    stay finite near ``±π/2``."""
+    return vertical_height / max(abs(math.cos(x_angle)), 1e-6)
+
+
+def length_and_height_from_extrusion(
+    extrusion_depth: float,
+    x_angle: float,
+    reference_line_x_extent: float,
+    unit_scale: float,
+) -> tuple[float, float]:
+    """SI ``(length, vertical_height)`` of a LAYER2 wall.
+
+    Height is the *vertical* projection of the slanted depth, not the
+    slanted depth itself."""
+    length = reference_line_x_extent * unit_scale
+    height = vertical_height_from_extrusion_depth(extrusion_depth * unit_scale, x_angle)
+    return length, height
+
+
 def are_axes_collinear(
     seg_a: tuple[tuple[float, float, float], tuple[float, float, float]],
     seg_b: tuple[tuple[float, float, float], tuple[float, float, float]],
-    parallel_threshold: float = 0.9994,
-    line_tolerance: float = 0.05,
+    parallel_threshold: float = PARALLEL_DOT_THRESHOLD,
+    line_tolerance: float = COLLINEAR_LINE_TOLERANCE,
 ) -> bool:
     """True if both axis segments lie on the same infinite line in plan.
 
-    Two conditions: directions must be (anti-)parallel within ``parallel_threshold``
-    (``cos(2°) ≈ 0.9994``), AND any endpoint of B must lie on A's infinite line
-    within ``line_tolerance``. Plan-only (Z ignored) — two parallel walls at
-    different elevations are still considered collinear because the merge operator
-    handles Z resolution itself.
-
-    Used by the wall-join gizmo's state machine: collinear pair → Merge icon at the
-    boundary, perpendicular pair → Join icon at the intersection."""
+    Two conditions: directions must be (anti-)parallel within ``parallel_threshold``,
+    AND any endpoint of B must lie on A's infinite line within ``line_tolerance``.
+    Plan-only (Z ignored)."""
     d1x, d1y = seg_a[1][0] - seg_a[0][0], seg_a[1][1] - seg_a[0][1]
     d2x, d2y = seg_b[1][0] - seg_b[0][0], seg_b[1][1] - seg_b[0][1]
     d1_len = (d1x * d1x + d1y * d1y) ** 0.5
@@ -300,11 +413,7 @@ def closest_endpoint_midpoint(
     seg_a: tuple[tuple[float, float, float], tuple[float, float, float]],
     seg_b: tuple[tuple[float, float, float], tuple[float, float, float]],
 ) -> tuple[float, float, float]:
-    """Midpoint of the closest pair of endpoints between two segments.
-
-    For walls that meet end-to-end this is the shared corner; for walls with a
-    small gap it's the midpoint of the gap. Either way it's the user-meaningful
-    "boundary" where a merge would graft the two segments together."""
+    """Midpoint of the closest endpoint pair between two segments."""
     endpoints_a = (seg_a[0], seg_a[1])
     endpoints_b = (seg_b[0], seg_b[1])
 
@@ -314,3 +423,209 @@ def closest_endpoint_midpoint(
     closest_pair = min(((a, b) for a in endpoints_a for b in endpoints_b), key=lambda pair: _distance_sq(*pair))
     a, b = closest_pair
     return ((a[0] + b[0]) / 2, (a[1] + b[1]) / 2, (a[2] + b[2]) / 2)
+
+
+def compute_path_connection_location(
+    seg_self: tuple[tuple[float, float, float], tuple[float, float, float]],
+    self_conn_type: str,
+    seg_other: tuple[tuple[float, float, float], tuple[float, float, float]],
+    other_conn_type: str,
+    parallel_threshold: float = PARALLEL_DOT_THRESHOLD,
+) -> tuple[float, float, float]:
+    """World-space location of a single ``IfcRelConnectsPathElements`` between
+    two wall axes.
+
+    Priority: ``self``'s ATSTART/ATEND endpoint → ``other``'s ATSTART/ATEND
+    endpoint → axis intersection → closest-endpoint midpoint fallback."""
+    if self_conn_type == "ATSTART":
+        return seg_self[0]
+    if self_conn_type == "ATEND":
+        return seg_self[1]
+    if other_conn_type == "ATSTART":
+        return seg_other[0]
+    if other_conn_type == "ATEND":
+        return seg_other[1]
+    intersection = project_axis_intersection(seg_self, seg_other, parallel_threshold)
+    if intersection is not None:
+        return intersection
+    return closest_endpoint_midpoint(seg_self, seg_other)
+
+
+def _vec_sub(a: tuple[float, float, float], b: tuple[float, float, float]) -> tuple[float, float, float]:
+    return (a[0] - b[0], a[1] - b[1], a[2] - b[2])
+
+
+def _vec_dot(a: tuple[float, float, float], b: tuple[float, float, float]) -> float:
+    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
+
+
+def _vec_cross(a: tuple[float, float, float], b: tuple[float, float, float]) -> tuple[float, float, float]:
+    return (a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0])
+
+
+def _vec_length(v: tuple[float, float, float]) -> float:
+    return (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) ** 0.5
+
+
+def _rotate_around_axis(
+    v: tuple[float, float, float],
+    axis: tuple[float, float, float],
+    angle: float,
+) -> tuple[float, float, float]:
+    """Rotate ``v`` around unit-length ``axis`` by ``angle`` radians."""
+    cos_a = math.cos(angle)
+    sin_a = math.sin(angle)
+    dot = _vec_dot(axis, v)
+    cross = _vec_cross(axis, v)
+    k = 1.0 - cos_a
+    return (
+        v[0] * cos_a + cross[0] * sin_a + axis[0] * dot * k,
+        v[1] * cos_a + cross[1] * sin_a + axis[1] * dot * k,
+        v[2] * cos_a + cross[2] * sin_a + axis[2] * dot * k,
+    )
+
+
+def compute_fillet_polylines(
+    seg_a: tuple[tuple[float, float, float], tuple[float, float, float]],
+    seg_b: tuple[tuple[float, float, float], tuple[float, float, float]],
+    radius: float,
+    arc_resolution: int = FILLET_DEFAULT_ARC_RESOLUTION,
+    parallel_threshold: float = PARALLEL_DOT_THRESHOLD,
+) -> dict:
+    """Preview polylines for a circular fillet at the junction of two axes.
+
+    Returns a dict with ``valid``, ``reason``, ``intersection``, ``tangent_a``
+    / ``tangent_b``, ``arc`` (``arc_resolution + 1`` samples), ``arc_center``,
+    ``arc_radius``, ``sweep_angle``, ``sweep_axis``, ``tangent_offset``,
+    ``wall_a_join_side`` / ``wall_b_join_side`` (ATSTART/ATEND/None),
+    ``invalid_radius`` (tangent overshoots — arc + tangents still populated
+    for warning rendering), and ``invalid_axes`` (set on parallel)."""
+    blank: dict = {
+        "valid": False,
+        "reason": None,
+        "intersection": None,
+        "tangent_a": None,
+        "tangent_b": None,
+        "arc": [],
+        "arc_center": None,
+        "arc_radius": radius,
+        "sweep_angle": 0.0,
+        "sweep_axis": None,
+        "tangent_offset": 0.0,
+        "wall_a_join_side": None,
+        "wall_b_join_side": None,
+        "invalid_radius": False,
+        "invalid_axes": None,
+    }
+
+    intersection = project_axis_intersection(seg_a, seg_b, parallel_threshold)
+    if intersection is None:
+        return {**blank, "reason": "parallel", "invalid_axes": [seg_a, seg_b]}
+
+    def _classify(seg, ipt):
+        d0 = (seg[0][0] - ipt[0]) ** 2 + (seg[0][1] - ipt[1]) ** 2 + (seg[0][2] - ipt[2]) ** 2
+        d1 = (seg[1][0] - ipt[0]) ** 2 + (seg[1][1] - ipt[1]) ** 2 + (seg[1][2] - ipt[2]) ** 2
+        if d0 <= d1:
+            return seg[0], seg[1], "ATSTART"
+        return seg[1], seg[0], "ATEND"
+
+    near_a, far_a, side_a = _classify(seg_a, intersection)
+    near_b, far_b, side_b = _classify(seg_b, intersection)
+
+    # Direction along each segment AWAY from the corner. ``far - intersection``
+    # handles both the shared-corner and extended-axes cases uniformly.
+    dir_a_raw = _vec_sub(far_a, intersection)
+    dir_b_raw = _vec_sub(far_b, intersection)
+    far_len_a = _vec_length(dir_a_raw)
+    far_len_b = _vec_length(dir_b_raw)
+    if far_len_a < 1e-9 or far_len_b < 1e-9:
+        return {**blank, "reason": "near_collinear", "intersection": intersection}
+    dir_a = (dir_a_raw[0] / far_len_a, dir_a_raw[1] / far_len_a, dir_a_raw[2] / far_len_a)
+    dir_b = (dir_b_raw[0] / far_len_b, dir_b_raw[1] / far_len_b, dir_b_raw[2] / far_len_b)
+
+    cos_angle = max(-1.0, min(1.0, _vec_dot(dir_a, dir_b)))
+    angle = math.acos(cos_angle)
+    sweep_angle = math.pi - angle
+    if sweep_angle < 1e-3 or sweep_angle > math.pi - 1e-3:
+        return {
+            **blank,
+            "reason": "near_collinear",
+            "intersection": intersection,
+            "sweep_angle": sweep_angle,
+            "wall_a_join_side": side_a,
+            "wall_b_join_side": side_b,
+        }
+
+    tangent_offset = radius * math.tan(sweep_angle / 2)
+    tangent_a = (
+        intersection[0] + dir_a[0] * tangent_offset,
+        intersection[1] + dir_a[1] * tangent_offset,
+        intersection[2] + dir_a[2] * tangent_offset,
+    )
+    tangent_b = (
+        intersection[0] + dir_b[0] * tangent_offset,
+        intersection[1] + dir_b[1] * tangent_offset,
+        intersection[2] + dir_b[2] * tangent_offset,
+    )
+
+    plane_normal_raw = _vec_cross(dir_a, dir_b)
+    pn_len = _vec_length(plane_normal_raw)
+    if pn_len < 1e-9:
+        return {**blank, "reason": "near_collinear", "intersection": intersection}
+    plane_normal = (
+        plane_normal_raw[0] / pn_len,
+        plane_normal_raw[1] / pn_len,
+        plane_normal_raw[2] / pn_len,
+    )
+
+    perp_a = _vec_cross(plane_normal, dir_a)
+    if _vec_dot(perp_a, dir_b) < 0:
+        perp_a = (-perp_a[0], -perp_a[1], -perp_a[2])
+
+    arc_center = (
+        tangent_a[0] + perp_a[0] * radius,
+        tangent_a[1] + perp_a[1] * radius,
+        tangent_a[2] + perp_a[2] * radius,
+    )
+
+    v_a = _vec_sub(tangent_a, arc_center)
+    v_b = _vec_sub(tangent_b, arc_center)
+    sweep_axis = plane_normal
+    if _vec_dot(_vec_cross(v_a, v_b), plane_normal) < 0:
+        sweep_axis = (-plane_normal[0], -plane_normal[1], -plane_normal[2])
+
+    arc_points: list[tuple[float, float, float]] = []
+    for i in range(arc_resolution + 1):
+        t = i / arc_resolution
+        rotated = _rotate_around_axis(v_a, sweep_axis, sweep_angle * t)
+        arc_points.append(
+            (
+                arc_center[0] + rotated[0],
+                arc_center[1] + rotated[1],
+                arc_center[2] + rotated[2],
+            )
+        )
+
+    # Overshoot check only for convex fillets (positive ``tangent_offset``);
+    # the inverted-fillet case puts tangents past the intersection.
+    invalid_radius = tangent_offset > 0 and (tangent_offset > far_len_a or tangent_offset > far_len_b)
+
+    return {
+        "valid": not invalid_radius,
+        "reason": "invalid_radius" if invalid_radius else None,
+        "intersection": intersection,
+        "tangent_a": tangent_a,
+        "tangent_b": tangent_b,
+        "arc": arc_points,
+        "arc_center": arc_center,
+        "arc_radius": radius,
+        "sweep_angle": sweep_angle,
+        "sweep_axis": sweep_axis,
+        "tangent_offset": tangent_offset,
+        "wall_a_join_side": side_a,
+        "wall_b_join_side": side_b,
+        "leg_a_available": far_len_a,
+        "leg_b_available": far_len_b,
+        "invalid_radius": invalid_radius,
+        "invalid_axes": None,
+    }
diff --git a/src/bonsai/bonsai/core/product.py b/src/bonsai/bonsai/core/product.py
new file mode 100644
index 0000000000..4eaddc833d
--- /dev/null
+++ b/src/bonsai/bonsai/core/product.py
@@ -0,0 +1,64 @@
+# Bonsai - OpenBIM Blender Add-on
+# Copyright (C) 2026
+#
+# This file is part of Bonsai.
+#
+# Bonsai is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# Bonsai is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Bonsai.  If not, see <http://www.gnu.org/licenses/>.
+#
+# This file was generated with the assistance of an AI coding tool.
+
+from __future__ import annotations
+
+import math
+from collections.abc import Iterable
+from typing import TYPE_CHECKING
+
+import bonsai.core.geometry
+
+if TYPE_CHECKING:
+    import bpy
+
+    import bonsai.tool as tool
+
+
+Z_ROTATION_ALIGNMENT_TOLERANCE = 1e-9
+
+
+def _z_rotation_diff(target_z: float, source_z: float) -> float:
+    """Signed Z-Euler difference wrapped to [-π, π]."""
+    return (target_z - source_z + math.pi) % (2 * math.pi) - math.pi
+
+
+def copy_z_rotation_to_selected(
+    ifc: type[tool.Ifc],
+    geometry: type[tool.Geometry],
+    surveyor: type[tool.Surveyor],
+    *,
+    active: bpy.types.Object,
+    targets: Iterable[bpy.types.Object],
+    flip: bool = False,
+) -> int:
+    """Apply ``active``'s Z-Euler rotation to each target."""
+    source_z = surveyor.get_z_rotation(active)
+    if flip:
+        source_z += math.pi
+    rotated = 0
+    for obj in targets:
+        if abs(_z_rotation_diff(surveyor.get_z_rotation(obj), source_z)) < Z_ROTATION_ALIGNMENT_TOLERANCE:
+            continue
+        surveyor.set_z_rotation(obj, source_z)
+        rotated += 1
+        if ifc.get_entity(obj) is not None:
+            bonsai.core.geometry.edit_object_placement(ifc, geometry, surveyor, obj=obj)
+    return rotated