carapace_core/analysis/
integrator.rs

1use nalgebra::DVector;
2use slotmap::Key;
3
4use crate::model::{Domain, ElementOps, NodeId, SparseMatrix};
5
6use super::{AnalysisError, ArcLength, SparseFactorization, SparseSolver};
7
8/// Static/pseudo-static integration strategy. Closed enum.
9/// `LoadControl`/`DisplacementControl` only change how this *step's* load
10/// factor is chosen (the "predictor") — the iteration that follows
11/// (see `Algorithm`) proceeds identically regardless of which one picked it.
12/// `ArcLength` is different in kind: it predicts displacement *and* load
13/// factor and corrects both together through a bordered system, so
14/// `Analysis::step` hands it to its own continuation driver instead of
15/// `predict`/`correct` (see `arclength.rs` and `docs/arclength.md`).
16///
17/// Generic over `NId` (the controlled node's ID type — `NodeId` or
18/// `Node3Id`), defaulted to the planar profile so bare `Integrator` keeps
19/// working, the same trick `Domain` uses. `predict`/`correct` take the rest
20/// of the profile (`NDIM`/`NDOF`/`ELEMENT_DOF`/`E`) as their own
21/// generic parameters instead of `Integrator`'s, since an integrator's
22/// identity only depends on which node/DOF it controls, not which `Domain`
23/// it's later used with.
24#[derive(Debug, Clone, PartialEq)]
25pub enum Integrator<NId = NodeId> {
26    LoadControl {
27        increment: f64,
28    },
29    /// Chooses this step's pseudo-time increment so that the controlled
30    /// node/DOF's displacement changes by exactly `increment`, via the
31    /// standard unit-load technique: solve the current tangent against the
32    /// *sensitivity* of the total applied load to a unit pseudo-time
33    /// increment (`Domain::assemble_reference_load_sensitivity` — every
34    /// non-frozen `LoadPattern`'s contribution, weighted by its series'
35    /// slope at the current pseudo-time; a `hold_pattern_constant`-frozen
36    /// or `LoadSeries::Constant` pattern contributes nothing, since it
37    /// doesn't respond to pseudo-time at all) once, then scale by whatever
38    /// factor makes the controlled DOF's response match `increment`.
39    DisplacementControl {
40        node: NId,
41        dof: usize,
42        increment: f64,
43    },
44    /// Arc-length continuation along the equilibrium path, through load
45    /// limits and snap-back — see `ArcLength`.
46    ArcLength(ArcLength<NId>),
47}
48
49impl<NId: Copy> Integrator<NId> {
50    /// Compute this step's new pseudo-time, given the domain's state as of
51    /// the end of the *previous* (converged) step, plus the displacement
52    /// predictor the caller must apply before iterating (`None` for
53    /// `LoadControl`). `DisplacementControl`'s predictor is the unit-load
54    /// response scaled by the pseudo-time increment, so the controlled DOF
55    /// starts the step exactly at its target and every Newton iteration
56    /// (the first included) measures a true correction. That is what lets
57    /// a linear step converge in one iteration instead of spending its
58    /// first on the predictor and a second on proving `du` is small.
59    pub(crate) fn predict<const NDIM: usize, const NDOF: usize, E>(
60        &self,
61        domain: &Domain<NDIM, NDOF, NId, E>,
62        solver: &SparseSolver,
63        current_pseudo_time: f64,
64    ) -> Result<(f64, Option<DVector<f64>>), AnalysisError>
65    where
66        NId: Key,
67        E: ElementOps<NDIM, NDOF, NId>,
68    {
69        match self {
70            Integrator::LoadControl { increment } => Ok((current_pseudo_time + increment, None)),
71            Integrator::DisplacementControl {
72                node,
73                dof,
74                increment,
75            } => {
76                let eq = domain
77                    .equation_of(*node, *dof)
78                    .ok_or(AnalysisError::InvalidConstraint)?;
79                let (k, _resistance) = domain.assemble_tangent_and_resistance(current_pseudo_time);
80                let sensitivity = domain.assemble_reference_load_sensitivity(current_pseudo_time);
81                let unit_response = solver.solve(&k, &sensitivity)?;
82                let delta_lambda = increment / unit_response[eq];
83                Ok((
84                    current_pseudo_time + delta_lambda,
85                    Some(unit_response * delta_lambda),
86                ))
87            }
88            Integrator::ArcLength(_) => {
89                unreachable!("Analysis::step routes ArcLength to its continuation driver")
90            }
91        }
92    }
93
94    /// Corrector for every Newton iteration of a step, applied after
95    /// `predict`'s displacement predictor.
96    /// `predict`'s pseudo-time only accounts for the tangent at the *start*
97    /// of the step — exact for `LoadControl` (whose pseudo-time never
98    /// changes mid-step, so `du_bar` unmodified is already the right
99    /// increment every iteration).
100    ///
101    /// For `DisplacementControl`, `predict` already moved the controlled DOF
102    /// to its target (OpenSees' `newStep` does the same), but the tangent it
103    /// used can go stale as the iteration proceeds (e.g. a step landing
104    /// exactly on a material's yield breakpoint sees the elastic tangent at
105    /// `predict` time, then the true, much softer, post-yield tangent) —
106    /// left uncorrected, ordinary fixed-load-factor Newton iteration would
107    /// then keep driving the controlled DOF *past* the target already
108    /// reached, chasing force equilibrium at the wrong load factor instead.
109    /// This is the standard Yang & Shieh consistent Displacement Control
110    /// corrector: solve the current tangent against the reference load
111    /// sensitivity again (`unit_response`, this iteration's version of
112    /// `predict`'s unit-load probe), then pick `delta_lambda` so that adding
113    /// `delta_lambda * unit_response` to `du_bar` cancels `du_bar`'s
114    /// contribution at the controlled DOF — i.e. this iteration leaves the
115    /// controlled DOF exactly where it already is, no matter how far the
116    /// tangent has drifted, while every other DOF still gets `du_bar`'s
117    /// residual-driven correction.
118    /// `current_factorization`, when supplied, must factor `k`. An older
119    /// Newton tangent must not be passed here: the load-sensitivity probe
120    /// continues to use the current tangent even under tangent reuse.
121    pub(crate) fn correct<const NDIM: usize, const NDOF: usize, E>(
122        &self,
123        domain: &Domain<NDIM, NDOF, NId, E>,
124        solver: &SparseSolver,
125        k: &SparseMatrix,
126        current_factorization: Option<&SparseFactorization>,
127        du_bar: DVector<f64>,
128        pseudo_time: f64,
129    ) -> Result<(f64, DVector<f64>), AnalysisError>
130    where
131        NId: Key,
132        E: ElementOps<NDIM, NDOF, NId>,
133    {
134        match self {
135            Integrator::LoadControl { .. } => Ok((0.0, du_bar)),
136            Integrator::DisplacementControl { node, dof, .. } => {
137                let eq = domain
138                    .equation_of(*node, *dof)
139                    .ok_or(AnalysisError::InvalidConstraint)?;
140                let sensitivity = domain.assemble_reference_load_sensitivity(pseudo_time);
141                let unit_response = match current_factorization {
142                    Some(factorization) => factorization.solve(&sensitivity),
143                    None => solver.solve(k, &sensitivity)?,
144                };
145                let delta_lambda = -du_bar[eq] / unit_response[eq];
146                let du = du_bar + unit_response * delta_lambda;
147                Ok((delta_lambda, du))
148            }
149            Integrator::ArcLength(_) => {
150                unreachable!("Analysis::step routes ArcLength to its continuation driver")
151            }
152        }
153    }
154}
carapace_core/analysis/integrator.rs

carapace_core/analysis/
integrator.rs
