carapace_core/model/
integration.rs

1/// Numerical integration scheme along a frame element's length, relating
2/// per-section response (`FiberSection`) to element-level response — needed
3/// by `DispBeamColumn`. Closed enum.
4///
5/// Point/weight tables are copied directly from OpenSees's own
6/// (`xara/SRC/quadrature/Frame/{Legendre,Lobatto}BeamIntegration.cpp`),
7/// not computed via a general-purpose quadrature crate — domain-specific
8/// numerics like this are copied from the reference implementation, the
9/// same way element/material formulations are, not sourced from a generic
10/// library (unlike heavily-optimized linear algebra, e.g. `nalgebra`/
11/// `faer`, where using a library is the right call). Supports 2-6 points,
12/// matching the range `DispBeamColumn` actually needs in practice — the
13/// C++ source supports up to 10; extend the tables if that's ever needed.
14#[derive(Debug, Clone, Copy, PartialEq)]
15pub enum BeamIntegration {
16    /// Gauss-Legendre: interior points only, `n`-point rule exact to
17    /// degree `2n-1`. `DispBeamColumn`'s usual default.
18    Legendre { points: usize },
19    /// Gauss-Lobatto: includes both end points, `n`-point rule exact to
20    /// degree `2n-3` — favored for `ForceBeamColumn` since end-section
21    /// behavior matters most for plastic hinges there, but not needed by
22    /// `DispBeamColumn` itself. Included now since both were asked for
23    /// together and the table-driven implementation cost is the
24    /// same either way.
25    Lobatto { points: usize },
26}
27
28impl BeamIntegration {
29    /// `(xi, weight)` pairs, `xi` in `[0,1]` along the element length,
30    /// weights summing to 1 — already in the same mapped-to-`[0,1]` form
31    /// `getSectionLocations`/`getSectionWeights` return in the C++ source
32    /// (`xi = 0.5*(xi+1)`, `wt *= 0.5`), not the raw `[-1,1]` table.
33    pub fn points(&self) -> Vec<(f64, f64)> {
34        let (raw_xi, raw_wt): (&[f64], &[f64]) = match self {
35            BeamIntegration::Legendre { points } => legendre_table(*points),
36            BeamIntegration::Lobatto { points } => lobatto_table(*points),
37        };
38        raw_xi
39            .iter()
40            .zip(raw_wt)
41            .map(|(&xi, &wt)| (0.5 * (xi + 1.0), 0.5 * wt))
42            .collect()
43    }
44}
45
46fn legendre_table(points: usize) -> (&'static [f64], &'static [f64]) {
47    match points {
48        2 => (&[-0.577350269189626, 0.577350269189626], &[1.0, 1.0]),
49        3 => (
50            &[-0.774596669241483, 0.0, 0.774596669241483],
51            &[0.555555555555556, 0.888888888888889, 0.555555555555556],
52        ),
53        4 => (
54            &[
55                -0.861136311594053,
56                -0.339981043584856,
57                0.339981043584856,
58                0.861136311594053,
59            ],
60            &[
61                0.347854845137454,
62                0.652145154862546,
63                0.652145154862546,
64                0.347854845137454,
65            ],
66        ),
67        5 => (
68            &[
69                -0.906179845938664,
70                -0.538469310105683,
71                0.0,
72                0.538469310105683,
73                0.906179845938664,
74            ],
75            &[
76                0.236926885056189,
77                0.478628670499366,
78                0.568888888888889,
79                0.478628670499366,
80                0.236926885056189,
81            ],
82        ),
83        6 => (
84            &[
85                -0.932469514203152,
86                -0.661209386466265,
87                -0.238619186083197,
88                0.238619186083197,
89                0.661209386466265,
90                0.932469514203152,
91            ],
92            &[
93                0.171324492379170,
94                0.360761573048139,
95                0.467913934572691,
96                0.467913934572691,
97                0.360761573048139,
98                0.171324492379170,
99            ],
100        ),
101        _ => panic!("BeamIntegration::Legendre supports 2-6 points, got {points}"),
102    }
103}
104
105fn lobatto_table(points: usize) -> (&'static [f64], &'static [f64]) {
106    match points {
107        2 => (&[-1.0, 1.0], &[1.0, 1.0]),
108        3 => (
109            &[-1.0, 0.0, 1.0],
110            &[0.333333333333333, 1.333333333333333, 0.333333333333333],
111        ),
112        4 => (
113            &[-1.0, -0.44721360, 0.44721360, 1.0],
114            &[
115                0.166666666666667,
116                0.833333333333333,
117                0.833333333333333,
118                0.166666666666667,
119            ],
120        ),
121        5 => (
122            &[-1.0, -0.65465367, 0.0, 0.65465367, 1.0],
123            &[0.1, 0.5444444444, 0.7111111111, 0.5444444444, 0.1],
124        ),
125        6 => (
126            &[
127                -1.0,
128                -0.7650553239,
129                -0.2852315164,
130                0.2852315164,
131                0.7650553239,
132                1.0,
133            ],
134            &[
135                0.06666666667,
136                0.3784749562,
137                0.5548583770,
138                0.5548583770,
139                0.3784749562,
140                0.06666666667,
141            ],
142        ),
143        _ => panic!("BeamIntegration::Lobatto supports 2-6 points, got {points}"),
144    }
145}
146
147#[cfg(test)]
148mod tests {
149    use super::*;
150
151    #[test]
152    fn points_are_symmetric_and_weights_sum_to_one() {
153        for points in 2..=6 {
154            for integration in [
155                BeamIntegration::Legendre { points },
156                BeamIntegration::Lobatto { points },
157            ] {
158                let pts = integration.points();
159                assert_eq!(pts.len(), points);
160                let sum: f64 = pts.iter().map(|(_, w)| w).sum();
161                // 1e-8, not tighter: OpenSees's own source literals for
162                // higher-order Lobatto weights are truncated to ~10
163                // significant digits (e.g. `0.06666666667`), not full f64
164                // precision.
165                assert!(
166                    (sum - 1.0).abs() < 1e-8,
167                    "{integration:?} weights should sum to 1, got {sum}"
168                );
169            }
170        }
171    }
172
173    #[test]
174    fn lobatto_includes_the_endpoints() {
175        let pts = BeamIntegration::Lobatto { points: 4 }.points();
176        assert!((pts.first().unwrap().0 - 0.0).abs() < 1e-12);
177        assert!((pts.last().unwrap().0 - 1.0).abs() < 1e-12);
178    }
179}
carapace_core/model/integration.rs

carapace_core/model/
integration.rs
