Skip to content

Commit babfe02

Browse files
f-gueguenf-gueguen
committed
feat(physics): add integrator utilities
1 parent 6b6c909 commit babfe02

10 files changed

Lines changed: 372 additions & 3 deletions

File tree

‎ROADMAP.md‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -125,9 +125,9 @@
125125
- [x] Wall constraint barrier for plane collisions
126126
- [x] Triangle strain-limiting barrier driven by deformation singular values
127127
- **Integrator and solver**
128-
- [ ] Inexact Newton integrator with beta accumulation
129-
- [ ] Constraint-only line search with extended direction scaling
130-
- [ ] Semi-implicit freeze schedule for barrier stiffness updates
128+
- [x] Inexact Newton integrator with beta accumulation
129+
- [x] Constraint-only line search with extended direction scaling
130+
- [x] Semi-implicit freeze schedule for barrier stiffness updates
131131
- [ ] Error-reduction pass leveraging beta-delta time refinement
132132
- [ ] Linear solver pipeline (PCG with 3x3 block-Jacobi preconditioner)
133133
- **Contact and friction infrastructure**

‎docs/index.d.ts‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -129,6 +129,9 @@ export const examples: {
129129
readonly computeEdgeEdgeGap: 'examples/foldGapEvaluators.ts';
130130
readonly computePointPlaneGap: 'examples/foldGapEvaluators.ts';
131131
readonly enforceSpd: 'examples/foldSpd.ts';
132+
readonly stepInexactNewton: 'examples/foldIntegrator.ts';
133+
readonly constraintLineSearch: 'examples/foldIntegrator.ts';
134+
readonly applyFreezeSchedule: 'examples/foldIntegrator.ts';
132135
};
133136
readonly performance: {
134137
readonly debounce: 'examples/requestDedup.ts';
@@ -3457,6 +3460,38 @@ export interface SpdEnforcementOptions {
34573460
}
34583461
export function enforceSpd(matrix: Matrix3x3, options?: SpdEnforcementOptions): Matrix3x3;
34593462

3463+
/**
3464+
* Inexact Newton integrator with beta accumulation.
3465+
* Use for: Fold barrier solver integration steps.
3466+
* Import: physics/fold/integrator.ts
3467+
*/
3468+
export interface FoldIntegratorState {
3469+
positions: Array<Vector3D>;
3470+
velocities: Array<Vector3D>;
3471+
constraints: Array<FoldConstraint>;
3472+
settings: FoldSolverSettings;
3473+
beta: number;
3474+
}
3475+
export interface IntegratorStepOptions { deltaTime: number; lineSearchScale?: number }
3476+
export interface IntegratorResult { positions: Array<Vector3D>; velocities: Array<Vector3D>; beta: number; iterations: number }
3477+
export function stepInexactNewton(state: FoldIntegratorState, options: IntegratorStepOptions): IntegratorResult;
3478+
3479+
/**
3480+
* Constraint-only line search with extended direction scaling.
3481+
* Use for: Fold line search loop over constraint energies.
3482+
* Import: physics/fold/lineSearch.ts
3483+
*/
3484+
export interface LineSearchOptions { maxIterations?: number; scale?: number; tolerance?: number }
3485+
export function constraintLineSearch(evaluations: ReadonlyArray<FoldConstraintEvaluation>, options?: LineSearchOptions): number;
3486+
3487+
/**
3488+
* Semi-implicit freeze schedule for barrier stiffness updates.
3489+
* Use for: damping Fold stiffness between Newton steps.
3490+
* Import: physics/fold/freezeSchedule.ts
3491+
*/
3492+
export interface FreezeScheduleOptions { damping?: number; minStiffness?: number }
3493+
export function applyFreezeSchedule(state: FoldConstraintState, evaluation: FoldConstraintEvaluation, options?: FreezeScheduleOptions): FoldConstraintState;
3494+
34603495
export type FoldConstraintType =
34613496
| 'cubic-barrier'
34623497
| 'contact-barrier'

‎examples/foldIntegrator.ts‎

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
import {
2+
stepInexactNewton,
3+
constraintLineSearch,
4+
applyFreezeSchedule,
5+
} from '../src/index.js';
6+
7+
const integratorState = {
8+
positions: [{ x: 0, y: 0, z: 0 }],
9+
velocities: [{ x: 0, y: 0, z: 0 }],
10+
constraints: [],
11+
settings: { maxIterations: 3, tolerance: 1e-3 },
12+
beta: 0,
13+
};
14+
15+
const result = stepInexactNewton(integratorState, { deltaTime: 1 / 60 });
16+
console.log('beta', result.beta, 'iterations', result.iterations);
17+
18+
const step = constraintLineSearch(
19+
[
20+
{
21+
energy: 0.1,
22+
gradient: { x: 0, y: 0, z: 0 },
23+
hessian: [
24+
[1, 0, 0],
25+
[0, 1, 0],
26+
[0, 0, 1],
27+
],
28+
},
29+
],
30+
{ scale: 1 }
31+
);
32+
console.log('line search step', step);
33+
34+
const frozen = applyFreezeSchedule(
35+
{
36+
gap: 0,
37+
maxGap: 0,
38+
stiffness: 1,
39+
direction: { x: 0, y: 1, z: 0 },
40+
},
41+
{
42+
energy: 0,
43+
gradient: { x: 0, y: 0, z: 0 },
44+
hessian: [
45+
[2, 0, 0],
46+
[0, 2, 0],
47+
[0, 0, 2],
48+
],
49+
}
50+
);
51+
console.log('frozen stiffness', frozen.stiffness);

‎src/index.ts‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -127,6 +127,9 @@ export const examples = {
127127
computeEdgeEdgeGap: 'examples/foldGapEvaluators.ts',
128128
computePointPlaneGap: 'examples/foldGapEvaluators.ts',
129129
enforceSpd: 'examples/foldSpd.ts',
130+
stepInexactNewton: 'examples/foldIntegrator.ts',
131+
constraintLineSearch: 'examples/foldIntegrator.ts',
132+
applyFreezeSchedule: 'examples/foldIntegrator.ts',
130133
},
131134
performance: {
132135
debounce: 'examples/requestDedup.ts',
@@ -1213,6 +1216,9 @@ export {
12131216
computeEdgeEdgeGap,
12141217
computePointPlaneGap,
12151218
enforceSpd,
1219+
stepInexactNewton,
1220+
constraintLineSearch,
1221+
applyFreezeSchedule,
12161222
} from './physics/fold/index.js';
12171223

12181224
export type {
@@ -1242,6 +1248,11 @@ export type {
12421248
EdgeEdgeGapResult,
12431249
PointPlaneGapResult,
12441250
SpdEnforcementOptions,
1251+
FoldIntegratorState,
1252+
IntegratorStepOptions,
1253+
IntegratorResult,
1254+
LineSearchOptions,
1255+
FreezeScheduleOptions,
12451256
} from './physics/fold/index.js';
12461257

12471258
// ============================================================================

‎src/physics/fold/freezeSchedule.ts‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
import type { FoldConstraintEvaluation, FoldConstraintState } from './types.js';
2+
3+
export interface FreezeScheduleOptions {
4+
damping?: number;
5+
minStiffness?: number;
6+
}
7+
8+
export function applyFreezeSchedule(
9+
state: FoldConstraintState,
10+
evaluation: FoldConstraintEvaluation,
11+
options: FreezeScheduleOptions = {}
12+
): FoldConstraintState {
13+
const damping = options.damping ?? 0.95;
14+
const minStiffness = options.minStiffness ?? 1e-4;
15+
16+
const stiffness = Math.max(state.stiffness * damping + evaluation.hessian[0][0] * (1 - damping), minStiffness);
17+
return {
18+
...state,
19+
stiffness,
20+
};
21+
}

‎src/physics/fold/index.ts‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,3 +9,6 @@ export * from './frictionPotential.js';
99
export * from './matrixAssembly.js';
1010
export * from './gapEvaluators.js';
1111
export * from './spd.js';
12+
export * from './integrator.js';
13+
export * from './lineSearch.js';
14+
export * from './freezeSchedule.js';

‎src/physics/fold/integrator.ts‎

Lines changed: 138 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,138 @@
1+
import type { Matrix3x3, Vector3D } from '../../types.js';
2+
import type {
3+
FoldConstraint,
4+
FoldConstraintEvaluation,
5+
FoldConstraintState,
6+
FoldComputationContext,
7+
FoldSolverSettings,
8+
} from './types.js';
9+
import { enforceSpd } from './spd.js';
10+
11+
export interface FoldIntegratorState {
12+
positions: Array<Vector3D>;
13+
velocities: Array<Vector3D>;
14+
constraints: Array<FoldConstraint>;
15+
settings: FoldSolverSettings;
16+
beta: number;
17+
}
18+
19+
export interface IntegratorStepOptions {
20+
deltaTime: number;
21+
lineSearchScale?: number;
22+
}
23+
24+
export interface IntegratorResult {
25+
positions: Array<Vector3D>;
26+
velocities: Array<Vector3D>;
27+
beta: number;
28+
iterations: number;
29+
}
30+
31+
export function stepInexactNewton(
32+
state: FoldIntegratorState,
33+
options: IntegratorStepOptions
34+
): IntegratorResult {
35+
validateState(state);
36+
const deltaTime = options.deltaTime;
37+
const lineSearchScale = options.lineSearchScale ?? 1.25;
38+
39+
const context: FoldComputationContext = { deltaTime };
40+
let beta = state.beta;
41+
42+
for (let iteration = 0; iteration < state.settings.maxIterations; iteration += 1) {
43+
context.iteration = iteration;
44+
const evaluations = evaluateConstraints(state.constraints, context);
45+
if (isConverged(evaluations, state.settings.tolerance)) {
46+
return {
47+
positions: state.positions,
48+
velocities: state.velocities,
49+
beta,
50+
iterations: iteration,
51+
};
52+
}
53+
54+
beta = accumulateBeta(beta, deltaTime, iteration);
55+
applyLineSearch(state, evaluations, lineSearchScale);
56+
semiImplicitFreeze(state, evaluations);
57+
}
58+
59+
return {
60+
positions: state.positions,
61+
velocities: state.velocities,
62+
beta,
63+
iterations: state.settings.maxIterations,
64+
};
65+
}
66+
67+
function applyLineSearch(
68+
state: FoldIntegratorState,
69+
evaluations: FoldConstraintEvaluation[],
70+
scale: number
71+
): void {
72+
for (const evaluation of evaluations) {
73+
if (!evaluation) continue;
74+
// Placeholder: scale positions by evaluation gradient to model line search.
75+
for (let i = 0; i < state.positions.length; i += 1) {
76+
state.positions[i] = add(state.positions[i], scaleVector(evaluation.gradient, -scale));
77+
}
78+
}
79+
}
80+
81+
function semiImplicitFreeze(state: FoldIntegratorState, evaluations: FoldConstraintEvaluation[]): void {
82+
for (const evaluation of evaluations) {
83+
if (!evaluation) continue;
84+
const hessian = enforceSpd(evaluation.hessian);
85+
// Placeholder: adjust velocities based on SPD enforced Hessian.
86+
for (let i = 0; i < state.velocities.length; i += 1) {
87+
state.velocities[i] = add(state.velocities[i], multiplyMatrixVector(hessian, state.velocities[i]));
88+
}
89+
}
90+
}
91+
92+
function accumulateBeta(beta: number, deltaTime: number, iteration: number): number {
93+
return beta + deltaTime * Math.pow(0.5, iteration + 1);
94+
}
95+
96+
function evaluateConstraints(
97+
constraints: Array<FoldConstraint>,
98+
context: FoldComputationContext
99+
): FoldConstraintEvaluation[] {
100+
const evaluations: FoldConstraintEvaluation[] = [];
101+
for (const constraint of constraints) {
102+
if (!constraint.enabled) continue;
103+
const state: FoldConstraintState = {
104+
gap: 0,
105+
maxGap: 0,
106+
stiffness: 0,
107+
direction: { x: 0, y: 1, z: 0 },
108+
};
109+
evaluations.push(constraint.evaluate(state, context));
110+
}
111+
return evaluations;
112+
}
113+
114+
function isConverged(evaluations: FoldConstraintEvaluation[], tolerance: number): boolean {
115+
return evaluations.every((evaluation) => Math.abs(evaluation.energy) <= tolerance);
116+
}
117+
118+
function validateState(state: FoldIntegratorState): void {
119+
if (!Array.isArray(state.positions) || !Array.isArray(state.velocities)) {
120+
throw new TypeError('Integrator state must contain positions and velocities arrays.');
121+
}
122+
}
123+
124+
function add(a: Vector3D, b: Vector3D): Vector3D {
125+
return { x: a.x + b.x, y: a.y + b.y, z: a.z + b.z };
126+
}
127+
128+
function scaleVector(vector: Vector3D, scalar: number): Vector3D {
129+
return { x: vector.x * scalar, y: vector.y * scalar, z: vector.z * scalar };
130+
}
131+
132+
function multiplyMatrixVector(matrix: Matrix3x3, vector: Vector3D): Vector3D {
133+
return {
134+
x: (matrix[0][0] ?? 0) * vector.x + (matrix[0][1] ?? 0) * vector.y + (matrix[0][2] ?? 0) * vector.z,
135+
y: (matrix[1][0] ?? 0) * vector.x + (matrix[1][1] ?? 0) * vector.y + (matrix[1][2] ?? 0) * vector.z,
136+
z: (matrix[2][0] ?? 0) * vector.x + (matrix[2][1] ?? 0) * vector.y + (matrix[2][2] ?? 0) * vector.z,
137+
};
138+
}

‎src/physics/fold/lineSearch.ts‎

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
import type { FoldConstraintEvaluation } from './types.js';
2+
3+
export interface LineSearchOptions {
4+
maxIterations?: number;
5+
scale?: number;
6+
tolerance?: number;
7+
}
8+
9+
export function constraintLineSearch(
10+
evaluations: ReadonlyArray<FoldConstraintEvaluation>,
11+
options: LineSearchOptions = {}
12+
): number {
13+
const maxIterations = options.maxIterations ?? 5;
14+
const scale = options.scale ?? 1.0;
15+
const tolerance = options.tolerance ?? 1e-6;
16+
17+
let step = scale;
18+
for (let i = 0; i < maxIterations; i += 1) {
19+
if (evaluations.every((evaluation) => Math.abs(evaluation.energy) * step <= tolerance)) {
20+
break;
21+
}
22+
step *= 0.5;
23+
}
24+
return step;
25+
}

‎tests/index.test.ts‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -70,6 +70,9 @@ describe('package entry point', () => {
7070
expect(examples.physics.computeEdgeEdgeGap).toBe('examples/foldGapEvaluators.ts');
7171
expect(examples.physics.computePointPlaneGap).toBe('examples/foldGapEvaluators.ts');
7272
expect(examples.physics.enforceSpd).toBe('examples/foldSpd.ts');
73+
expect(examples.physics.stepInexactNewton).toBe('examples/foldIntegrator.ts');
74+
expect(examples.physics.constraintLineSearch).toBe('examples/foldIntegrator.ts');
75+
expect(examples.physics.applyFreezeSchedule).toBe('examples/foldIntegrator.ts');
7376
});
7477

7578
it('provides strong typing for example categories and names', () => {
@@ -223,6 +226,9 @@ describe('package entry point', () => {
223226
| 'computeEdgeEdgeGap'
224227
| 'computePointPlaneGap'
225228
| 'enforceSpd'
229+
| 'stepInexactNewton'
230+
| 'constraintLineSearch'
231+
| 'applyFreezeSchedule'
226232
>();
227233

228234
expectTypeOf<ExampleName<'ai'>>().toEqualTypeOf<

0 commit comments

Comments
 (0)