runmat_analysis_fea/assembly/elements/
shell.rs

1use serde::{Deserialize, Serialize};
2use thiserror::Error;
3
4pub const SHELL_NODE_DOF_COUNT: usize = 6;
5pub const SHELL_ELEMENT_NODE_COUNT: usize = 3;
6pub const SHELL_ELEMENT_DOF_COUNT: usize = SHELL_NODE_DOF_COUNT * SHELL_ELEMENT_NODE_COUNT;
7
8pub type ShellMatrix18 = [[f64; SHELL_ELEMENT_DOF_COUNT]; SHELL_ELEMENT_DOF_COUNT];
9
10#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
11pub struct ShellSection {
12    pub thickness_m: f64,
13    pub shear_correction: f64,
14    pub drilling_stiffness_scale: f64,
15}
16
17#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
18pub struct ShellMaterial {
19    pub youngs_modulus_pa: f64,
20    pub poisson_ratio: f64,
21    pub shear_modulus_pa: f64,
22}
23
24#[derive(Debug, Clone, Copy, PartialEq)]
25pub struct ShellElementGeometry {
26    pub nodes_m: [[f64; 3]; SHELL_ELEMENT_NODE_COUNT],
27    pub reference_axis: [f64; 3],
28}
29
30#[derive(Debug, Clone, Copy, PartialEq)]
31pub struct ShellLocalFrame {
32    pub x: [f64; 3],
33    pub y: [f64; 3],
34    pub z: [f64; 3],
35    pub area_m2: f64,
36}
37
38#[derive(Debug, Error, Clone, PartialEq)]
39pub enum ShellElementError {
40    #[error("shell element area must be positive and finite")]
41    DegenerateArea,
42    #[error("shell reference axis must be finite and non-parallel to the shell normal")]
43    DegenerateReferenceAxis,
44    #[error("shell thickness, shear correction, and drilling scale must be positive and finite")]
45    InvalidSection,
46    #[error("shell material properties must be positive and finite")]
47    InvalidMaterial,
48}
49
50impl ShellSection {
51    pub fn validate(self) -> Result<(), ShellElementError> {
52        if positive_finite(self.thickness_m)
53            && positive_finite(self.shear_correction)
54            && positive_finite(self.drilling_stiffness_scale)
55        {
56            Ok(())
57        } else {
58            Err(ShellElementError::InvalidSection)
59        }
60    }
61}
62
63impl ShellMaterial {
64    pub fn validate(self) -> Result<(), ShellElementError> {
65        if positive_finite(self.youngs_modulus_pa)
66            && positive_finite(self.shear_modulus_pa)
67            && self.poisson_ratio.is_finite()
68            && self.poisson_ratio > -1.0
69            && self.poisson_ratio < 0.5
70        {
71            Ok(())
72        } else {
73            Err(ShellElementError::InvalidMaterial)
74        }
75    }
76}
77
78impl ShellElementGeometry {
79    pub fn local_frame(self) -> Result<ShellLocalFrame, ShellElementError> {
80        let edge_a = sub(self.nodes_m[1], self.nodes_m[0]);
81        let edge_b = sub(self.nodes_m[2], self.nodes_m[0]);
82        let normal_area = cross(edge_a, edge_b);
83        let normal_norm = norm(normal_area);
84        if !positive_finite(normal_norm) {
85            return Err(ShellElementError::DegenerateArea);
86        }
87        let z = scale(normal_area, 1.0 / normal_norm);
88        if !self.reference_axis.iter().all(|value| value.is_finite()) {
89            return Err(ShellElementError::DegenerateReferenceAxis);
90        }
91        let reference_projection = sub(self.reference_axis, scale(z, dot(self.reference_axis, z)));
92        let projection_norm = norm(reference_projection);
93        let x = if projection_norm > 1.0e-12 && projection_norm.is_finite() {
94            scale(reference_projection, 1.0 / projection_norm)
95        } else {
96            let edge_norm = norm(edge_a);
97            if !positive_finite(edge_norm) {
98                return Err(ShellElementError::DegenerateReferenceAxis);
99            }
100            scale(edge_a, 1.0 / edge_norm)
101        };
102        let y = cross(z, x);
103        Ok(ShellLocalFrame {
104            x,
105            y,
106            z,
107            area_m2: 0.5 * normal_norm,
108        })
109    }
110}
111
112pub fn global_stiffness_matrix(
113    section: ShellSection,
114    material: ShellMaterial,
115    geometry: ShellElementGeometry,
116) -> Result<ShellMatrix18, ShellElementError> {
117    section.validate()?;
118    material.validate()?;
119    let frame = geometry.local_frame()?;
120    let mut k = [[0.0; SHELL_ELEMENT_DOF_COUNT]; SHELL_ELEMENT_DOF_COUNT];
121
122    let membrane = material.youngs_modulus_pa * section.thickness_m * frame.area_m2;
123    let bending = material.youngs_modulus_pa * section.thickness_m.powi(3) * frame.area_m2
124        / (12.0 * (1.0 - material.poisson_ratio.powi(2)).max(1.0e-9));
125    let shear =
126        material.shear_modulus_pa * section.thickness_m * section.shear_correction * frame.area_m2;
127    let drilling = membrane * section.drilling_stiffness_scale;
128
129    for (a, b) in [(0usize, 1usize), (1, 2), (2, 0)] {
130        let length = distance(geometry.nodes_m[a], geometry.nodes_m[b]).max(1.0e-12);
131        let membrane_k = membrane / length.powi(2);
132        let shear_k = shear / length.powi(2);
133        let bending_k = bending / length.powi(2);
134        let drilling_k = drilling / length.powi(2);
135        add_pair_spring(&mut k, a, b, 0, membrane_k);
136        add_pair_spring(&mut k, a, b, 1, membrane_k);
137        add_pair_spring(&mut k, a, b, 2, shear_k);
138        add_pair_spring(&mut k, a, b, 3, bending_k);
139        add_pair_spring(&mut k, a, b, 4, bending_k);
140        add_pair_spring(&mut k, a, b, 5, drilling_k);
141    }
142
143    let stabilization = (membrane + shear + bending + drilling).max(1.0) * 1.0e-12;
144    for (index, row) in k.iter_mut().enumerate() {
145        row[index] += stabilization;
146    }
147    Ok(k)
148}
149
150fn add_pair_spring(
151    matrix: &mut ShellMatrix18,
152    node_a: usize,
153    node_b: usize,
154    component: usize,
155    stiffness: f64,
156) {
157    let row = node_a * SHELL_NODE_DOF_COUNT + component;
158    let col = node_b * SHELL_NODE_DOF_COUNT + component;
159    matrix[row][row] += stiffness;
160    matrix[col][col] += stiffness;
161    matrix[row][col] -= stiffness;
162    matrix[col][row] -= stiffness;
163}
164
165fn positive_finite(value: f64) -> bool {
166    value.is_finite() && value > 0.0
167}
168
169fn sub(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
170    [a[0] - b[0], a[1] - b[1], a[2] - b[2]]
171}
172
173fn scale(a: [f64; 3], factor: f64) -> [f64; 3] {
174    [a[0] * factor, a[1] * factor, a[2] * factor]
175}
176
177fn dot(a: [f64; 3], b: [f64; 3]) -> f64 {
178    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
179}
180
181fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
182    [
183        a[1] * b[2] - a[2] * b[1],
184        a[2] * b[0] - a[0] * b[2],
185        a[0] * b[1] - a[1] * b[0],
186    ]
187}
188
189fn norm(a: [f64; 3]) -> f64 {
190    dot(a, a).sqrt()
191}
192
193fn distance(a: [f64; 3], b: [f64; 3]) -> f64 {
194    norm(sub(a, b))
195}
196
197#[cfg(test)]
198mod tests {
199    use super::*;
200
201    fn section() -> ShellSection {
202        ShellSection {
203            thickness_m: 0.002,
204            shear_correction: 5.0 / 6.0,
205            drilling_stiffness_scale: 1.0e-4,
206        }
207    }
208
209    fn material() -> ShellMaterial {
210        ShellMaterial {
211            youngs_modulus_pa: 200.0e9,
212            poisson_ratio: 0.3,
213            shear_modulus_pa: 76.9e9,
214        }
215    }
216
217    fn geometry() -> ShellElementGeometry {
218        ShellElementGeometry {
219            nodes_m: [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]],
220            reference_axis: [1.0, 0.0, 0.0],
221        }
222    }
223
224    #[test]
225    fn shell_local_frame_is_orthonormal() {
226        let frame = geometry().local_frame().expect("frame should build");
227        assert_close(frame.area_m2, 0.5, 1.0e-12);
228        assert_close(dot(frame.x, frame.y), 0.0, 1.0e-12);
229        assert_close(dot(frame.x, frame.z), 0.0, 1.0e-12);
230        assert_close(dot(frame.y, frame.z), 0.0, 1.0e-12);
231        assert_close(norm(frame.x), 1.0, 1.0e-12);
232        assert_close(norm(frame.y), 1.0, 1.0e-12);
233        assert_close(norm(frame.z), 1.0, 1.0e-12);
234    }
235
236    #[test]
237    fn shell_global_stiffness_is_symmetric() {
238        let k = global_stiffness_matrix(section(), material(), geometry())
239            .expect("stiffness should build");
240        for row in 0..SHELL_ELEMENT_DOF_COUNT {
241            for col in 0..SHELL_ELEMENT_DOF_COUNT {
242                assert_close(k[row][col], k[col][row], 1.0e-6);
243            }
244        }
245    }
246
247    #[test]
248    fn shell_frame_rejects_degenerate_inputs() {
249        assert_eq!(
250            ShellElementGeometry {
251                nodes_m: [[0.0, 0.0, 0.0]; 3],
252                reference_axis: [1.0, 0.0, 0.0],
253            }
254            .local_frame()
255            .expect_err("zero-area shell should fail"),
256            ShellElementError::DegenerateArea
257        );
258    }
259
260    fn assert_close(actual: f64, expected: f64, tolerance: f64) {
261        assert!(
262            (actual - expected).abs() <= tolerance,
263            "actual={actual} expected={expected} tolerance={tolerance}",
264        );
265    }
266}
runmat_analysis_fea/assembly/elements/shell.rs

runmat_analysis_fea/assembly/elements/
shell.rs
