Merge pull request 'feat(assembly): fixed reference planes (Top/Front/Right) + solver docs' (#307) from feat/assembly-origin-planes into main

Reviewed-on: #307
2026-02-21 15:09:55 +00:00
parent 148bed59f6 acc255972d
commit 0f8fa0be86
14 changed files with 1259 additions and 9 deletions
--- a/docs/src/solver/assembly-integration.md
+++ b/docs/src/solver/assembly-integration.md
@@ -0,0 +1,165 @@
+# Assembly Integration
+
+The Kindred solver integrates with FreeCAD's Assembly workbench through the KCSolve pluggable solver framework. This page describes the bridge layer, preference system, and interactive drag protocol.
+
+## KindredSolver class
+
+**Source:** `mods/solver/kindred_solver/solver.py`
+
+`KindredSolver` subclasses `kcsolve.IKCSolver` and implements the solver interface:
+
+```python
+class KindredSolver(kcsolve.IKCSolver):
+    def name(self):
+        return "Kindred (Newton-Raphson)"
+
+    def supported_joints(self):
+        return list(_SUPPORTED)  # 20 of 24 BaseJointKind values
+
+    def solve(self, ctx):        # Static solve
+    def diagnose(self, ctx):     # Constraint analysis
+    def pre_drag(self, ctx, drag_parts):   # Begin drag session
+    def drag_step(self, drag_placements):  # Mouse move during drag
+    def post_drag(self):                   # End drag session
+    def is_deterministic(self):  # Returns True
+```
+
+### Registration
+
+The solver is registered at addon load time via `Init.py`:
+
+```python
+import kcsolve
+from kindred_solver import KindredSolver
+kcsolve.register_solver("kindred", KindredSolver)
+```
+
+The `mods/solver/` directory is a FreeCAD addon discovered by the addon loader through its `package.xml` manifest.
+
+### Supported joints
+
+The Kindred solver handles 20 of the 24 `BaseJointKind` values. The remaining 4 are stubs that produce no residuals:
+
+| Supported | Stub (no residuals) |
+|-----------|-------------------|
+| Coincident, PointOnLine, PointInPlane, Concentric, Tangent, Planar, LineInPlane, Parallel, Perpendicular, Angle, Fixed, Revolute, Cylindrical, Slider, Ball, Screw, Universal, Gear, RackPinion, DistancePointPoint | Cam, Slot, DistanceCylSph, Custom |
+
+### Joint limits
+
+Joint travel limits (`Constraint.limits`) are accepted but not enforced. The solver logs a warning once per instance when limits are encountered. Enforcing inequality constraints requires active-set or barrier-method extensions beyond the current Newton-Raphson formulation.
+
+## Solver selection
+
+### C++ preference
+
+`AssemblyObject::getOrCreateSolver()` reads the user preference to select the solver backend:
+
+```cpp
+ParameterGrp::handle hGrp = App::GetApplication().GetParameterGroupByPath(
+    "User parameter:BaseApp/Preferences/Mod/Assembly");
+std::string solverName = hGrp->GetASCII("Solver", "");
+solver_ = KCSolve::SolverRegistry::instance().get(solverName);
+```
+
+An empty string (`""`) returns the registry default (the first solver registered, which is OndselSolver). Setting `"kindred"` selects the Kindred solver.
+
+`resetSolver()` clears the cached solver instance so the next solve picks up preference changes.
+
+### Preferences UI
+
+The Assembly preferences page (`Edit > Preferences > Assembly`) includes a "Solver backend" combo box populated from the registry at load time:
+
+- **Default** -- empty string, uses the registry default (OndselSolver)
+- **OndselSolver (Lagrangian)** -- `"ondsel"`
+- **Kindred (Newton-Raphson)** -- `"kindred"` (available when the solver addon is loaded)
+
+The preference is stored as `Mod/Assembly/Solver` in the FreeCAD parameter system.
+
+### Programmatic switching
+
+From the Python console:
+
+```python
+import FreeCAD
+pref = FreeCAD.ParamGet("User parameter:BaseApp/Preferences/Mod/Assembly")
+
+# Switch to Kindred
+pref.SetString("Solver", "kindred")
+
+# Switch back to default
+pref.SetString("Solver", "")
+
+# Force the active assembly to pick up the change
+if hasattr(FreeCADGui, "ActiveDocument"):
+    for obj in FreeCAD.ActiveDocument.Objects:
+        if hasattr(obj, "resetSolver"):
+            obj.resetSolver()
+```
+
+## Interactive drag protocol
+
+The drag protocol provides real-time constraint solving during viewport part dragging. It is a three-phase protocol:
+
+### pre_drag(ctx, drag_parts)
+
+Called when the user begins dragging. Stores the context and dragged part IDs, then runs a full solve to establish the starting state.
+
+```python
+def pre_drag(self, ctx, drag_parts):
+    self._drag_ctx = ctx
+    self._drag_parts = set(drag_parts)
+    return self.solve(ctx)
+```
+
+### drag_step(drag_placements)
+
+Called on each mouse move. Updates the dragged parts' placements in the stored context, then re-solves. Since the parts moved only slightly from the previous position, Newton-Raphson converges in 1-2 iterations.
+
+```python
+def drag_step(self, drag_placements):
+    ctx = self._drag_ctx
+    for pr in drag_placements:
+        for part in ctx.parts:
+            if part.id == pr.id:
+                part.placement = pr.placement
+                break
+    return self.solve(ctx)
+```
+
+### post_drag()
+
+Called when the drag ends. Clears the stored state.
+
+```python
+def post_drag(self):
+    self._drag_ctx = None
+    self._drag_parts = None
+```
+
+### Performance notes
+
+The current implementation re-solves from scratch on each drag step, using the updated placements as the initial guess. This is correct and simple. For assemblies with fewer than ~50 parts, interactive frame rates are maintained because:
+
+- Newton-Raphson converges in 1-2 iterations from a nearby initial guess
+- Pre-passes eliminate fixed parameters before the iterative loop
+- The symbolic Jacobian is recomputed each step (no caching yet)
+
+For larger assemblies, cached incremental solving (reusing the decomposition and Jacobian structure across drag steps) is planned as a future optimization.
+
+## Diagnostics integration
+
+`diagnose(ctx)` builds the constraint system and runs overconstrained detection, returning a list of `kcsolve.ConstraintDiagnostic` objects. The Assembly module calls this to populate the constraint diagnostics panel.
+
+```python
+def diagnose(self, ctx):
+    system = _build_system(ctx)
+    residuals = substitution_pass(system.all_residuals, system.params)
+    return _run_diagnostics(residuals, system.params, system.residual_ranges, ctx)
+```
+
+## Not yet implemented
+
+- **Kinematic simulation** (`run_kinematic`, `num_frames`, `update_for_frame`) -- the base class defaults return `Failed`. Requires time-stepping integration with motion driver expression evaluation.
+- **Joint limit enforcement** -- inequality constraints need active-set or barrier solver extensions.
+- **Fixed-joint bundling** (`supports_bundle_fixed()` returns `False`) -- the solver receives unbundled parts; the Assembly module pre-bundles when needed.
+- **Native export** (`export_native()`) -- no solver-native debug format defined.
--- a/docs/src/solver/constraints.md
+++ b/docs/src/solver/constraints.md
@@ -0,0 +1,116 @@
+# Constraints
+
+Each constraint type maps to a class that produces residual expressions. The residuals equal zero when the constraint is satisfied. The number of residuals equals the number of degrees of freedom removed.
+
+**Source:** `mods/solver/kindred_solver/constraints.py`, `mods/solver/kindred_solver/geometry.py`
+
+## Constraint vocabulary
+
+### Point constraints
+
+| Type | DOF removed | Residuals |
+|------|-------------|-----------|
+| **Coincident** | 3 | `p_i - p_j` (world-frame marker origins coincide) |
+| **PointOnLine** | 2 | Two components of `(p_i - p_j) x z_j` (point lies on line through `p_j` along `z_j`) |
+| **PointInPlane** | 1 | `(p_i - p_j) . z_j - offset` (signed distance to plane) |
+
+### Orientation constraints
+
+| Type | DOF removed | Residuals |
+|------|-------------|-----------|
+| **Parallel** | 2 | Two components of `z_i x z_j` (cross product of Z-axes is zero) |
+| **Perpendicular** | 1 | `z_i . z_j` (dot product of Z-axes is zero) |
+| **Angle** | 1 | `z_i . z_j - cos(angle)` |
+
+### Axis/surface constraints
+
+| Type | DOF removed | Residuals |
+|------|-------------|-----------|
+| **Concentric** | 4 | Parallel Z-axes (2) + point-on-line (2) |
+| **Tangent** | 1 | `(p_i - p_j) . z_j` (signed distance along normal) |
+| **Planar** | 3 | Parallel normals (2) + point-in-plane (1) |
+| **LineInPlane** | 2 | Point-in-plane (1) + `z_i . n_j` (line direction perpendicular to normal) (1) |
+
+### Kinematic joints
+
+| Type | DOF removed | DOF remaining | Residuals |
+|------|-------------|---------------|-----------|
+| **Fixed** | 6 | 0 | Coincident origins (3) + quaternion error imaginary parts (3) |
+| **Ball** | 3 | 3 | Coincident origins (same as Coincident) |
+| **Revolute** | 5 | 1 (rotation about Z) | Coincident origins (3) + parallel Z-axes (2) |
+| **Cylindrical** | 4 | 2 (rotation + slide) | Parallel Z-axes (2) + point-on-line (2) |
+| **Slider** | 5 | 1 (slide along Z) | Parallel Z-axes (2) + point-on-line (2) + twist lock: `x_i . y_j` (1) |
+| **Screw** | 5 | 1 (helical) | Cylindrical (4) + pitch coupling: `axial - pitch * qz_rel / pi` (1) |
+| **Universal** | 4 | 2 (rotation about each Z) | Coincident origins (3) + perpendicular Z-axes (1) |
+
+### Mechanical elements
+
+| Type | DOF removed | Residuals |
+|------|-------------|-----------|
+| **Gear** | 1 | `r_i * qz_i + r_j * qz_j` (coupled rotation via quaternion Z-components) |
+| **RackPinion** | 1 | `translation - 2 * pitch_radius * qz_i` (rotation-translation coupling) |
+| **Cam** | 0 | Stub (no residuals) |
+| **Slot** | 0 | Stub (no residuals) |
+
+### Distance constraints
+
+| Type | DOF removed | Residuals |
+|------|-------------|-----------|
+| **DistancePointPoint** | 1 | `\|p_i - p_j\|^2 - d^2` (squared form avoids sqrt in Jacobian) |
+| **DistanceCylSph** | 0 | Stub (geometry classification dependent) |
+
+## Marker convention
+
+Every constraint references two parts (`body_i`, `body_j`) with local coordinate frames called markers. Each marker has a position (attachment point on the part) and a quaternion (orientation).
+
+The marker Z-axis defines the constraint direction:
+- **Revolute:** Z-axis = hinge axis
+- **Planar:** Z-axis = face normal
+- **PointOnLine:** Z-axis = line direction
+- **Slider:** Z-axis = slide direction
+
+The solver computes world-frame marker axes by composing the body quaternion with the marker quaternion: `q_world = q_body * q_marker`, then rotating unit vectors through the result.
+
+## Fixed constraint orientation
+
+The Fixed constraint locks all 6 DOF using a quaternion error formulation:
+
+1. Compute total orientation: `q_i = q_body_i * q_marker_i`, `q_j = q_body_j * q_marker_j`
+2. Compute relative quaternion: `q_err = conj(q_i) * q_j`
+3. When orientations match, `q_err` is the identity quaternion `(1, 0, 0, 0)`
+4. Residuals are the three imaginary components of `q_err` (should be zero)
+
+The quaternion normalization constraint on each body provides the fourth equation needed to fully determine the quaternion.
+
+## Rotation proxies for mechanical constraints
+
+Gear, RackPinion, and Screw constraints need to measure rotation angles. Rather than extracting Euler angles (which would introduce transcendentals), they use the Z-component of a relative quaternion as a proxy:
+
+```
+q_local = conj(q_marker) * q_body * q_marker
+angle ~ 2 * qz_local  (for small angles)
+```
+
+This is exact at the solution and has correct gradient direction, which is sufficient for Newton-Raphson convergence from a nearby initial guess.
+
+## Geometry helpers
+
+The `geometry.py` module provides Expr-level vector operations used by constraint classes:
+
+- `marker_z_axis(body, marker_quat)` -- world-frame Z-axis via `quat_rotate(q_body * q_marker, [0,0,1])`
+- `marker_x_axis(body, marker_quat)` -- world-frame X-axis (used by Slider twist lock)
+- `marker_y_axis(body, marker_quat)` -- world-frame Y-axis (used by Slider twist lock)
+- `dot3(a, b)` -- dot product of Expr triples
+- `cross3(a, b)` -- cross product of Expr triples
+- `point_plane_distance(point, origin, normal)` -- signed distance
+- `point_line_perp_components(point, origin, dir)` -- two perpendicular distance components
+
+## Writing a new constraint
+
+To add a constraint type:
+
+1. Subclass `ConstraintBase` in `constraints.py`
+2. Implement `residuals()` returning a list of `Expr` nodes
+3. Add a case in `solver.py:_build_constraint()` to instantiate it from `BaseJointKind`
+4. Add the `BaseJointKind` value to `_SUPPORTED` in `solver.py`
+5. Add the residual count to the tables in `decompose.py`
--- a/docs/src/solver/diagnostics.md
+++ b/docs/src/solver/diagnostics.md
@@ -0,0 +1,117 @@
+# Diagnostics
+
+The solver provides three levels of constraint analysis: system-wide DOF counting, per-entity DOF decomposition, and overconstrained/conflicting constraint detection.
+
+## DOF counting
+
+**Source:** `mods/solver/kindred_solver/dof.py`
+
+Degrees of freedom are computed from the Jacobian rank:
+
+```
+DOF = n_free_params - rank(J)
+```
+
+Where `n_free_params` is the number of non-fixed parameters and `rank(J)` is the numerical rank of the Jacobian evaluated at current parameter values (SVD with tolerance `1e-8`).
+
+A well-constrained assembly has `DOF = 0` (exactly enough constraints to determine all positions). Positive DOF means underconstrained (parts can still move). Negative DOF is not possible with this formulation -- instead, rank deficiency in an overdetermined system indicates redundant constraints.
+
+The DOF value is reported in `SolveResult.dof` after every solve.
+
+## Per-entity DOF
+
+**Source:** `mods/solver/kindred_solver/diagnostics.py`
+
+`per_entity_dof()` breaks down the DOF count per body, identifying which motions remain free for each part:
+
+1. Build the full Jacobian
+2. For each non-grounded body, extract the 7 columns corresponding to its parameters
+3. Compute SVD of the sub-matrix; rank = number of constrained directions
+4. `remaining_dof = 7 - rank` (includes the quaternion normalization constraint counted in the rank)
+5. Classify null-space vectors as free motions by analyzing their translation vs. rotation components:
+   - Pure translation: >80% of the null vector's energy is in `tx, ty, tz` components
+   - Pure rotation: >80% of the energy is in `qw, qx, qy, qz` components
+   - Helical: mixed
+
+Returns a list of `EntityDOF` dataclasses:
+
+```python
+@dataclass
+class EntityDOF:
+    entity_id: str
+    remaining_dof: int
+    free_motions: list[str]  # e.g., ["rotation about Z", "translation along X"]
+```
+
+## Overconstrained detection
+
+**Source:** `mods/solver/kindred_solver/diagnostics.py`
+
+`find_overconstrained()` identifies redundant and conflicting constraints when the system is overconstrained (Jacobian is rank-deficient). It runs automatically when `solve()` fails to converge.
+
+### Algorithm
+
+Following the approach used by SolvSpace:
+
+1. **Check rank.** Build the full Jacobian `J`, compute its rank via SVD. If `rank == n_residuals`, the system is not overconstrained -- return empty.
+
+2. **Find redundant constraints.** For each constraint, temporarily remove its rows from J and re-check rank. If the rank is preserved, the constraint is **redundant** (removing it doesn't change the system's effective equations).
+
+3. **Distinguish conflicting from merely redundant.** Compute the left null space of J (columns of U beyond the rank). Project the residual vector onto this null space:
+   ```
+   null_residual = U_null^T @ r
+   residual_projection = U_null @ null_residual
+   ```
+   If a redundant constraint's residuals have significant projection onto the null space, it is **conflicting** -- it's both redundant and unsatisfied, meaning it contradicts other constraints.
+
+### Diagnostic output
+
+Returns `ConstraintDiag` dataclasses:
+
+```python
+@dataclass
+class ConstraintDiag:
+    constraint_index: int
+    kind: str           # "redundant" or "conflicting"
+    detail: str         # Human-readable explanation
+```
+
+These are converted to `kcsolve.ConstraintDiagnostic` objects in the IKCSolver bridge:
+
+| ConstraintDiag.kind | kcsolve.DiagnosticKind |
+|---------------------|----------------------|
+| `"redundant"` | `Redundant` |
+| `"conflicting"` | `Conflicting` |
+
+### Example
+
+Two Fixed joints between the same pair of parts:
+
+- Joint A: 6 residuals (3 position + 3 orientation)
+- Joint B: 6 residuals (same as Joint A)
+
+Jacobian rank = 6 (Joint B's rows are linearly dependent on Joint A's). Both joints are detected as redundant. If the joints specify different relative positions, both are also flagged as conflicting.
+
+## Solution preferences
+
+**Source:** `mods/solver/kindred_solver/preference.py`
+
+Solution preferences guide the solver toward physically intuitive solutions when multiple valid configurations exist.
+
+### Minimum-movement weighting
+
+The weight vector scales the Newton step to prefer solutions near the initial configuration. Translation parameters get weight `1.0`, quaternion parameters get weight `(180/pi)^2 ~ 3283`. This makes a 1-radian rotation equally "expensive" as a ~57-unit translation.
+
+The weighted minimum-norm step is:
+
+```
+J_scaled = J @ diag(W^{-1/2})
+dx_scaled = lstsq(J_scaled, -r)
+dx = dx_scaled * W^{-1/2}
+```
+
+This produces the minimum-norm solution in the weighted parameter space, biasing toward small movements.
+
+### Half-space tracking
+
+Described in detail in [Solving Algorithms: Half-space tracking](solving.md#half-space-tracking). Preserves the initial configuration's "branch" for constraints with multiple valid solutions by detecting and correcting branch crossings during iteration.
--- a/docs/src/solver/expression-dag.md
+++ b/docs/src/solver/expression-dag.md
@@ -0,0 +1,96 @@
+# Expression DAG
+
+The expression DAG is the foundation of the Kindred solver. All constraint equations, Jacobian entries, and residuals are built as immutable trees of `Expr` nodes. This lets the solver compute exact symbolic derivatives and simplify constant sub-expressions before the iterative solve loop.
+
+**Source:** `mods/solver/kindred_solver/expr.py`
+
+## Node types
+
+Every node is a subclass of `Expr` and implements three methods:
+
+- `eval(env)` -- evaluate the expression given a name-to-value dictionary
+- `diff(var)` -- return a new Expr tree for the partial derivative with respect to `var`
+- `simplify()` -- return an algebraically simplified copy
+
+### Leaf nodes
+
+| Node | Description | diff(x) |
+|------|-------------|---------|
+| `Const(v)` | Literal floating-point value | 0 |
+| `Var(name)` | Named parameter (from `ParamTable`) | 1 if name matches, else 0 |
+
+### Unary nodes
+
+| Node | Description | diff(x) |
+|------|-------------|---------|
+| `Neg(f)` | Negation: `-f` | `-f'` |
+| `Sin(f)` | Sine: `sin(f)` | `cos(f) * f'` |
+| `Cos(f)` | Cosine: `cos(f)` | `-sin(f) * f'` |
+| `Sqrt(f)` | Square root: `sqrt(f)` | `f' / (2 * sqrt(f))` |
+
+### Binary nodes
+
+| Node | Description | diff(x) |
+|------|-------------|---------|
+| `Add(a, b)` | Sum: `a + b` | `a' + b'` |
+| `Sub(a, b)` | Difference: `a - b` | `a' - b'` |
+| `Mul(a, b)` | Product: `a * b` | `a'b + ab'` (product rule) |
+| `Div(a, b)` | Quotient: `a / b` | `(a'b - ab') / b^2` (quotient rule) |
+| `Pow(a, n)` | Power: `a^n` (constant exponent only) | `n * a^(n-1) * a'` |
+
+### Sentinels
+
+`ZERO = Const(0.0)` and `ONE = Const(1.0)` are pre-allocated constants used by `diff()` to avoid allocating trivial nodes.
+
+## Operator overloading
+
+Python's arithmetic operators are overloaded on `Expr`, so constraints can be written in natural notation:
+
+```python
+from kindred_solver.expr import Var, Const
+
+x = Var("x")
+y = Var("y")
+
+# Build the expression: x^2 + 2*x*y - 1
+expr = x**2 + 2*x*y - Const(1.0)
+
+# Evaluate at x=3, y=4
+expr.eval({"x": 3.0, "y": 4.0})  # 32.0
+
+# Symbolic derivative w.r.t. x
+dx = expr.diff("x").simplify()  # 2*x + 2*y
+dx.eval({"x": 3.0, "y": 4.0})   # 14.0
+```
+
+The `_wrap()` helper coerces plain `int` and `float` values to `Const` nodes automatically, so `2 * x` works without wrapping the `2`.
+
+## Simplification
+
+`simplify()` applies algebraic identities bottom-up:
+
+- Constant folding: `Const(2) + Const(3)` becomes `Const(5)`
+- Identity elimination: `x + 0 = x`, `x * 1 = x`, `x^0 = 1`, `x^1 = x`
+- Zero propagation: `0 * x = 0`
+- Negation collapse: `-(-x) = x`
+- Power expansion: `x^2` becomes `x * x` (avoids `pow()` in evaluation)
+
+Simplification is applied once to each Jacobian entry after symbolic differentiation, before the solve loop begins. This reduces the expression tree size and speeds up repeated evaluation.
+
+## How the solver uses expressions
+
+1. **Parameter registration.** `ParamTable.add("Part001/tx", 10.0)` creates a `Var("Part001/tx")` node and records its current value.
+
+2. **Constraint building.** Constraint classes compose `Var` nodes with arithmetic to produce residual `Expr` trees. For example, `CoincidentConstraint` builds `body_i.world_point() - body_j.world_point()`, producing 3 residual expressions.
+
+3. **Jacobian construction.** Newton-Raphson calls `r.diff(name).simplify()` for every (residual, free parameter) pair to build the symbolic Jacobian. This happens once before the solve loop.
+
+4. **Evaluation.** Each Newton iteration calls `expr.eval(env)` on every residual and Jacobian entry using the current parameter snapshot. `eval()` is a simple recursive tree walk with dictionary lookups.
+
+## Design notes
+
+**Why not numpy directly?** Symbolic expressions give exact derivatives without finite-difference approximations, and enable pre-passes (substitution, single-equation solve) that can eliminate variables before the iterative solver runs. The overhead of tree evaluation is acceptable for the problem sizes encountered in assembly solving (typically tens to hundreds of variables).
+
+**Why immutable?** Immutability means `diff()` can safely share sub-tree references between the original and derivative expressions. It also simplifies the substitution pass, which rebuilds trees with `Const` nodes replacing fixed `Var` nodes.
+
+**Limitations.** `Pow` differentiation only supports constant exponents. Variable exponents would require logarithmic differentiation (`d/dx f^g = f^g * (g' * ln(f) + g * f'/f)`), which hasn't been needed for assembly constraints.
--- a/docs/src/solver/overview.md
+++ b/docs/src/solver/overview.md
@@ -0,0 +1,92 @@
+# Kindred Solver Overview
+
+The Kindred solver is an expression-based Newton-Raphson constraint solver for the Assembly workbench. It is a pure-Python implementation that registers as a pluggable backend through the [KCSolve framework](../architecture/ondsel-solver.md), providing an alternative to the built-in OndselSolver (Lagrangian) backend.
+
+## Architecture
+
+```
+                            Assembly Module
+                                  │
+                      ┌───────────┴───────────┐
+                      │    SolverRegistry      │
+                      │  get("kindred")        │
+                      └───────────┬───────────┘
+                                  │
+                      ┌───────────┴───────────┐
+                      │    KindredSolver       │
+                      │  (kcsolve.IKCSolver)   │
+                      └───────────┬───────────┘
+                                  │
+              ┌───────────────────┼───────────────────┐
+              │                   │                   │
+     ┌────────┴────────┐  ┌──────┴──────┐  ┌────────┴────────┐
+     │  _build_system  │  │  Solve      │  │  Diagnostics    │
+     │  ────────────── │  │  ─────      │  │  ───────────    │
+     │  ParamTable     │  │  pre-passes │  │  DOF counting   │
+     │  RigidBody      │  │  Newton-R   │  │  overconstrained│
+     │  Constraints    │  │  BFGS       │  │  per-entity DOF │
+     │  Residuals      │  │  decompose  │  │  half-spaces    │
+     └─────────────────┘  └─────────────┘  └─────────────────┘
+```
+
+## Design principles
+
+**Symbolic differentiation.** All constraint equations are built as immutable expression DAGs (`Expr` trees). The Jacobian is computed symbolically via `expr.diff()` rather than finite differences. This gives exact derivatives, avoids numerical step-size tuning, and allows pre-passes to simplify or eliminate trivial equations before the iterative solver runs.
+
+**Residual-based formulation.** Each constraint produces a list of residual expressions that should evaluate to zero when satisfied. A Coincident constraint produces 3 residuals (dx, dy, dz), a Revolute produces 5 (3 position + 2 axis alignment), and so on. The solver minimizes the residual vector norm.
+
+**Unit quaternions for rotation.** Orientation is parameterized as a unit quaternion (w, x, y, z) rather than Euler angles, avoiding gimbal lock. A quaternion normalization residual (qw^2 + qx^2 + qy^2 + qz^2 - 1 = 0) is added for each free body, and quaternions are re-projected onto the unit sphere after each Newton step.
+
+**Current placements as initial guess.** The solver uses the parts' current positions as the initial guess, so it naturally converges to the nearest solution. Combined with half-space tracking, this produces physically intuitive results without branch-switching surprises.
+
+## Solve pipeline
+
+When `KindredSolver.solve(ctx)` is called with a `SolveContext`:
+
+1. **Build system** (`_build_system`) -- Create a `ParamTable` with 7 parameters per part (tx, ty, tz, qw, qx, qy, qz). Grounded parts have all parameters fixed. Build constraint objects from the context, collect their residual expressions, and add quaternion normalization residuals for free bodies.
+
+2. **Solution preferences** -- Compute half-space trackers for branching constraints (Distance, Parallel, Angle, Perpendicular) and build a minimum-movement weight vector that penalizes quaternion changes more than translation changes.
+
+3. **Pre-passes** -- Run the substitution pass (replace fixed parameters with constants) and the single-equation pass (analytically solve residuals with only one free variable).
+
+4. **Solve** -- For assemblies with 8+ free bodies, decompose the constraint graph into biconnected components and solve each cluster independently. For smaller assemblies, solve the full system monolithically. In both cases, use Newton-Raphson first, falling back to L-BFGS-B if Newton doesn't converge.
+
+5. **Post-process** -- Count degrees of freedom via Jacobian SVD rank. On failure, run overconstrained detection to identify redundant or conflicting constraints. Extract solved placements from the parameter table.
+
+## Module map
+
+| Module | Purpose |
+|--------|---------|
+| `solver.py` | `KindredSolver` class: IKCSolver bridge, solve/diagnose/drag entry points |
+| `expr.py` | Immutable expression DAG with eval, diff, simplify |
+| `params.py` | Parameter table: named variables with fixed/free tracking |
+| `entities.py` | `RigidBody`: 7-DOF entity owning solver parameters |
+| `quat.py` | Quaternion rotation as polynomial Expr trees |
+| `geometry.py` | Marker axis extraction, vector ops (dot, cross, point-plane, point-line) |
+| `constraints.py` | 24 constraint classes producing residual expressions |
+| `newton.py` | Newton-Raphson with symbolic Jacobian, quaternion renormalization |
+| `bfgs.py` | L-BFGS-B fallback via scipy |
+| `prepass.py` | Substitution pass and single-equation analytical solve |
+| `decompose.py` | Biconnected component graph decomposition and cluster-by-cluster solving |
+| `dof.py` | DOF counting via Jacobian SVD rank |
+| `diagnostics.py` | Overconstrained detection, per-entity DOF classification |
+| `preference.py` | Half-space tracking and minimum-movement weighting |
+
+## File locations
+
+- **Solver addon:** `mods/solver/` (git submodule)
+- **KCSolve C++ framework:** `src/Mod/Assembly/Solver/`
+- **Python bindings:** `src/Mod/Assembly/Solver/bindings/`
+- **Integration tests:** `src/Mod/Assembly/AssemblyTests/TestKindredSolverIntegration.py`
+- **Unit tests:** `mods/solver/tests/`
+
+## Related
+
+- [Expression DAG](expression-dag.md) -- the Expr type system
+- [Constraints](constraints.md) -- constraint vocabulary and residuals
+- [Solving algorithms](solving.md) -- Newton-Raphson, BFGS, decomposition
+- [Diagnostics](diagnostics.md) -- DOF counting, overconstrained detection
+- [Assembly integration](assembly-integration.md) -- IKCSolver bridge, preferences, drag
+- [Writing a custom solver](writing-a-solver.md) -- tutorial
+- [KCSolve architecture](../architecture/ondsel-solver.md) -- pluggable solver framework
+- [KCSolve Python API](../reference/kcsolve-python.md) -- kcsolve module reference
--- a/docs/src/solver/solving.md
+++ b/docs/src/solver/solving.md
@@ -0,0 +1,128 @@
+# Solving Algorithms
+
+The Kindred solver uses a multi-stage pipeline: pre-passes reduce the system, Newton-Raphson iterates toward a solution, and L-BFGS-B provides a fallback. For large assemblies, graph decomposition splits the system into independent clusters solved in sequence.
+
+## Pre-passes
+
+**Source:** `mods/solver/kindred_solver/prepass.py`
+
+Pre-passes run before the iterative solver and can eliminate variables analytically, reducing the problem size and improving convergence.
+
+### Substitution pass
+
+Replaces all fixed-parameter `Var` nodes with `Const` nodes carrying their current values, then simplifies. This compiles grounded-body parameters and previously-solved variables out of the expression trees.
+
+After substitution, residuals involving only fixed parameters simplify to constants (typically zero), and Jacobian entries for those parameters become exactly zero. This reduces the effective system size without changing the linear algebra.
+
+### Single-equation pass
+
+Scans residuals for any that depend on exactly one free variable. If the residual is linear in that variable (`a*x + b = 0`), it solves `x = -b/a` analytically, fixes the variable, and re-substitutes.
+
+The pass repeats until no more single-variable residuals can be solved. This handles cascading dependencies: solving one variable may reduce another residual to single-variable form.
+
+## Newton-Raphson
+
+**Source:** `mods/solver/kindred_solver/newton.py`
+
+The primary iterative solver. Each iteration:
+
+1. Evaluate the residual vector `r` and check convergence (`||r|| < tol`)
+2. Evaluate the Jacobian matrix `J` by calling `expr.eval()` on pre-computed symbolic derivatives
+3. Solve `J @ dx = -r` via `numpy.linalg.lstsq` (handles rank-deficient systems)
+4. Update parameters: `x += dx`
+5. Apply half-space correction (if configured)
+6. Re-normalize quaternions to unit length
+
+### Symbolic Jacobian
+
+The Jacobian is built once before the solve loop by calling `r.diff(name).simplify()` for every (residual, free parameter) pair. The resulting `Expr` trees are stored and re-evaluated at the current parameter values each iteration. This gives exact derivatives with no step-size tuning.
+
+### Weighted minimum-norm
+
+When a weight vector is provided, the step is column-scaled to produce the weighted minimum-norm solution. The solver scales J by W^{-1/2}, solves the scaled system, then unscales the step. This biases the solver toward solutions requiring smaller parameter changes in high-weight dimensions.
+
+The default weight vector assigns `1.0` to translation parameters and `~3283` to quaternion parameters (the square of 180/pi), making a 1-radian rotation equivalent to a ~57-unit translation. This produces physically intuitive solutions that prefer translating over rotating.
+
+### Quaternion renormalization
+
+After each Newton step, quaternion parameter groups `(qw, qx, qy, qz)` are re-projected onto the unit sphere by dividing by their norm. This prevents the quaternion from drifting away from unit length during iteration (the quaternion normalization residual only enforces this at convergence, not during intermediate steps).
+
+If a quaternion degenerates to near-zero norm, it is reset to the identity quaternion `(1, 0, 0, 0)`.
+
+### Convergence
+
+Newton-Raphson runs for up to 100 iterations with tolerance `1e-10` on the residual norm. For well-conditioned systems near the solution, convergence is typically quadratic (3-5 iterations). Interactive drag from a nearby position typically converges in 1-2 iterations.
+
+## L-BFGS-B fallback
+
+**Source:** `mods/solver/kindred_solver/bfgs.py`
+
+If Newton-Raphson fails to converge, L-BFGS-B minimizes the sum of squared residuals: `f(x) = 0.5 * sum(r_i^2)`. This is a quasi-Newton method that approximates the Hessian from gradient history, with bounded memory usage.
+
+The gradient is computed analytically from the same symbolic Jacobian: `grad = J^T @ r`. This is passed directly to `scipy.optimize.minimize` via the `jac=True` interface to avoid redundant function evaluations.
+
+L-BFGS-B is more robust for ill-conditioned systems where the Jacobian is nearly singular, but converges more slowly (superlinear rather than quadratic). It runs for up to 200 iterations.
+
+If scipy is not available, the fallback is skipped gracefully.
+
+## Graph decomposition
+
+**Source:** `mods/solver/kindred_solver/decompose.py`
+
+For assemblies with 8 or more free bodies, the solver decomposes the constraint graph into clusters and solves them independently. This improves performance for large assemblies by reducing the Jacobian size from O(n^2) to the sum of smaller cluster Jacobians.
+
+### Algorithm
+
+1. **Build constraint graph.** Bodies are nodes, constraints are edges weighted by their residual count (DOF removed). Grounded bodies are tagged.
+
+2. **Find biconnected components.** Using `networkx.biconnected_components()`, decompose the graph into rigid clusters. Articulation points (bodies shared between clusters) are identified.
+
+3. **Build block-cut tree.** A bipartite graph of clusters and articulation points, rooted at a grounded cluster.
+
+4. **BFS ordering.** Traverse the block-cut tree root-to-leaf, producing a solve order where grounded clusters come first and boundary conditions propagate outward.
+
+5. **Solve each cluster.** For each cluster in order:
+   - Fix boundary bodies that were already solved by previous clusters (their parameters become constants)
+   - Collect the cluster's residuals and quaternion normalization equations
+   - Run substitution pass (compiles fixed boundary values to constants)
+   - Newton-Raphson + BFGS fallback on the reduced system
+   - Mark the cluster's bodies as solved
+   - Unfix boundary parameters for downstream clusters
+
+### Example
+
+Consider a chain of 4 bodies: `Ground -- A -- B -- C` with joints at each connection. This decomposes into two biconnected components (if the joints create articulation points):
+
+- Cluster 1: {Ground, A} -- solved first (grounded)
+- Cluster 2: {A, B, C} -- solved second with A's parameters fixed to Cluster 1's result
+
+The 21-variable monolithic system (3 free bodies x 7 params) becomes two smaller systems solved in sequence.
+
+### Disconnected sub-assemblies
+
+The decomposition also handles disconnected components. Each connected component of the constraint graph is decomposed independently. Components without a grounded body will fail to solve (returning `NoGroundedParts`).
+
+### Pebble game integration
+
+The `classify_cluster_rigidity()` function uses the pebble game algorithm from `GNN/solver/datagen/` to classify clusters as well-constrained, underconstrained, overconstrained, or mixed. This provides fast O(n) rigidity analysis without running the full solver.
+
+## Half-space tracking
+
+**Source:** `mods/solver/kindred_solver/preference.py`
+
+Many constraints have multiple valid solutions (branches). A distance constraint between two points can be satisfied with the points on either side of each other. Parallel axes can point in the same or opposite directions.
+
+Half-space tracking preserves the initial configuration branch:
+
+1. **At setup:** Evaluate an indicator function for each branching constraint. Record its sign as the reference branch.
+
+2. **After each Newton step:** Re-evaluate the indicator. If the sign flipped, apply a correction to push the solution back to the reference branch.
+
+Tracked constraint types:
+
+| Constraint | Indicator | Correction |
+|-----------|-----------|------------|
+| DistancePointPoint (d > 0) | Dot product of displacement with reference direction | Reflect the moving body's position |
+| Parallel | `z_i . z_j` (same vs. opposite direction) | None (tracked for monitoring) |
+| Angle | Dominant cross product component | None (tracked for monitoring) |
+| Perpendicular | Dominant cross product component | None (tracked for monitoring) |
--- a/docs/src/solver/writing-a-solver.md
+++ b/docs/src/solver/writing-a-solver.md
@@ -0,0 +1,256 @@
+# Writing a Custom Solver
+
+The KCSolve framework lets you implement a solver backend in pure Python, register it at runtime, and select it through the Assembly preferences. This tutorial walks through building a minimal solver and then extending it.
+
+## Minimal solver
+
+A solver must subclass `kcsolve.IKCSolver` and implement three methods:
+
+```python
+import kcsolve
+
+class MySolver(kcsolve.IKCSolver):
+    def __init__(self):
+        super().__init__()  # required for pybind11 trampoline
+
+    def name(self):
+        return "My Custom Solver"
+
+    def supported_joints(self):
+        return [
+            kcsolve.BaseJointKind.Fixed,
+            kcsolve.BaseJointKind.Revolute,
+        ]
+
+    def solve(self, ctx):
+        result = kcsolve.SolveResult()
+
+        # Find grounded parts
+        grounded = {p.id for p in ctx.parts if p.grounded}
+        if not grounded:
+            result.status = kcsolve.SolveStatus.NoGroundedParts
+            return result
+
+        # Your solving logic here...
+        # For each non-grounded part, compute its solved placement
+        for part in ctx.parts:
+            if part.grounded:
+                continue
+            pr = kcsolve.SolveResult.PartResult()
+            pr.id = part.id
+            pr.placement = part.placement  # use current placement as placeholder
+            result.placements = result.placements + [pr]
+
+        result.status = kcsolve.SolveStatus.Success
+        result.dof = 0
+        return result
+```
+
+Register it:
+
+```python
+kcsolve.register_solver("my_solver", MySolver)
+```
+
+Test it from the FreeCAD console:
+
+```python
+solver = kcsolve.load("my_solver")
+print(solver.name())  # "My Custom Solver"
+
+ctx = kcsolve.SolveContext()
+# ... build context ...
+result = solver.solve(ctx)
+print(result.status)  # SolveStatus.Success
+```
+
+## Addon packaging
+
+To make your solver load automatically, create a FreeCAD addon:
+
+```
+my_solver_addon/
+  package.xml          # Addon manifest
+  Init.py              # Registration entry point
+  my_solver/
+    __init__.py
+    solver.py          # MySolver class
+```
+
+**package.xml:**
+```xml
+<?xml version="1.0" encoding="UTF-8"?>
+<package format="1">
+  <name>MyCustomSolver</name>
+  <description>Custom assembly constraint solver</description>
+  <version>0.1.0</version>
+  <content>
+    <preferencepack>
+      <name>MySolver</name>
+    </preferencepack>
+  </content>
+</package>
+```
+
+**Init.py:**
+```python
+import kcsolve
+from my_solver.solver import MySolver
+kcsolve.register_solver("my_solver", MySolver)
+```
+
+Place the addon in the FreeCAD Mod directory or as a git submodule in `mods/`.
+
+## Working with SolveContext
+
+The `SolveContext` contains everything the solver needs:
+
+### Parts
+
+```python
+for part in ctx.parts:
+    print(f"{part.id}: grounded={part.grounded}")
+    print(f"  position: {list(part.placement.position)}")
+    print(f"  quaternion: {list(part.placement.quaternion)}")
+    print(f"  mass: {part.mass}")
+```
+
+Each part has 7 degrees of freedom: 3 translation (x, y, z) and 4 quaternion components (w, x, y, z) with a unit-norm constraint reducing the rotational DOF to 3.
+
+**Quaternion convention:** `(w, x, y, z)` where `w` is the scalar part. This differs from FreeCAD's `Base.Rotation(x, y, z, w)`. The adapter layer handles the swap.
+
+### Constraints
+
+```python
+for c in ctx.constraints:
+    if not c.activated:
+        continue
+    print(f"{c.id}: {c.type} between {c.part_i} and {c.part_j}")
+    print(f"  marker_i: pos={list(c.marker_i.position)}, "
+          f"quat={list(c.marker_i.quaternion)}")
+    print(f"  params: {list(c.params)}")
+    print(f"  limits: {len(c.limits)}")
+```
+
+The marker transforms define local coordinate frames on each part. The constraint type determines what geometric relationship is enforced between these frames.
+
+### Returning results
+
+```python
+result = kcsolve.SolveResult()
+result.status = kcsolve.SolveStatus.Success
+result.dof = computed_dof
+
+placements = []
+for part_id, pos, quat in solved_parts:
+    pr = kcsolve.SolveResult.PartResult()
+    pr.id = part_id
+    pr.placement = kcsolve.Transform()
+    pr.placement.position = list(pos)
+    pr.placement.quaternion = list(quat)
+    placements.append(pr)
+result.placements = placements
+
+return result
+```
+
+**Important:** pybind11 list fields return copies. Use `result.placements = [...]` (whole-list assignment), not `result.placements.append(...)`.
+
+## Adding optional capabilities
+
+### Diagnostics
+
+Override `diagnose()` to detect overconstrained or malformed assemblies:
+
+```python
+def diagnose(self, ctx):
+    diagnostics = []
+    # ... analyze constraints ...
+    d = kcsolve.ConstraintDiagnostic()
+    d.constraint_id = "Joint001"
+    d.kind = kcsolve.DiagnosticKind.Redundant
+    d.detail = "This joint duplicates Joint002"
+    diagnostics.append(d)
+    return diagnostics
+```
+
+### Interactive drag
+
+Override the three drag methods for real-time viewport dragging:
+
+```python
+def pre_drag(self, ctx, drag_parts):
+    self._ctx = ctx
+    self._dragging = set(drag_parts)
+    return self.solve(ctx)
+
+def drag_step(self, drag_placements):
+    # Update dragged parts in stored context
+    for pr in drag_placements:
+        for part in self._ctx.parts:
+            if part.id == pr.id:
+                part.placement = pr.placement
+                break
+    return self.solve(self._ctx)
+
+def post_drag(self):
+    self._ctx = None
+    self._dragging = None
+```
+
+For responsive dragging, the solver should converge quickly from a nearby initial guess. Use warm-starting (current placements as initial guess) and consider caching internal state across drag steps.
+
+### Incremental update
+
+Override `update()` for the case where only constraint parameters changed (not topology):
+
+```python
+def update(self, ctx):
+    # Reuse cached factorization, only re-evaluate changed residuals
+    return self.solve(ctx)  # default: just re-solve
+```
+
+## Testing
+
+### Unit tests (without FreeCAD)
+
+Test your solver logic with hand-built `SolveContext` objects:
+
+```python
+import kcsolve
+
+def test_fixed_joint():
+    ctx = kcsolve.SolveContext()
+
+    base = kcsolve.Part()
+    base.id = "base"
+    base.grounded = True
+
+    arm = kcsolve.Part()
+    arm.id = "arm"
+    arm.placement.position = [100.0, 0.0, 0.0]
+
+    joint = kcsolve.Constraint()
+    joint.id = "Joint001"
+    joint.part_i = "base"
+    joint.part_j = "arm"
+    joint.type = kcsolve.BaseJointKind.Fixed
+
+    ctx.parts = [base, arm]
+    ctx.constraints = [joint]
+
+    solver = MySolver()
+    result = solver.solve(ctx)
+    assert result.status == kcsolve.SolveStatus.Success
+```
+
+### Integration tests (with FreeCAD)
+
+For integration testing within FreeCAD, follow the pattern in `TestKindredSolverIntegration.py`: set the solver preference in `setUp()`, create document objects, and verify solve results.
+
+## Reference
+
+- [KCSolve Python API](../reference/kcsolve-python.md) -- complete type and function reference
+- [KCSolve Architecture](../architecture/ondsel-solver.md) -- C++ framework details
+- [Constraints](constraints.md) -- constraint types and residual counts
+- [Kindred Solver Overview](overview.md) -- how the built-in Kindred solver works