carapace_wasm/input_v1/
session.rs

1//! `Session`/`StageRunner`/`advance` — `Session` is chosen once at decode time from the
2//! header's `ndm` discriminant and never branches on it again; stepping
3//! is driven by the caller via a step budget rather than run to completion
4//! inside one call, the mechanism cooperative cancellation is built on.
5//!
6//! `ModelSession<NDIM, NDOF, NId, E>` is generic over the same
7//! kinematic profile `core`'s own `Domain`/`Analysis` are generic over
8//! (`AnalysisBuilder::build` already infers that whole profile from its
9//! `domain` argument — see its doc comment) — the 2D/3D split only
10//! actually differs in element physics and DTO shapes (`decode/elements_2d.rs`/
11//! `elements_3d.rs`), never in this stepping/
12//! recording bookkeeping, so `Session2`/`Session3` are just two
13//! instantiations of one implementation rather than two hand-duplicated
14//! ones.
15
16use carapace_core::analysis::{
17    modal_analysis, Algorithm, Analysis, AnalysisBuilder, AnalysisError, ArcFailure, ArcStepInfo,
18    ConstraintHandler, ConvergenceTest, GroundMotion, Integrator, Mode, RayleighDamping,
19    StopReason, TransientAnalysis,
20};
21use carapace_core::model::{
22    Domain, Element, Element3, ElementOps, LoadPatternId, Node3Id, NodeId, NDF, PLANAR_NDIM,
23    SPATIAL_NDF, SPATIAL_NDIM,
24};
25use serde::Serialize;
26use slotmap::Key;
27
28use super::error::ModelErrorDetail;
29use super::sequence::{FiberResponseKind, GaussQuantity, ModalDampingAnchors};
30
31/// Shared by `record_sample`'s `Static`/`Transient` arms: resolves a
32/// `ResolvedRecorder::Fiber` against whatever `element`'s current fiber
33/// state is, or `None` if that element has no fiber section at all
34/// (`Domain::element_fiber_responses`'s doc comment) or `point`/`fiber`
35/// is out of range for it — the same graceful-skip handling `ModeShape`'s
36/// out-of-range `mode` gets, since neither can be checked before the
37/// analysis actually runs.
38fn fiber_value<const NDIM: usize, const NDOF: usize, NId, E>(
39    domain: &Domain<NDIM, NDOF, NId, E>,
40    element: E::Id,
41    point: u32,
42    fiber: u32,
43    response: FiberResponseKind,
44) -> Option<f64>
45where
46    NId: Key,
47    E: ElementOps<NDIM, NDOF, NId>,
48{
49    let responses = domain.element_fiber_responses(element)?;
50    let &(strain, stress) = responses.get(point as usize)?.get(fiber as usize)?;
51    Some(match response {
52        FiberResponseKind::Strain => strain,
53        FiberResponseKind::Stress => stress,
54    })
55}
56
57fn gauss_value<const NDIM: usize, const NDOF: usize, NId, E>(
58    domain: &Domain<NDIM, NDOF, NId, E>,
59    element: E::Id,
60    point: u32,
61    quantity: GaussQuantity,
62    component: u8,
63) -> Option<f64>
64where
65    NId: Key,
66    E: ElementOps<NDIM, NDOF, NId>,
67{
68    if let Some(responses) = domain.element_shell_responses(element) {
69        let response = responses.get(point as usize)?;
70        let values = match quantity {
71            GaussQuantity::Strain => response.strain,
72            GaussQuantity::Stress => response.resultant,
73        };
74        return values.get(component as usize).copied();
75    }
76    let response = *domain
77        .element_gauss_responses(element)?
78        .get(point as usize)?;
79    let values = match quantity {
80        GaussQuantity::Strain => response.strain,
81        GaussQuantity::Stress => response.stress,
82    };
83    values.get(component as usize).copied()
84}
85
86/// `AnalysisError`'s fields, restated so `advance`'s result doesn't need to
87/// name `carapace_core`'s error type directly — kept in the same tagged-
88/// variant style.
89#[derive(Debug, Clone, Copy, PartialEq, Serialize, tsify::Tsify)]
90#[serde(
91    tag = "kind",
92    rename_all = "camelCase",
93    rename_all_fields = "camelCase"
94)]
95pub enum AnalysisErrorDetail {
96    FailedToConverge {
97        step: usize,
98    },
99    SingularSystem,
100    InvalidConstraint,
101    InvalidModeCount {
102        requested: usize,
103        free_dofs: usize,
104    },
105    /// A configuration value is out of range or unsupported in combination
106    /// (`field` is its wire spelling). Decode rejects these up front; this
107    /// covers anything only detectable once the model is known (e.g. a
108    /// missing rotation scale for a model with rotational DOFs).
109    InvalidOption {
110        field: &'static str,
111    },
112    /// Arc length found an unfrozen `path` load series.
113    UnsupportedLoadSeries,
114    /// Arc length has no reference load to continue.
115    ZeroLoadSensitivity,
116    /// Arc length must start from equilibrium; `measure` is the force-test
117    /// measure at entry.
118    InitialStateNotInEquilibrium {
119        measure: f64,
120    },
121    /// The first arc-length tangent is singular and no seed was given.
122    MissingSeedDirection,
123    /// Every arc-length attempt for `step` failed.
124    CutbacksExhausted {
125        step: usize,
126        attempts: usize,
127        radius: f64,
128        last_failure: ArcFailureDetail,
129    },
130    /// An arc-length stop criterion already ended the stage.
131    ContinuationComplete,
132    /// The model failed `core`'s `Domain::validate` when the stage started,
133    /// for example a nodal load registered by this stage on a DOF nothing
134    /// uses.
135    InvalidModel {
136        error: ModelErrorDetail,
137    },
138}
139
140/// `core::ArcFailure`, restated for the wire.
141#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, tsify::Tsify)]
142#[serde(rename_all = "camelCase")]
143pub enum ArcFailureDetail {
144    NotConverged,
145    Stagnated,
146    NoDescent,
147    NonFinite,
148    SingularSystem,
149    InaccurateSolve,
150    BranchOrientation,
151    DegenerateConstraint,
152}
153
154impl From<ArcFailure> for ArcFailureDetail {
155    fn from(failure: ArcFailure) -> Self {
156        match failure {
157            ArcFailure::NotConverged => Self::NotConverged,
158            ArcFailure::Stagnated => Self::Stagnated,
159            ArcFailure::NoDescent => Self::NoDescent,
160            ArcFailure::NonFinite => Self::NonFinite,
161            ArcFailure::SingularSystem => Self::SingularSystem,
162            ArcFailure::InaccurateSolve => Self::InaccurateSolve,
163            ArcFailure::BranchOrientation => Self::BranchOrientation,
164            ArcFailure::DegenerateConstraint => Self::DegenerateConstraint,
165        }
166    }
167}
168
169/// `core::ArcStepInfo` for the last accepted arc-length step of an
170/// `advance()` call (`StepOutcome::continuation`).
171#[derive(Debug, Clone, Copy, PartialEq, Serialize, tsify::Tsify)]
172#[serde(rename_all = "camelCase")]
173pub struct ContinuationDetail {
174    pub radius: f64,
175    pub next_radius: f64,
176    pub retries: usize,
177    pub corrector_solves: usize,
178    pub factorizations: usize,
179    pub force_measure: f64,
180    pub arc_error: f64,
181    pub chord_length: f64,
182    pub det_sign: Option<i8>,
183    pub bifurcation_suspected: bool,
184    pub displacement_scale: f64,
185    pub rotation_scale: Option<f64>,
186    pub stop: Option<ContinuationStopDetail>,
187}
188
189#[derive(Debug, Clone, Copy, PartialEq, Serialize, tsify::Tsify)]
190#[serde(rename_all = "camelCase")]
191pub struct ContinuationStopDetail {
192    pub reason: StopReasonDetail,
193    pub landed_exactly: bool,
194    pub overshoot: f64,
195}
196
197#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, tsify::Tsify)]
198#[serde(rename_all = "camelCase")]
199pub enum StopReasonDetail {
200    DisplacementTarget,
201    LoadFactorTarget,
202    LoadFactorZeroCrossing,
203    ChordLength,
204    StepCount,
205}
206
207impl From<ArcStepInfo> for ContinuationDetail {
208    fn from(info: ArcStepInfo) -> Self {
209        Self {
210            radius: info.radius,
211            next_radius: info.next_radius,
212            retries: info.retries,
213            corrector_solves: info.corrector_solves,
214            factorizations: info.factorizations,
215            force_measure: info.force_measure,
216            arc_error: info.arc_error,
217            chord_length: info.chord_length,
218            det_sign: info.det_sign,
219            bifurcation_suspected: info.bifurcation_suspected,
220            displacement_scale: info.displacement_scale,
221            rotation_scale: info.rotation_scale,
222            stop: info.stop.map(|stop| ContinuationStopDetail {
223                reason: match stop.reason {
224                    StopReason::DisplacementTarget => StopReasonDetail::DisplacementTarget,
225                    StopReason::LoadFactorTarget => StopReasonDetail::LoadFactorTarget,
226                    StopReason::LoadFactorZeroCrossing => StopReasonDetail::LoadFactorZeroCrossing,
227                    StopReason::ChordLength => StopReasonDetail::ChordLength,
228                    StopReason::StepCount => StopReasonDetail::StepCount,
229                },
230                landed_exactly: stop.landed_exactly,
231                overshoot: stop.overshoot,
232            }),
233        }
234    }
235}
236
237impl From<AnalysisError> for AnalysisErrorDetail {
238    fn from(error: AnalysisError) -> Self {
239        match error {
240            AnalysisError::FailedToConverge { step } => {
241                AnalysisErrorDetail::FailedToConverge { step }
242            }
243            AnalysisError::SingularSystem => AnalysisErrorDetail::SingularSystem,
244            AnalysisError::InvalidConstraint => AnalysisErrorDetail::InvalidConstraint,
245            AnalysisError::InvalidModeCount {
246                requested,
247                free_dofs,
248            } => AnalysisErrorDetail::InvalidModeCount {
249                requested,
250                free_dofs,
251            },
252            AnalysisError::InvalidOption { field } => AnalysisErrorDetail::InvalidOption { field },
253            AnalysisError::UnsupportedLoadSeries => AnalysisErrorDetail::UnsupportedLoadSeries,
254            AnalysisError::ZeroLoadSensitivity => AnalysisErrorDetail::ZeroLoadSensitivity,
255            AnalysisError::InitialStateNotInEquilibrium { measure } => {
256                AnalysisErrorDetail::InitialStateNotInEquilibrium { measure }
257            }
258            AnalysisError::MissingSeedDirection => AnalysisErrorDetail::MissingSeedDirection,
259            AnalysisError::CutbacksExhausted {
260                step,
261                attempts,
262                radius,
263                last_failure,
264            } => AnalysisErrorDetail::CutbacksExhausted {
265                step,
266                attempts,
267                radius,
268                last_failure: last_failure.into(),
269            },
270            AnalysisError::ContinuationComplete => AnalysisErrorDetail::ContinuationComplete,
271            AnalysisError::InvalidModel(error) => AnalysisErrorDetail::InvalidModel {
272                error: error.into(),
273            },
274        }
275    }
276}
277
278/// A `RecorderSpec` resolved against a decoded `Domain` —
279/// real node/element ID handles instead of wire-format table indices.
280/// Generic over the profile's ID types so one definition serves both
281/// `Session2` (`NodeId`/`ElementId`) and `Session3`
282/// (`Node3Id`/`Element3Id`).
283#[derive(Debug, Clone, Copy, PartialEq)]
284pub enum ResolvedRecorder<NId, EId> {
285    NodeDisp {
286        node: NId,
287        dof: u8,
288    },
289    /// Only meaningful during a `Transient` stage (`RecorderSpec::NodeVel`'s
290    /// doc comment) — a `Static`/`Modal` stage's `record_sample` skips it.
291    NodeVel {
292        node: NId,
293        dof: u8,
294    },
295    NodeAccel {
296        node: NId,
297        dof: u8,
298    },
299    ElementForce {
300        element: EId,
301        component: u8,
302    },
303    /// See `RecorderSpec::ElementLoad`; skipped during `Modal`.
304    ElementLoad {
305        element: EId,
306        component: u8,
307    },
308    /// Only meaningful during a `Modal` stage (`RecorderSpec::ModeShape`'s
309    /// doc comment).
310    ModeShape {
311        mode: u32,
312        node: NId,
313        dof: u8,
314    },
315    /// Valid during `Static`/`Transient` (`RecorderSpec::Reaction`'s doc
316    /// comment); skipped during `Modal`.
317    Reaction {
318        node: NId,
319        dof: u8,
320    },
321    /// Valid during `Static`/`Transient` (`RecorderSpec::Fiber`'s doc
322    /// comment); skipped during `Modal`.
323    Fiber {
324        element: EId,
325        point: u32,
326        fiber: u32,
327        response: super::sequence::FiberResponseKind,
328    },
329    /// Valid during `Static`/`Transient` (`RecorderSpec::GaussPoint`'s doc
330    /// comment); skipped during `Modal`.
331    GaussPoint {
332        element: EId,
333        point: u32,
334        quantity: super::sequence::GaussQuantity,
335        component: u8,
336    },
337}
338
339pub(super) struct CompiledStage<NId, EId, Load> {
340    pub id: String,
341    pub steps: u32,
342    pub kind: CompiledStageKind<NId>,
343    /// Registered on the `Domain` only when this stage starts — see
344    /// `tables::NodalLoadTable::stage`.
345    pub pending_nodal_loads: Vec<(LoadPatternId, NId, usize, f64)>,
346    pub pending_element_loads: Vec<(LoadPatternId, EId, Load)>,
347}
348
349pub(super) enum CompiledStageKind<NId> {
350    Static {
351        integrator: Integrator<NId>,
352        algorithm: Algorithm,
353        convergence: ConvergenceTest,
354        hold_patterns_after: Vec<LoadPatternId>,
355    },
356    /// A single eigensolve, not an iterative loop — `CompiledStage::steps`
357    /// is always `1` for this kind (`decode.rs`'s `compile_stages`).
358    Modal { num_modes: usize },
359    /// Swaps the domain for its as-built copy (`ModelSession::pristine`); `steps` is `0`.
360    Reset,
361    Transient {
362        algorithm: Algorithm,
363        convergence: ConvergenceTest,
364        damping: RayleighDamping,
365        modal_anchors: Option<ModalDampingAnchors>,
366        dt: f64,
367        ground_motions: Vec<GroundMotion>,
368    },
369}
370
371enum StageRunner<const NDIM: usize, const NDOF: usize, NId, E>
372where
373    NId: Key + Copy,
374    E: ElementOps<NDIM, NDOF, NId> + Clone,
375    E::Load: Clone,
376{
377    Static {
378        analysis: Analysis<NDIM, NDOF, NId, E>,
379        steps_remaining: u32,
380    },
381    /// `modal_analysis` already ran (in `start_current_stage`) by the time
382    /// this variant exists — there's no `core` type to iteratively step
383    /// the way `Analysis`/`TransientAnalysis` are, so this just holds the
384    /// domain (handed back to `finish_current_stage`, exactly like
385    /// `Analysis::into_domain` does for the other two kinds) alongside the
386    /// already-computed modes.
387    Modal {
388        domain: Domain<NDIM, NDOF, NId, E>,
389        modes: Vec<Mode>,
390        steps_remaining: u32,
391    },
392    Transient {
393        analysis: TransientAnalysis<NDIM, NDOF, NId, E>,
394        steps_remaining: u32,
395    },
396    /// The restored domain, handed to the next stage by `finish_current_stage`.
397    Reset {
398        domain: Domain<NDIM, NDOF, NId, E>,
399        steps_remaining: u32,
400    },
401}
402
403/// `advance`'s result — `{ done, stageComplete, stepsTaken,
404/// progressSnapshot, recorderBatch? }`, minus the parts that
405/// are the worker/JS boundary's job (`progressSnapshot`'s throttling,
406/// `recorderBatch`'s `response_blocks` byte layout): `load_factor` here is
407/// the raw signal those would be built from — for a `Static` stage the
408/// integrator's load factor, for a `Transient` stage the elapsed time
409/// (`TransientAnalysis::time`), and for a `Modal` stage always `0.0` (a
410/// single eigensolve has no comparable incremental progress scalar; read
411/// `Mode::frequency` from a `ModeShape` recorder's batch instead). `error`,
412/// once set, is sticky — later stages are not attempted, matching "a
413/// stage's `AnalysisError` stops the sequence".
414#[derive(Debug, Clone, PartialEq, Serialize, tsify::Tsify)]
415#[serde(rename_all = "camelCase")]
416pub struct StepOutcome {
417    pub done: bool,
418    pub stage_complete: bool,
419    pub steps_taken: u32,
420    pub load_factor: f64,
421    pub error: Option<AnalysisErrorDetail>,
422    /// Diagnostics of the last accepted arc-length step this call took
423    /// (radius, retries, solves, chord length, determinant sign, stop
424    /// criterion), or `None` when it took none. With arc length the load
425    /// factor can decrease and repeat, so samples are ordered by their
426    /// sample index, never by load factor.
427    pub continuation: Option<ContinuationDetail>,
428    /// Only the samples *this* `advance()` call produced, one entry per recorder that recorded
429    /// at least one sample this call (every recorder samples every step today, so in practice
430    /// this is either empty — no step taken — or has one entry per recorder). This is the
431    /// results-storage plan's (docs/results-storage-indexeddb.md) `recorderBatch`, restated in
432    /// field names that match pysees's frozen `ResultBlock` protocol type
433    /// (`src/app/types/resultsStorage.ts`); the caller assigns `blockIndex` and packs `data`,
434    /// since neither concept exists on this side of the wasm boundary.
435    pub recorder_batches: Vec<RecorderBatch>,
436}
437
438/// One computed mode of a finished `Modal` stage. `shape` is node-major over every node in the
439/// input's node-table order (`shape[node * ndf + dof]`, fixed DOFs `0.0`), M-normalized and
440/// sign-fixed so its largest-magnitude entry is positive. `participation[d]` is the modal
441/// participation factor `phi^T M r_d` for global translation direction `d`, and `mass_ratio[d]`
442/// the effective modal mass `participation^2` as a fraction of `ModalStageResult::total_mass[d]`.
443#[derive(Debug, Clone, PartialEq, Serialize, tsify::Tsify)]
444#[serde(rename_all = "camelCase")]
445pub struct ModeResult {
446    /// Natural circular frequency, rad/time.
447    pub frequency: f64,
448    pub shape: Vec<f64>,
449    pub participation: Vec<f64>,
450    pub mass_ratio: Vec<f64>,
451}
452
453/// Everything a finished `Modal` stage computed.
454#[derive(Debug, Clone, PartialEq, Serialize, tsify::Tsify)]
455#[serde(rename_all = "camelCase")]
456pub struct ModalStageResult {
457    pub stage_index: usize,
458    pub stage_id: String,
459    /// DOFs per node (the stride of `ModeResult::shape`).
460    pub ndf: usize,
461    pub modes: Vec<ModeResult>,
462    /// `r_d^T M r_d` per translation direction `d`: the mass the mass ratios are relative to.
463    pub total_mass: Vec<f64>,
464}
465
466/// `Session::modal_results`'s payload: every `Modal` stage that has finished so far, in order.
467#[derive(Debug, Clone, PartialEq, Serialize, tsify::Tsify)]
468#[serde(rename_all = "camelCase")]
469pub struct ModalResultsReport {
470    pub stages: Vec<ModalStageResult>,
471}
472
473/// One recorder's new samples from a single `advance()` call. `first_sample` plus
474/// `samples.len()` gives the sample range this batch covers, so a caller can retry a failed
475/// persist without renumbering — the "each `advance()` batch deterministic and identifies its
476/// first sample/count" contract from results-storage-indexeddb.md's "Failure and cancellation
477/// contract". Scalar-only (`(pseudo_time, value)` pairs, one channel) until a vector recorder
478/// needs more.
479#[derive(Debug, Clone, PartialEq, Serialize, tsify::Tsify)]
480#[serde(rename_all = "camelCase")]
481pub struct RecorderBatch {
482    pub recorder_index: usize,
483    pub stage_index: usize,
484    pub first_sample: u32,
485    pub samples: Vec<(f64, f64)>,
486}
487
488/// One opaque handle exposed to the caller, chosen once from `ndm` at
489/// decode time.
490#[derive(Debug)]
491pub enum Session {
492    D2(Session2),
493    D3(Session3),
494}
495
496impl Session {
497    pub fn advance(&mut self, step_budget: u32) -> StepOutcome {
498        match self {
499            Session::D2(session) => session.advance(step_budget),
500            Session::D3(session) => session.advance(step_budget),
501        }
502    }
503
504    /// Every `Modal` stage finished so far (see `ModalStageResult`).
505    pub fn modal_results(&self) -> ModalResultsReport {
506        match self {
507            Session::D2(session) => session.modal_results(),
508            Session::D3(session) => session.modal_results(),
509        }
510    }
511
512    /// The `AnalysisSequence` stage id `advance` is currently in (or, once
513    /// `done`, the id of whichever stage stopped it) — for surfacing which
514    /// stage an `error` belongs to, so a failed run can be marked with the error's
515    /// structured detail.
516    pub fn current_stage_id(&self) -> Option<&str> {
517        match self {
518            Session::D2(session) => session.current_stage_id(),
519            Session::D3(session) => session.current_stage_id(),
520        }
521    }
522}
523
524pub type Session2 = ModelSession<PLANAR_NDIM, NDF, NodeId, Element>;
525pub type Session3 = ModelSession<SPATIAL_NDIM, SPATIAL_NDF, Node3Id, Element3>;
526
527pub struct ModelSession<const NDIM: usize, const NDOF: usize, NId, E>
528where
529    NId: Key + Copy,
530    E: ElementOps<NDIM, NDOF, NId> + Clone,
531    E::Load: Clone,
532{
533    domain: Option<Domain<NDIM, NDOF, NId, E>>,
534    // Node handles in input node-table order, for `ModeResult::shape`.
535    node_ids: Vec<NId>,
536    modal_results: Vec<ModalStageResult>,
537    // The domain as built, kept only when a `Reset` stage exists.
538    pristine: Option<Domain<NDIM, NDOF, NId, E>>,
539    // Every `holdPatternsAfter` freeze applied so far, in order, so a `Reset` can re-apply them.
540    holds: Vec<(LoadPatternId, f64)>,
541    record_initial: bool,
542    stages: Vec<CompiledStage<NId, E::Id, E::Load>>,
543    current_stage: usize,
544    runner: Option<StageRunner<NDIM, NDOF, NId, E>>,
545    load_factor: f64,
546    error: Option<AnalysisErrorDetail>,
547    recorders: Vec<ResolvedRecorder<NId, E::Id>>,
548    // Per-advance only: the last arc-length step's diagnostics.
549    continuation: Option<ContinuationDetail>,
550    // Per-advance only: cleared at the top of every `advance()` call, drained into that call's
551    // `StepOutcome::recorder_batches` at the end. Nothing here survives across calls except the
552    // running counts below — this is results-storage-indexeddb.md's "Bound memory" step:
553    // `recorder_history`'s old whole-run `Vec<Vec<(f64, f64)>>` is gone.
554    current_batch: Vec<Vec<(f64, f64)>>,
555    // Total samples emitted so far per recorder, so the next batch's `first_sample` is correct
556    // without retaining the samples themselves.
557    recorder_sample_counts: Vec<u32>,
558}
559
560impl<const NDIM: usize, const NDOF: usize, NId, E> std::fmt::Debug
561    for ModelSession<NDIM, NDOF, NId, E>
562where
563    NId: Key + Copy,
564    E: ElementOps<NDIM, NDOF, NId> + Clone,
565    E::Load: Clone,
566{
567    // `Analysis` (inside `StageRunner`, held via `runner`/`domain`) doesn't
568    // derive `Debug`, so this reports the fields useful for a test/error
569    // message rather than the full solver state.
570    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
571        f.debug_struct("ModelSession")
572            .field("current_stage", &self.current_stage)
573            .field("load_factor", &self.load_factor)
574            .field("error", &self.error)
575            .finish()
576    }
577}
578
579impl<const NDIM: usize, const NDOF: usize, NId, E> ModelSession<NDIM, NDOF, NId, E>
580where
581    NId: Key + Copy,
582    E: ElementOps<NDIM, NDOF, NId> + Clone,
583    E::Load: Clone,
584{
585    pub(super) fn new(
586        domain: Domain<NDIM, NDOF, NId, E>,
587        node_ids: Vec<NId>,
588        stages: Vec<CompiledStage<NId, E::Id, E::Load>>,
589        recorders: Vec<ResolvedRecorder<NId, E::Id>>,
590        record_initial: bool,
591    ) -> Self {
592        let pristine = stages
593            .iter()
594            .any(|stage| matches!(stage.kind, CompiledStageKind::Reset))
595            .then(|| domain.clone());
596        let current_batch = vec![Vec::new(); recorders.len()];
597        let recorder_sample_counts = vec![0; recorders.len()];
598        Self {
599            domain: Some(domain),
600            node_ids,
601            modal_results: Vec::new(),
602            pristine,
603            holds: Vec::new(),
604            record_initial,
605            stages,
606            current_stage: 0,
607            runner: None,
608            load_factor: 0.0,
609            error: None,
610            recorders,
611            continuation: None,
612            current_batch,
613            recorder_sample_counts,
614        }
615    }
616
617    pub fn modal_results(&self) -> ModalResultsReport {
618        ModalResultsReport {
619            stages: self.modal_results.clone(),
620        }
621    }
622
623    pub fn current_stage_id(&self) -> Option<&str> {
624        self.stages
625            .get(self.current_stage)
626            .map(|stage| stage.id.as_str())
627    }
628
629    /// Steps up to `step_budget` times, stopping early if the current stage
630    /// completes or a step fails. Each call only ever executes up to
631    /// `step_budget` steps before returning — the mechanism a worker's
632    /// cooperative-cancellation loop is built on (checking a pending cancel
633    /// between calls), even though there is no such loop calling this yet.
634    pub fn advance(&mut self, step_budget: u32) -> StepOutcome {
635        for batch in &mut self.current_batch {
636            batch.clear();
637        }
638        self.continuation = None;
639
640        if self.error.is_none() && self.runner.is_none() && self.current_stage < self.stages.len() {
641            self.start_current_stage();
642        }
643        if self.error.is_some() || self.current_stage >= self.stages.len() {
644            return StepOutcome {
645                done: true,
646                stage_complete: false,
647                steps_taken: 0,
648                load_factor: self.load_factor,
649                error: self.error,
650                continuation: None,
651                recorder_batches: Vec::new(),
652            };
653        }
654
655        // Captured before stepping: `finish_current_stage` below advances `current_stage`, but
656        // every sample this call records (the stepping loop never starts a new stage mid-call)
657        // belongs to the stage active when the call began.
658        let stage_index_for_batch = self.current_stage;
659
660        let mut steps_taken = 0;
661        let mut stage_complete = false;
662        while steps_taken < step_budget {
663            if self.steps_remaining() == 0 {
664                stage_complete = true;
665                break;
666            }
667            match self.step_once() {
668                Ok(()) => {
669                    steps_taken += 1;
670                    if self.steps_remaining() == 0 {
671                        stage_complete = true;
672                        break;
673                    }
674                }
675                Err(error) => {
676                    self.error = Some(error.into());
677                    break;
678                }
679            }
680        }
681
682        if stage_complete {
683            self.finish_current_stage();
684        }
685
686        let recorder_batches = self
687            .current_batch
688            .iter()
689            .zip(self.recorder_sample_counts.iter_mut())
690            .enumerate()
691            .filter_map(|(recorder_index, (samples, sample_count))| {
692                if samples.is_empty() {
693                    return None;
694                }
695                let first_sample = *sample_count;
696                *sample_count += samples.len() as u32;
697                Some(RecorderBatch {
698                    recorder_index,
699                    stage_index: stage_index_for_batch,
700                    first_sample,
701                    samples: samples.clone(),
702                })
703            })
704            .collect();
705
706        StepOutcome {
707            done: self.error.is_some() || self.current_stage >= self.stages.len(),
708            stage_complete,
709            steps_taken,
710            load_factor: self.load_factor,
711            error: self.error,
712            continuation: self.continuation,
713            recorder_batches,
714        }
715    }
716
717    fn start_current_stage(&mut self) {
718        let stage = &self.stages[self.current_stage];
719        let mut domain = self.domain.take().expect("domain held between stages");
720        if matches!(stage.kind, CompiledStageKind::Reset) {
721            // As-built domain plus every load an earlier stage had registered (loads are
722            // registered when their stage starts, so the pristine copy has none), and the freezes
723            // (`loadConst`) already applied: like OpenSees' `reset`, this reverts the state but not
724            // which patterns are held constant, so a held pattern is applied in full from the
725            // first step after the reset.
726            domain = self.pristine.clone().expect("kept for Reset stages");
727            for earlier in &self.stages[..self.current_stage] {
728                for &(pattern, node, dof, value) in &earlier.pending_nodal_loads {
729                    domain.add_nodal_load(pattern, node, dof, value);
730                }
731                for (pattern, element, load) in &earlier.pending_element_loads {
732                    domain.add_element_load(*pattern, *element, load.clone());
733                }
734            }
735            for &(pattern, factor) in &self.holds {
736                domain.hold_pattern_constant(pattern, factor);
737            }
738            self.runner = Some(StageRunner::Reset {
739                domain,
740                steps_remaining: 0,
741            });
742            self.load_factor = 0.0;
743            return;
744        }
745        for &(pattern, node, dof, value) in &stage.pending_nodal_loads {
746            domain.add_nodal_load(pattern, node, dof, value);
747        }
748        for (pattern, element, load) in &stage.pending_element_loads {
749            domain.add_element_load(*pattern, *element, load.clone());
750        }
751        // `Transformation` resolves `equal_dof`/`rigid_diaphragm` ties;
752        // `Plain` is a strictly cheaper no-op for a domain with none (and
753        // `AnalysisBuilder::build` would reject `Plain` outright against a
754        // domain that does have them — see its own doc comment) — so this
755        // just asks the domain which one it actually needs rather than the
756        // caller ever choosing.
757        let constraint_handler = if domain.has_mp_constraints() {
758            ConstraintHandler::Transformation
759        } else {
760            ConstraintHandler::Plain
761        };
762        match &stage.kind {
763            CompiledStageKind::Static {
764                integrator,
765                algorithm,
766                convergence,
767                ..
768            } => {
769                match AnalysisBuilder::new()
770                    .constraint_handler(constraint_handler)
771                    .integrator(integrator.clone())
772                    .algorithm(*algorithm)
773                    .test(*convergence)
774                    .try_build(domain)
775                {
776                    Ok(analysis) => {
777                        self.runner = Some(StageRunner::Static {
778                            analysis,
779                            steps_remaining: stage.steps,
780                        });
781                    }
782                    Err(error) => self.error = Some(AnalysisError::InvalidModel(error).into()),
783                }
784            }
785            // The eigensolve runs immediately, here, rather than lazily on
786            // the first `step_once` call — there's no cheaper "start" step
787            // to defer it to, unlike `Static`/`Transient`'s first Newton/
788            // Newmark solve, which naturally happens on the first `step()`.
789            CompiledStageKind::Modal { num_modes } => match modal_analysis(&mut domain, *num_modes)
790            {
791                Ok(modes) => {
792                    self.runner = Some(StageRunner::Modal {
793                        domain,
794                        modes,
795                        steps_remaining: stage.steps,
796                    });
797                }
798                Err(error) => self.error = Some(error.into()),
799            },
800            CompiledStageKind::Reset => unreachable!("handled above"),
801            CompiledStageKind::Transient {
802                algorithm,
803                convergence,
804                damping,
805                modal_anchors,
806                dt,
807                ground_motions,
808            } => {
809                let resolved = (|| -> Result<RayleighDamping, AnalysisError> {
810                    let Some(anchors) = modal_anchors else {
811                        return Ok(*damping);
812                    };
813                    let count = anchors.mode1.max(anchors.mode2) as usize;
814                    let modes = modal_analysis(&mut domain, count)?;
815                    let w1 = modes.get(anchors.mode1 as usize - 1).map(|m| m.frequency);
816                    let w2 = modes.get(anchors.mode2 as usize - 1).map(|m| m.frequency);
817                    let (Some(w1), Some(w2)) = (w1, w2) else {
818                        return Err(AnalysisError::InvalidModeCount {
819                            requested: count,
820                            free_dofs: modes.len(),
821                        });
822                    };
823                    if !w1.is_finite() || !w2.is_finite() || w1 <= 0.0 || w2 <= 0.0 {
824                        return Err(AnalysisError::InvalidOption {
825                            field: "damping.modalAnchors",
826                        });
827                    }
828                    Ok(RayleighDamping::new(
829                        2.0 * anchors.ratio * w1 * w2 / (w1 + w2),
830                        2.0 * anchors.ratio / (w1 + w2),
831                    ))
832                })();
833                match resolved.and_then(|damping| TransientAnalysis::new(domain, damping, *dt)) {
834                    Ok(mut analysis) => {
835                        analysis = analysis.with_algorithm(*algorithm, *convergence);
836                        for motion in ground_motions.iter().cloned() {
837                            analysis = analysis.with_ground_motion(motion);
838                        }
839                        self.runner = Some(StageRunner::Transient {
840                            analysis,
841                            steps_remaining: stage.steps,
842                        });
843                    }
844                    Err(error) => self.error = Some(error.into()),
845                }
846            }
847        }
848        self.load_factor = 0.0;
849        if self.record_initial
850            && matches!(
851                self.stages[self.current_stage].kind,
852                CompiledStageKind::Static { .. } | CompiledStageKind::Transient { .. }
853            )
854        {
855            self.record_sample();
856        }
857    }
858
859    fn steps_remaining(&self) -> u32 {
860        match self.runner.as_ref().expect("stage started above") {
861            StageRunner::Static {
862                steps_remaining, ..
863            }
864            | StageRunner::Modal {
865                steps_remaining, ..
866            }
867            | StageRunner::Reset {
868                steps_remaining, ..
869            }
870            | StageRunner::Transient {
871                steps_remaining, ..
872            } => *steps_remaining,
873        }
874    }
875
876    /// Advances whichever kind of stage is currently running by one unit —
877    /// a Newton step (`Static`), a Newmark step (`Transient`), or (for
878    /// `Modal`, whose eigensolve already ran in `start_current_stage`) just
879    /// consuming its one always-available "step" — then records a sample
880    /// for every recorder this runner kind supports.
881    fn step_once(&mut self) -> Result<(), AnalysisError> {
882        match self.runner.as_mut().expect("stage started above") {
883            StageRunner::Static {
884                analysis,
885                steps_remaining,
886            } => {
887                let result = analysis.step()?;
888                self.load_factor = result.load_factor;
889                *steps_remaining -= 1;
890                if let Some(info) = result.arc {
891                    self.continuation = Some(info.into());
892                    // A met stop criterion ends the stage early.
893                    if info.stop.is_some() {
894                        *steps_remaining = 0;
895                    }
896                }
897            }
898            StageRunner::Transient {
899                analysis,
900                steps_remaining,
901            } => {
902                let result = analysis.step()?;
903                self.load_factor = result.time;
904                *steps_remaining -= 1;
905            }
906            StageRunner::Modal {
907                steps_remaining, ..
908            } => {
909                *steps_remaining -= 1;
910            }
911            StageRunner::Reset { .. } => {}
912        }
913        self.record_sample();
914        Ok(())
915    }
916
917    /// `Analysis::into_domain`/fresh-`AnalysisBuilder` multi-phase pattern
918    /// (as in the two-phase tests in
919    /// `core/tests/elements/force_beam.rs`): freeze this stage's held patterns (a `Static`-only concept —
920    /// `Modal`/`Transient` stages hold none) at its final load factor, then
921    /// hand the same `Domain` (committed material state and all) to the
922    /// next stage's fresh `Analysis`/`TransientAnalysis`/`modal_analysis`
923    /// call.
924    fn finish_current_stage(&mut self) {
925        let hold_patterns_after = match &self.stages[self.current_stage].kind {
926            CompiledStageKind::Static {
927                hold_patterns_after,
928                ..
929            } => hold_patterns_after.clone(),
930            CompiledStageKind::Modal { .. }
931            | CompiledStageKind::Transient { .. }
932            | CompiledStageKind::Reset => Vec::new(),
933        };
934        let load_factor = self.load_factor;
935
936        let mut domain = match self.runner.take().expect("stage running") {
937            StageRunner::Static { analysis, .. } => analysis.into_domain(),
938            StageRunner::Modal {
939                mut domain, modes, ..
940            } => {
941                let result = self.modal_stage_result(&mut domain, &modes);
942                self.modal_results.push(result);
943                domain
944            }
945            StageRunner::Transient { analysis, .. } => analysis.into_domain(),
946            StageRunner::Reset { domain, .. } => domain,
947        };
948        for pattern in hold_patterns_after {
949            domain.hold_pattern_constant(pattern, load_factor);
950            self.holds.push((pattern, load_factor));
951        }
952        self.domain = Some(domain);
953        self.current_stage += 1;
954    }
955
956    fn modal_stage_result(
957        &self,
958        domain: &mut Domain<NDIM, NDOF, NId, E>,
959        modes: &[Mode],
960    ) -> ModalStageResult {
961        let mass = domain.assemble_mass_diagonal();
962        // Influence vector per global translation direction: a rigid translation of the whole
963        // model in reduced coordinates (1 on every free unknown of that direction), which every
964        // rigid-body constraint (diaphragm, rigid link) satisfies exactly.
965        let influence: Vec<Vec<f64>> = (0..NDIM)
966            .map(|direction| {
967                domain
968                    .direction_incidence(direction)
969                    .iter()
970                    .copied()
971                    .collect()
972            })
973            .collect();
974        let total_mass: Vec<f64> = influence
975            .iter()
976            .map(|r| r.iter().zip(mass.iter()).map(|(r, m)| r * r * m).sum())
977            .collect();
978        let modes = modes
979            .iter()
980            .map(|mode| {
981                let sign = mode
982                    .shape
983                    .iter()
984                    .copied()
985                    .fold(
986                        0.0_f64,
987                        |best, v| if v.abs() > best.abs() { v } else { best },
988                    )
989                    .signum();
990                let sign = if sign == 0.0 { 1.0 } else { sign };
991                let mut shape = Vec::with_capacity(self.node_ids.len() * NDOF);
992                for &node in &self.node_ids {
993                    for dof in 0..NDOF {
994                        shape.push(sign * domain.value_at(&mode.shape, node, dof));
995                    }
996                }
997                let participation: Vec<f64> = influence
998                    .iter()
999                    .map(|r| {
1000                        sign * (0..mass.len())
1001                            .map(|i| mode.shape[i] * mass[i] * r[i])
1002                            .sum::<f64>()
1003                    })
1004                    .collect();
1005                let mass_ratio = participation
1006                    .iter()
1007                    .zip(&total_mass)
1008                    .map(|(g, total)| if *total > 0.0 { g * g / total } else { 0.0 })
1009                    .collect();
1010                ModeResult {
1011                    frequency: mode.frequency,
1012                    shape,
1013                    participation,
1014                    mass_ratio,
1015                }
1016            })
1017            .collect();
1018        ModalStageResult {
1019            stage_index: self.current_stage,
1020            stage_id: self.stages[self.current_stage].id.clone(),
1021            ndf: NDOF,
1022            modes,
1023            total_mass,
1024        }
1025    }
1026
1027    fn record_sample(&mut self) {
1028        let Some(runner) = &self.runner else {
1029            return;
1030        };
1031        match runner {
1032            StageRunner::Static { analysis, .. } => {
1033                let progress = self.load_factor;
1034                for (recorder, batch) in self.recorders.iter().zip(self.current_batch.iter_mut()) {
1035                    let value = match *recorder {
1036                        ResolvedRecorder::NodeDisp { node, dof } => {
1037                            Some(analysis.domain().node(node).displacement[dof as usize])
1038                        }
1039                        ResolvedRecorder::ElementForce { element, component } => Some(
1040                            analysis.domain().element_end_force(element, progress)
1041                                [component as usize],
1042                        ),
1043                        ResolvedRecorder::ElementLoad { element, component } => {
1044                            Some(analysis.domain().element_load_component(
1045                                element,
1046                                progress,
1047                                component as usize,
1048                            ))
1049                        }
1050                        ResolvedRecorder::Reaction { node, dof } => {
1051                            Some(analysis.domain().reaction(node, dof as usize, progress))
1052                        }
1053                        ResolvedRecorder::Fiber {
1054                            element,
1055                            point,
1056                            fiber,
1057                            response,
1058                        } => fiber_value(analysis.domain(), element, point, fiber, response),
1059                        ResolvedRecorder::GaussPoint {
1060                            element,
1061                            point,
1062                            quantity,
1063                            component,
1064                        } => gauss_value(analysis.domain(), element, point, quantity, component),
1065                        ResolvedRecorder::NodeVel { .. }
1066                        | ResolvedRecorder::NodeAccel { .. }
1067                        | ResolvedRecorder::ModeShape { .. } => None,
1068                    };
1069                    if let Some(value) = value {
1070                        batch.push((progress, value));
1071                    }
1072                }
1073            }
1074            StageRunner::Transient { analysis, .. } => {
1075                let progress = self.load_factor; // == analysis.time(), set by step_once
1076                for (recorder, batch) in self.recorders.iter().zip(self.current_batch.iter_mut()) {
1077                    let value = match *recorder {
1078                        ResolvedRecorder::NodeDisp { node, dof } => {
1079                            Some(analysis.domain().node(node).displacement[dof as usize])
1080                        }
1081                        ResolvedRecorder::NodeVel { node, dof } => {
1082                            Some(analysis.domain().node(node).velocity[dof as usize])
1083                        }
1084                        ResolvedRecorder::NodeAccel { node, dof } => {
1085                            Some(analysis.domain().node(node).acceleration[dof as usize])
1086                        }
1087                        ResolvedRecorder::ElementForce { element, component } => Some(
1088                            analysis.domain().element_end_force(element, progress)
1089                                [component as usize],
1090                        ),
1091                        ResolvedRecorder::ElementLoad { element, component } => {
1092                            Some(analysis.domain().element_load_component(
1093                                element,
1094                                progress,
1095                                component as usize,
1096                            ))
1097                        }
1098                        ResolvedRecorder::Reaction { node, dof } => {
1099                            Some(analysis.reaction(node, dof as usize))
1100                        }
1101                        ResolvedRecorder::Fiber {
1102                            element,
1103                            point,
1104                            fiber,
1105                            response,
1106                        } => fiber_value(analysis.domain(), element, point, fiber, response),
1107                        ResolvedRecorder::GaussPoint {
1108                            element,
1109                            point,
1110                            quantity,
1111                            component,
1112                        } => gauss_value(analysis.domain(), element, point, quantity, component),
1113                        ResolvedRecorder::ModeShape { .. } => None,
1114                    };
1115                    if let Some(value) = value {
1116                        batch.push((progress, value));
1117                    }
1118                }
1119            }
1120            StageRunner::Reset { .. } => {}
1121            StageRunner::Modal { domain, modes, .. } => {
1122                for (recorder, batch) in self.recorders.iter().zip(self.current_batch.iter_mut()) {
1123                    let ResolvedRecorder::ModeShape { mode, node, dof } = *recorder else {
1124                        continue;
1125                    };
1126                    let Some(computed) = modes.get(mode as usize) else {
1127                        continue;
1128                    };
1129                    // `value_at` reads a constrained DOF (a diaphragm or
1130                    // rigid-link slave) through its constraint and is zero
1131                    // for a fixed or inactive DOF, which can't participate
1132                    // in any mode.
1133                    let value = domain.value_at(&computed.shape, node, dof as usize);
1134                    batch.push((computed.frequency, value));
1135                }
1136            }
1137        }
1138    }
1139}
carapace_wasm/input_v1/session.rs

carapace_wasm/input_v1/
session.rs
