carapace_core/model/elements/
truss.rs

1use nalgebra::{SMatrix, SVector};
2
3use super::super::{
4    Material, Node, Node3, Node3Id, NodeId, SpatialDof, ELEMENT_DOF, SPATIAL_ELEMENT_DOF,
5    SPATIAL_NDF, SPATIAL_NDIM,
6};
7
8/// Fixed-size spatial two-node element tangent/resistance pair types.
9pub type SpatialElementMatrix = SMatrix<f64, SPATIAL_ELEMENT_DOF, SPATIAL_ELEMENT_DOF>;
10pub type SpatialElementVector = SVector<f64, SPATIAL_ELEMENT_DOF>;
11
12/// A 2-node axial truss: fixed-size element-local linear algebra, no
13/// heap allocation in the hot path. Only ever populates the translational
14/// DOF entries of the (3-DOF-per-node) local/global matrices; the rotation
15/// entries stay zero. Its `dof_mask` is therefore translations only, so a
16/// node connected only to trusses has no rotation unknown and the model
17/// need not fix it.
18#[derive(Debug, Clone)]
19pub struct Truss {
20    pub node_i: NodeId,
21    pub node_j: NodeId,
22    pub area: f64,
23    pub material: Material,
24    /// Mass per unit volume. Zero (the default) means massless — existing
25    /// models are unaffected unless they opt in via `with_density`.
26    pub density: f64,
27}
28
29impl Truss {
30    pub fn new(node_i: NodeId, node_j: NodeId, area: f64, material: Material) -> Self {
31        Truss {
32            node_i,
33            node_j,
34            area,
35            material,
36            density: 0.0,
37        }
38    }
39
40    pub fn with_density(mut self, density: f64) -> Self {
41        self.density = density;
42        self
43    }
44
45    fn geometry(&self, node_i: &Node, node_j: &Node) -> (f64, f64, f64) {
46        let dx = node_j.coords[0] - node_i.coords[0];
47        let dy = node_j.coords[1] - node_i.coords[1];
48        let length = (dx * dx + dy * dy).sqrt();
49        (length, dx / length, dy / length)
50    }
51
52    /// Axial elongation from current nodal displacements, projected onto
53    /// the (undeformed) bar axis — small-displacement theory. Shared by
54    /// `form_tangent_and_resistance` (trial) and `commit`.
55    fn strain(&self, node_i: &Node, node_j: &Node) -> f64 {
56        let (length, cx, cy) = self.geometry(node_i, node_j);
57        ((node_j.displacement[0] - node_i.displacement[0]) * cx
58            + (node_j.displacement[1] - node_i.displacement[1]) * cy)
59            / length
60    }
61
62    pub(super) fn form_tangent_and_resistance(
63        &self,
64        node_i: &Node,
65        node_j: &Node,
66    ) -> (
67        SMatrix<f64, ELEMENT_DOF, ELEMENT_DOF>,
68        SVector<f64, ELEMENT_DOF>,
69    ) {
70        let (length, cx, cy) = self.geometry(node_i, node_j);
71        let strain = self.strain(node_i, node_j);
72        let (stress, tangent_modulus) = self.material.trial_stress_tangent(strain);
73
74        // Local DOF order [ux_i, uy_i, rz_i, ux_j, uy_j, rz_j]; rotation
75        // entries (2, 5) stay zero — a truss carries no moment.
76        let b = SVector::<f64, ELEMENT_DOF>::from_column_slice(&[-cx, -cy, 0.0, cx, cy, 0.0]);
77
78        let k = (tangent_modulus * self.area / length) * (b * b.transpose());
79        let resistance = (stress * self.area) * b;
80
81        (k, resistance)
82    }
83
84    pub(super) fn commit(&mut self, node_i: &Node, node_j: &Node) {
85        let strain = self.strain(node_i, node_j);
86        self.material = self.material.commit(strain);
87    }
88
89    /// A truss has no distinct local-axis frame to rotate into (its
90    /// direction cosines are already baked into `form_tangent_and_resistance`'s
91    /// `b`, unlike a beam-column's fixed `t`/current-chord transform) — its
92    /// resistance vector already *is* the local nodal force, so this is
93    /// that same value, recomputed fresh from the current committed strain
94    /// (cheap: one material evaluation, no history to cache).
95    pub(super) fn local_force(&self, node_i: &Node, node_j: &Node) -> SVector<f64, ELEMENT_DOF> {
96        self.form_tangent_and_resistance(node_i, node_j).1
97    }
98
99    /// Lumped mass: half the element's total mass (`density * area * length`)
100    /// at each node, split equally between that node's translational DOFs —
101    /// a point mass has no directional preference. Zero rotational
102    /// contribution (index 2, 5) — a truss carries no moment either.
103    pub(super) fn form_mass(&self, node_i: &Node, node_j: &Node) -> SVector<f64, ELEMENT_DOF> {
104        let (length, _cx, _cy) = self.geometry(node_i, node_j);
105        let half = self.density * self.area * length / 2.0;
106        SVector::<f64, ELEMENT_DOF>::from_column_slice(&[half, half, 0.0, half, half, 0.0])
107    }
108}
109
110/// A small-displacement axial 3D truss.
111///
112/// The bar has three translational components at each end; its rotational
113/// rows/columns remain zero. Its `dof_mask` is translations only, so a node
114/// connected only to `Truss3` elements has no rotation unknowns and the model
115/// need not fix them.
116#[derive(Debug, Clone)]
117pub struct Truss3 {
118    pub node_i: Node3Id,
119    pub node_j: Node3Id,
120    pub area: f64,
121    pub material: Material,
122    /// Mass per unit volume; zero means massless.
123    pub density: f64,
124}
125
126impl Truss3 {
127    pub fn new(node_i: Node3Id, node_j: Node3Id, area: f64, material: Material) -> Self {
128        Truss3 {
129            node_i,
130            node_j,
131            area,
132            material,
133            density: 0.0,
134        }
135    }
136
137    pub fn with_density(mut self, density: f64) -> Self {
138        self.density = density;
139        self
140    }
141
142    /// Initial length and global direction cosines. A zero-length bar is a
143    /// model construction error, not a valid spatial element.
144    fn geometry(&self, node_i: &Node3, node_j: &Node3) -> (f64, [f64; SPATIAL_NDIM]) {
145        let dx = node_j.coords[0] - node_i.coords[0];
146        let dy = node_j.coords[1] - node_i.coords[1];
147        let dz = node_j.coords[2] - node_i.coords[2];
148        let length = (dx * dx + dy * dy + dz * dz).sqrt();
149        assert!(length > 0.0, "Truss3 endpoints must not coincide");
150        (length, [dx / length, dy / length, dz / length])
151    }
152
153    fn axial_strain(&self, node_i: &Node3, node_j: &Node3) -> f64 {
154        let (length, [cx, cy, cz]) = self.geometry(node_i, node_j);
155        let relative = (node_j.displacement[SpatialDof::Ux as usize]
156            - node_i.displacement[SpatialDof::Ux as usize])
157            * cx
158            + (node_j.displacement[SpatialDof::Uy as usize]
159                - node_i.displacement[SpatialDof::Uy as usize])
160                * cy
161            + (node_j.displacement[SpatialDof::Uz as usize]
162                - node_i.displacement[SpatialDof::Uz as usize])
163                * cz;
164        relative / length
165    }
166
167    /// Element tangent and internal resistance in global DOF order
168    /// `[ux_i, uy_i, uz_i, rx_i, ry_i, rz_i, ux_j, ..., rz_j]`.
169    pub fn form_tangent_and_resistance(
170        &self,
171        node_i: &Node3,
172        node_j: &Node3,
173    ) -> (SpatialElementMatrix, SpatialElementVector) {
174        let (length, [cx, cy, cz]) = self.geometry(node_i, node_j);
175        let strain = self.axial_strain(node_i, node_j);
176        let (stress, tangent_modulus) = self.material.trial_stress_tangent(strain);
177
178        let mut b = SpatialElementVector::zeros();
179        b[SpatialDof::Ux as usize] = -cx;
180        b[SpatialDof::Uy as usize] = -cy;
181        b[SpatialDof::Uz as usize] = -cz;
182        b[SPATIAL_NDF + SpatialDof::Ux as usize] = cx;
183        b[SPATIAL_NDF + SpatialDof::Uy as usize] = cy;
184        b[SPATIAL_NDF + SpatialDof::Uz as usize] = cz;
185
186        let tangent = (tangent_modulus * self.area / length) * (b.clone() * b.transpose());
187        let resistance = (stress * self.area) * b;
188        (tangent, resistance)
189    }
190
191    pub fn commit(&mut self, node_i: &Node3, node_j: &Node3) {
192        self.material = self.material.commit(self.axial_strain(node_i, node_j));
193    }
194
195    /// See `Truss::local_force`'s doc comment — same reasoning, no separate
196    /// local frame to rotate into.
197    pub fn local_force(&self, node_i: &Node3, node_j: &Node3) -> SpatialElementVector {
198        self.form_tangent_and_resistance(node_i, node_j).1
199    }
200
201    /// Diagonal lumped mass: half of the member mass at each node, applied to
202    /// each of its three translational DOFs. Rotational inertia remains zero
203    /// until a spatial beam/mass formulation supplies it explicitly.
204    pub fn form_mass(&self, node_i: &Node3, node_j: &Node3) -> SpatialElementVector {
205        let (length, _) = self.geometry(node_i, node_j);
206        let half = self.density * self.area * length / 2.0;
207        let mut mass = SpatialElementVector::zeros();
208        for dof in [SpatialDof::Ux, SpatialDof::Uy, SpatialDof::Uz] {
209            mass[dof as usize] = half;
210            mass[SPATIAL_NDF + dof as usize] = half;
211        }
212        mass
213    }
214}
215
216#[cfg(test)]
217mod spatial_tests {
218    use super::*;
219
220    #[test]
221    fn skew_truss_matches_direction_cosine_stiffness_and_resistance() {
222        let node_i = Node3::new([0.0, 0.0, 0.0]);
223        let mut node_j = Node3::new([2.0, 3.0, 6.0]);
224        node_j.displacement[SpatialDof::Ux as usize] = 0.004;
225        node_j.displacement[SpatialDof::Uy as usize] = 0.006;
226        node_j.displacement[SpatialDof::Uz as usize] = 0.012;
227
228        let truss = Truss3::new(
229            Node3Id::default(),
230            Node3Id::default(),
231            2.0,
232            Material::Elastic { e: 30_000.0 },
233        );
234        let (k, r) = truss.form_tangent_and_resistance(&node_i, &node_j);
235
236        assert!((r[6] - 240.0 / 7.0).abs() < 1e-12);
237        assert!((r[7] - 360.0 / 7.0).abs() < 1e-12);
238        assert!((r[8] - 720.0 / 7.0).abs() < 1e-12);
239        assert_eq!(r[3], 0.0);
240        assert_eq!(r[9], 0.0);
241
242        let mut u = SpatialElementVector::zeros();
243        u[6] = 0.004;
244        u[7] = 0.006;
245        u[8] = 0.012;
246        assert!((k * u - r).norm() < 1e-12);
247    }
248
249    #[test]
250    fn mass_is_lumped_over_three_translations_only() {
251        let node_i = Node3::new([0.0, 0.0, 0.0]);
252        let node_j = Node3::new([0.0, 0.0, 4.0]);
253        let truss = Truss3::new(
254            Node3Id::default(),
255            Node3Id::default(),
256            3.0,
257            Material::Elastic { e: 1.0 },
258        )
259        .with_density(2.0);
260        let mass = truss.form_mass(&node_i, &node_j);
261
262        for i in [0, 1, 2, 6, 7, 8] {
263            assert_eq!(mass[i], 12.0);
264        }
265        for i in [3, 4, 5, 9, 10, 11] {
266            assert_eq!(mass[i], 0.0);
267        }
268    }
269}
carapace_core/model/elements/truss.rs

carapace_core/model/elements/
truss.rs
