carapace_core/model/
load_pattern.rs

1use std::collections::HashMap;
2use std::hash::Hash;
3
4use slotmap::new_key_type;
5
6use super::{ElementId, NodeId, NDF};
7
8new_key_type! {
9    /// Generational index into `Domain`'s load-pattern store.
10    pub struct LoadPatternId;
11}
12
13/// The pseudo-time-to-load-factor function a `LoadPattern` scales its
14/// reference loads by — Xara/OpenSees's `TimeSeries` family, cut down to
15/// exactly the three shapes Carapace's scope actually needs (`Trig`/
16/// `Rectangular`/`Pulse` etc. are for shaping a *dynamic* excitation
17/// pattern, not a quasi-static load/displacement protocol — out of scope
18/// until something needs them). Closed enum.
19#[derive(Debug, Clone)]
20pub enum LoadSeries {
21    /// `factor(t) = 1.0` for every `t` — the pattern is always fully
22    /// applied, independent of pseudo-time.
23    Constant,
24    /// `factor(t) = slope * t`. `slope: 1.0` reproduces the single-pattern,
25    /// `load_factor`-scales-everything behavior every `Analysis` had before
26    /// multi-pattern support existed.
27    Linear { slope: f64 },
28    /// Piecewise-linear interpolation between `(times[i], factors[i])`
29    /// control points — Xara's `PathSeries`/`PathTimeSeries`. This is the
30    /// load-path / quasi-static cyclic-protocol vehicle: drive it with
31    /// `Integrator::LoadControl` (stepping pseudo-time along `times`) for a
32    /// load-controlled protocol, or read an accelerogram's acceleration
33    /// values into `factors` for `GroundMotion`. Holds the last point's
34    /// value beyond the last control point (Xara's `useLast` default),
35    /// and the first point's value before the first. `times` must be
36    /// sorted ascending and the same length as `factors` — not enforced by
37    /// the type (a malformed path is a model-construction error, not a
38    /// runtime failure, same class as `AnalysisBuilder::build`'s
39    /// `ConstraintHandler` assertion).
40    Path { times: Vec<f64>, factors: Vec<f64> },
41    /// OpenSees accelerogram bounds: zero before the first sample and after the last,
42    /// unless use_last holds the last value.
43    BoundedPath {
44        times: Vec<f64>,
45        factors: Vec<f64>,
46        use_last: bool,
47    },
48}
49
50impl LoadSeries {
51    pub(crate) fn factor(&self, pseudo_time: f64) -> f64 {
52        match self {
53            LoadSeries::Constant => 1.0,
54            LoadSeries::Linear { slope } => slope * pseudo_time,
55            LoadSeries::Path { times, factors } => path_interpolate(times, factors, pseudo_time),
56            LoadSeries::BoundedPath {
57                times,
58                factors,
59                use_last,
60            } => {
61                if pseudo_time < times[0] || (!use_last && pseudo_time > times[times.len() - 1]) {
62                    0.0
63                } else {
64                    path_interpolate(times, factors, pseudo_time)
65                }
66            }
67        }
68    }
69
70    /// `d(factor)/d(pseudo_time)` at `pseudo_time` — what `Integrator::
71    /// DisplacementControl`'s unit-load technique needs (it must know how
72    /// the load *responds* to a unit pseudo-time increment, not its
73    /// current value). Zero for `Constant` (doesn't respond to pseudo-time
74    /// at all — same reason a `hold_pattern_constant`-frozen pattern is
75    /// excluded from the sensitivity sum entirely, see `Domain::
76    /// assemble_reference_load_sensitivity`), `slope` for `Linear`, the
77    /// active segment's slope for `Path`.
78    pub(crate) fn slope(&self, pseudo_time: f64) -> f64 {
79        match self {
80            LoadSeries::Constant => 0.0,
81            LoadSeries::Linear { slope } => *slope,
82            LoadSeries::Path { times, factors }
83            | LoadSeries::BoundedPath { times, factors, .. } => {
84                path_slope(times, factors, pseudo_time)
85            }
86        }
87    }
88}
89
90fn path_segment(times: &[f64], pseudo_time: f64) -> usize {
91    // Index `i` such that `pseudo_time` falls in `[times[i], times[i+1])`,
92    // clamped to the path's first/last segment outside its range (Xara's
93    // `useLast`/pre-first-point behavior). `times` has at least 2 entries
94    // here (the `times.len() == 1` case is handled by callers before this
95    // is reached).
96    let n = times.partition_point(|&t| t <= pseudo_time);
97    n.saturating_sub(1).min(times.len() - 2)
98}
99
100fn path_interpolate(times: &[f64], factors: &[f64], pseudo_time: f64) -> f64 {
101    if times.len() == 1 {
102        return factors[0];
103    }
104    let i = path_segment(times, pseudo_time);
105    let (t0, t1) = (times[i], times[i + 1]);
106    let (f0, f1) = (factors[i], factors[i + 1]);
107    if pseudo_time <= t0 {
108        f0
109    } else if pseudo_time >= t1 {
110        f1
111    } else {
112        f0 + (f1 - f0) * (pseudo_time - t0) / (t1 - t0)
113    }
114}
115
116fn path_slope(times: &[f64], factors: &[f64], pseudo_time: f64) -> f64 {
117    if times.len() == 1 {
118        return 0.0;
119    }
120    let i = path_segment(times, pseudo_time);
121    let (t0, t1) = (times[i], times[i + 1]);
122    if pseudo_time < t0 || pseudo_time > t1 {
123        return 0.0; // held constant outside the path's range — no response to pseudo-time there
124    }
125    (factors[i + 1] - factors[i]) / (t1 - t0)
126}
127
128/// The load on a 2D element: an additive accumulator with one field group per
129/// load kind, so that the sum of every pattern's load on an element (each
130/// scaled by its factor) is just a field-wise sum. `Uniform + Body` is simply
131/// both fields set. Which kinds an element accepts is decided by
132/// `ElementOps::accepts_load`, checked by `Domain::validate` before any load
133/// is accumulated.
134///
135/// Component frames differ by kind: `beam_uniform` is in the element's
136/// *local* axes (local `x` is the member axis), as in OpenSees's
137/// `-beamUniform`; `body`, `edge_traction` and `edge_pressure` are in *global*
138/// axes (pressure along the inward normal).
139#[derive(Debug, Clone, Copy, PartialEq, Default)]
140pub struct ElementLoad {
141    /// Beam-columns: uniform load (force/length) `[wx, wy]`, `wx` along the
142    /// member axis and `wy` transverse in local +y.
143    pub beam_uniform: [f64; 2],
144    /// Continuum elements: body force per unit volume `[bx, by]`.
145    pub body: [f64; 2],
146    /// Continuum elements: traction `[tx, ty]` per unit edge length on edge
147    /// `k`, which joins local node `k` to `k + 1`.
148    pub edge_traction: [[f64; 2]; 4],
149    /// Continuum elements: pressure on edge `k`, positive acting into the
150    /// element along the inward normal.
151    pub edge_pressure: [f64; 4],
152}
153
154impl ElementLoad {
155    /// A uniform beam load (the former `ElementLoad::Uniform { wx, wy }`).
156    pub fn uniform(wx: f64, wy: f64) -> Self {
157        ElementLoad {
158            beam_uniform: [wx, wy],
159            ..Default::default()
160        }
161    }
162
163    /// A body force per unit volume in global axes.
164    pub fn body(bx: f64, by: f64) -> Self {
165        ElementLoad {
166            body: [bx, by],
167            ..Default::default()
168        }
169    }
170
171    /// A global-axes traction per unit length on edge `edge` (`0..4`).
172    pub fn edge_traction(edge: usize, tx: f64, ty: f64) -> Self {
173        let mut load = ElementLoad::default();
174        load.edge_traction[edge] = [tx, ty];
175        load
176    }
177
178    /// A pressure on edge `edge` (`0..4`), positive into the element.
179    pub fn edge_pressure(edge: usize, pressure: f64) -> Self {
180        let mut load = ElementLoad::default();
181        load.edge_pressure[edge] = pressure;
182        load
183    }
184
185    /// Whether any beam-load field is nonzero.
186    pub fn has_beam(&self) -> bool {
187        self.beam_uniform != [0.0; 2]
188    }
189
190    /// Whether the body force is nonzero.
191    pub fn has_body(&self) -> bool {
192        self.body != [0.0; 2]
193    }
194
195    /// Whether any edge `>= edges` carries a traction or pressure.
196    pub fn has_edge_from(&self, edges: usize) -> bool {
197        (edges..4).any(|k| self.edge_traction[k] != [0.0; 2] || self.edge_pressure[k] != 0.0)
198    }
199
200    /// Whether any edge carries a traction or pressure.
201    pub fn has_edge(&self) -> bool {
202        self.has_edge_from(0)
203    }
204}
205
206impl std::ops::Add for ElementLoad {
207    type Output = ElementLoad;
208
209    /// Field-wise sum — several loads on one element accumulate, as in OpenSees.
210    fn add(self, other: ElementLoad) -> ElementLoad {
211        ElementLoad {
212            beam_uniform: std::array::from_fn(|i| self.beam_uniform[i] + other.beam_uniform[i]),
213            body: std::array::from_fn(|i| self.body[i] + other.body[i]),
214            edge_traction: std::array::from_fn(|k| {
215                std::array::from_fn(|i| self.edge_traction[k][i] + other.edge_traction[k][i])
216            }),
217            edge_pressure: std::array::from_fn(|k| self.edge_pressure[k] + other.edge_pressure[k]),
218        }
219    }
220}
221
222/// Scalar view of an element load's components (local axes), for recording
223/// the load an element actually carries at some pseudo-time.
224pub trait ElementLoadComponents {
225    /// Number of components (`ElementLoad`: 2, `ElementLoad3`: 3); the bound a
226    /// recorder's `component` index is checked against.
227    const COUNT: usize;
228
229    /// Component `index` (`ElementLoad`: `0` = `wx`, `1` = `wy`;
230    /// `ElementLoad3`: `0..3` = `wx`, `wy`, `wz`); `0.0` past the last.
231    fn component(&self, index: usize) -> f64;
232}
233
234/// Component layout: `[wx, wy, bx, by, t0x, t0y, p0, t1x, t1y, p1, t2x, t2y, p2, t3x, t3y, p3]`
235/// (beam uniform, body force, then traction and pressure per edge).
236impl ElementLoadComponents for ElementLoad {
237    const COUNT: usize = 16;
238
239    fn component(&self, index: usize) -> f64 {
240        match index {
241            0 | 1 => self.beam_uniform[index],
242            2 | 3 => self.body[index - 2],
243            4..=15 => {
244                let (edge, slot) = ((index - 4) / 3, (index - 4) % 3);
245                match slot {
246                    0 | 1 => self.edge_traction[edge][slot],
247                    _ => self.edge_pressure[edge],
248                }
249            }
250            _ => 0.0,
251        }
252    }
253}
254
255/// Component layout: `[wx, wy, wz, bx, by, bz, pressure]` (beam uniform, shell body, shell pressure).
256impl ElementLoadComponents for ElementLoad3 {
257    const COUNT: usize = 7;
258
259    fn component(&self, index: usize) -> f64 {
260        match index {
261            0..=2 => self.uniform[index],
262            3..=5 => self.body[index - 3],
263            6 => self.pressure,
264            _ => 0.0,
265        }
266    }
267}
268
269impl std::ops::Mul<f64> for ElementLoad {
270    type Output = ElementLoad;
271
272    fn mul(self, factor: f64) -> ElementLoad {
273        ElementLoad {
274            beam_uniform: self.beam_uniform.map(|v| v * factor),
275            body: self.body.map(|v| v * factor),
276            edge_traction: self.edge_traction.map(|t| t.map(|v| v * factor)),
277            edge_pressure: self.edge_pressure.map(|v| v * factor),
278        }
279    }
280}
281
282/// `ElementLoad`'s 3D counterpart: an additive accumulator, one field group per load kind.
283/// `ElementOps::Load` for the 3D profile.
284#[derive(Debug, Clone, Copy, PartialEq, Default)]
285pub struct ElementLoad3 {
286    /// Beam-columns: uniform load (force/length) in the member's local axes: `wx` along the
287    /// member axis, `wy`/`wz` transverse in local +y/+z (`vec_xz` defines `y`/`z`).
288    /// Xara/OpenSees's `Beam3dUniformLoad`.
289    pub uniform: [f64; 3],
290    /// Shells: body acceleration in global axes (a gravity vector `g` points down); the force
291    /// per unit area is `+rho h b`. OpenSees' `ShellMITC4` applies `-rho h b` for the same
292    /// `-selfWeight` factors, so exporters negate.
293    pub body: [f64; 3],
294    /// Shells: pressure per unit area, positive along the local normal `e3` (right-hand rule
295    /// from the node order).
296    pub pressure: f64,
297}
298
299impl ElementLoad3 {
300    /// A uniform beam load.
301    pub const fn uniform(wx: f64, wy: f64, wz: f64) -> Self {
302        ElementLoad3 {
303            uniform: [wx, wy, wz],
304            body: [0.0; 3],
305            pressure: 0.0,
306        }
307    }
308
309    /// A shell self-weight (acceleration factors in global axes).
310    pub fn body(bx: f64, by: f64, bz: f64) -> Self {
311        ElementLoad3 {
312            body: [bx, by, bz],
313            ..Default::default()
314        }
315    }
316
317    /// A shell pressure along the local normal.
318    pub fn pressure(pressure: f64) -> Self {
319        ElementLoad3 {
320            pressure,
321            ..Default::default()
322        }
323    }
324
325    pub fn has_uniform(&self) -> bool {
326        self.uniform != [0.0; 3]
327    }
328
329    pub fn has_body(&self) -> bool {
330        self.body != [0.0; 3]
331    }
332
333    pub fn has_pressure(&self) -> bool {
334        self.pressure != 0.0
335    }
336}
337
338impl std::ops::Add for ElementLoad3 {
339    type Output = ElementLoad3;
340
341    /// See `ElementLoad`'s `Add`.
342    fn add(self, other: ElementLoad3) -> ElementLoad3 {
343        ElementLoad3 {
344            uniform: std::array::from_fn(|i| self.uniform[i] + other.uniform[i]),
345            body: std::array::from_fn(|i| self.body[i] + other.body[i]),
346            pressure: self.pressure + other.pressure,
347        }
348    }
349}
350
351impl std::ops::Mul<f64> for ElementLoad3 {
352    type Output = ElementLoad3;
353
354    fn mul(self, factor: f64) -> ElementLoad3 {
355        ElementLoad3 {
356            uniform: self.uniform.map(|v| v * factor),
357            body: self.body.map(|v| v * factor),
358            pressure: self.pressure * factor,
359        }
360    }
361}
362
363/// A named collection of reference loads (nodal and element), scaled by a
364/// `LoadSeries` function of pseudo-time — Xara/OpenSees's `LoadPattern`.
365/// `Domain` owns any number of these; multiple patterns coexist and sum
366/// (`Domain::assemble_reference_load`), each independently scaled, so e.g.
367/// gravity can be ramped to full scale and then held constant
368/// (`Domain::hold_pattern_constant`) while a separate lateral pattern
369/// continues ramping in a later analysis phase.
370///
371/// Generic over `NDOF` (nodal-load array size), the node/element ID types,
372/// and `EL` (the element-load kind — `ElementLoad` for the planar profile,
373/// `ElementLoad3` for the spatial one, matching `ElementOps::Load`). `NDOF` defaults to the planar profile so bare `LoadPattern`
374/// keeps working unchanged, the same trick `Node`'s defaults use.
375#[derive(Debug, Clone)]
376pub(crate) struct LoadPattern<
377    const NDOF: usize = NDF,
378    NId = NodeId,
379    EId = ElementId,
380    EL = ElementLoad,
381> where
382    NId: Eq + Hash + Copy,
383    EId: Eq + Hash + Copy,
384{
385    series: LoadSeries,
386    scale_factor: f64,
387    /// `Some(frozen_factor)` once `Domain::hold_pattern_constant` has been
388    /// called — computed once, at freeze time, from the exact pseudo-time
389    /// the caller was at (not lazily cached on every assemble call, which
390    /// would require `assemble_reference_load` to take `&mut self`). `None` means
391    /// "still driven by `series`".
392    frozen_factor: Option<f64>,
393    nodal_loads: HashMap<NId, [f64; NDOF]>,
394    element_loads: HashMap<EId, EL>,
395}
396
397impl<const NDOF: usize, NId, EId, EL> LoadPattern<NDOF, NId, EId, EL>
398where
399    NId: Eq + Hash + Copy,
400    EId: Eq + Hash + Copy,
401    EL: Copy + std::ops::Add<Output = EL>,
402{
403    pub(crate) fn has_element_loads(&self) -> bool {
404        !self.element_loads.is_empty()
405    }
406
407    pub(crate) fn new(series: LoadSeries) -> Self {
408        LoadPattern {
409            series,
410            scale_factor: 1.0,
411            frozen_factor: None,
412            nodal_loads: HashMap::new(),
413            element_loads: HashMap::new(),
414        }
415    }
416
417    pub(crate) fn with_scale_factor(mut self, scale_factor: f64) -> Self {
418        self.scale_factor = scale_factor;
419        self
420    }
421
422    pub(crate) fn factor(&self, pseudo_time: f64) -> f64 {
423        self.frozen_factor
424            .unwrap_or_else(|| self.series.factor(pseudo_time) * self.scale_factor)
425    }
426
427    /// Zero once frozen (see `LoadSeries::slope`) — a frozen
428    /// pattern doesn't respond to further pseudo-time change.
429    pub(crate) fn sensitivity(&self, pseudo_time: f64) -> f64 {
430        if self.frozen_factor.is_some() {
431            0.0
432        } else {
433            self.series.slope(pseudo_time) * self.scale_factor
434        }
435    }
436
437    /// Whether this pattern still follows a `LoadSeries::Path` — arc-length
438    /// continuation rejects these, since a path's slope isn't valid across
439    /// a load reversal or a breakpoint.
440    pub(crate) fn is_unfrozen_path(&self) -> bool {
441        self.frozen_factor.is_none()
442            && matches!(
443                self.series,
444                LoadSeries::Path { .. } | LoadSeries::BoundedPath { .. }
445            )
446    }
447
448    pub(crate) fn hold_constant(&mut self, pseudo_time: f64) {
449        self.frozen_factor = Some(self.factor(pseudo_time));
450    }
451
452    pub(crate) fn add_nodal_load(&mut self, node: NId, dof: usize, value: f64) {
453        self.nodal_loads.entry(node).or_insert([0.0; NDOF])[dof] = value;
454    }
455
456    /// Adds `load` to whatever this pattern already applies to `element`
457    /// (loads accumulate rather than replace).
458    pub(crate) fn add_element_load(&mut self, element: EId, load: EL) {
459        self.element_loads
460            .entry(element)
461            .and_modify(|existing| *existing = *existing + load)
462            .or_insert(load);
463    }
464
465    pub(crate) fn nodal_load(&self, node: NId) -> Option<&[f64; NDOF]> {
466        self.nodal_loads.get(&node)
467    }
468
469    pub(crate) fn element_loads(&self) -> impl Iterator<Item = (&EId, &EL)> {
470        self.element_loads.iter()
471    }
472
473    pub(crate) fn element_load(&self, element: EId) -> Option<&EL> {
474        self.element_loads.get(&element)
475    }
476}
477
478/// Patterns that currently contribute element loads, each with its factor
479/// at `pseudo_time` — computed once per assembly so the per-element lookup
480/// below doesn't re-evaluate every pattern's time series per element.
481pub(crate) fn active_element_patterns<'a, const NDOF: usize, NId, EId, EL>(
482    patterns: impl Iterator<Item = &'a LoadPattern<NDOF, NId, EId, EL>>,
483    pseudo_time: f64,
484) -> Vec<(f64, &'a LoadPattern<NDOF, NId, EId, EL>)>
485where
486    NId: Eq + Hash + Copy + 'a,
487    EId: Eq + Hash + Copy + 'a,
488    EL: Copy + std::ops::Add<Output = EL> + 'a,
489{
490    patterns
491        .filter(|pattern| pattern.has_element_loads())
492        .map(|pattern| (pattern.factor(pseudo_time), pattern))
493        .filter(|(factor, _)| *factor != 0.0)
494        .collect()
495}
496
497/// An element's *effective* load: every active pattern's load on it, each
498/// scaled by that pattern's factor, summed — `None` if no active pattern
499/// loads it. What state-dependent elements (`ForceBeamColumn`) need inside
500/// their own state determination.
501pub(crate) fn effective_element_load<const NDOF: usize, NId, EId, EL>(
502    active: &[(f64, &LoadPattern<NDOF, NId, EId, EL>)],
503    element: EId,
504) -> Option<EL>
505where
506    NId: Eq + Hash + Copy,
507    EId: Eq + Hash + Copy,
508    EL: Copy + std::ops::Add<Output = EL> + std::ops::Mul<f64, Output = EL>,
509{
510    let mut total: Option<EL> = None;
511    for (factor, pattern) in active {
512        if let Some(load) = pattern.element_load(element) {
513            let scaled = *load * *factor;
514            total = Some(match total {
515                Some(sum) => sum + scaled,
516                None => scaled,
517            });
518        }
519    }
520    total
521}
522
523#[cfg(test)]
524mod bounded_path_tests {
525    use super::LoadSeries;
526
527    /// OpenSees-style accelerogram bounds, interpolation and useLast; exact arithmetic.
528    #[test]
529    fn bounded_accelerogram_is_zero_outside_its_samples() {
530        for use_last in [false, true] {
531            let path = LoadSeries::BoundedPath {
532                times: vec![1.0, 2.0],
533                factors: vec![2.0, 4.0],
534                use_last,
535            };
536            assert_eq!(path.factor(0.0), 0.0);
537            assert_eq!(path.factor(1.0), 2.0);
538            assert_eq!(path.factor(1.5), 3.0);
539            assert_eq!(path.factor(2.0), 4.0);
540            assert_eq!(path.factor(3.0), if use_last { 4.0 } else { 0.0 });
541            assert_eq!(path.slope(0.0), 0.0);
542            assert_eq!(path.slope(1.5), 2.0);
543            assert_eq!(path.slope(3.0), 0.0);
544        }
545    }
546}
carapace_core/model/load_pattern.rs

carapace_core/model/
load_pattern.rs
