feat(physics): implement cubic barrier potential

ROADMAP.md

@@ -118,7 +118,7 @@
 ## Milestone 0.6.0 – Fold Barrier Physics Suite (Planned)
 - [x] Align Fold barrier scope with the paper and define shared constraint interfaces in `src/physics/fold`
 - **Barrier primitives** (each item: runtime module + `docs/index.d.ts` entry + Vitest coverage + runnable example when feasible)
-  - [ ] Cubic barrier potential (energy, gradient, Hessian evaluation)
+  - [x] Cubic barrier potential (energy, gradient, Hessian evaluation)
   - [ ] Stiffness design principle for frozen barrier stiffness
   - [ ] Contact barrier with extended direction handling
   - [ ] Pin constraint barrier using cubic barrier formulation

docs/index.d.ts

@@ -117,6 +117,7 @@ export const examples: {
   };
   readonly physics: {
     readonly createFoldConstraintRegistry: 'examples/foldSetup.ts';
+    readonly createCubicBarrier: 'examples/foldCubicBarrier.ts';
   };
   readonly performance: {
     readonly debounce: 'examples/requestDedup.ts';
@@ -3299,6 +3300,19 @@ export function closestPair(points: ReadonlyArray<Point>): ClosestPairResult;
  */
 export function createFoldConstraintRegistry(): FoldConstraintRegistry;
 
+/**
+ * Factory for the cubic barrier potential described in Fold.
+ * Use for: enforcing inequality constraints with C^2 continuity.
+ * Import: physics/fold/cubicBarrier.ts
+ */
+export interface CubicBarrierOptions {
+  id?: string;
+  stiffness?: number;
+  maxGap?: number;
+  direction?: Vector3D;
+}
+export function createCubicBarrier(options?: CubicBarrierOptions): FoldConstraint;
+
 export type FoldConstraintType =
   | 'cubic-barrier'
   | 'contact-barrier'

examples/foldCubicBarrier.ts

@@ -0,0 +1,16 @@
+import { createCubicBarrier } from '../src/index.js';
+
+const cubicBarrier = createCubicBarrier({ stiffness: 50 });
+
+const evaluation = cubicBarrier.evaluate(
+  {
+    gap: -0.05,
+    maxGap: 0,
+    stiffness: 30,
+    direction: { x: 0, y: 0, z: 1 },
+  },
+  { deltaTime: 1 / 60 }
+);
+
+console.log('energy', evaluation.energy);
+console.log('gradient', evaluation.gradient);

src/index.ts

@@ -115,6 +115,7 @@ export const examples = {
   },
   physics: {
     createFoldConstraintRegistry: 'examples/foldSetup.ts',
+    createCubicBarrier: 'examples/foldCubicBarrier.ts',
   },
   performance: {
     debounce: 'examples/requestDedup.ts',
@@ -1187,9 +1188,10 @@ export type { ClosestPairResult } from './geometry/closestPair.js';
 // ⚙️ PHYSICS & FOLD BARRIERS
 // ============================================================================
 
-export { createFoldConstraintRegistry } from './physics/fold/index.js';
+export { createFoldConstraintRegistry, createCubicBarrier } from './physics/fold/index.js';
 
 export type {
+  CubicBarrierOptions,
   FoldConstraintType,
   FoldConstraintState,
   FoldComputationContext,

src/physics/fold/cubicBarrier.ts

@@ -0,0 +1,84 @@
+import type { Matrix3x3, Vector3D } from '../../types.js';
+import type {
+  FoldComputationContext,
+  FoldConstraint,
+  FoldConstraintEvaluation,
+  FoldConstraintState,
+} from './types.js';
+
+export interface CubicBarrierOptions {
+  id?: string;
+  stiffness?: number;
+  maxGap?: number;
+  direction?: Vector3D;
+}
+
+const ZERO_GRADIENT: Vector3D = { x: 0, y: 0, z: 0 };
+const ZERO_HESSIAN: Matrix3x3 = [
+  [0, 0, 0],
+  [0, 0, 0],
+  [0, 0, 0],
+];
+
+export function createCubicBarrier(options: CubicBarrierOptions = {}): FoldConstraint {
+  return {
+    type: 'cubic-barrier',
+    id: options.id,
+    enabled: true,
+    evaluate(state: FoldConstraintState, context: FoldComputationContext): FoldConstraintEvaluation {
+      void context;
+      const maxGap = options.maxGap ?? state.maxGap;
+      const stiffness = options.stiffness ?? state.stiffness;
+      const direction = options.direction ?? state.direction;
+
+      if (stiffness <= 0) {
+        return { energy: 0, gradient: ZERO_GRADIENT, hessian: ZERO_HESSIAN };
+      }
+
+      const violation = Math.max(0, maxGap - state.gap);
+      if (violation <= 0) {
+        return { energy: 0, gradient: ZERO_GRADIENT, hessian: ZERO_HESSIAN };
+      }
+
+      const unit = normalise(direction);
+      if (!unit) {
+        return { energy: 0, gradient: ZERO_GRADIENT, hessian: ZERO_HESSIAN };
+      }
+
+      const energy = stiffness * (violation ** 3) / 3;
+      const gradientMagnitude = stiffness * (violation ** 2);
+      const gradient = scaleVector(unit, gradientMagnitude);
+      const hessianMagnitude = stiffness * 2 * violation;
+      const hessian = outerProduct(unit, hessianMagnitude);
+
+      return { energy, gradient, hessian };
+    },
+  };
+}
+
+function normalise(vector: Vector3D | undefined): Vector3D | null {
+  if (!vector) return null;
+  const length = Math.hypot(vector.x, vector.y, vector.z);
+  if (length === 0) return null;
+  return {
+    x: vector.x / length,
+    y: vector.y / length,
+    z: vector.z / length,
+  };
+}
+
+function scaleVector(vector: Vector3D, scalar: number): Vector3D {
+  return {
+    x: vector.x * scalar,
+    y: vector.y * scalar,
+    z: vector.z * scalar,
+  };
+}
+
+function outerProduct(vector: Vector3D, scalar: number): Matrix3x3 {
+  return [
+    [vector.x * vector.x * scalar, vector.x * vector.y * scalar, vector.x * vector.z * scalar],
+    [vector.y * vector.x * scalar, vector.y * vector.y * scalar, vector.y * vector.z * scalar],
+    [vector.z * vector.x * scalar, vector.z * vector.y * scalar, vector.z * vector.z * scalar],
+  ];
+}

src/physics/fold/index.ts

@@ -1 +1,2 @@
 export * from './types.js';
+export * from './cubicBarrier.js';

tests/cubicBarrier.test.ts

@@ -0,0 +1,43 @@
+import { describe, expect, it } from 'vitest';
+
+import { createCubicBarrier } from '../src/physics/fold/cubicBarrier.js';
+
+describe('cubic barrier potential', () => {
+  it('returns zero output when constraint is satisfied', () => {
+    const barrier = createCubicBarrier({ stiffness: 10 });
+    const evaluation = barrier.evaluate(
+      {
+        gap: 0.5,
+        maxGap: 0.5,
+        stiffness: 10,
+        direction: { x: 0, y: 0, z: 1 },
+      },
+      { deltaTime: 1 }
+    );
+
+    expect(evaluation.energy).toBe(0);
+    expect(evaluation.gradient).toEqual({ x: 0, y: 0, z: 0 });
+    expect(evaluation.hessian).toEqual([
+      [0, 0, 0],
+      [0, 0, 0],
+      [0, 0, 0],
+    ]);
+  });
+
+  it('produces cubic energy growth for penetration', () => {
+    const barrier = createCubicBarrier();
+    const evaluation = barrier.evaluate(
+      {
+        gap: -0.1,
+        maxGap: 0,
+        stiffness: 20,
+        direction: { x: 0, y: 0, z: 1 },
+      },
+      { deltaTime: 1 / 60 }
+    );
+
+    expect(evaluation.energy).toBeCloseTo(20 * (0.1 ** 3) / 3, 6);
+    expect(evaluation.gradient.z).toBeGreaterThan(0);
+    expect(evaluation.hessian[2]?.[2]).toBeGreaterThan(0);
+  });
+});

tests/index.test.ts

@@ -58,6 +58,7 @@ describe('package entry point', () => {
     expect(examples.spatial.queryBvh).toBe('examples/bvh.ts');
     expect(examples.spatial.raycastBvh).toBe('examples/bvh.ts');
     expect(examples.physics.createFoldConstraintRegistry).toBe('examples/foldSetup.ts');
+    expect(examples.physics.createCubicBarrier).toBe('examples/foldCubicBarrier.ts');
   });
 
   it('provides strong typing for example categories and names', () => {
@@ -198,7 +199,8 @@ describe('package entry point', () => {
     >();
 
     expectTypeOf<ExampleName<'physics'>>().toEqualTypeOf<
-      'createFoldConstraintRegistry'
+      | 'createFoldConstraintRegistry'
+      | 'createCubicBarrier'
     >();
 
     expectTypeOf<ExampleName<'ai'>>().toEqualTypeOf<