carapace_core/model/elements/
elastic_beam_column.rs

1use nalgebra::{SMatrix, SVector};
2
3use super::super::transform::Corotational2d;
4use super::super::{GeomTransf, GeomTransf3, Node, Node3, Node3Id, NodeId};
5use super::truss::{SpatialElementMatrix, SpatialElementVector};
6
7/// A 2-node, prismatic, linear-elastic 2D beam-column (Euler-Bernoulli, no
8/// shear deformation): closed-form stiffness, no iteration — unlike
9/// `Truss`/`ZeroLength`, there's no `Material` dispatch here, since the
10/// element's response is fully determined by `e`/`a`/`iz` with no nonlinear
11/// stress-strain law to evaluate.
12#[derive(Debug, Clone)]
13pub struct ElasticBeamColumn {
14    pub node_i: NodeId,
15    pub node_j: NodeId,
16    pub e: f64,
17    pub a: f64,
18    pub iz: f64,
19    pub transform: GeomTransf,
20    /// Mass per unit volume. Zero (the default) means massless — existing
21    /// models are unaffected unless they opt in via `with_density`.
22    pub density: f64,
23}
24
25impl ElasticBeamColumn {
26    pub fn new(
27        node_i: NodeId,
28        node_j: NodeId,
29        e: f64,
30        a: f64,
31        iz: f64,
32        transform: GeomTransf,
33    ) -> Self {
34        ElasticBeamColumn {
35            node_i,
36            node_j,
37            e,
38            a,
39            iz,
40            transform,
41            density: 0.0,
42        }
43    }
44
45    pub fn with_density(mut self, density: f64) -> Self {
46        self.density = density;
47        self
48    }
49
50    fn geometry(&self, node_i: &Node, node_j: &Node) -> (f64, f64, f64) {
51        let dx = node_j.coords[0] - node_i.coords[0];
52        let dy = node_j.coords[1] - node_i.coords[1];
53        let length = (dx * dx + dy * dy).sqrt();
54        (length, dx / length, dy / length)
55    }
56
57    /// Local-DOF-order `[u1, v1, θ1, u2, v2, θ2]` elastic stiffness —
58    /// standard Euler-Bernoulli prismatic beam (Przemieniecki / Cook §2).
59    fn local_elastic_stiffness(&self, length: f64) -> SMatrix<f64, 6, 6> {
60        let (e, a, i) = (self.e, self.a, self.iz);
61        let l = length;
62        let ea_l = e * a / l;
63        let ei = e * i;
64
65        #[rustfmt::skip]
66        let k = SMatrix::<f64, 6, 6>::from_row_slice(&[
67             ea_l,             0.0,             0.0,  -ea_l,             0.0,             0.0,
68              0.0,  12.0*ei/l.powi(3),  6.0*ei/l.powi(2),   0.0, -12.0*ei/l.powi(3),  6.0*ei/l.powi(2),
69              0.0,   6.0*ei/l.powi(2),      4.0*ei/l,       0.0,  -6.0*ei/l.powi(2),      2.0*ei/l,
70            -ea_l,             0.0,             0.0,   ea_l,             0.0,             0.0,
71              0.0, -12.0*ei/l.powi(3), -6.0*ei/l.powi(2),   0.0,  12.0*ei/l.powi(3), -6.0*ei/l.powi(2),
72              0.0,   6.0*ei/l.powi(2),      2.0*ei/l,       0.0,  -6.0*ei/l.powi(2),      4.0*ei/l,
73        ]);
74        k
75    }
76
77    /// First-order (linearized) geometric-stiffness correction from the
78    /// current axial force `p` (tension positive), added to the elastic
79    /// bending block for `GeomTransf::PDelta` — the standard consistent
80    /// geometric stiffness matrix (e.g. Przemieniecki §2), *not* a full
81    /// corotational update. Only the tangent used to solve for the next
82    /// displacement increment is corrected this way; the resisting-force
83    /// recovery (`form_tangent_and_resistance`'s second return value)
84    /// stays purely elastic. Consistent path-following P-Delta (where the
85    /// resisting force itself must reflect the same correction) needs
86    /// Newton iteration to track correctly and is not implemented.
87    fn geometric_stiffness(&self, p: f64, length: f64) -> SMatrix<f64, 6, 6> {
88        let l = length;
89        let mut kg = SMatrix::<f64, 6, 6>::zeros();
90        #[rustfmt::skip]
91        let block = [
92            [ 6.0/5.0,       l/10.0,      -6.0/5.0,       l/10.0],
93            [   l/10.0,  2.0*l*l/15.0,      -l/10.0,   -l*l/30.0],
94            [-6.0/5.0,      -l/10.0,       6.0/5.0,      -l/10.0],
95            [   l/10.0,    -l*l/30.0,      -l/10.0,  2.0*l*l/15.0],
96        ];
97        let indices = [1, 2, 4, 5];
98        for (bi, &gi) in indices.iter().enumerate() {
99            for (bj, &gj) in indices.iter().enumerate() {
100                kg[(gi, gj)] = (p / l) * block[bi][bj];
101            }
102        }
103        kg
104    }
105
106    /// Block-diagonal rotation from global to local DOFs: `d_local = T * d_global`.
107    fn transformation(&self, cx: f64, cy: f64) -> SMatrix<f64, 6, 6> {
108        #[rustfmt::skip]
109        let block = SMatrix::<f64, 3, 3>::new(
110             cx,  cy, 0.0,
111            -cy,  cx, 0.0,
112            0.0, 0.0, 1.0,
113        );
114        let mut t = SMatrix::<f64, 6, 6>::zeros();
115        t.fixed_view_mut::<3, 3>(0, 0).copy_from(&block);
116        t.fixed_view_mut::<3, 3>(3, 3).copy_from(&block);
117        t
118    }
119
120    pub(super) fn form_tangent_and_resistance(
121        &self,
122        node_i: &Node,
123        node_j: &Node,
124    ) -> (SMatrix<f64, 6, 6>, SVector<f64, 6>) {
125        if self.transform == GeomTransf::Corotational {
126            let transform = Corotational2d::new(node_i, node_j);
127            let l = transform.initial_length();
128            let ea_l = self.e * self.a / l;
129            let eiz_l = self.e * self.iz / l;
130            let k_basic = SMatrix::<f64, 3, 3>::new(
131                ea_l,
132                0.0,
133                0.0,
134                0.0,
135                4.0 * eiz_l,
136                2.0 * eiz_l,
137                0.0,
138                2.0 * eiz_l,
139                4.0 * eiz_l,
140            );
141            let q = k_basic * transform.basic_deformation();
142            return (
143                transform.global_tangent(&k_basic, &q),
144                transform.global_resistance(&q),
145            );
146        }
147
148        let (length, cx, cy) = self.geometry(node_i, node_j);
149        let t = self.transformation(cx, cy);
150        let k_local = self.local_elastic_stiffness(length);
151
152        let d_global = SVector::<f64, 6>::from_row_slice(&[
153            node_i.displacement[0],
154            node_i.displacement[1],
155            node_i.displacement[2],
156            node_j.displacement[0],
157            node_j.displacement[1],
158            node_j.displacement[2],
159        ]);
160        let d_local = t * d_global;
161
162        let resistance_local = k_local * d_local;
163        let resistance = t.transpose() * resistance_local;
164
165        let k_total_local = match self.transform {
166            GeomTransf::Linear => k_local,
167            GeomTransf::PDelta => {
168                let axial_force = resistance_local[3]; // tension-positive axial force at node j
169                k_local + self.geometric_stiffness(axial_force, length)
170            }
171            GeomTransf::Corotational => unreachable!("handled above"),
172        };
173        let k = t.transpose() * k_total_local * t;
174
175        (k, resistance)
176    }
177
178    /// Local nodal force — `resistance_local`/`Corotational2d::
179    /// local_resistance`, the pre-transform value `form_tangent_and_
180    /// resistance` already computes and discards, recomputed fresh here
181    /// (stateless/elastic, so this is cheap and always exactly consistent
182    /// with the current committed displacement — no cache to keep in
183    /// sync).
184    pub(super) fn local_force(&self, node_i: &Node, node_j: &Node) -> SVector<f64, 6> {
185        if self.transform == GeomTransf::Corotational {
186            let transform = Corotational2d::new(node_i, node_j);
187            let l = transform.initial_length();
188            let ea_l = self.e * self.a / l;
189            let eiz_l = self.e * self.iz / l;
190            let k_basic = SMatrix::<f64, 3, 3>::new(
191                ea_l,
192                0.0,
193                0.0,
194                0.0,
195                4.0 * eiz_l,
196                2.0 * eiz_l,
197                0.0,
198                2.0 * eiz_l,
199                4.0 * eiz_l,
200            );
201            let q = k_basic * transform.basic_deformation();
202            return transform.local_resistance(&q);
203        }
204
205        let (length, cx, cy) = self.geometry(node_i, node_j);
206        let t = self.transformation(cx, cy);
207        let k_local = self.local_elastic_stiffness(length);
208        let d_global = SVector::<f64, 6>::from_row_slice(&[
209            node_i.displacement[0],
210            node_i.displacement[1],
211            node_i.displacement[2],
212            node_j.displacement[0],
213            node_j.displacement[1],
214            node_j.displacement[2],
215        ]);
216        k_local * (t * d_global)
217    }
218
219    /// Equivalent nodal load (global coordinates) from a uniform local-axis
220    /// load (`wx` axial, `wy` transverse in local +y; force/length) applied
221    /// to this element by whichever `LoadPattern` is currently being
222    /// assembled (`Domain::assemble_reference_load` passes the load in
223    /// — a pattern-scoped load can't live on `Element`, same reasoning as
224    /// `Node`'s load leaving `Node`), via consistent (virtual-work) Hermite
225    /// cubic / linear-axial shape-function integration.
226    pub(super) fn form_load_vector(
227        &self,
228        node_i: &Node,
229        node_j: &Node,
230        wx: f64,
231        wy: f64,
232    ) -> SVector<f64, 6> {
233        if wx == 0.0 && wy == 0.0 {
234            return SVector::<f64, 6>::zeros();
235        }
236        if self.transform == GeomTransf::Corotational {
237            let transform = Corotational2d::new(node_i, node_j);
238            return transform.global_to_local().transpose()
239                * Self::uniform_load_local(transform.initial_length(), wx, wy);
240        }
241        let (length, cx, cy) = self.geometry(node_i, node_j);
242        let t = self.transformation(cx, cy);
243        t.transpose() * Self::uniform_load_local(length, wx, wy)
244    }
245
246    /// `form_load_vector`'s equivalent nodal load in the element-local frame
247    /// `local_force` reports in (the current chord frame for `Corotational`).
248    pub(super) fn form_local_load_vector(
249        &self,
250        node_i: &Node,
251        node_j: &Node,
252        wx: f64,
253        wy: f64,
254    ) -> SVector<f64, 6> {
255        if self.transform == GeomTransf::Corotational {
256            let length = Corotational2d::new(node_i, node_j).initial_length();
257            return Self::uniform_load_local(length, wx, wy);
258        }
259        let (length, _cx, _cy) = self.geometry(node_i, node_j);
260        Self::uniform_load_local(length, wx, wy)
261    }
262
263    fn uniform_load_local(l: f64, wx: f64, wy: f64) -> SVector<f64, 6> {
264        super::uniform_load::planar_local(l, wx, wy)
265    }
266
267    /// Lumped mass: half the element's total mass (`density * a * length`)
268    /// at each node's translational DOFs, zero rotational contribution —
269    /// the simplest standard lumped-mass model;
270    /// a consistent (non-diagonal) mass matrix or a nonzero rotational
271    /// lumped inertia is a further refinement, not built until needed.
272    pub(super) fn form_mass(&self, node_i: &Node, node_j: &Node) -> SVector<f64, 6> {
273        let (length, _cx, _cy) = self.geometry(node_i, node_j);
274        let half = self.density * self.a * length / 2.0;
275        SVector::<f64, 6>::from_column_slice(&[half, half, 0.0, half, half, 0.0])
276    }
277}
278
279/// `ElasticBeamColumn`'s spatial counterpart: a 2-node, prismatic,
280/// linear-elastic 3D Euler-Bernoulli beam-column (`E, G, A, J, Iy, Iz`),
281/// closed-form, `GeomTransf3::Linear3` only — see this struct's
282/// module-level context for why 3D P-Delta isn't implemented yet: it needs
283/// force/tangent recovery already consistent in the linear case first.
284///
285/// Local DOF order per node `[u, v, w, rx, ry, rz]` (axial, then bending
286/// about local z causing local-y displacement `v`/`rz`, then bending about
287/// local y causing local-z displacement `w`/`ry`, torsion `rx`) — matches
288/// `SpatialDof`'s `[Ux, Uy, Uz, Rx, Ry, Rz]` order exactly, so no per-node
289/// DOF reordering is needed between local and global, only `transform`'s
290/// rotation. Verified against Xara/OpenSees's `ElasticBeam3d` (basic-
291/// stiffness terms, local DOF order).
292#[derive(Debug, Clone)]
293pub struct ElasticBeamColumn3 {
294    pub node_i: Node3Id,
295    pub node_j: Node3Id,
296    pub e: f64,
297    pub g: f64,
298    pub a: f64,
299    /// Torsion constant (St. Venant), not a polar moment of inertia for
300    /// non-circular sections.
301    pub j: f64,
302    /// Second moment of area about local y — governs bending that
303    /// displaces local z (`w`/`ry`).
304    pub iy: f64,
305    /// Second moment of area about local z — governs bending that
306    /// displaces local y (`v`/`rz`), the plane a 2D `ElasticBeamColumn`
307    /// implicitly bends in.
308    pub iz: f64,
309    pub transform: GeomTransf3,
310    /// Mass per unit volume; zero means massless.
311    pub density: f64,
312}
313
314impl ElasticBeamColumn3 {
315    #[allow(clippy::too_many_arguments)]
316    pub fn new(
317        node_i: Node3Id,
318        node_j: Node3Id,
319        e: f64,
320        g: f64,
321        a: f64,
322        j: f64,
323        iy: f64,
324        iz: f64,
325        transform: GeomTransf3,
326    ) -> Self {
327        ElasticBeamColumn3 {
328            node_i,
329            node_j,
330            e,
331            g,
332            a,
333            j,
334            iy,
335            iz,
336            transform,
337            density: 0.0,
338        }
339    }
340
341    pub fn with_density(mut self, density: f64) -> Self {
342        self.density = density;
343        self
344    }
345
346    /// Local-DOF-order 3D Euler-Bernoulli stiffness (no shear deformation).
347    /// Axial and torsion are simple 2×2 blocks; the two bending planes are
348    /// each the same 2D beam 4×4 block `ElasticBeamColumn::local_elastic_
349    /// stiffness` uses, keyed to `[v, rz]` (local z bending, `Iz`) or `[w,
350    /// ry]` (local y bending, `Iy`) — the `Iy` block's off-diagonal
351    /// (shear/moment coupling) terms carry the opposite sign from `Iz`'s,
352    /// the standard result of `y`/`z` forming a right-handed triad with
353    /// `x`: positive rotation about `+y` produces a *negative* `w` slope,
354    /// where positive rotation about `+z` produces a *positive* `v` slope.
355    fn local_elastic_stiffness(&self, length: f64) -> SpatialElementMatrix {
356        let (e, g, a, j, iy, iz) = (self.e, self.g, self.a, self.j, self.iy, self.iz);
357        let l = length;
358        let mut k = SpatialElementMatrix::zeros();
359
360        let ea_l = e * a / l;
361        k[(0, 0)] = ea_l;
362        k[(6, 6)] = ea_l;
363        k[(0, 6)] = -ea_l;
364        k[(6, 0)] = -ea_l;
365
366        let gj_l = g * j / l;
367        k[(3, 3)] = gj_l;
368        k[(9, 9)] = gj_l;
369        k[(3, 9)] = -gj_l;
370        k[(9, 3)] = -gj_l;
371
372        // Bending about local z: [v1, rz1, v2, rz2] = indices [1, 5, 7, 11].
373        let eiz = e * iz;
374        let z_dof = [1, 5, 7, 11];
375        #[rustfmt::skip]
376        let z_block = [
377            [ 12.0*eiz/l.powi(3),   6.0*eiz/l.powi(2), -12.0*eiz/l.powi(3),   6.0*eiz/l.powi(2)],
378            [  6.0*eiz/l.powi(2),       4.0*eiz/l,       -6.0*eiz/l.powi(2),      2.0*eiz/l],
379            [-12.0*eiz/l.powi(3),  -6.0*eiz/l.powi(2),  12.0*eiz/l.powi(3),  -6.0*eiz/l.powi(2)],
380            [  6.0*eiz/l.powi(2),       2.0*eiz/l,       -6.0*eiz/l.powi(2),      4.0*eiz/l],
381        ];
382        for (bi, &gi) in z_dof.iter().enumerate() {
383            for (bj, &gj) in z_dof.iter().enumerate() {
384                k[(gi, gj)] = z_block[bi][bj];
385            }
386        }
387
388        // Bending about local y: [w1, ry1, w2, ry2] = indices [2, 4, 8, 10];
389        // off-diagonal signs flipped relative to the z block (see doc comment).
390        let eiy = e * iy;
391        let y_dof = [2, 4, 8, 10];
392        #[rustfmt::skip]
393        let y_block = [
394            [ 12.0*eiy/l.powi(3),  -6.0*eiy/l.powi(2), -12.0*eiy/l.powi(3),  -6.0*eiy/l.powi(2)],
395            [ -6.0*eiy/l.powi(2),       4.0*eiy/l,        6.0*eiy/l.powi(2),      2.0*eiy/l],
396            [-12.0*eiy/l.powi(3),   6.0*eiy/l.powi(2),  12.0*eiy/l.powi(3),   6.0*eiy/l.powi(2)],
397            [ -6.0*eiy/l.powi(2),       2.0*eiy/l,        6.0*eiy/l.powi(2),      4.0*eiy/l],
398        ];
399        for (bi, &gi) in y_dof.iter().enumerate() {
400            for (bj, &gj) in y_dof.iter().enumerate() {
401                k[(gi, gj)] = y_block[bi][bj];
402            }
403        }
404
405        k
406    }
407
408    /// First-order (linearized) geometric-stiffness correction from the
409    /// current axial force `p` (tension positive), added to the elastic
410    /// bending blocks for `GeomTransf3::PDelta3` — the spatial analogue of
411    /// planar `ElasticBeamColumn::geometric_stiffness` (see its doc comment
412    /// for the "linearized, not full corotational" caveat, which applies
413    /// identically here). Corrects *both* bending planes: the standard
414    /// consistent geometric-stiffness block depends only on the transverse
415    /// shape functions and axial force, not on `EI`, so it's the exact same
416    /// numeric block in each plane — but expressed in `[w, ry]` rather than
417    /// `[v, rz]`, its off-diagonal (shear/rotation coupling) entries flip
418    /// sign for the same reason `local_elastic_stiffness`'s `y_block` does
419    /// relative to `z_block` (see that method's doc comment): `ry = -dw/dx`
420    /// while `rz = +dv/dx`, and this block's off-diagonal terms are each
421    /// linear in exactly one rotation, so substituting the sign-flipped
422    /// rotation flips exactly those terms, leaving the pure-translation and
423    /// pure-rotation entries unchanged.
424    fn geometric_stiffness(&self, p: f64, length: f64) -> SpatialElementMatrix {
425        let l = length;
426        let mut kg = SpatialElementMatrix::zeros();
427
428        // Bending about local z: [v1, rz1, v2, rz2] = indices [1, 5, 7, 11].
429        let z_dof = [1, 5, 7, 11];
430        #[rustfmt::skip]
431        let z_block = [
432            [ 6.0/5.0,       l/10.0,      -6.0/5.0,       l/10.0],
433            [   l/10.0,  2.0*l*l/15.0,      -l/10.0,   -l*l/30.0],
434            [-6.0/5.0,      -l/10.0,       6.0/5.0,      -l/10.0],
435            [   l/10.0,    -l*l/30.0,      -l/10.0,  2.0*l*l/15.0],
436        ];
437        for (bi, &gi) in z_dof.iter().enumerate() {
438            for (bj, &gj) in z_dof.iter().enumerate() {
439                kg[(gi, gj)] = (p / l) * z_block[bi][bj];
440            }
441        }
442
443        // Bending about local y: [w1, ry1, w2, ry2] = indices [2, 4, 8, 10];
444        // off-diagonal (w-ry coupling) signs flipped relative to the z
445        // block — see this method's doc comment.
446        let y_dof = [2, 4, 8, 10];
447        #[rustfmt::skip]
448        let y_block = [
449            [ 6.0/5.0,      -l/10.0,      -6.0/5.0,      -l/10.0],
450            [  -l/10.0,  2.0*l*l/15.0,       l/10.0,   -l*l/30.0],
451            [-6.0/5.0,       l/10.0,       6.0/5.0,       l/10.0],
452            [  -l/10.0,    -l*l/30.0,       l/10.0,  2.0*l*l/15.0],
453        ];
454        for (bi, &gi) in y_dof.iter().enumerate() {
455            for (bj, &gj) in y_dof.iter().enumerate() {
456                kg[(gi, gj)] = (p / l) * y_block[bi][bj];
457            }
458        }
459
460        kg
461    }
462
463    pub(super) fn form_tangent_and_resistance(
464        &self,
465        node_i: &Node3,
466        node_j: &Node3,
467    ) -> (SpatialElementMatrix, SpatialElementVector) {
468        let (length, r) = self.transform.local_axes(node_i, node_j);
469        let t = GeomTransf3::rotation_matrix(&r);
470        let k_local = self.local_elastic_stiffness(length);
471
472        let d_global = SpatialElementVector::from_iterator(
473            node_i
474                .displacement
475                .iter()
476                .chain(node_j.displacement.iter())
477                .copied(),
478        );
479        let d_local = t * d_global;
480
481        let resistance_local = k_local * d_local;
482        let resistance = t.transpose() * resistance_local;
483
484        let k_total_local = match self.transform {
485            GeomTransf3::Linear3 { .. } => k_local,
486            GeomTransf3::PDelta3 { .. } => {
487                let axial_force = resistance_local[6]; // tension-positive axial force at node j
488                k_local + self.geometric_stiffness(axial_force, length)
489            }
490        };
491        let k = t.transpose() * k_total_local * t;
492
493        (k, resistance)
494    }
495
496    /// Local nodal force, spatial counterpart of the planar `local_force`
497    /// above — same `resistance_local` value `form_tangent_and_resistance`
498    /// already computes and discards, recomputed fresh (cheap, stateless).
499    pub(super) fn local_force(&self, node_i: &Node3, node_j: &Node3) -> SpatialElementVector {
500        let (length, r) = self.transform.local_axes(node_i, node_j);
501        let t = GeomTransf3::rotation_matrix(&r);
502        let k_local = self.local_elastic_stiffness(length);
503        let d_global = SpatialElementVector::from_iterator(
504            node_i
505                .displacement
506                .iter()
507                .chain(node_j.displacement.iter())
508                .copied(),
509        );
510        k_local * (t * d_global)
511    }
512
513    /// Equivalent nodal load (global coordinates) from a uniform local-axis
514    /// load: axial `wx` plus transverse `wy`/`wz` (force/length, local
515    /// `+y`/`+z` directions). The transverse part is the
516    /// direct spatial generalization of planar `ElasticBeamColumn::
517    /// form_load_vector`'s consistent (virtual-work) Hermite-cubic
518    /// equivalent load, applied independently in each bending plane. The
519    /// `wy`-induced terms (`[v, rz]`, indices `[1,5,7,11]`) are numerically
520    /// identical to the planar formula; the `wz`-induced terms (`[w, ry]`,
521    /// indices `[2,4,8,10]`) carry the same rotation-coefficient sign flip
522    /// as `local_elastic_stiffness`'s `y_block` relative to `z_block` (`ry
523    /// = -dw/dx` vs `rz = +dv/dx`) — verified against Xara/OpenSees's
524    /// `ElasticBeam3d::addLoad`'s `Beam3dUniformLoad` fixed-end-force
525    /// derivation (`q0`: `MI=-Mz, MJ=Mz` for the `wy`-driven `z`-bending
526    /// pair and `MI=My, MJ=-My` for the `wz`-driven `y`-bending pair, each
527    /// negated here since a fixed-end restraint force is the negative of
528    /// the equivalent nodal load it corresponds to), not just asserted by
529    /// analogy to the planar case.
530    pub(super) fn form_load_vector(
531        &self,
532        node_i: &Node3,
533        node_j: &Node3,
534        wx: f64,
535        wy: f64,
536        wz: f64,
537    ) -> SpatialElementVector {
538        if wx == 0.0 && wy == 0.0 && wz == 0.0 {
539            return SpatialElementVector::zeros();
540        }
541        let (length, r) = self.transform.local_axes(node_i, node_j);
542        let t = GeomTransf3::rotation_matrix(&r);
543        t.transpose() * Self::uniform_load_local(length, wx, wy, wz)
544    }
545
546    /// `form_load_vector`'s equivalent nodal load in the element-local
547    /// frame `local_force` reports in.
548    pub(super) fn form_local_load_vector(
549        &self,
550        node_i: &Node3,
551        node_j: &Node3,
552        wx: f64,
553        wy: f64,
554        wz: f64,
555    ) -> SpatialElementVector {
556        let (length, _r) = self.transform.local_axes(node_i, node_j);
557        Self::uniform_load_local(length, wx, wy, wz)
558    }
559
560    fn uniform_load_local(l: f64, wx: f64, wy: f64, wz: f64) -> SpatialElementVector {
561        super::uniform_load::spatial_local(l, wx, wy, wz)
562    }
563
564    /// Lumped mass: half the element's total mass at each node's three
565    /// translational DOFs, zero rotational contribution — same
566    /// simplification as planar `ElasticBeamColumn::form_mass`.
567    pub(super) fn form_mass(&self, node_i: &Node3, node_j: &Node3) -> SpatialElementVector {
568        let (length, _r) = self.transform.local_axes(node_i, node_j);
569        let half = self.density * self.a * length / 2.0;
570        let mut mass = SpatialElementVector::zeros();
571        for dof in [0, 1, 2, 6, 7, 8] {
572            mass[dof] = half;
573        }
574        mass
575    }
576}
577
578#[cfg(test)]
579mod spatial_tests {
580    use nalgebra::Matrix3;
581
582    use super::*;
583
584    /// A cantilever along global x with `vec_xz = [0, 0, 1]` (global z),
585    /// which makes local `[x, y, z] = global [x, y, z]` exactly (`R` is the
586    /// identity) — isolates the local stiffness formula itself from the
587    /// `GeomTransf3` rotation, which `transform.rs`'s and `truss.rs`'s own
588    /// tests already cover for the transform-only case. With `node_i` fully
589    /// fixed, the free-DOF system reduces to `local_elastic_stiffness`'s
590    /// `[6..12, 6..12]` block applied directly (no rotation needed), so tip
591    /// displacements can be checked against textbook cantilever closed
592    /// forms for axial, biaxial bending, and torsion simultaneously (all
593    /// four modes are decoupled in a prismatic doubly-symmetric section).
594    #[test]
595    fn axis_aligned_cantilever_matches_closed_form_axial_biaxial_bending_and_torsion() {
596        let (e, g, a, j, iy, iz) = (30_000.0, 12_000.0, 20.0, 5.0, 400.0, 800.0);
597        let length = 100.0;
598        let transform = GeomTransf3::linear([0.0, 0.0, 1.0]);
599        let beam = ElasticBeamColumn3::new(
600            Node3Id::default(),
601            Node3Id::default(),
602            e,
603            g,
604            a,
605            j,
606            iy,
607            iz,
608            transform,
609        );
610
611        let node_i = Node3::new([0.0, 0.0, 0.0]);
612        let node_j = Node3::new([length, 0.0, 0.0]);
613        let (geom_length, r) = beam.transform.local_axes(&node_i, &node_j);
614        assert_eq!(geom_length, length);
615        assert!(
616            (r - Matrix3::identity()).norm() < 1e-12,
617            "vec_xz=[0,0,1] with a member along x must give R = identity"
618        );
619
620        let k_local = beam.local_elastic_stiffness(length);
621
622        // node_i fixed at zero, so the free (node_j) system is exactly
623        // k_local's bottom-right 6x6 block.
624        let k_free = k_local.fixed_view::<6, 6>(6, 6).into_owned();
625        let (axial_force, fy, fz, torque) = (100.0, 50.0, 30.0, 20.0);
626        // [ux, uy, uz, rx, ry, rz] at node_j; no applied end moments (my, mz).
627        let f = SVector::<f64, 6>::from_row_slice(&[axial_force, fy, fz, torque, 0.0, 0.0]);
628        let d = k_free
629            .lu()
630            .solve(&f)
631            .expect("k_free is nonsingular for a real prismatic section");
632
633        let expected_ux = axial_force * length / (e * a);
634        let expected_uy = fy * length.powi(3) / (3.0 * e * iz);
635        let expected_rz = fy * length.powi(2) / (2.0 * e * iz);
636        let expected_uz = fz * length.powi(3) / (3.0 * e * iy);
637        let expected_ry = -fz * length.powi(2) / (2.0 * e * iy);
638        let expected_rx = torque * length / (g * j);
639
640        let tol = 1e-9;
641        assert!(
642            (d[0] - expected_ux).abs() < tol,
643            "axial: {} vs {}",
644            d[0],
645            expected_ux
646        );
647        assert!(
648            (d[1] - expected_uy).abs() < tol,
649            "z-bending deflection: {} vs {}",
650            d[1],
651            expected_uy
652        );
653        assert!(
654            (d[5] - expected_rz).abs() < tol,
655            "z-bending rotation: {} vs {}",
656            d[5],
657            expected_rz
658        );
659        assert!(
660            (d[2] - expected_uz).abs() < tol,
661            "y-bending deflection: {} vs {}",
662            d[2],
663            expected_uz
664        );
665        assert!(
666            (d[4] - expected_ry).abs() < tol,
667            "y-bending rotation: {} vs {}",
668            d[4],
669            expected_ry
670        );
671        assert!(
672            (d[3] - expected_rx).abs() < tol,
673            "torsion: {} vs {}",
674            d[3],
675            expected_rx
676        );
677    }
678
679    /// The same cantilever driven through the full `form_tangent_and_
680    /// resistance` (transform included, not just the local block) with a
681    /// skew member and a non-trivial orientation, checked by round-tripping:
682    /// applying the closed-form local displacement (rotated into global
683    /// coordinates) must reproduce the applied global load exactly, proving
684    /// `GeomTransf3`/the local stiffness compose correctly rather than
685    /// being separately-plausible but inconsistent.
686    #[test]
687    fn skew_cantilever_transform_round_trips_a_known_local_solution() {
688        let (e, g, a, j, iy, iz) = (30_000.0, 12_000.0, 20.0, 5.0, 400.0, 800.0);
689        let node_i = Node3::new([0.0, 0.0, 0.0]);
690        let node_j = Node3::new([3.0, 4.0, 12.0]); // length 13, skew in all three axes
691        let transform = GeomTransf3::linear([1.0, 0.0, 0.0]);
692        let beam = ElasticBeamColumn3::new(
693            Node3Id::default(),
694            Node3Id::default(),
695            e,
696            g,
697            a,
698            j,
699            iy,
700            iz,
701            transform,
702        );
703
704        let (length, r) = beam.transform.local_axes(&node_i, &node_j);
705        assert!((length - 13.0).abs() < 1e-12);
706
707        // A purely local axial elongation at node_j: local d = [delta,0,0,0,0,0].
708        let delta = 0.005;
709        let mut local_d = SVector::<f64, 6>::zeros();
710        local_d[0] = delta;
711        let global_d_j = r.transpose() * local_d.fixed_rows::<3>(0).into_owned();
712        let global_rot_j = r.transpose() * local_d.fixed_rows::<3>(3).into_owned();
713
714        let mut node_j_displaced = node_j.clone();
715        node_j_displaced.displacement[0] = global_d_j[0];
716        node_j_displaced.displacement[1] = global_d_j[1];
717        node_j_displaced.displacement[2] = global_d_j[2];
718        node_j_displaced.displacement[3] = global_rot_j[0];
719        node_j_displaced.displacement[4] = global_rot_j[1];
720        node_j_displaced.displacement[5] = global_rot_j[2];
721
722        let (_, resistance) = beam.form_tangent_and_resistance(&node_i, &node_j_displaced);
723
724        // Pure axial elongation should produce a pure axial force pair,
725        // recoverable via the same projection `Truss3`'s skew test uses.
726        let expected_axial_force = e * a * delta / length;
727        let axial_disp_of_force = resistance
728            .fixed_rows::<3>(6)
729            .into_owned()
730            .dot(&r.row(0).transpose());
731        assert!((axial_disp_of_force - expected_axial_force).abs() < 1e-9);
732    }
733}
carapace_core/model/elements/elastic_beam_column.rs

carapace_core/model/elements/
elastic_beam_column.rs
