docs(3d): Plan BVH, Physics, and WGSL modularization (Task #49, #50)

- Added detailed design for dynamic BVH and GPU traversal in 3D.md.
- Outlined SDF-based collision detection and semi-implicit Euler physics loop.
- Updated TODO.md with atomic tasks for physics implementation and granular WGSL refactoring.

1 files changed, 136 insertions, 0 deletions

diff --git a/doc/3D.md b/doc/3D.md
index 5c5d1bd..12c7b15 100644
--- a/doc/3D.md
+++ b/doc/3D.md
@@ -38,6 +38,142 @@ removed with STRIP_ALL) that:
 This must be captured and tracked as **Task #39: Visual Debugging System**.
 
 
+## Global Acceleration Structure (BVH)
+
+To support efficient global queries (e.g., for soft shadows via raymarching, or physics), we will implement a dynamic Bounding Volume Hierarchy (BVH).
+
+### Data Structures
+
+**CPU/GPU Common Node Layout:**
+The node structure is designed to be 32-byte aligned for efficient GPU access.
+
+```cpp
+struct BVHNode {
+    // 16 bytes
+    float min_x, min_y, min_z;
+    int32_t left_idx;   // If < 0, this is a leaf node.
+
+    // 16 bytes
+    float max_x, max_y, max_z;
+    int32_t right_idx;  // If leaf, this holds the object_index.
+};
+```
+
+**WGSL Representation:**
+```wgsl
+struct BVHNode {
+    min: vec3<f32>,
+    left_idx: i32,
+    max: vec3<f32>,
+    obj_idx_or_right: i32,
+};
+@group(0) @binding(2) var<storage, read> bvh_nodes: array<BVHNode>;
+```
+
+### Construction (CPU)
+
+Since the scene is small (< 100 objects) and dynamic, we will **rebuild** the BVH from scratch every frame. This avoids the complexity of refitting and balancing.
+
+**Algorithm: Recursive Midpoint Split**
+1.  **Calculate Centroids**: Compute the center of the AABB for all active objects.
+2.  **Split**:
+    *   Find the axis with the largest variance in centroid positions.
+    *   Sort objects along that axis (or just partition based on the median).
+    *   Create a node, assign `left` and `right` children recursively.
+    *   Stop when a node has 1 object (Leaf).
+3.  **Linearize**: Flatten the tree into a `std::vector<BVHNode>` for GPU upload.
+
+### Traversal (GPU/WGSL)
+
+Since WGSL does not support recursion, we use a fixed-size stack.
+
+```wgsl
+fn map_scene(pos: vec3<f32>) -> f32 {
+    var min_dist = 10000.0;
+    var stack: array<i32, 32>;
+    var stack_ptr = 0;
+    
+    // Push root
+    stack[stack_ptr] = 0;
+    stack_ptr++;
+
+    while (stack_ptr > 0) {
+        stack_ptr--;
+        let node_idx = stack[stack_ptr];
+        let node = bvh_nodes[node_idx];
+
+        // Check AABB intersection
+        if (aabb_sdf(pos, node.min, node.max) < min_dist) {
+            if (node.left_idx < 0) {
+                // Leaf: Evaluate Object SDF
+                let obj_idx = node.obj_idx_or_right;
+                let d = object_sdf(pos, obj_idx);
+                min_dist = min(min_dist, d);
+            } else {
+                // Internal: Push children
+                // Optimization: Push furthest child first so we pop closest first
+                stack[stack_ptr] = node.left_idx;
+                stack_ptr++;
+                stack[stack_ptr] = node.obj_idx_or_right;
+                stack_ptr++;
+            }
+        }
+    }
+    return min_dist;
+}
+```
+
+## Physics & Collision System
+
+We need a lightweight physics engine that leverages our SDF representation.
+
+### Components
+We will extend `Object3D` or add a `PhysicsBody` struct:
+*   **Mass/InverseMass**: For dynamic response.
+*   **Velocity**: Linear velocity `vec3`.
+*   **Restitution**: Bounciness (0.0 - 1.0).
+*   **Friction**: Damping factor.
+
+### Broad Phase: BVH Re-use
+We reuse the same BVH constructed for rendering to accelerate CPU-side collision detection.
+*   **Query**: `QueryBVH(AABB query_box)` returns a list of potential candidates.
+*   This avoids $O(N^2)$ checks.
+
+### Narrow Phase: SDF-Based Collision
+Since we don't have explicit meshes for collision, we use the analytical SDFs.
+
+**Algorithm: Proxy Point Probing**
+Instead of full SDF-on-SDF intersection (which is expensive), we approximate dynamic objects as a set of "Probe Points" (e.g., bounding box corners, center).
+
+1.  **Transform**: For a dynamic object $A$ and candidate neighbor $B$:
+    *   Transform $A$'s probe points into $B$'s local space: $P_{local} = InverseModel_B * P_{world}$.
+2.  **Evaluate**: Calculate $d = SDF_B(P_{local})$.
+3.  **Collision**: If $d < 0$:
+    *   **Penetration Depth**: $p = -d$.
+    *   **Normal**: Compute gradient $N = \nabla SDF_B(P_{local})$ via central differences. Transform $N$ back to world space.
+    *   **Contact Point**: $P_{contact} = P_{world} - N * p$.
+
+### Solver Loop (Semi-Implicit Euler)
+1.  **Integrate Velocity**: `vel += gravity * dt`.
+2.  **Broad Phase**: Find pairs using BVH.
+3.  **Narrow Phase & Resolution**:
+    *   If collision detected (depth $p > 0$):
+        *   **Positional Correction**: Move object by $N * p$ (to resolve penetration).
+        *   **Velocity Response**: `vel = vel - (1 + restitution) * dot(vel, N) * N`.
+        *   **Friction**: Apply tangential damping.
+4.  **Integrate Position**: `pos += vel * dt`.
+
+### Code Integration Plan
+
+1.  **CPU-Side SDF Library**:
+    *   Create `src/3d/sdf_cpu.h` implementing the same primitives (`sdSphere`, `sdBox`, `sdTorus`) as the WGSL shaders, using `src/util/mini_math.h`.
+2.  **Physics System Class**:
+    *   Create `src/3d/physics.h/cc`.
+    *   `PhysicsSystem::Update(Scene& scene, float dt)`.
+    *   `PhysicsSystem::ResolveCollision(Object3D& a, Object3D& b)`.
+3.  **Main Loop Integration**:
+    *   Call `PhysicsSystem::Update` before `Renderer3D::render`.
+
 ## future step
 
 Have an exporter from Blender modelling software. That would be a converter