carapace_core/model/elements/
tri3.rs

1//! `Tri3`: the 3-node constant-strain triangle. One Gauss point, 6 DOFs
2//! `[ux, uy]` per node, plane stress or strain through its `PlaneMaterial`.
3
4use nalgebra::{SMatrix, SVector};
5
6use super::plane_common::{
7    b_matrix, check_section, load_vector, plane_coords, plane_displacement, plane_dofs, PlaneView,
8};
9use super::{DofMask, ElementForce, GaussResponse, TangentSink, VectorSink};
10use crate::model::continuum::{physical_derivatives, tri3_dshape, validate_triangle, Jacobian};
11use crate::model::{ElementLoad, NodeId, PlaneMaterial, PlaneVector};
12
13#[derive(Debug, Clone)]
14pub struct Tri3 {
15    pub nodes: [NodeId; 3],
16    pub thickness: f64,
17    pub density: f64,
18    /// The one Gauss point's material (an independent copy of the prototype).
19    material: PlaneMaterial,
20    cache: Option<Geometry>,
21}
22
23/// `B` and the area, which depend only on the (fixed) coordinates.
24#[derive(Debug, Clone, Copy, PartialEq)]
25struct Geometry {
26    b: SMatrix<f64, 3, 6>,
27    area: f64,
28}
29
30impl Tri3 {
31    pub fn new(
32        n1: NodeId,
33        n2: NodeId,
34        n3: NodeId,
35        thickness: f64,
36        material: PlaneMaterial,
37    ) -> Self {
38        Tri3 {
39            nodes: [n1, n2, n3],
40            thickness,
41            density: 0.0,
42            material,
43            cache: None,
44        }
45    }
46
47    pub fn with_density(mut self, density: f64) -> Self {
48        self.density = density;
49        self
50    }
51
52    pub fn material(&self) -> &PlaneMaterial {
53        &self.material
54    }
55
56    fn compute_geometry(&self, view: &PlaneView<'_>) -> Geometry {
57        let coords = plane_coords(&self.nodes, view);
58        let dn = tri3_dshape();
59        let jacobian = Jacobian::new(&coords, &dn);
60        Geometry {
61            b: b_matrix(&physical_derivatives(&jacobian, &dn)),
62            area: 0.5 * jacobian.det,
63        }
64    }
65
66    fn geometry(&self, view: &PlaneView<'_>) -> Geometry {
67        self.cache.unwrap_or_else(|| self.compute_geometry(view))
68    }
69
70    pub(super) fn dof_mask(&self) -> DofMask {
71        DofMask::none().with(0).with(1)
72    }
73
74    pub(super) fn validate(&self, view: &PlaneView<'_>) -> Result<(), &'static str> {
75        check_section(self.thickness, self.density, &self.material)?;
76        validate_triangle(&plane_coords(&self.nodes, view))
77    }
78
79    pub(super) fn prepare(&mut self, view: &PlaneView<'_>) {
80        self.cache = Some(self.compute_geometry(view));
81    }
82
83    fn strain(&self, g: &Geometry, view: &PlaneView<'_>) -> PlaneVector {
84        g.b * plane_displacement::<3, 6>(&self.nodes, view)
85    }
86
87    pub(super) fn assemble_tangent<S: TangentSink<NodeId>>(
88        &self,
89        view: &PlaneView<'_>,
90        sink: &mut S,
91    ) {
92        let g = self.geometry(view);
93        let (stress, d) = self.material.trial_stress_tangent(&self.strain(&g, view));
94        let scale = self.thickness * g.area;
95        let k = scale * g.b.transpose() * d * g.b;
96        let r = scale * g.b.transpose() * stress;
97        sink.add(&plane_dofs::<3, 6>(&self.nodes), &k, &r);
98    }
99
100    /// Exact lumped mass `rho t A / 3` on each translational DOF of each node.
101    pub(super) fn assemble_mass<S: VectorSink<NodeId>>(&self, view: &PlaneView<'_>, sink: &mut S) {
102        let m = self.density * self.thickness * self.geometry(view).area / 3.0;
103        sink.add(
104            &plane_dofs::<3, 6>(&self.nodes),
105            &SVector::<f64, 6>::repeat(m),
106        );
107    }
108
109    /// Consistent nodal forces of `load` (global DOF order): body force
110    /// `b t A / 3` per node, plus edge tractions and pressures.
111    pub(super) fn load_vector(&self, view: &PlaneView<'_>, load: &ElementLoad) -> SVector<f64, 6> {
112        let share = self.thickness * self.geometry(view).area / 3.0;
113        load_vector(
114            &plane_coords(&self.nodes, view),
115            &[share; 3],
116            self.thickness,
117            load,
118        )
119    }
120
121    pub(super) fn assemble_load<S: VectorSink<NodeId>>(
122        &self,
123        view: &PlaneView<'_>,
124        load: &ElementLoad,
125        sink: &mut S,
126    ) {
127        sink.add(
128            &plane_dofs::<3, 6>(&self.nodes),
129            &self.load_vector(view, load),
130        );
131    }
132
133    pub(super) fn commit(&mut self, view: &PlaneView<'_>) {
134        let g = self.geometry(view);
135        self.material = self.material.commit(&self.strain(&g, view));
136    }
137
138    /// The nodal resisting forces in global DOF order.
139    pub(super) fn local_force(&self, view: &PlaneView<'_>) -> ElementForce {
140        let g = self.geometry(view);
141        let (stress, _) = self.material.trial_stress_tangent(&self.strain(&g, view));
142        (self.thickness * g.area * g.b.transpose() * stress).into()
143    }
144
145    pub(super) fn gauss_responses(&self, view: &PlaneView<'_>) -> Vec<GaussResponse> {
146        let g = self.geometry(view);
147        let strain = self.strain(&g, view);
148        let (stress, _) = self.material.trial_stress_tangent(&strain);
149        vec![GaussResponse {
150            strain: strain.into(),
151            stress: stress.into(),
152        }]
153    }
154
155    #[cfg(test)]
156    pub(crate) fn cached(&self) -> bool {
157        self.cache.is_some()
158    }
159}
160
161#[cfg(test)]
162mod tests {
163    use nalgebra::SymmetricEigen;
164    use slotmap::SlotMap;
165
166    use super::*;
167    use crate::model::{DofRef, ElementLoad, Node};
168
169    struct Collect {
170        k: SMatrix<f64, 6, 6>,
171        r: SVector<f64, 6>,
172    }
173
174    impl TangentSink<NodeId> for Collect {
175        fn add<const N: usize>(
176            &mut self,
177            _dofs: &[DofRef<NodeId>; N],
178            k: &SMatrix<f64, N, N>,
179            r: &SVector<f64, N>,
180        ) {
181            self.k = SMatrix::from_fn(|i, j| k[(i, j)]);
182            self.r = SVector::from_fn(|i, _| r[i]);
183        }
184    }
185
186    struct Mass(SVector<f64, 6>);
187
188    impl VectorSink<NodeId> for Mass {
189        fn add<const N: usize>(&mut self, _dofs: &[DofRef<NodeId>; N], v: &SVector<f64, N>) {
190            self.0 = SVector::from_fn(|i, _| v[i]);
191        }
192    }
193
194    fn setup(coords: [[f64; 2]; 3]) -> (SlotMap<NodeId, Node>, [NodeId; 3]) {
195        let mut nodes = SlotMap::with_key();
196        let ids = coords.map(|c| nodes.insert(Node::new(c)));
197        (nodes, ids)
198    }
199
200    fn assemble(el: &Tri3, nodes: &SlotMap<NodeId, Node>) -> Collect {
201        let mut sink = Collect {
202            k: SMatrix::zeros(),
203            r: SVector::zeros(),
204        };
205        el.assemble_tangent(&PlaneView::new(nodes), &mut sink);
206        sink
207    }
208
209    fn unit_triangle(e: f64, nu: f64, t: f64) -> (Tri3, SlotMap<NodeId, Node>) {
210        let (nodes, [a, b, c]) = setup([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]);
211        (
212            Tri3::new(a, b, c, t, PlaneMaterial::plane_stress(e, nu).unwrap()),
213            nodes,
214        )
215    }
216
217    #[test]
218    fn stiffness_matches_the_closed_form_for_the_unit_right_triangle() {
219        // E = 1, nu = 0, t = 1: K = A B^T D B with D = diag(1, 1, 1/2) and A = 1/2; rows sum to zero
220        // (rigid translation) and e.g. K[ux1, ux1] = A * 1 = 1/2, K[ux0, ux0] = A * (1 + 1/2) = 3/4.
221        let (el, nodes) = unit_triangle(1.0, 0.0, 1.0);
222        let k = assemble(&el, &nodes).k;
223        #[rustfmt::skip]
224        let expected = SMatrix::<f64, 6, 6>::from_row_slice(&[
225             0.75,  0.25, -0.5,  -0.25, -0.25,  0.0,
226             0.25,  0.75,  0.0,  -0.25, -0.25, -0.5,
227            -0.5,   0.0,   0.5,   0.0,   0.0,   0.0,
228            -0.25, -0.25,  0.0,   0.25,  0.25,  0.0,
229            -0.25, -0.25,  0.0,   0.25,  0.25,  0.0,
230             0.0,  -0.5,   0.0,   0.0,   0.0,   0.5,
231        ]);
232        assert!((k - expected).abs().max() < 1e-14, "{k}");
233    }
234
235    #[test]
236    fn stiffness_is_symmetric_with_exactly_three_rigid_body_modes() {
237        let (nodes, [a, b, c]) = setup([[0.3, 0.1], [2.1, 0.4], [0.9, 1.7]]);
238        let el = Tri3::new(
239            a,
240            b,
241            c,
242            0.4,
243            PlaneMaterial::plane_strain(30e3, 0.2).unwrap(),
244        );
245        let k = assemble(&el, &nodes).k;
246        assert!((k - k.transpose()).abs().max() < 1e-9 * k.abs().max());
247        let eigen = SymmetricEigen::new(k).eigenvalues;
248        let largest = eigen.max();
249        let zero = eigen.iter().filter(|&&l| l.abs() < 1e-9 * largest).count();
250        assert_eq!(zero, 3, "{eigen}");
251        assert!(
252            eigen.iter().all(|&l| l > -1e-9 * largest),
253            "not positive semidefinite: {eigen}"
254        );
255    }
256
257    #[test]
258    fn a_rigid_motion_is_stress_free() {
259        let (mut nodes, ids) = setup([[0.3, 0.1], [2.1, 0.4], [0.9, 1.7]]);
260        let el = Tri3::new(
261            ids[0],
262            ids[1],
263            ids[2],
264            0.4,
265            PlaneMaterial::plane_stress(30e3, 0.3).unwrap(),
266        );
267        // Small rotation about the origin plus a translation: u = t + theta * [-y, x].
268        let (t, theta) = ([0.01, -0.02], 1e-5);
269        for id in ids {
270            let [x, y] = nodes[id].coords;
271            nodes[id].displacement = [t[0] - theta * y, t[1] + theta * x, 0.0];
272        }
273        let r = assemble(&el, &nodes).r;
274        assert!(r.abs().max() < 1e-9, "{r}");
275    }
276
277    #[test]
278    fn lumped_mass_sums_to_rho_t_a_in_each_direction() {
279        let (nodes, [a, b, c]) = setup([[0.0, 0.0], [4.0, 0.0], [1.0, 3.0]]);
280        let (rho, t) = (2.5, 0.2);
281        let el =
282            Tri3::new(a, b, c, t, PlaneMaterial::plane_stress(1.0, 0.0).unwrap()).with_density(rho);
283        let mut sink = Mass(SVector::zeros());
284        el.assemble_mass(&PlaneView::new(&nodes), &mut sink);
285        let area = 6.0;
286        let (sum_x, sum_y): (f64, f64) = (
287            (0..3).map(|i| sink.0[2 * i]).sum(),
288            (0..3).map(|i| sink.0[2 * i + 1]).sum(),
289        );
290        assert!((sum_x - rho * t * area).abs() < 1e-12 && (sum_y - rho * t * area).abs() < 1e-12);
291        assert!(sink
292            .0
293            .iter()
294            .all(|&m| (m - rho * t * area / 3.0).abs() < 1e-12));
295    }
296
297    #[test]
298    fn prepare_is_idempotent_and_the_cache_equals_recomputation() {
299        let (nodes, [a, b, c]) = setup([[0.3, 0.1], [2.1, 0.4], [0.9, 1.7]]);
300        let mut el = Tri3::new(
301            a,
302            b,
303            c,
304            0.4,
305            PlaneMaterial::plane_stress(30e3, 0.3).unwrap(),
306        );
307        let before = assemble(&el, &nodes).k;
308        assert!(!el.cached());
309        let view = PlaneView::new(&nodes);
310        el.prepare(&view);
311        el.prepare(&view);
312        assert!(el.cached());
313        assert_eq!(before, assemble(&el, &nodes).k);
314        assert_eq!(el.cache, Some(el.compute_geometry(&view)));
315    }
316
317    #[test]
318    fn invalid_geometry_and_section_are_reported() {
319        let (nodes, [a, b, c]) = setup([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]);
320        let view = PlaneView::new(&nodes);
321        let material = PlaneMaterial::plane_stress(1.0, 0.0).unwrap();
322        assert!(Tri3::new(a, b, c, 1.0, material.clone())
323            .validate(&view)
324            .is_ok());
325        assert!(Tri3::new(a, c, b, 1.0, material.clone())
326            .validate(&view)
327            .unwrap_err()
328            .contains("clockwise"));
329        assert!(Tri3::new(a, b, c, 0.0, material.clone())
330            .validate(&view)
331            .is_err());
332        assert!(Tri3::new(a, b, c, 1.0, material.clone())
333            .with_density(-1.0)
334            .validate(&view)
335            .is_err());
336        let bad = PlaneMaterial::PlaneStress { e: -1.0, nu: 0.0 };
337        assert!(Tri3::new(a, b, c, 1.0, bad).validate(&view).is_err());
338    }
339
340    #[test]
341    fn gauss_response_reports_strain_and_the_stress_it_implies() {
342        let (mut nodes, ids) = setup([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]);
343        let material = PlaneMaterial::plane_stress(100.0, 0.25).unwrap();
344        let el = Tri3::new(ids[0], ids[1], ids[2], 1.0, material.clone());
345        // u = (a x + b y, c x + d y): constant strain [a, d, b + c].
346        let (a, b, c, d) = (1e-3, 2e-3, -3e-4, 5e-4);
347        for id in ids {
348            let [x, y] = nodes[id].coords;
349            nodes[id].displacement = [a * x + b * y, c * x + d * y, 0.0];
350        }
351        let response = el.gauss_responses(&PlaneView::new(&nodes));
352        assert_eq!(response.len(), 1);
353        let strain = PlaneVector::new(a, d, b + c);
354        assert!((PlaneVector::from(response[0].strain) - strain).norm() < 1e-15);
355        assert!(
356            (PlaneVector::from(response[0].stress) - material.d_matrix() * strain).norm() < 1e-12
357        );
358    }
359
360    #[test]
361    fn body_force_and_edge_loads_total_their_resultants() {
362        let (nodes, [a, b, c]) = setup([[0.0, 0.0], [4.0, 0.0], [1.0, 3.0]]);
363        let el = Tri3::new(a, b, c, 0.2, PlaneMaterial::plane_stress(1.0, 0.0).unwrap());
364        let view = PlaneView::new(&nodes);
365        let sums = |v: SVector<f64, 6>| {
366            [
367                (0..3).map(|i| v[2 * i]).sum::<f64>(),
368                (0..3).map(|i| v[2 * i + 1]).sum::<f64>(),
369            ]
370        };
371        let body = sums(el.load_vector(&view, &ElementLoad::body(1.5, -4.0)));
372        assert!(
373            (body[0] - 1.5 * 0.2 * 6.0).abs() < 1e-14 && (body[1] + 4.0 * 0.2 * 6.0).abs() < 1e-14
374        );
375        // Edge 0 is the base (length 4), edge 1 runs (4,0) -> (1,3) (length sqrt(18)).
376        let base = sums(el.load_vector(&view, &ElementLoad::edge_traction(0, 2.0, 3.0)));
377        assert!(
378            (base[0] - 2.0 * 0.2 * 4.0).abs() < 1e-13 && (base[1] - 3.0 * 0.2 * 4.0).abs() < 1e-13
379        );
380        // Pressure on the base pushes up (+y, the inward normal of a counter-clockwise triangle).
381        let pressure = sums(el.load_vector(&view, &ElementLoad::edge_pressure(0, 5.0)));
382        assert!(pressure[0].abs() < 1e-13 && (pressure[1] - 5.0 * 0.2 * 4.0).abs() < 1e-13);
383        // Edge 1's pressure resultant is its length times the unit inward normal, i.e. (-3, -3) * p t / ... .
384        let slanted = sums(el.load_vector(&view, &ElementLoad::edge_pressure(1, 1.0)));
385        // Tangent (-3, 3), inward (left) normal (-3, -3)/|.|; resultant = p t |edge| n = p t (-3, -3).
386        assert!((slanted[0] + 0.6).abs() < 1e-13 && (slanted[1] + 0.6).abs() < 1e-13);
387    }
388}
carapace_core/model/elements/tri3.rs

carapace_core/model/elements/
tri3.rs
