carapace_core/model/elements/
disp_beam_column.rs

1use nalgebra::{SMatrix, SVector};
2
3use super::super::transform::Corotational2d;
4use super::super::{
5    BeamIntegration, Fiber, Fiber3, FiberSection, FiberSection3, GeomTransf, GeomTransf3, Node,
6    Node3, Node3Id, NodeId,
7};
8use super::truss::{SpatialElementMatrix, SpatialElementVector};
9use super::uniform_load;
10
11/// A 2-node, displacement-based, fiber-discretized 2D beam-column:
12/// nodal displacements directly give the strain/curvature field along the
13/// length (cubic Hermite transverse + linear axial shape functions, same
14/// as `ElasticBeamColumn`'s), each `BeamIntegration` point's `FiberSection`
15/// converts that to stress-resultants, and one length-integral gives the
16/// element's resisting force/tangent — a single, direct evaluation, unlike
17/// `ForceBeamColumn`, which needs its own internal equilibrium
18/// iteration per element.
19///
20/// Prismatic-member assumption: every integration point starts from the
21/// same fiber layout (`fibers`, replicated `points` times) but evolves
22/// independently thereafter, since each point's material history depends
23/// on that point's own strain path. A tapered/non-prismatic member (a
24/// different section per point) is a natural extension, not built until
25/// needed.
26///
27/// Small-displacement by default, with an opt-in `GeomTransf::Corotational`
28/// formulation that integrates fiber response from the three objective
29/// basic deformations. `PDelta` remains unsupported for this element (a
30/// fiber section's force state does not reduce to the scalar axial force
31/// used by `ElasticBeamColumn`'s closed-form correction).
32/// Uniform element loads use the same consistent equivalent nodal load as
33/// `ElasticBeamColumn` (`form_load_vector`), which is independent of the
34/// section model and integration rule. Note the field is still the
35/// Hermite-interpolated linear curvature: a distributed load's parabolic
36/// moment is only approximated within one element, so nonlinear response under
37/// distributed load wants the member subdivided.
38#[derive(Debug, Clone)]
39pub struct DispBeamColumn {
40    pub node_i: NodeId,
41    pub node_j: NodeId,
42    integration: BeamIntegration,
43    sections: Vec<FiberSection>,
44    transform: GeomTransf,
45    /// Mass per unit volume, applied uniformly over the section's total
46    /// fiber area. Zero (the default) means massless.
47    pub density: f64,
48}
49
50impl DispBeamColumn {
51    pub fn new(
52        node_i: NodeId,
53        node_j: NodeId,
54        fibers: Vec<Fiber>,
55        integration: BeamIntegration,
56    ) -> Self {
57        let n_points = integration.points().len();
58        let sections = (0..n_points)
59            .map(|_| FiberSection::new(fibers.clone()))
60            .collect();
61        DispBeamColumn {
62            node_i,
63            node_j,
64            integration,
65            sections,
66            transform: GeomTransf::Linear,
67            density: 0.0,
68        }
69    }
70
71    pub fn with_density(mut self, density: f64) -> Self {
72        self.density = density;
73        self
74    }
75
76    /// Use finite-rotation, chord-following kinematics for this member.
77    pub fn with_corotational(mut self) -> Self {
78        self.transform = GeomTransf::Corotational;
79        self
80    }
81
82    fn geometry(&self, node_i: &Node, node_j: &Node) -> (f64, f64, f64) {
83        let dx = node_j.coords[0] - node_i.coords[0];
84        let dy = node_j.coords[1] - node_i.coords[1];
85        let length = (dx * dx + dy * dy).sqrt();
86        (length, dx / length, dy / length)
87    }
88
89    fn transformation(&self, cx: f64, cy: f64) -> SMatrix<f64, 6, 6> {
90        #[rustfmt::skip]
91        let block = SMatrix::<f64, 3, 3>::new(
92             cx,  cy, 0.0,
93            -cy,  cx, 0.0,
94            0.0, 0.0, 1.0,
95        );
96        let mut t = SMatrix::<f64, 6, 6>::zeros();
97        t.fixed_view_mut::<3, 3>(0, 0).copy_from(&block);
98        t.fixed_view_mut::<3, 3>(3, 3).copy_from(&block);
99        t
100    }
101
102    /// Strain-displacement vectors `(b_eps0, b_kappa)` at `xi` in `[0,1]`
103    /// along the length, relating local nodal displacements
104    /// `[u1,v1,th1,u2,v2,th2]` to `[eps0, kappa]` — `b_eps0` from the
105    /// linear axial shape functions (constant), `b_kappa` the second
106    /// derivative of the same cubic Hermite transverse shape functions
107    /// `ElasticBeamColumn` uses, derived directly (not copied from a
108    /// reference — this is standard Euler-Bernoulli/Hermite beam theory,
109    /// verified against `ElasticBeamColumn`'s closed-form stiffness in
110    /// `core/tests/elements/disp_beam.rs`: for constant `EA`/`EI` a rule
111    /// of 2 or more points must reproduce it exactly, since `b_kappa` is
112    /// linear in `xi` and its square is then only degree 2).
113    fn strain_displacement(xi: f64, length: f64) -> (SVector<f64, 6>, SVector<f64, 6>) {
114        let l = length;
115        let b_eps0 = SVector::<f64, 6>::from_column_slice(&[-1.0 / l, 0.0, 0.0, 1.0 / l, 0.0, 0.0]);
116        #[rustfmt::skip]
117        let b_kappa = SVector::<f64, 6>::from_column_slice(&[
118            0.0,
119            (-6.0 + 12.0 * xi) / (l * l),
120            (-4.0 + 6.0 * xi) / l,
121            0.0,
122            (6.0 - 12.0 * xi) / (l * l),
123            (-2.0 + 6.0 * xi) / l,
124        ]);
125        (b_eps0, b_kappa)
126    }
127
128    fn local_displacement(
129        &self,
130        node_i: &Node,
131        node_j: &Node,
132    ) -> (f64, SMatrix<f64, 6, 6>, SVector<f64, 6>) {
133        let (length, cx, cy) = self.geometry(node_i, node_j);
134        let t = self.transformation(cx, cy);
135        let d_global = SVector::<f64, 6>::from_row_slice(&[
136            node_i.displacement[0],
137            node_i.displacement[1],
138            node_i.displacement[2],
139            node_j.displacement[0],
140            node_j.displacement[1],
141            node_j.displacement[2],
142        ]);
143        (length, t, t * d_global)
144    }
145
146    pub(super) fn form_tangent_and_resistance(
147        &self,
148        node_i: &Node,
149        node_j: &Node,
150    ) -> (SMatrix<f64, 6, 6>, SVector<f64, 6>) {
151        let corotational = (self.transform == GeomTransf::Corotational)
152            .then(|| Corotational2d::new(node_i, node_j));
153        let (length, t, d_local) = if let Some(state) = &corotational {
154            (
155                state.initial_length(),
156                SMatrix::<f64, 6, 6>::identity(),
157                Corotational2d::local_basic_modes() * state.basic_deformation(),
158            )
159        } else {
160            self.local_displacement(node_i, node_j)
161        };
162
163        let mut k_local = SMatrix::<f64, 6, 6>::zeros();
164        let mut r_local = SVector::<f64, 6>::zeros();
165
166        for ((xi, w), section) in self.integration.points().iter().zip(&self.sections) {
167            let (b_eps0, b_kappa) = Self::strain_displacement(*xi, length);
168            let eps0 = b_eps0.dot(&d_local);
169            let kappa = b_kappa.dot(&d_local);
170            let (n, m, k_section) = section.trial(eps0, kappa);
171
172            let scale = w * length;
173            r_local += scale * (b_eps0 * n + b_kappa * m);
174            k_local += scale
175                * (k_section[0][0] * (b_eps0 * b_eps0.transpose())
176                    + k_section[0][1] * (b_eps0 * b_kappa.transpose())
177                    + k_section[1][0] * (b_kappa * b_eps0.transpose())
178                    + k_section[1][1] * (b_kappa * b_kappa.transpose()));
179        }
180
181        if let Some(state) = &corotational {
182            let modes = Corotational2d::local_basic_modes();
183            let basic_tangent = modes.transpose() * k_local * modes;
184            let basic_resistance = modes.transpose() * r_local;
185            (
186                state.global_tangent(&basic_tangent, &basic_resistance),
187                state.global_resistance(&basic_resistance),
188            )
189        } else {
190            (t.transpose() * k_local * t, t.transpose() * r_local)
191        }
192    }
193
194    pub(super) fn commit(&mut self, node_i: &Node, node_j: &Node) {
195        let (length, d_local) = if self.transform == GeomTransf::Corotational {
196            let state = Corotational2d::new(node_i, node_j);
197            (
198                state.initial_length(),
199                Corotational2d::local_basic_modes() * state.basic_deformation(),
200            )
201        } else {
202            let (length, _t, d_local) = self.local_displacement(node_i, node_j);
203            (length, d_local)
204        };
205        let points = self.integration.points();
206        for ((xi, _w), section) in points.iter().zip(&mut self.sections) {
207            let (b_eps0, b_kappa) = Self::strain_displacement(*xi, length);
208            let eps0 = b_eps0.dot(&d_local);
209            let kappa = b_kappa.dot(&d_local);
210            section.commit(eps0, kappa);
211        }
212    }
213
214    /// Local nodal force — the same `r_local`/`basic_resistance` value
215    /// `form_tangent_and_resistance` already computes and discards,
216    /// recomputed fresh here (cheaper than that method, in fact: no
217    /// tangent accumulation needed). Always exactly consistent with the
218    /// current committed section state — no cache to keep in sync, since
219    /// this element has no Newton loop between a target and the committed
220    /// answer the way `ForceBeamColumn` does.
221    pub(super) fn local_force(&self, node_i: &Node, node_j: &Node) -> SVector<f64, 6> {
222        let corotational = (self.transform == GeomTransf::Corotational)
223            .then(|| Corotational2d::new(node_i, node_j));
224        let (length, d_local) = if let Some(state) = &corotational {
225            (
226                state.initial_length(),
227                Corotational2d::local_basic_modes() * state.basic_deformation(),
228            )
229        } else {
230            let (length, _t, d_local) = self.local_displacement(node_i, node_j);
231            (length, d_local)
232        };
233
234        let mut r_local = SVector::<f64, 6>::zeros();
235        for ((xi, w), section) in self.integration.points().iter().zip(&self.sections) {
236            let (b_eps0, b_kappa) = Self::strain_displacement(*xi, length);
237            let eps0 = b_eps0.dot(&d_local);
238            let kappa = b_kappa.dot(&d_local);
239            let (n, m, _k_section) = section.trial(eps0, kappa);
240            r_local += (w * length) * (b_eps0 * n + b_kappa * m);
241        }
242
243        // Non-corotational: `r_local` is already the local nodal force
244        // (`d_local` came from the fixed-orientation `t`, same as
245        // `form_tangent_and_resistance`'s pre-transform value) — no further
246        // rotation needed, unlike that method's `t.transpose() * r_local`
247        // step, which converts *to* global.
248        match &corotational {
249            Some(state) => {
250                state.local_resistance(&(Corotational2d::local_basic_modes().transpose() * r_local))
251            }
252            None => r_local,
253        }
254    }
255
256    /// Equivalent nodal load (global coordinates) from a uniform local-axis
257    /// load (`wx` axial, `wy` transverse). This is the same consistent
258    /// virtual-work vector `ElasticBeamColumn` uses — it depends only on the
259    /// load and member length, not on the section model or integration rule
260    /// (OpenSees's `DispBeamColumn2d` fixed-end forces `p0` are likewise
261    /// `wL/2`, `wL²/12`) — and, like the elastic beam, a corotational member
262    /// treats it as a follower load over the undeformed length.
263    pub(super) fn form_load_vector(
264        &self,
265        node_i: &Node,
266        node_j: &Node,
267        wx: f64,
268        wy: f64,
269    ) -> SVector<f64, 6> {
270        if wx == 0.0 && wy == 0.0 {
271            return SVector::<f64, 6>::zeros();
272        }
273        if self.transform == GeomTransf::Corotational {
274            let state = Corotational2d::new(node_i, node_j);
275            return state.global_to_local().transpose()
276                * uniform_load::planar_local(state.initial_length(), wx, wy);
277        }
278        let (length, cx, cy) = self.geometry(node_i, node_j);
279        self.transformation(cx, cy).transpose() * uniform_load::planar_local(length, wx, wy)
280    }
281
282    /// `form_load_vector` in the element-local frame `local_force` reports in.
283    pub(super) fn form_local_load_vector(
284        &self,
285        node_i: &Node,
286        node_j: &Node,
287        wx: f64,
288        wy: f64,
289    ) -> SVector<f64, 6> {
290        let length = if self.transform == GeomTransf::Corotational {
291            Corotational2d::new(node_i, node_j).initial_length()
292        } else {
293            self.geometry(node_i, node_j).0
294        };
295        uniform_load::planar_local(length, wx, wy)
296    }
297
298    pub(super) fn form_mass(&self, node_i: &Node, node_j: &Node) -> SVector<f64, 6> {
299        let (length, _cx, _cy) = self.geometry(node_i, node_j);
300        // Prismatic assumption (see the type doc comment): every section
301        // has the same fiber layout, so any one of them gives the total
302        // area.
303        let total_area = self.sections[0].total_area();
304        let half = self.density * total_area * length / 2.0;
305        SVector::<f64, 6>::from_column_slice(&[half, half, 0.0, half, half, 0.0])
306    }
307
308    /// Every integration point's per-fiber `(strain, stress)`, freshly
309    /// recomputed from `node_i`/`node_j`'s current committed displacement —
310    /// same "recompute, don't cache" reasoning as `local_force` (this
311    /// element has no Newton loop between a target and the committed
312    /// answer, unlike `ForceBeamColumn`). Outer `Vec` is one entry per
313    /// integration point (`BeamIntegration::points()` order), inner `Vec`
314    /// one entry per fiber (the order originally passed to `new`).
315    pub fn fiber_responses(&self, node_i: &Node, node_j: &Node) -> Vec<Vec<(f64, f64)>> {
316        let corotational = (self.transform == GeomTransf::Corotational)
317            .then(|| Corotational2d::new(node_i, node_j));
318        let (length, d_local) = if let Some(state) = &corotational {
319            (
320                state.initial_length(),
321                Corotational2d::local_basic_modes() * state.basic_deformation(),
322            )
323        } else {
324            let (length, _t, d_local) = self.local_displacement(node_i, node_j);
325            (length, d_local)
326        };
327
328        self.integration
329            .points()
330            .iter()
331            .zip(&self.sections)
332            .map(|((xi, _w), section)| {
333                let (b_eps0, b_kappa) = Self::strain_displacement(*xi, length);
334                let eps0 = b_eps0.dot(&d_local);
335                let kappa = b_kappa.dot(&d_local);
336                section.fiber_responses(eps0, kappa)
337            })
338            .collect()
339    }
340}
341
342/// `DispBeamColumn`'s spatial (biaxial) counterpart: a 2-node,
343/// displacement-based, fiber-discretized 3D beam-column. Bending in both
344/// planes comes directly from `FiberSection3` the same way `DispBeamColumn`
345/// gets its single bending plane from `FiberSection`; torsion is *not*
346/// fiber-derived (see `FiberSection3`'s doc comment) and is instead a
347/// simple decoupled elastic `G*J/L` term, folded directly into the local
348/// tangent/resistance alongside the fiber-integrated axial/biaxial-bending
349/// terms — the same decoupling `ElasticBeamColumn3` uses for torsion.
350///
351/// `vec_xz` plus `Node3`'s coordinates give the local axes via the same
352/// `GeomTransf3::local_axes`/`rotation_matrix` machinery `ElasticBeamColumn3`
353/// uses, but only ever in its `Linear3` sense — no `PDelta3` option is
354/// exposed, matching planar `DispBeamColumn`'s "`GeomTransf::Linear` only"
355/// scope (see its doc comment: a fiber section's state isn't a single
356/// scalar axial force, so it doesn't plug into the closed-form geometric-
357/// stiffness formula directly, and building the right generalization is
358/// deferred, not silently wrong).
359///
360/// Local DOF order per node `[u, v, w, rx, ry, rz]`, matching
361/// `ElasticBeamColumn3`/`SpatialDof` exactly.
362#[derive(Debug, Clone)]
363pub struct DispBeamColumn3 {
364    pub node_i: Node3Id,
365    pub node_j: Node3Id,
366    pub g: f64,
367    pub j: f64,
368    vec_xz: [f64; 3],
369    integration: BeamIntegration,
370    sections: Vec<FiberSection3>,
371    /// Mass per unit volume, same convention as `DispBeamColumn::density`.
372    pub density: f64,
373}
374
375impl DispBeamColumn3 {
376    pub fn new(
377        node_i: Node3Id,
378        node_j: Node3Id,
379        g: f64,
380        j: f64,
381        vec_xz: [f64; 3],
382        fibers: Vec<Fiber3>,
383        integration: BeamIntegration,
384    ) -> Self {
385        let n_points = integration.points().len();
386        let sections = (0..n_points)
387            .map(|_| FiberSection3::new(fibers.clone()))
388            .collect();
389        DispBeamColumn3 {
390            node_i,
391            node_j,
392            g,
393            j,
394            vec_xz,
395            integration,
396            sections,
397            density: 0.0,
398        }
399    }
400
401    pub fn with_density(mut self, density: f64) -> Self {
402        self.density = density;
403        self
404    }
405
406    /// Strain-displacement vectors `(b_eps0, b_kappa_z, b_kappa_y)` at `xi`
407    /// in `[0,1]`, relating local nodal displacements
408    /// `[u1,v1,w1,rx1,ry1,rz1,u2,v2,w2,rx2,ry2,rz2]` to `[eps0, kappa_z,
409    /// kappa_y]`. `b_eps0`/`b_kappa_z` are `DispBeamColumn::strain_
410    /// displacement`'s `b_eps0`/`b_kappa` re-indexed onto `[u,v,rz]`'s
411    /// spatial slots (`0,1,5` / `6,7,11`); `b_kappa_y` is the same Hermite
412    /// curvature operator re-derived for `[w, ry]`'s slots (`2,4` / `8,10`)
413    /// with its rotation-column signs flipped, from substituting
414    /// `theta_equiv = -ry` (the standard-Hermite formula assumes
415    /// `theta = dv/dx`, but this element's convention is `ry = -dw/dx`, see
416    /// `ElasticBeamColumn3::local_elastic_stiffness`'s doc comment) into the
417    /// unmodified `b_kappa_z` formula — the same substitution
418    /// `ElasticBeamColumn3::geometric_stiffness`'s `y_block` uses relative
419    /// to `z_block`. Verified, not just plausible-by-analogy, in
420    /// `core/tests/elements/disp_beam.rs` against `ElasticBeamColumn3`'s
421    /// exact closed-form biaxial-bending stiffness.
422    fn strain_displacement(
423        xi: f64,
424        length: f64,
425    ) -> (
426        SpatialElementVector,
427        SpatialElementVector,
428        SpatialElementVector,
429    ) {
430        let l = length;
431        let mut b_eps0 = SpatialElementVector::zeros();
432        b_eps0[0] = -1.0 / l;
433        b_eps0[6] = 1.0 / l;
434
435        let mut b_kappa_z = SpatialElementVector::zeros();
436        b_kappa_z[1] = (-6.0 + 12.0 * xi) / (l * l);
437        b_kappa_z[5] = (-4.0 + 6.0 * xi) / l;
438        b_kappa_z[7] = (6.0 - 12.0 * xi) / (l * l);
439        b_kappa_z[11] = (-2.0 + 6.0 * xi) / l;
440
441        let mut b_kappa_y = SpatialElementVector::zeros();
442        b_kappa_y[2] = (-6.0 + 12.0 * xi) / (l * l);
443        b_kappa_y[4] = (4.0 - 6.0 * xi) / l;
444        b_kappa_y[8] = (6.0 - 12.0 * xi) / (l * l);
445        b_kappa_y[10] = (2.0 - 6.0 * xi) / l;
446
447        (b_eps0, b_kappa_z, b_kappa_y)
448    }
449
450    fn local_displacement(
451        &self,
452        node_i: &Node3,
453        node_j: &Node3,
454    ) -> (f64, SpatialElementMatrix, SpatialElementVector) {
455        let transform = GeomTransf3::linear(self.vec_xz);
456        let (length, r) = transform.local_axes(node_i, node_j);
457        let t = GeomTransf3::rotation_matrix(&r);
458        let d_global = SpatialElementVector::from_iterator(
459            node_i
460                .displacement
461                .iter()
462                .chain(node_j.displacement.iter())
463                .copied(),
464        );
465        (length, t, t * d_global)
466    }
467
468    /// Decoupled elastic torsion contribution `[rx1, rx2]` — see the type
469    /// doc comment for why torsion isn't fiber-derived.
470    fn torsion_stiffness(&self, length: f64) -> SpatialElementMatrix {
471        let gj_l = self.g * self.j / length;
472        let mut k = SpatialElementMatrix::zeros();
473        k[(3, 3)] = gj_l;
474        k[(9, 9)] = gj_l;
475        k[(3, 9)] = -gj_l;
476        k[(9, 3)] = -gj_l;
477        k
478    }
479
480    pub(super) fn form_tangent_and_resistance(
481        &self,
482        node_i: &Node3,
483        node_j: &Node3,
484    ) -> (SpatialElementMatrix, SpatialElementVector) {
485        let (length, t, d_local) = self.local_displacement(node_i, node_j);
486
487        let mut k_local = self.torsion_stiffness(length);
488        let mut r_local = k_local * d_local;
489
490        for ((xi, w), section) in self.integration.points().iter().zip(&self.sections) {
491            let (b_eps0, b_kappa_z, b_kappa_y) = Self::strain_displacement(*xi, length);
492            let eps0 = b_eps0.dot(&d_local);
493            let kappa_z = b_kappa_z.dot(&d_local);
494            let kappa_y = b_kappa_y.dot(&d_local);
495            let (n, mz, my, k_section) = section.trial(eps0, kappa_z, kappa_y);
496
497            let scale = w * length;
498            r_local += scale * (b_eps0 * n + b_kappa_z * mz + b_kappa_y * my);
499
500            let b = [&b_eps0, &b_kappa_z, &b_kappa_y];
501            for (bi, row_i) in b.iter().enumerate() {
502                for (bj, row_j) in b.iter().enumerate() {
503                    k_local += scale * k_section[bi][bj] * (*row_i * row_j.transpose());
504                }
505            }
506        }
507
508        (t.transpose() * k_local * t, t.transpose() * r_local)
509    }
510
511    /// Local nodal force (including the decoupled torsion term) — the same
512    /// pre-transform `r_local` value above, recomputed fresh (cheap: no
513    /// tangent accumulation needed, and no Newton loop to avoid re-running).
514    pub(super) fn local_force(&self, node_i: &Node3, node_j: &Node3) -> SpatialElementVector {
515        let (length, _t, d_local) = self.local_displacement(node_i, node_j);
516        let mut r_local = self.torsion_stiffness(length) * d_local;
517        for ((xi, w), section) in self.integration.points().iter().zip(&self.sections) {
518            let (b_eps0, b_kappa_z, b_kappa_y) = Self::strain_displacement(*xi, length);
519            let eps0 = b_eps0.dot(&d_local);
520            let kappa_z = b_kappa_z.dot(&d_local);
521            let kappa_y = b_kappa_y.dot(&d_local);
522            let (n, mz, my, _k_section) = section.trial(eps0, kappa_z, kappa_y);
523            r_local += (w * length) * (b_eps0 * n + b_kappa_z * mz + b_kappa_y * my);
524        }
525        r_local
526    }
527
528    pub(super) fn commit(&mut self, node_i: &Node3, node_j: &Node3) {
529        let (length, _t, d_local) = self.local_displacement(node_i, node_j);
530        let points = self.integration.points();
531        for ((xi, _w), section) in points.iter().zip(&mut self.sections) {
532            let (b_eps0, b_kappa_z, b_kappa_y) = Self::strain_displacement(*xi, length);
533            let eps0 = b_eps0.dot(&d_local);
534            let kappa_z = b_kappa_z.dot(&d_local);
535            let kappa_y = b_kappa_y.dot(&d_local);
536            section.commit(eps0, kappa_z, kappa_y);
537        }
538    }
539
540    /// Spatial counterpart of `DispBeamColumn::form_load_vector`: equivalent
541    /// nodal load (global coordinates) from a uniform local-axis load —
542    /// axial `wx` and transverse `wy`/`wz`.
543    pub(super) fn form_load_vector(
544        &self,
545        node_i: &Node3,
546        node_j: &Node3,
547        wx: f64,
548        wy: f64,
549        wz: f64,
550    ) -> SpatialElementVector {
551        if wx == 0.0 && wy == 0.0 && wz == 0.0 {
552            return SpatialElementVector::zeros();
553        }
554        let (length, r) = GeomTransf3::linear(self.vec_xz).local_axes(node_i, node_j);
555        let t = GeomTransf3::rotation_matrix(&r);
556        t.transpose() * uniform_load::spatial_local(length, wx, wy, wz)
557    }
558
559    /// `form_load_vector` in the element-local frame `local_force` reports in.
560    pub(super) fn form_local_load_vector(
561        &self,
562        node_i: &Node3,
563        node_j: &Node3,
564        wx: f64,
565        wy: f64,
566        wz: f64,
567    ) -> SpatialElementVector {
568        let (length, _r) = GeomTransf3::linear(self.vec_xz).local_axes(node_i, node_j);
569        uniform_load::spatial_local(length, wx, wy, wz)
570    }
571
572    pub(super) fn form_mass(&self, node_i: &Node3, node_j: &Node3) -> SpatialElementVector {
573        let transform = GeomTransf3::linear(self.vec_xz);
574        let (length, _r) = transform.local_axes(node_i, node_j);
575        let total_area = self.sections[0].total_area();
576        let half = self.density * total_area * length / 2.0;
577        let mut mass = SpatialElementVector::zeros();
578        for dof in [0, 1, 2, 6, 7, 8] {
579            mass[dof] = half;
580        }
581        mass
582    }
583
584    /// See `DispBeamColumn::fiber_responses`'s doc comment — same
585    /// contract, biaxial strain field.
586    pub fn fiber_responses(&self, node_i: &Node3, node_j: &Node3) -> Vec<Vec<(f64, f64)>> {
587        let (length, _t, d_local) = self.local_displacement(node_i, node_j);
588        self.integration
589            .points()
590            .iter()
591            .zip(&self.sections)
592            .map(|((xi, _w), section)| {
593                let (b_eps0, b_kappa_z, b_kappa_y) = Self::strain_displacement(*xi, length);
594                let eps0 = b_eps0.dot(&d_local);
595                let kappa_z = b_kappa_z.dot(&d_local);
596                let kappa_y = b_kappa_y.dot(&d_local);
597                section.fiber_responses(eps0, kappa_z, kappa_y)
598            })
599            .collect()
600    }
601}
carapace_core/model/elements/disp_beam_column.rs

carapace_core/model/elements/
disp_beam_column.rs
