1use serde::{Deserialize, Serialize};
2use thiserror::Error;
3
4pub const BEAM_NODE_DOF_COUNT: usize = 6;
5pub const BEAM_ELEMENT_DOF_COUNT: usize = 12;
6
7pub type BeamMatrix12 = [[f64; BEAM_ELEMENT_DOF_COUNT]; BEAM_ELEMENT_DOF_COUNT];
8pub type BeamTransform12 = [[f64; BEAM_ELEMENT_DOF_COUNT]; BEAM_ELEMENT_DOF_COUNT];
9
10#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
11pub struct BeamSection {
12 pub area_m2: f64,
13 pub iy_m4: f64,
14 pub iz_m4: f64,
15 pub torsion_j_m4: f64,
16 #[serde(default)]
17 pub outer_fiber_y_m: f64,
18 #[serde(default)]
19 pub outer_fiber_z_m: f64,
20 #[serde(default)]
21 pub torsion_outer_radius_m: f64,
22}
23
24#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
25pub struct BeamMaterial {
26 pub youngs_modulus_pa: f64,
27 pub shear_modulus_pa: f64,
28}
29
30#[derive(Debug, Clone, Copy, PartialEq)]
31pub struct BeamElementGeometry {
32 pub node_i_m: [f64; 3],
33 pub node_j_m: [f64; 3],
34 pub reference_axis: [f64; 3],
35}
36
37#[derive(Debug, Clone, Copy, PartialEq)]
38pub struct BeamLocalFrame {
39 pub x: [f64; 3],
40 pub y: [f64; 3],
41 pub z: [f64; 3],
42 pub length_m: f64,
43}
44
45#[derive(Debug, Error, Clone, PartialEq)]
46pub enum BeamElementError {
47 #[error("beam element length must be positive and finite")]
48 DegenerateLength,
49 #[error("beam reference axis must be finite and non-parallel to the beam axis")]
50 DegenerateReferenceAxis,
51 #[error("beam section area, second moments, and torsion constant must be positive and finite")]
52 InvalidSection,
53 #[error("beam Young's modulus and shear modulus must be positive and finite")]
54 InvalidMaterial,
55}
56
57impl BeamSection {
58 pub fn validate(self) -> Result<(), BeamElementError> {
59 if positive_finite(self.area_m2)
60 && positive_finite(self.iy_m4)
61 && positive_finite(self.iz_m4)
62 && positive_finite(self.torsion_j_m4)
63 {
64 Ok(())
65 } else {
66 Err(BeamElementError::InvalidSection)
67 }
68 }
69}
70
71impl BeamMaterial {
72 pub fn validate(self) -> Result<(), BeamElementError> {
73 if positive_finite(self.youngs_modulus_pa) && positive_finite(self.shear_modulus_pa) {
74 Ok(())
75 } else {
76 Err(BeamElementError::InvalidMaterial)
77 }
78 }
79}
80
81impl BeamElementGeometry {
82 pub fn local_frame(self) -> Result<BeamLocalFrame, BeamElementError> {
83 let axis = sub(self.node_j_m, self.node_i_m);
84 let length_m = norm(axis);
85 if !positive_finite(length_m) {
86 return Err(BeamElementError::DegenerateLength);
87 }
88 let x = scale(axis, 1.0 / length_m);
89 if !self.reference_axis.iter().all(|value| value.is_finite()) {
90 return Err(BeamElementError::DegenerateReferenceAxis);
91 }
92 let reference_projection = sub(self.reference_axis, scale(x, dot(self.reference_axis, x)));
93 let projection_norm = norm(reference_projection);
94 if projection_norm <= 1.0e-12 || !projection_norm.is_finite() {
95 return Err(BeamElementError::DegenerateReferenceAxis);
96 }
97 let y = scale(reference_projection, 1.0 / projection_norm);
98 let z = cross(x, y);
99 Ok(BeamLocalFrame { x, y, z, length_m })
100 }
101}
102
103pub fn local_stiffness_matrix(
104 section: BeamSection,
105 material: BeamMaterial,
106 length_m: f64,
107) -> Result<BeamMatrix12, BeamElementError> {
108 section.validate()?;
109 material.validate()?;
110 if !positive_finite(length_m) {
111 return Err(BeamElementError::DegenerateLength);
112 }
113
114 let l = length_m;
115 let l2 = l * l;
116 let l3 = l2 * l;
117 let ea_l = material.youngs_modulus_pa * section.area_m2 / l;
118 let gj_l = material.shear_modulus_pa * section.torsion_j_m4 / l;
119 let eiy = material.youngs_modulus_pa * section.iy_m4;
120 let eiz = material.youngs_modulus_pa * section.iz_m4;
121
122 let mut k = [[0.0; BEAM_ELEMENT_DOF_COUNT]; BEAM_ELEMENT_DOF_COUNT];
123
124 add_symmetric(&mut k, 0, 0, ea_l);
125 add_symmetric(&mut k, 0, 6, -ea_l);
126 add_symmetric(&mut k, 6, 6, ea_l);
127
128 add_symmetric(&mut k, 3, 3, gj_l);
129 add_symmetric(&mut k, 3, 9, -gj_l);
130 add_symmetric(&mut k, 9, 9, gj_l);
131
132 add_bending_z(&mut k, eiz, l, l2, l3);
133 add_bending_y(&mut k, eiy, l, l2, l3);
134
135 Ok(k)
136}
137
138pub fn transformation_matrix(frame: BeamLocalFrame) -> BeamTransform12 {
139 let rotation = [frame.x, frame.y, frame.z];
140 let mut transform = [[0.0; BEAM_ELEMENT_DOF_COUNT]; BEAM_ELEMENT_DOF_COUNT];
141 for block in 0..4 {
142 let offset = block * 3;
143 for row in 0..3 {
144 for col in 0..3 {
145 transform[offset + row][offset + col] = rotation[row][col];
146 }
147 }
148 }
149 transform
150}
151
152pub fn global_stiffness_matrix(
153 section: BeamSection,
154 material: BeamMaterial,
155 geometry: BeamElementGeometry,
156) -> Result<BeamMatrix12, BeamElementError> {
157 let frame = geometry.local_frame()?;
158 let local = local_stiffness_matrix(section, material, frame.length_m)?;
159 let transform = transformation_matrix(frame);
160 Ok(transform_transpose_multiply(&transform, &local))
161}
162
163fn add_bending_z(k: &mut BeamMatrix12, ei: f64, l: f64, l2: f64, l3: f64) {
164 let c12 = 12.0 * ei / l3;
165 let c6 = 6.0 * ei / l2;
166 let c4 = 4.0 * ei / l;
167 let c2 = 2.0 * ei / l;
168 add_symmetric(k, 1, 1, c12);
169 add_symmetric(k, 1, 5, c6);
170 add_symmetric(k, 1, 7, -c12);
171 add_symmetric(k, 1, 11, c6);
172 add_symmetric(k, 5, 5, c4);
173 add_symmetric(k, 5, 7, -c6);
174 add_symmetric(k, 5, 11, c2);
175 add_symmetric(k, 7, 7, c12);
176 add_symmetric(k, 7, 11, -c6);
177 add_symmetric(k, 11, 11, c4);
178}
179
180fn add_bending_y(k: &mut BeamMatrix12, ei: f64, l: f64, l2: f64, l3: f64) {
181 let c12 = 12.0 * ei / l3;
182 let c6 = 6.0 * ei / l2;
183 let c4 = 4.0 * ei / l;
184 let c2 = 2.0 * ei / l;
185 add_symmetric(k, 2, 2, c12);
186 add_symmetric(k, 2, 4, -c6);
187 add_symmetric(k, 2, 8, -c12);
188 add_symmetric(k, 2, 10, -c6);
189 add_symmetric(k, 4, 4, c4);
190 add_symmetric(k, 4, 8, c6);
191 add_symmetric(k, 4, 10, c2);
192 add_symmetric(k, 8, 8, c12);
193 add_symmetric(k, 8, 10, c6);
194 add_symmetric(k, 10, 10, c4);
195}
196
197fn add_symmetric(matrix: &mut BeamMatrix12, row: usize, col: usize, value: f64) {
198 matrix[row][col] += value;
199 if row != col {
200 matrix[col][row] += value;
201 }
202}
203
204fn transform_transpose_multiply(transform: &BeamTransform12, local: &BeamMatrix12) -> BeamMatrix12 {
205 let mut temp = [[0.0; BEAM_ELEMENT_DOF_COUNT]; BEAM_ELEMENT_DOF_COUNT];
206 for (row, temp_row) in temp.iter_mut().enumerate() {
207 for (col, value) in temp_row.iter_mut().enumerate() {
208 *value = (0..BEAM_ELEMENT_DOF_COUNT)
209 .map(|idx| local[row][idx] * transform[idx][col])
210 .sum();
211 }
212 }
213
214 let mut global = [[0.0; BEAM_ELEMENT_DOF_COUNT]; BEAM_ELEMENT_DOF_COUNT];
215 for (row, global_row) in global.iter_mut().enumerate() {
216 for (col, value) in global_row.iter_mut().enumerate() {
217 *value = (0..BEAM_ELEMENT_DOF_COUNT)
218 .map(|idx| transform[idx][row] * temp[idx][col])
219 .sum();
220 }
221 }
222 global
223}
224
225fn positive_finite(value: f64) -> bool {
226 value.is_finite() && value > 0.0
227}
228
229fn sub(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
230 [a[0] - b[0], a[1] - b[1], a[2] - b[2]]
231}
232
233fn scale(a: [f64; 3], factor: f64) -> [f64; 3] {
234 [a[0] * factor, a[1] * factor, a[2] * factor]
235}
236
237fn dot(a: [f64; 3], b: [f64; 3]) -> f64 {
238 a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
239}
240
241fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
242 [
243 a[1] * b[2] - a[2] * b[1],
244 a[2] * b[0] - a[0] * b[2],
245 a[0] * b[1] - a[1] * b[0],
246 ]
247}
248
249fn norm(a: [f64; 3]) -> f64 {
250 dot(a, a).sqrt()
251}
252
253#[cfg(test)]
254mod tests {
255 use super::*;
256
257 fn section() -> BeamSection {
258 BeamSection {
259 area_m2: 2.0e-4,
260 iy_m4: 1.6e-9,
261 iz_m4: 6.4e-9,
262 torsion_j_m4: 2.4e-9,
263 outer_fiber_y_m: 0.01,
264 outer_fiber_z_m: 0.005,
265 torsion_outer_radius_m: 0.011_180_339_887_498_949,
266 }
267 }
268
269 fn material() -> BeamMaterial {
270 BeamMaterial {
271 youngs_modulus_pa: 200.0e9,
272 shear_modulus_pa: 79.3e9,
273 }
274 }
275
276 #[test]
277 fn beam_local_frame_is_orthonormal() {
278 let frame = BeamElementGeometry {
279 node_i_m: [0.0, 0.0, 0.0],
280 node_j_m: [2.0, 0.0, 0.0],
281 reference_axis: [0.0, 0.0, 1.0],
282 }
283 .local_frame()
284 .expect("frame should build");
285
286 assert_close(frame.length_m, 2.0, 1.0e-12);
287 assert_close(dot(frame.x, frame.y), 0.0, 1.0e-12);
288 assert_close(dot(frame.x, frame.z), 0.0, 1.0e-12);
289 assert_close(dot(frame.y, frame.z), 0.0, 1.0e-12);
290 assert_close(norm(frame.x), 1.0, 1.0e-12);
291 assert_close(norm(frame.y), 1.0, 1.0e-12);
292 assert_close(norm(frame.z), 1.0, 1.0e-12);
293 }
294
295 #[test]
296 fn beam_local_stiffness_matches_closed_form_terms() {
297 let l = 2.5;
298 let k = local_stiffness_matrix(section(), material(), l).expect("matrix should build");
299 let ea_l = material().youngs_modulus_pa * section().area_m2 / l;
300 let gj_l = material().shear_modulus_pa * section().torsion_j_m4 / l;
301 let eiy = material().youngs_modulus_pa * section().iy_m4;
302 let eiz = material().youngs_modulus_pa * section().iz_m4;
303
304 assert_close(k[0][0], ea_l, 1.0e-6);
305 assert_close(k[0][6], -ea_l, 1.0e-6);
306 assert_close(k[3][3], gj_l, 1.0e-9);
307 assert_close(k[3][9], -gj_l, 1.0e-9);
308 assert_close(k[1][1], 12.0 * eiz / l.powi(3), 1.0e-6);
309 assert_close(k[1][5], 6.0 * eiz / l.powi(2), 1.0e-6);
310 assert_close(k[5][11], 2.0 * eiz / l, 1.0e-6);
311 assert_close(k[2][2], 12.0 * eiy / l.powi(3), 1.0e-6);
312 assert_close(k[2][4], -6.0 * eiy / l.powi(2), 1.0e-6);
313 assert_close(k[4][10], 2.0 * eiy / l, 1.0e-6);
314 }
315
316 #[test]
317 fn beam_local_stiffness_is_symmetric() {
318 let k = local_stiffness_matrix(section(), material(), 3.0).expect("matrix should build");
319 for row in 0..BEAM_ELEMENT_DOF_COUNT {
320 for col in 0..BEAM_ELEMENT_DOF_COUNT {
321 assert_close(k[row][col], k[col][row], 1.0e-9);
322 }
323 }
324 }
325
326 #[test]
327 fn beam_transformation_is_block_orthonormal() {
328 let frame = BeamElementGeometry {
329 node_i_m: [0.0, 0.0, 0.0],
330 node_j_m: [1.0, 1.0, 0.0],
331 reference_axis: [0.0, 0.0, 1.0],
332 }
333 .local_frame()
334 .expect("frame should build");
335 let transform = transformation_matrix(frame);
336
337 for row in 0..3 {
338 for col in 0..3 {
339 let value: f64 = (0..3)
340 .map(|idx| transform[row][idx] * transform[col][idx])
341 .sum();
342 assert_close(value, if row == col { 1.0 } else { 0.0 }, 1.0e-12);
343 }
344 }
345 }
346
347 #[test]
348 fn beam_global_stiffness_is_symmetric() {
349 let k = global_stiffness_matrix(
350 section(),
351 material(),
352 BeamElementGeometry {
353 node_i_m: [0.0, 0.0, 0.0],
354 node_j_m: [1.0, 1.0, 0.25],
355 reference_axis: [0.0, 0.0, 1.0],
356 },
357 )
358 .expect("global matrix should build");
359
360 for row in 0..BEAM_ELEMENT_DOF_COUNT {
361 for col in 0..BEAM_ELEMENT_DOF_COUNT {
362 assert_close(k[row][col], k[col][row], 1.0e-6);
363 }
364 }
365 }
366
367 #[test]
368 fn beam_frame_rejects_degenerate_inputs() {
369 assert_eq!(
370 BeamElementGeometry {
371 node_i_m: [0.0, 0.0, 0.0],
372 node_j_m: [0.0, 0.0, 0.0],
373 reference_axis: [0.0, 1.0, 0.0],
374 }
375 .local_frame()
376 .expect_err("zero-length beam should fail"),
377 BeamElementError::DegenerateLength
378 );
379 assert_eq!(
380 BeamElementGeometry {
381 node_i_m: [0.0, 0.0, 0.0],
382 node_j_m: [1.0, 0.0, 0.0],
383 reference_axis: [2.0, 0.0, 0.0],
384 }
385 .local_frame()
386 .expect_err("parallel reference should fail"),
387 BeamElementError::DegenerateReferenceAxis
388 );
389 }
390
391 fn assert_close(actual: f64, expected: f64, tolerance: f64) {
392 assert!(
393 (actual - expected).abs() <= tolerance,
394 "actual={actual} expected={expected} tolerance={tolerance}",
395 );
396 }
397}