carapace_core/model/elements/
mod.rs

1use nalgebra::{SMatrix, SVector};
2use slotmap::new_key_type;
3
4use super::{
5    ElementLoad, ElementLoad3, Node, Node3, Node3Id, NodeId, ELEMENT_DOF, NDF, PLANAR_NDIM,
6    SPATIAL_ELEMENT_DOF, SPATIAL_NDF, SPATIAL_NDIM,
7};
8
9mod ops;
10pub use ops::{
11    two_node_dofs, DofMask, DofRef, ElementForce, ElementOps, GaussResponse, NodeList, NodeView,
12    ShellResponse, TangentSink, VectorSink, MAX_ELEMENT_NODES,
13};
14
15mod disp_beam_column;
16mod elastic_beam_column;
17mod force_beam_column;
18mod plane_common;
19mod quad4;
20mod shell3;
21mod shell4;
22mod tri3;
23mod truss;
24mod uniform_load;
25mod zero_length;
26
27pub use disp_beam_column::{DispBeamColumn, DispBeamColumn3};
28pub use elastic_beam_column::{ElasticBeamColumn, ElasticBeamColumn3};
29pub use force_beam_column::{ForceBeamColumn, ForceBeamColumn3};
30pub use quad4::{Quad4, Quad4Formulation};
31pub use shell3::Shell3;
32pub use shell4::Shell4;
33pub use tri3::Tri3;
34pub use truss::{SpatialElementMatrix, SpatialElementVector, Truss, Truss3};
35pub use zero_length::{
36    Friction, Friction3, Orientation, OrientationError, ZeroLength, ZeroLength3, ZeroLengthSection,
37    ZeroLengthSection3,
38};
39
40new_key_type! {
41    /// Generational index into `Domain`'s element store.
42    pub struct ElementId;
43
44    /// Generational index into `Domain3`'s element store. Distinct from
45    /// `ElementId` for the same reason `Node3Id` is distinct from `NodeId`
46    /// (see its doc comment).
47    pub struct Element3Id;
48}
49
50/// Element catalog. Closed enum, `match`-based dispatch, no `Box<dyn Trait>`
51///.
52#[derive(Debug, Clone)]
53pub enum Element {
54    Truss(Truss),
55    ZeroLength(ZeroLength),
56    ZeroLengthSection(ZeroLengthSection),
57    ElasticBeamColumn(ElasticBeamColumn),
58    DispBeamColumn(DispBeamColumn),
59    ForceBeamColumn(ForceBeamColumn),
60    Tri3(Tri3),
61    Quad4(Quad4),
62}
63
64const CONTINUUM: &str = "continuum elements assemble through ElementOps, not the two-node path";
65
66impl Element {
67    /// The nodes this element connects, in local order.
68    pub fn nodes(&self) -> NodeList<NodeId> {
69        match self {
70            Element::Tri3(t) => t.nodes.into_iter().collect(),
71            Element::Quad4(q) => q.nodes.into_iter().collect(),
72            _ => self.pair().into_iter().collect(),
73        }
74    }
75
76    /// The two end nodes of a two-node element.
77    fn pair(&self) -> [NodeId; 2] {
78        match self {
79            Element::Tri3(_) | Element::Quad4(_) => unreachable!("{CONTINUUM}"),
80            Element::Truss(t) => [t.node_i, t.node_j],
81            Element::ZeroLength(z) => [z.node_i, z.node_j],
82            Element::ZeroLengthSection(z) => [z.node_i, z.node_j],
83            Element::ElasticBeamColumn(b) => [b.node_i, b.node_j],
84            Element::DispBeamColumn(b) => [b.node_i, b.node_j],
85            Element::ForceBeamColumn(b) => [b.node_i, b.node_j],
86        }
87    }
88
89    /// Form this element's contribution to the global tangent stiffness and
90    /// internal resisting force, given its two nodes' current state.
91    /// Returned in element-local DOF order
92    /// `[ux_i, uy_i, rz_i, ux_j, uy_j, rz_j]`; the caller
93    /// (`Domain::form_tangent_and_residual`) scatters these into the global
94    /// system using each node's equation numbers.
95    fn form_tangent_and_resistance(
96        &self,
97        node_i: &Node,
98        node_j: &Node,
99        load: Option<&ElementLoad>,
100    ) -> (
101        SMatrix<f64, ELEMENT_DOF, ELEMENT_DOF>,
102        SVector<f64, ELEMENT_DOF>,
103    ) {
104        match self {
105            Element::Truss(t) => t.form_tangent_and_resistance(node_i, node_j),
106            Element::ZeroLength(z) => z.form_tangent_and_resistance(node_i, node_j),
107            Element::ZeroLengthSection(z) => z.form_tangent_and_resistance(node_i, node_j),
108            Element::ElasticBeamColumn(b) => b.form_tangent_and_resistance(node_i, node_j),
109            Element::DispBeamColumn(b) => b.form_tangent_and_resistance(node_i, node_j),
110            Element::ForceBeamColumn(b) => b.form_tangent_and_resistance(node_i, node_j, load),
111            Element::Tri3(_) | Element::Quad4(_) => unreachable!("{CONTINUUM}"),
112        }
113    }
114
115    /// This element's equivalent nodal load vector (global coordinates,
116    /// same DOF order as above) from `load` — the `ElementLoad` (if any)
117    /// that whichever `LoadPattern` is currently being assembled has on
118    /// this element — e.g. a beam-column's distributed transverse
119    /// load. Zero when `load` is `None`, and for elements with no
120    /// element-load support at all (`Truss`, `ZeroLength`).
121    fn form_load_vector(
122        &self,
123        node_i: &Node,
124        node_j: &Node,
125        load: Option<&ElementLoad>,
126    ) -> SVector<f64, ELEMENT_DOF> {
127        match (self, load) {
128            (Element::ElasticBeamColumn(b), Some(load)) => {
129                b.form_load_vector(node_i, node_j, load.beam_uniform[0], load.beam_uniform[1])
130            }
131            (Element::DispBeamColumn(b), Some(load)) => {
132                b.form_load_vector(node_i, node_j, load.beam_uniform[0], load.beam_uniform[1])
133            }
134            (Element::ForceBeamColumn(b), Some(load)) => {
135                b.form_load_vector(node_i, node_j, load.beam_uniform[0], load.beam_uniform[1])
136            }
137            _ => SVector::<f64, ELEMENT_DOF>::zeros(),
138        }
139    }
140
141    /// See `ElementOps::form_local_load_vector`'s doc comment.
142    fn form_local_load_vector(
143        &self,
144        node_i: &Node,
145        node_j: &Node,
146        load: Option<&ElementLoad>,
147    ) -> SVector<f64, ELEMENT_DOF> {
148        match (self, load) {
149            (Element::ElasticBeamColumn(b), Some(load)) => {
150                b.form_local_load_vector(node_i, node_j, load.beam_uniform[0], load.beam_uniform[1])
151            }
152            (Element::DispBeamColumn(b), Some(load)) => {
153                b.form_local_load_vector(node_i, node_j, load.beam_uniform[0], load.beam_uniform[1])
154            }
155            (Element::ForceBeamColumn(b), Some(load)) => {
156                b.form_local_load_vector(node_i, node_j, load.beam_uniform[0], load.beam_uniform[1])
157            }
158            _ => SVector::<f64, ELEMENT_DOF>::zeros(),
159        }
160    }
161
162    /// This element's lumped-mass contribution (diagonal only, lumped) in the same local DOF
163    /// order, geometry-dependent (`length`) so it needs both nodes. Zero for `ZeroLength` (a spring/
164    /// connector, not a mass-bearing member) and for any element with the
165    /// default `density = 0.0`.
166    fn form_mass(&self, node_i: &Node, node_j: &Node) -> SVector<f64, ELEMENT_DOF> {
167        match self {
168            Element::Truss(t) => t.form_mass(node_i, node_j),
169            Element::ZeroLength(_) => SVector::<f64, ELEMENT_DOF>::zeros(),
170            Element::ZeroLengthSection(_) => SVector::<f64, ELEMENT_DOF>::zeros(),
171            Element::ElasticBeamColumn(b) => b.form_mass(node_i, node_j),
172            Element::DispBeamColumn(b) => b.form_mass(node_i, node_j),
173            Element::ForceBeamColumn(b) => b.form_mass(node_i, node_j),
174            Element::Tri3(_) | Element::Quad4(_) => unreachable!("{CONTINUUM}"),
175        }
176    }
177
178    /// Commit this element's material(s) at the given (final, converged)
179    /// node state — see `Material`'s doc comment for why this is the only
180    /// place a `Material` ever mutates. A no-op for `ElasticBeamColumn`
181    /// (no `Material` — its response is closed-form) and for any
182    /// `ZeroLength` direction with no material assigned.
183    fn commit(&mut self, node_i: &Node, node_j: &Node, load: Option<&ElementLoad>) {
184        match self {
185            Element::Truss(t) => t.commit(node_i, node_j),
186            Element::ZeroLength(z) => z.commit(node_i, node_j),
187            Element::ZeroLengthSection(z) => z.commit(node_i, node_j),
188            Element::ElasticBeamColumn(_) => {}
189            Element::DispBeamColumn(b) => b.commit(node_i, node_j),
190            Element::ForceBeamColumn(b) => b.commit(node_i, node_j, load),
191            Element::Tri3(_) | Element::Quad4(_) => unreachable!("{CONTINUUM}"),
192        }
193    }
194
195    /// See `ElementOps::local_force`'s doc comment. `ForceBeamColumn`'s
196    /// local force is cached at `commit` time (its own `local_force()`
197    /// takes no node arguments — reading the cache, not recomputing);
198    /// every other variant recomputes fresh from `node_i`/`node_j`.
199    fn local_force(&self, node_i: &Node, node_j: &Node) -> SVector<f64, ELEMENT_DOF> {
200        match self {
201            Element::Truss(t) => t.local_force(node_i, node_j),
202            Element::ZeroLength(z) => z.local_force(node_i, node_j),
203            Element::ZeroLengthSection(z) => z.local_force(node_i, node_j),
204            Element::ElasticBeamColumn(b) => b.local_force(node_i, node_j),
205            Element::DispBeamColumn(b) => b.local_force(node_i, node_j),
206            Element::ForceBeamColumn(b) => b.local_force(),
207            Element::Tri3(_) | Element::Quad4(_) => unreachable!("{CONTINUUM}"),
208        }
209    }
210
211    /// See `ElementOps::fiber_responses`'s doc comment.
212    fn fiber_responses(&self, node_i: &Node, node_j: &Node) -> Option<Vec<Vec<(f64, f64)>>> {
213        match self {
214            Element::DispBeamColumn(b) => Some(b.fiber_responses(node_i, node_j)),
215            Element::ForceBeamColumn(b) => Some(b.fiber_responses()),
216            Element::Truss(_)
217            | Element::ZeroLength(_)
218            | Element::ZeroLengthSection(_)
219            | Element::ElasticBeamColumn(_)
220            | Element::Tri3(_)
221            | Element::Quad4(_) => None,
222        }
223    }
224}
225
226impl ElementOps<PLANAR_NDIM, NDF, NodeId> for Element {
227    type Load = ElementLoad;
228    type Id = ElementId;
229
230    fn nodes(&self) -> NodeList<NodeId> {
231        Element::nodes(self)
232    }
233
234    fn dof_mask(&self) -> DofMask {
235        match self {
236            // A truss carries no moment: translations only.
237            Element::Truss(_) => DofMask::none().with(0).with(1),
238            Element::ZeroLength(z) => z.dof_mask(),
239            Element::ZeroLengthSection(z) => z.dof_mask(),
240            Element::ElasticBeamColumn(_)
241            | Element::DispBeamColumn(_)
242            | Element::ForceBeamColumn(_) => DofMask::all(NDF),
243            Element::Tri3(t) => t.dof_mask(),
244            Element::Quad4(q) => q.dof_mask(),
245        }
246    }
247
248    fn validate(&self, nodes: &NodeView<'_, PLANAR_NDIM, NDF, NodeId>) -> Result<(), &'static str> {
249        match self {
250            Element::Tri3(t) => t.validate(nodes),
251            Element::Quad4(q) => q.validate(nodes),
252            _ => Ok(()),
253        }
254    }
255
256    fn accepts_load(&self, load: &ElementLoad) -> bool {
257        match self {
258            Element::ElasticBeamColumn(_)
259            | Element::DispBeamColumn(_)
260            | Element::ForceBeamColumn(_) => !load.has_body() && !load.has_edge(),
261            Element::Tri3(_) => !load.has_beam() && !load.has_edge_from(3),
262            Element::Quad4(_) => !load.has_beam(),
263            Element::Truss(_) | Element::ZeroLength(_) | Element::ZeroLengthSection(_) => {
264                !load.has_beam() && !load.has_body() && !load.has_edge()
265            }
266        }
267    }
268
269    fn prepare(&mut self, nodes: &NodeView<'_, PLANAR_NDIM, NDF, NodeId>) {
270        match self {
271            Element::Tri3(t) => t.prepare(nodes),
272            Element::Quad4(q) => q.prepare(nodes),
273            _ => {}
274        }
275    }
276
277    fn gauss_point_count(&self) -> usize {
278        match self {
279            Element::Tri3(_) => 1,
280            Element::Quad4(_) => 4,
281            _ => 0,
282        }
283    }
284
285    fn gauss_responses(
286        &self,
287        nodes: &NodeView<'_, PLANAR_NDIM, NDF, NodeId>,
288    ) -> Option<Vec<GaussResponse>> {
289        match self {
290            Element::Tri3(t) => Some(t.gauss_responses(nodes)),
291            Element::Quad4(q) => Some(q.gauss_responses(nodes)),
292            _ => None,
293        }
294    }
295
296    fn assemble_tangent<S: TangentSink<NodeId>>(
297        &self,
298        nodes: &NodeView<'_, PLANAR_NDIM, NDF, NodeId>,
299        load: Option<&ElementLoad>,
300        sink: &mut S,
301    ) {
302        match self {
303            Element::Tri3(t) => return t.assemble_tangent(nodes, sink),
304            Element::Quad4(q) => return q.assemble_tangent(nodes, sink),
305            _ => {}
306        }
307        let [i, j] = Element::pair(self);
308        let (k, r) = Element::form_tangent_and_resistance(self, nodes.get(i), nodes.get(j), load);
309        sink.add(&two_node_dofs::<_, ELEMENT_DOF>(i, j, NDF), &k, &r);
310    }
311
312    fn assemble_load<S: VectorSink<NodeId>>(
313        &self,
314        nodes: &NodeView<'_, PLANAR_NDIM, NDF, NodeId>,
315        load: Option<&ElementLoad>,
316        sink: &mut S,
317    ) {
318        match (self, load) {
319            (Element::Tri3(t), load) => {
320                return load.map_or((), |load| t.assemble_load(nodes, load, sink))
321            }
322            (Element::Quad4(q), load) => {
323                return load.map_or((), |load| q.assemble_load(nodes, load, sink))
324            }
325            _ => {}
326        }
327        let [i, j] = Element::pair(self);
328        let v = Element::form_load_vector(self, nodes.get(i), nodes.get(j), load);
329        sink.add(&two_node_dofs::<_, ELEMENT_DOF>(i, j, NDF), &v);
330    }
331
332    fn assemble_mass<S: VectorSink<NodeId>>(
333        &self,
334        nodes: &NodeView<'_, PLANAR_NDIM, NDF, NodeId>,
335        sink: &mut S,
336    ) {
337        match self {
338            Element::Tri3(t) => return t.assemble_mass(nodes, sink),
339            Element::Quad4(q) => return q.assemble_mass(nodes, sink),
340            _ => {}
341        }
342        let [i, j] = Element::pair(self);
343        let v = Element::form_mass(self, nodes.get(i), nodes.get(j));
344        sink.add(&two_node_dofs::<_, ELEMENT_DOF>(i, j, NDF), &v);
345    }
346
347    fn commit(
348        &mut self,
349        nodes: &NodeView<'_, PLANAR_NDIM, NDF, NodeId>,
350        load: Option<&ElementLoad>,
351    ) {
352        match self {
353            Element::Tri3(t) => return t.commit(nodes),
354            Element::Quad4(q) => return q.commit(nodes),
355            _ => {}
356        }
357        let [i, j] = Element::pair(self);
358        Element::commit(self, nodes.get(i), nodes.get(j), load)
359    }
360
361    fn local_force_width(&self) -> usize {
362        match self {
363            Element::Tri3(_) => 6,
364            Element::Quad4(_) => 8,
365            _ => ELEMENT_DOF,
366        }
367    }
368
369    fn local_force(&self, nodes: &NodeView<'_, PLANAR_NDIM, NDF, NodeId>) -> ElementForce {
370        match self {
371            Element::Tri3(t) => return t.local_force(nodes),
372            Element::Quad4(q) => return q.local_force(nodes),
373            _ => {}
374        }
375        let [i, j] = Element::pair(self);
376        Element::local_force(self, nodes.get(i), nodes.get(j)).into()
377    }
378
379    fn local_load_force(
380        &self,
381        nodes: &NodeView<'_, PLANAR_NDIM, NDF, NodeId>,
382        load: Option<&ElementLoad>,
383    ) -> ElementForce {
384        match self {
385            Element::Tri3(t) => {
386                return load
387                    .map_or(SVector::<f64, 6>::zeros(), |l| t.load_vector(nodes, l))
388                    .into()
389            }
390            Element::Quad4(q) => {
391                return load
392                    .map_or(SVector::<f64, 8>::zeros(), |l| q.load_vector(nodes, l))
393                    .into()
394            }
395            _ => {}
396        }
397        let [i, j] = Element::pair(self);
398        Element::form_local_load_vector(self, nodes.get(i), nodes.get(j), load).into()
399    }
400
401    fn fiber_responses(
402        &self,
403        nodes: &NodeView<'_, PLANAR_NDIM, NDF, NodeId>,
404    ) -> Option<Vec<Vec<(f64, f64)>>> {
405        if let Element::Tri3(_) | Element::Quad4(_) = self {
406            return None;
407        }
408        let [i, j] = Element::pair(self);
409        Element::fiber_responses(self, nodes.get(i), nodes.get(j))
410    }
411}
412/// 3D element catalog — `Domain3`'s counterpart to `Element`. Closed
413/// enum, `match`-based dispatch, same reasoning as `Element`.
414#[derive(Debug, Clone)]
415pub enum Element3 {
416    Truss3(Truss3),
417    ZeroLength3(ZeroLength3),
418    ZeroLengthSection3(ZeroLengthSection3),
419    ElasticBeamColumn3(ElasticBeamColumn3),
420    DispBeamColumn3(DispBeamColumn3),
421    ForceBeamColumn3(ForceBeamColumn3),
422    Shell3(Shell3),
423    Shell4(Shell4),
424}
425
426const SHELL: &str = "shell elements assemble through ElementOps, not the two-node path";
427
428impl Element3 {
429    /// The nodes this element connects, in local order.
430    pub fn nodes(&self) -> NodeList<Node3Id> {
431        match self {
432            Element3::Shell3(s) => s.nodes.into_iter().collect(),
433            Element3::Shell4(s) => s.nodes.into_iter().collect(),
434            _ => self.pair().into_iter().collect(),
435        }
436    }
437
438    /// The two end nodes of a two-node element.
439    fn pair(&self) -> [Node3Id; 2] {
440        match self {
441            Element3::Shell3(_) | Element3::Shell4(_) => unreachable!("{SHELL}"),
442            Element3::Truss3(t) => [t.node_i, t.node_j],
443            Element3::ZeroLength3(z) => [z.node_i, z.node_j],
444            Element3::ZeroLengthSection3(z) => [z.node_i, z.node_j],
445            Element3::ElasticBeamColumn3(b) => [b.node_i, b.node_j],
446            Element3::DispBeamColumn3(b) => [b.node_i, b.node_j],
447            Element3::ForceBeamColumn3(b) => [b.node_i, b.node_j],
448        }
449    }
450
451    fn form_tangent_and_resistance(
452        &self,
453        node_i: &Node3,
454        node_j: &Node3,
455        load: Option<&ElementLoad3>,
456    ) -> (
457        SMatrix<f64, SPATIAL_ELEMENT_DOF, SPATIAL_ELEMENT_DOF>,
458        SVector<f64, SPATIAL_ELEMENT_DOF>,
459    ) {
460        match self {
461            Element3::Truss3(t) => t.form_tangent_and_resistance(node_i, node_j),
462            Element3::ZeroLength3(z) => z.form_tangent_and_resistance(node_i, node_j),
463            Element3::ZeroLengthSection3(z) => z.form_tangent_and_resistance(node_i, node_j),
464            Element3::ElasticBeamColumn3(b) => b.form_tangent_and_resistance(node_i, node_j),
465            Element3::DispBeamColumn3(b) => b.form_tangent_and_resistance(node_i, node_j),
466            Element3::ForceBeamColumn3(b) => b.form_tangent_and_resistance(node_i, node_j, load),
467            Element3::Shell3(_) | Element3::Shell4(_) => unreachable!("{SHELL}"),
468        }
469    }
470
471    /// This element's equivalent nodal load vector from `load` — the
472    /// 3D counterpart of `Element::form_load_vector`. Zero when
473    /// `load` is `None`, and for elements with no element-load support
474    /// (`Truss3`, `ZeroLength3`).
475    fn form_load_vector(
476        &self,
477        node_i: &Node3,
478        node_j: &Node3,
479        load: Option<&ElementLoad3>,
480    ) -> SVector<f64, SPATIAL_ELEMENT_DOF> {
481        match (self, load) {
482            (Element3::ElasticBeamColumn3(b), Some(l)) => {
483                let [wx, wy, wz] = l.uniform;
484                b.form_load_vector(node_i, node_j, wx, wy, wz)
485            }
486            (Element3::DispBeamColumn3(b), Some(l)) => {
487                let [wx, wy, wz] = l.uniform;
488                b.form_load_vector(node_i, node_j, wx, wy, wz)
489            }
490            (Element3::ForceBeamColumn3(b), Some(l)) => {
491                let [wx, wy, wz] = l.uniform;
492                b.form_load_vector(node_i, node_j, wx, wy, wz)
493            }
494            _ => SVector::<f64, SPATIAL_ELEMENT_DOF>::zeros(),
495        }
496    }
497
498    /// See `ElementOps::form_local_load_vector`'s doc comment.
499    fn form_local_load_vector(
500        &self,
501        node_i: &Node3,
502        node_j: &Node3,
503        load: Option<&ElementLoad3>,
504    ) -> SVector<f64, SPATIAL_ELEMENT_DOF> {
505        match (self, load) {
506            (Element3::ElasticBeamColumn3(b), Some(l)) => {
507                let [wx, wy, wz] = l.uniform;
508                b.form_local_load_vector(node_i, node_j, wx, wy, wz)
509            }
510            (Element3::DispBeamColumn3(b), Some(l)) => {
511                let [wx, wy, wz] = l.uniform;
512                b.form_local_load_vector(node_i, node_j, wx, wy, wz)
513            }
514            (Element3::ForceBeamColumn3(b), Some(l)) => {
515                let [wx, wy, wz] = l.uniform;
516                b.form_local_load_vector(node_i, node_j, wx, wy, wz)
517            }
518            _ => SVector::<f64, SPATIAL_ELEMENT_DOF>::zeros(),
519        }
520    }
521
522    /// This element's lumped-mass contribution — see `Element::form_mass`.
523    fn form_mass(&self, node_i: &Node3, node_j: &Node3) -> SVector<f64, SPATIAL_ELEMENT_DOF> {
524        match self {
525            Element3::Truss3(t) => t.form_mass(node_i, node_j),
526            Element3::ZeroLength3(_) => SVector::<f64, SPATIAL_ELEMENT_DOF>::zeros(),
527            Element3::ZeroLengthSection3(_) => SVector::<f64, SPATIAL_ELEMENT_DOF>::zeros(),
528            Element3::ElasticBeamColumn3(b) => b.form_mass(node_i, node_j),
529            Element3::DispBeamColumn3(b) => b.form_mass(node_i, node_j),
530            Element3::ForceBeamColumn3(b) => b.form_mass(node_i, node_j),
531            Element3::Shell3(_) | Element3::Shell4(_) => unreachable!("{SHELL}"),
532        }
533    }
534
535    fn commit(&mut self, node_i: &Node3, node_j: &Node3, load: Option<&ElementLoad3>) {
536        match self {
537            Element3::Truss3(t) => t.commit(node_i, node_j),
538            Element3::ZeroLength3(z) => z.commit(node_i, node_j),
539            Element3::ZeroLengthSection3(z) => z.commit(node_i, node_j),
540            Element3::ElasticBeamColumn3(_) => {}
541            Element3::DispBeamColumn3(b) => b.commit(node_i, node_j),
542            Element3::ForceBeamColumn3(b) => b.commit(node_i, node_j, load),
543            Element3::Shell3(_) | Element3::Shell4(_) => unreachable!("{SHELL}"),
544        }
545    }
546
547    /// See `Element::local_force`'s doc comment.
548    fn local_force(&self, node_i: &Node3, node_j: &Node3) -> SVector<f64, SPATIAL_ELEMENT_DOF> {
549        match self {
550            Element3::Truss3(t) => t.local_force(node_i, node_j),
551            Element3::ZeroLength3(z) => z.local_force(node_i, node_j),
552            Element3::ZeroLengthSection3(z) => z.local_force(node_i, node_j),
553            Element3::ElasticBeamColumn3(b) => b.local_force(node_i, node_j),
554            Element3::DispBeamColumn3(b) => b.local_force(node_i, node_j),
555            Element3::ForceBeamColumn3(b) => b.local_force(),
556            Element3::Shell3(_) | Element3::Shell4(_) => unreachable!("{SHELL}"),
557        }
558    }
559
560    /// See `ElementOps::fiber_responses`'s doc comment.
561    fn fiber_responses(&self, node_i: &Node3, node_j: &Node3) -> Option<Vec<Vec<(f64, f64)>>> {
562        match self {
563            Element3::DispBeamColumn3(b) => Some(b.fiber_responses(node_i, node_j)),
564            Element3::ForceBeamColumn3(b) => Some(b.fiber_responses()),
565            Element3::Truss3(_)
566            | Element3::ZeroLength3(_)
567            | Element3::ZeroLengthSection3(_)
568            | Element3::ElasticBeamColumn3(_)
569            | Element3::Shell3(_)
570            | Element3::Shell4(_) => None,
571        }
572    }
573}
574
575impl ElementOps<SPATIAL_NDIM, SPATIAL_NDF, Node3Id> for Element3 {
576    type Load = ElementLoad3;
577    type Id = Element3Id;
578
579    fn nodes(&self) -> NodeList<Node3Id> {
580        Element3::nodes(self)
581    }
582
583    fn dof_mask(&self) -> DofMask {
584        match self {
585            Element3::Truss3(_) => DofMask::none().with(0).with(1).with(2),
586            Element3::ZeroLength3(z) => z.dof_mask(),
587            Element3::ZeroLengthSection3(z) => z.dof_mask(),
588            Element3::ElasticBeamColumn3(_)
589            | Element3::DispBeamColumn3(_)
590            | Element3::ForceBeamColumn3(_) => DofMask::all(SPATIAL_NDF),
591            Element3::Shell3(s) => s.dof_mask(),
592            Element3::Shell4(s) => s.dof_mask(),
593        }
594    }
595
596    fn validate(
597        &self,
598        nodes: &NodeView<'_, SPATIAL_NDIM, SPATIAL_NDF, Node3Id>,
599    ) -> Result<(), &'static str> {
600        match self {
601            Element3::Shell3(s) => s.validate(nodes),
602            Element3::Shell4(s) => s.validate(nodes),
603            _ => Ok(()),
604        }
605    }
606
607    fn accepts_load(&self, load: &ElementLoad3) -> bool {
608        match self {
609            Element3::Shell3(_) | Element3::Shell4(_) => !load.has_uniform(),
610            _ => !load.has_body() && !load.has_pressure(),
611        }
612    }
613
614    fn prepare(&mut self, nodes: &NodeView<'_, SPATIAL_NDIM, SPATIAL_NDF, Node3Id>) {
615        match self {
616            Element3::Shell3(s) => s.prepare(nodes),
617            Element3::Shell4(s) => s.prepare(nodes),
618            _ => {}
619        }
620    }
621
622    fn gauss_point_count(&self) -> usize {
623        match self {
624            Element3::Shell3(_) | Element3::Shell4(_) => 4,
625            _ => 0,
626        }
627    }
628
629    fn gauss_component_count(&self) -> usize {
630        match self {
631            Element3::Shell3(_) | Element3::Shell4(_) => 8,
632            _ => 3,
633        }
634    }
635
636    fn shell_responses(
637        &self,
638        nodes: &NodeView<'_, SPATIAL_NDIM, SPATIAL_NDF, Node3Id>,
639    ) -> Option<Vec<ShellResponse>> {
640        match self {
641            Element3::Shell3(s) => Some(s.shell_responses(nodes)),
642            Element3::Shell4(s) => Some(s.shell_responses(nodes)),
643            _ => None,
644        }
645    }
646
647    fn assemble_tangent<S: TangentSink<Node3Id>>(
648        &self,
649        nodes: &NodeView<'_, SPATIAL_NDIM, SPATIAL_NDF, Node3Id>,
650        load: Option<&ElementLoad3>,
651        sink: &mut S,
652    ) {
653        match self {
654            Element3::Shell3(s) => return s.assemble_tangent(nodes, sink),
655            Element3::Shell4(s) => return s.assemble_tangent(nodes, sink),
656            _ => {}
657        }
658        let [i, j] = Element3::pair(self);
659        let (k, r) = Element3::form_tangent_and_resistance(self, nodes.get(i), nodes.get(j), load);
660        sink.add(
661            &two_node_dofs::<_, SPATIAL_ELEMENT_DOF>(i, j, SPATIAL_NDF),
662            &k,
663            &r,
664        );
665    }
666
667    fn assemble_load<S: VectorSink<Node3Id>>(
668        &self,
669        nodes: &NodeView<'_, SPATIAL_NDIM, SPATIAL_NDF, Node3Id>,
670        load: Option<&ElementLoad3>,
671        sink: &mut S,
672    ) {
673        match self {
674            Element3::Shell3(s) => {
675                return load.map_or((), |load| s.assemble_load(nodes, load, sink))
676            }
677            Element3::Shell4(s) => {
678                return load.map_or((), |load| s.assemble_load(nodes, load, sink))
679            }
680            _ => {}
681        }
682        let [i, j] = Element3::pair(self);
683        let v = Element3::form_load_vector(self, nodes.get(i), nodes.get(j), load);
684        sink.add(
685            &two_node_dofs::<_, SPATIAL_ELEMENT_DOF>(i, j, SPATIAL_NDF),
686            &v,
687        );
688    }
689
690    fn assemble_mass<S: VectorSink<Node3Id>>(
691        &self,
692        nodes: &NodeView<'_, SPATIAL_NDIM, SPATIAL_NDF, Node3Id>,
693        sink: &mut S,
694    ) {
695        match self {
696            Element3::Shell3(s) => return s.assemble_mass(nodes, sink),
697            Element3::Shell4(s) => return s.assemble_mass(nodes, sink),
698            _ => {}
699        }
700        let [i, j] = Element3::pair(self);
701        let v = Element3::form_mass(self, nodes.get(i), nodes.get(j));
702        sink.add(
703            &two_node_dofs::<_, SPATIAL_ELEMENT_DOF>(i, j, SPATIAL_NDF),
704            &v,
705        );
706    }
707
708    fn commit(
709        &mut self,
710        nodes: &NodeView<'_, SPATIAL_NDIM, SPATIAL_NDF, Node3Id>,
711        load: Option<&ElementLoad3>,
712    ) {
713        match self {
714            Element3::Shell3(s) => return s.commit(nodes),
715            Element3::Shell4(s) => return s.commit(nodes),
716            _ => {}
717        }
718        let [i, j] = Element3::pair(self);
719        Element3::commit(self, nodes.get(i), nodes.get(j), load)
720    }
721
722    fn local_force_width(&self) -> usize {
723        match self {
724            Element3::Shell3(_) => 18,
725            Element3::Shell4(_) => 24,
726            _ => SPATIAL_ELEMENT_DOF,
727        }
728    }
729
730    fn local_force(
731        &self,
732        nodes: &NodeView<'_, SPATIAL_NDIM, SPATIAL_NDF, Node3Id>,
733    ) -> ElementForce {
734        match self {
735            Element3::Shell3(s) => return s.local_force(nodes),
736            Element3::Shell4(s) => return s.local_force(nodes),
737            _ => {}
738        }
739        let [i, j] = Element3::pair(self);
740        Element3::local_force(self, nodes.get(i), nodes.get(j)).into()
741    }
742
743    fn local_load_force(
744        &self,
745        nodes: &NodeView<'_, SPATIAL_NDIM, SPATIAL_NDF, Node3Id>,
746        load: Option<&ElementLoad3>,
747    ) -> ElementForce {
748        match self {
749            Element3::Shell3(s) => {
750                return load
751                    .map_or(SVector::<f64, 18>::zeros(), |l| s.load_vector(nodes, l))
752                    .into()
753            }
754            Element3::Shell4(s) => {
755                return load
756                    .map_or(SVector::<f64, 24>::zeros(), |l| s.load_vector(nodes, l))
757                    .into()
758            }
759            _ => {}
760        }
761        let [i, j] = Element3::pair(self);
762        Element3::form_local_load_vector(self, nodes.get(i), nodes.get(j), load).into()
763    }
764
765    fn fiber_responses(
766        &self,
767        nodes: &NodeView<'_, SPATIAL_NDIM, SPATIAL_NDF, Node3Id>,
768    ) -> Option<Vec<Vec<(f64, f64)>>> {
769        if let Element3::Shell3(_) | Element3::Shell4(_) = self {
770            return None;
771        }
772        let [i, j] = Element3::pair(self);
773        Element3::fiber_responses(self, nodes.get(i), nodes.get(j))
774    }
775}
776
777/// `ElementOps::dof_mask` must cover every node DOF slot on which an element
778/// produces stiffness or resistance: a slot outside the mask that nothing
779/// else uses is not an equation, so the contribution would be silently
780/// dropped. This drives every element kind at a generic nonzero state and
781/// compares the slots it actually stiffens with its declared mask. Any new
782/// element kind needs an entry here.
783#[cfg(test)]
784mod mask_conformance {
785    use std::collections::BTreeSet;
786
787    use slotmap::{Key, SlotMap};
788
789    use super::*;
790    use crate::model::{
791        BeamIntegration, ElasticBeamColumn, ElasticBeamColumn3, Fiber, Fiber3, FiberSection,
792        FiberSection3, ForceBeamColumn, ForceBeamColumn3, Friction, Friction3, GeomTransf,
793        GeomTransf3, Material, Orientation, Truss, Truss3, ZeroLength, ZeroLength3,
794        ZeroLengthSection, ZeroLengthSection3,
795    };
796
797    /// Records which `(node, slot)` an element's tangent or resistance is
798    /// nonzero on, relative to the element's own largest entry.
799    struct Touched<NId> {
800        slots: BTreeSet<u8>,
801        _marker: std::marker::PhantomData<NId>,
802    }
803
804    impl<NId: Key> TangentSink<NId> for Touched<NId> {
805        fn add<const N: usize>(
806            &mut self,
807            dofs: &[DofRef<NId>; N],
808            k: &SMatrix<f64, N, N>,
809            r: &SVector<f64, N>,
810        ) {
811            let k_scale = (0..N).map(|i| k[(i, i)].abs()).fold(0.0, f64::max);
812            let r_scale = r.iter().map(|v| v.abs()).fold(0.0, f64::max);
813            for (a, &(_, slot)) in dofs.iter().enumerate() {
814                if k[(a, a)].abs() > 1e-10 * k_scale || r[a].abs() > 1e-10 * r_scale {
815                    self.slots.insert(slot);
816                }
817            }
818        }
819    }
820
821    fn touched<const NDIM: usize, const NDOF: usize, NId: Key, E: ElementOps<NDIM, NDOF, NId>>(
822        element: &E,
823        nodes: &SlotMap<NId, Node<NDIM, NDOF>>,
824    ) -> BTreeSet<u8> {
825        let mut sink = Touched::<NId> {
826            slots: BTreeSet::new(),
827            _marker: std::marker::PhantomData,
828        };
829        element.assemble_tangent(&NodeView::new(nodes), None, &mut sink);
830        sink.slots
831    }
832
833    fn mask_slots<const NDOF: usize>(mask: DofMask) -> BTreeSet<u8> {
834        (0..NDOF)
835            .filter(|&s| mask.contains(s))
836            .map(|s| s as u8)
837            .collect()
838    }
839
840    /// `exact`: the mask must equal the touched slots (no needless
841    /// activation); otherwise it only has to cover them.
842    fn check<const NDIM: usize, const NDOF: usize, NId: Key, E: ElementOps<NDIM, NDOF, NId>>(
843        name: &str,
844        element: &E,
845        nodes: &SlotMap<NId, Node<NDIM, NDOF>>,
846        exact: bool,
847    ) {
848        let touched = touched::<NDIM, NDOF, NId, E>(element, nodes);
849        let mask = mask_slots::<NDOF>(element.dof_mask());
850        assert!(
851            touched.is_subset(&mask),
852            "{name}: stiffens slots {touched:?} but declares only {mask:?}"
853        );
854        if exact {
855            assert_eq!(
856                touched, mask,
857                "{name}: declared mask is wider than what it stiffens"
858            );
859        }
860    }
861
862    /// Two 2D nodes with a generic nonzero state on every DOF, skew to the axes
863    /// (so a truss uses both translations).
864    fn nodes_2d() -> (SlotMap<NodeId, Node>, NodeId, NodeId) {
865        let mut nodes = SlotMap::with_key();
866        let a = nodes.insert(Node::new([0.0, 0.0]));
867        let mut j = Node::new([300.0, 200.0]);
868        j.displacement = [0.4, -0.7, 0.003];
869        let b = nodes.insert(j);
870        (nodes, a, b)
871    }
872
873    fn nodes_3d() -> (SlotMap<Node3Id, Node3>, Node3Id, Node3Id) {
874        let mut nodes = SlotMap::with_key();
875        let a = nodes.insert(Node3::new([0.0, 0.0, 0.0]));
876        let mut j = Node3::new([300.0, 200.0, 150.0]);
877        j.displacement = [0.4, -0.7, 0.3, 0.002, -0.003, 0.004];
878        let b = nodes.insert(j);
879        (nodes, a, b)
880    }
881
882    fn steel() -> Material {
883        Material::Elastic { e: 200_000.0 }
884    }
885
886    /// A skew four-node patch with a generic nonzero displacement, for the continuum elements.
887    fn nodes_plane() -> (SlotMap<NodeId, Node>, [NodeId; 4]) {
888        let mut nodes = SlotMap::with_key();
889        let coords = [[0.0, 0.0], [2.0, 0.2], [2.3, 1.7], [-0.1, 1.4]];
890        let ids = coords.map(|c| {
891            let mut node = Node::new(c);
892            node.displacement = [0.01 * (c[1] + 1.0), -0.02 * (c[0] + 0.5), 0.003];
893            nodes.insert(node)
894        });
895        (nodes, ids)
896    }
897
898    #[test]
899    fn continuum_elements_declare_exactly_the_translations_they_stiffen() {
900        let (nodes, [a, b, c, d]) = nodes_plane();
901        let material = || crate::model::PlaneMaterial::plane_stress(30e3, 0.2).unwrap();
902        check(
903            "Tri3",
904            &Element::Tri3(Tri3::new(a, b, c, 0.5, material())),
905            &nodes,
906            true,
907        );
908        check(
909            "Quad4",
910            &Element::Quad4(Quad4::new([a, b, c, d], 0.5, material())),
911            &nodes,
912            true,
913        );
914        check(
915            "Quad4 enhanced",
916            &Element::Quad4(
917                Quad4::new([a, b, c, d], 0.5, material())
918                    .with_formulation(crate::model::Quad4Formulation::Enhanced),
919            ),
920            &nodes,
921            true,
922        );
923    }
924
925    #[test]
926    fn shell_declares_all_six_slots_it_stiffens() {
927        let mut nodes = SlotMap::with_key();
928        let coords = [
929            [0.0, 0.0, 0.0],
930            [2.0, 0.2, 0.3],
931            [2.3, 1.7, 0.9],
932            [-0.1, 1.4, 0.2],
933        ];
934        let ids = coords.map(|c| {
935            let mut node = Node3::new(c);
936            node.displacement = [0.01, -0.02, 0.03, 0.002, -0.003, 0.004];
937            nodes.insert(node)
938        });
939        let section =
940            crate::model::ShellSection::elastic_membrane_plate(3e4, 0.25, 0.4, 0.0).unwrap();
941        check(
942            "Shell4",
943            &Element3::Shell4(Shell4::new(ids, section.clone())),
944            &nodes,
945            true,
946        );
947        check(
948            "Shell3",
949            &Element3::Shell3(Shell3::new([ids[0], ids[1], ids[2]], section)),
950            &nodes,
951            true,
952        );
953    }
954
955    #[test]
956    fn planar_elements_declare_the_slots_they_stiffen() {
957        let (nodes, a, b) = nodes_2d();
958        let fibers = || {
959            vec![
960                Fiber::new(10.0, 100.0, steel()),
961                Fiber::new(-10.0, 100.0, steel()),
962            ]
963        };
964        let rot = 30.0_f64.to_radians();
965        let oriented =
966            || Orientation::new([rot.cos(), rot.sin(), 0.0], [-rot.sin(), rot.cos(), 0.0]).unwrap();
967
968        check(
969            "Truss",
970            &Element::Truss(Truss::new(a, b, 10.0, steel())),
971            &nodes,
972            true,
973        );
974        check(
975            "ZeroLength ux only",
976            &Element::ZeroLength(ZeroLength::new(a, b).with_material(0, steel())),
977            &nodes,
978            true,
979        );
980        check(
981            "ZeroLength rz only",
982            &Element::ZeroLength(ZeroLength::new(a, b).with_material(2, steel())),
983            &nodes,
984            true,
985        );
986        check(
987            "ZeroLength oriented",
988            &Element::ZeroLength(
989                ZeroLength::new(a, b)
990                    .with_material(0, steel())
991                    .with_material(1, steel())
992                    .with_orientation(oriented())
993                    .unwrap(),
994            ),
995            &nodes,
996            true,
997        );
998        check(
999            "ZeroLength friction",
1000            &Element::ZeroLength(
1001                ZeroLength::new(a, b)
1002                    .with_material(0, steel())
1003                    .with_friction(Friction::new(0, 1, 0.3, 500.0, 0.01)),
1004            ),
1005            &nodes,
1006            true,
1007        );
1008        check(
1009            "ZeroLengthSection",
1010            &Element::ZeroLengthSection(ZeroLengthSection::new(a, b, FiberSection::new(fibers()))),
1011            &nodes,
1012            true,
1013        );
1014        check(
1015            "ZeroLengthSection + shear spring",
1016            &Element::ZeroLengthSection(
1017                ZeroLengthSection::new(a, b, FiberSection::new(fibers())).with_material(1, steel()),
1018            ),
1019            &nodes,
1020            true,
1021        );
1022        check(
1023            "ElasticBeamColumn",
1024            &Element::ElasticBeamColumn(ElasticBeamColumn::new(
1025                a,
1026                b,
1027                200_000.0,
1028                1e4,
1029                1e8,
1030                GeomTransf::PDelta,
1031            )),
1032            &nodes,
1033            true,
1034        );
1035        check(
1036            "DispBeamColumn",
1037            &Element::DispBeamColumn(DispBeamColumn::new(
1038                a,
1039                b,
1040                fibers(),
1041                BeamIntegration::Lobatto { points: 3 },
1042            )),
1043            &nodes,
1044            true,
1045        );
1046        check(
1047            "ForceBeamColumn",
1048            &Element::ForceBeamColumn(ForceBeamColumn::new(
1049                a,
1050                b,
1051                fibers(),
1052                BeamIntegration::Lobatto { points: 3 },
1053            )),
1054            &nodes,
1055            true,
1056        );
1057    }
1058
1059    #[test]
1060    fn spatial_elements_declare_the_slots_they_stiffen() {
1061        let (nodes, a, b) = nodes_3d();
1062        let fibers = || {
1063            vec![
1064                Fiber3::new(10.0, 10.0, 50.0, steel()),
1065                Fiber3::new(10.0, -10.0, 50.0, steel()),
1066                Fiber3::new(-10.0, 10.0, 50.0, steel()),
1067                Fiber3::new(-10.0, -10.0, 50.0, steel()),
1068            ]
1069        };
1070        let tilt = Orientation::new([0.0, 0.0, 1.0], [1.0, 0.0, 0.0]).unwrap();
1071
1072        check(
1073            "Truss3",
1074            &Element3::Truss3(Truss3::new(a, b, 10.0, steel())),
1075            &nodes,
1076            true,
1077        );
1078        check(
1079            "ZeroLength3 translations",
1080            &Element3::ZeroLength3(
1081                ZeroLength3::new(a, b)
1082                    .with_material(0, steel())
1083                    .with_material(1, steel())
1084                    .with_material(2, steel()),
1085            ),
1086            &nodes,
1087            true,
1088        );
1089        check(
1090            "ZeroLength3 torsion only",
1091            &Element3::ZeroLength3(ZeroLength3::new(a, b).with_material(3, steel())),
1092            &nodes,
1093            true,
1094        );
1095        check(
1096            "ZeroLength3 oriented",
1097            &Element3::ZeroLength3(
1098                ZeroLength3::new(a, b)
1099                    .with_material(0, steel())
1100                    .with_orientation(tilt.clone()),
1101            ),
1102            &nodes,
1103            false,
1104        );
1105        check(
1106            "ZeroLength3 friction",
1107            &Element3::ZeroLength3(
1108                ZeroLength3::new(a, b)
1109                    .with_material(0, steel())
1110                    .with_friction(Friction3::new(0, [1, 2], 0.3, 500.0, 0.01)),
1111            ),
1112            &nodes,
1113            true,
1114        );
1115        check(
1116            "ZeroLengthSection3",
1117            &Element3::ZeroLengthSection3(ZeroLengthSection3::new(
1118                a,
1119                b,
1120                FiberSection3::new(fibers()),
1121            )),
1122            &nodes,
1123            true,
1124        );
1125        check(
1126            "ZeroLengthSection3 + torsion spring",
1127            &Element3::ZeroLengthSection3(
1128                ZeroLengthSection3::new(a, b, FiberSection3::new(fibers()))
1129                    .with_material(3, steel()),
1130            ),
1131            &nodes,
1132            true,
1133        );
1134        check(
1135            "ElasticBeamColumn3",
1136            &Element3::ElasticBeamColumn3(ElasticBeamColumn3::new(
1137                a,
1138                b,
1139                200_000.0,
1140                80_000.0,
1141                1e4,
1142                1e7,
1143                5e7,
1144                1e8,
1145                GeomTransf3::Linear3 {
1146                    vec_xz: [0.0, 0.0, 1.0],
1147                },
1148            )),
1149            &nodes,
1150            true,
1151        );
1152        check(
1153            "DispBeamColumn3",
1154            &Element3::DispBeamColumn3(DispBeamColumn3::new(
1155                a,
1156                b,
1157                80_000.0,
1158                1e7,
1159                [0.0, 0.0, 1.0],
1160                fibers(),
1161                BeamIntegration::Lobatto { points: 3 },
1162            )),
1163            &nodes,
1164            true,
1165        );
1166        check(
1167            "ForceBeamColumn3",
1168            &Element3::ForceBeamColumn3(ForceBeamColumn3::new(
1169                a,
1170                b,
1171                80_000.0,
1172                1e7,
1173                [0.0, 0.0, 1.0],
1174                fibers(),
1175                BeamIntegration::Lobatto { points: 3 },
1176            )),
1177            &nodes,
1178            true,
1179        );
1180    }
1181}
carapace_core/model/elements/mod.rs

carapace_core/model/elements/
mod.rs
