1use std::collections::BTreeMap;
2
3use runmat_meshing_core::{AnalysisMeshArtifact, VolumeElementKind};
4
5use crate::operator::CsrMatrix;
6
7use super::elements::solid::{
8 global_stiffness_matrix as tetrahedron4_global_stiffness_matrix, SolidMaterial,
9 Tetrahedron4ElementGeometry, TETRAHEDRON4_ELEMENT_DOF_COUNT, TETRAHEDRON4_NODE_DOF_COUNT,
10};
11
12#[derive(Debug, Clone, Copy, PartialEq, Eq)]
13pub struct SolidAssemblyTopology {
14 pub dof_count: usize,
15 pub node_count: usize,
16 pub volume_element_count: usize,
17}
18
19#[derive(Debug, Clone, PartialEq, Eq)]
20pub enum SolidAssemblyError {
21 EmptyNodes,
22 EmptyVolumeElements,
23 UnsupportedVolumeElementKind { element_id: String },
24 UnknownElementNode { element_id: String, node_id: u32 },
25 InvalidElementNodeCount { element_id: String, actual: usize },
26 ElementStiffness { element_id: String, message: String },
27}
28
29pub fn solid_topology_from_analysis_mesh(
30 mesh: &AnalysisMeshArtifact,
31 base_dof_count: usize,
32) -> Result<SolidAssemblyTopology, SolidAssemblyError> {
33 if mesh.nodes.is_empty() {
34 return Err(SolidAssemblyError::EmptyNodes);
35 }
36 if mesh.volume_elements.is_empty() {
37 return Err(SolidAssemblyError::EmptyVolumeElements);
38 }
39 for element in &mesh.volume_elements {
40 if !matches!(element.kind, VolumeElementKind::Tetrahedron4) {
41 return Err(SolidAssemblyError::UnsupportedVolumeElementKind {
42 element_id: element.element_id.clone(),
43 });
44 }
45 }
46 Ok(SolidAssemblyTopology {
47 dof_count: mesh.nodes.len().saturating_mul(3).max(base_dof_count),
48 node_count: mesh.nodes.len(),
49 volume_element_count: mesh.volume_elements.len(),
50 })
51}
52
53pub fn assemble_solid_stiffness_dense(
54 mesh: &AnalysisMeshArtifact,
55 material: SolidMaterial,
56 base_dof_count: usize,
57) -> Result<Vec<f64>, SolidAssemblyError> {
58 let topology = solid_topology_from_analysis_mesh(mesh, base_dof_count)?;
59 let mut node_offsets = BTreeMap::<u32, usize>::new();
60 for (index, node) in mesh.nodes.iter().enumerate() {
61 node_offsets.insert(node.node_id, index * TETRAHEDRON4_NODE_DOF_COUNT);
62 }
63
64 let mut dense = vec![0.0_f64; topology.dof_count * topology.dof_count];
65 for element in &mesh.volume_elements {
66 if element.node_ids.len() != 4 {
67 return Err(SolidAssemblyError::InvalidElementNodeCount {
68 element_id: element.element_id.clone(),
69 actual: element.node_ids.len(),
70 });
71 }
72 let mut nodes_m = [[0.0_f64; 3]; 4];
73 let mut dof_offsets = [0_usize; 4];
74 for (local_index, node_id) in element.node_ids.iter().copied().enumerate() {
75 let node_index = mesh
76 .nodes
77 .iter()
78 .position(|node| node.node_id == node_id)
79 .ok_or_else(|| SolidAssemblyError::UnknownElementNode {
80 element_id: element.element_id.clone(),
81 node_id,
82 })?;
83 nodes_m[local_index] = mesh.nodes[node_index].coordinates_m;
84 dof_offsets[local_index] = *node_offsets.get(&node_id).ok_or_else(|| {
85 SolidAssemblyError::UnknownElementNode {
86 element_id: element.element_id.clone(),
87 node_id,
88 }
89 })?;
90 }
91 let element_stiffness =
92 tetrahedron4_global_stiffness_matrix(material, Tetrahedron4ElementGeometry { nodes_m })
93 .map_err(|err| SolidAssemblyError::ElementStiffness {
94 element_id: element.element_id.clone(),
95 message: err.to_string(),
96 })?;
97 scatter_tetrahedron4(
98 &mut dense,
99 topology.dof_count,
100 dof_offsets,
101 &element_stiffness,
102 );
103 }
104 Ok(dense)
105}
106
107pub fn assemble_solid_stiffness_csr(
108 mesh: &AnalysisMeshArtifact,
109 material: SolidMaterial,
110 base_dof_count: usize,
111) -> Result<CsrMatrix, SolidAssemblyError> {
112 assemble_solid_stiffness_csr_with_materials(mesh, material, &BTreeMap::new(), base_dof_count)
113}
114
115pub fn assemble_solid_stiffness_csr_with_materials(
116 mesh: &AnalysisMeshArtifact,
117 default_material: SolidMaterial,
118 materials_by_region: &BTreeMap<String, SolidMaterial>,
119 base_dof_count: usize,
120) -> Result<CsrMatrix, SolidAssemblyError> {
121 let topology = solid_topology_from_analysis_mesh(mesh, base_dof_count)?;
122 let mut node_offsets = BTreeMap::<u32, usize>::new();
123 for (index, node) in mesh.nodes.iter().enumerate() {
124 node_offsets.insert(node.node_id, index * TETRAHEDRON4_NODE_DOF_COUNT);
125 }
126
127 let mut rows = (0..topology.dof_count)
128 .map(|_| BTreeMap::<usize, f64>::new())
129 .collect::<Vec<_>>();
130 for element in &mesh.volume_elements {
131 if element.node_ids.len() != 4 {
132 return Err(SolidAssemblyError::InvalidElementNodeCount {
133 element_id: element.element_id.clone(),
134 actual: element.node_ids.len(),
135 });
136 }
137 let mut nodes_m = [[0.0_f64; 3]; 4];
138 let mut dof_offsets = [0_usize; 4];
139 for (local_index, node_id) in element.node_ids.iter().copied().enumerate() {
140 let node_index = mesh
141 .nodes
142 .iter()
143 .position(|node| node.node_id == node_id)
144 .ok_or_else(|| SolidAssemblyError::UnknownElementNode {
145 element_id: element.element_id.clone(),
146 node_id,
147 })?;
148 nodes_m[local_index] = mesh.nodes[node_index].coordinates_m;
149 dof_offsets[local_index] = *node_offsets.get(&node_id).ok_or_else(|| {
150 SolidAssemblyError::UnknownElementNode {
151 element_id: element.element_id.clone(),
152 node_id,
153 }
154 })?;
155 }
156 let material = materials_by_region
157 .get(element.material_region_id.as_str())
158 .copied()
159 .unwrap_or(default_material);
160 let element_stiffness =
161 tetrahedron4_global_stiffness_matrix(material, Tetrahedron4ElementGeometry { nodes_m })
162 .map_err(|err| SolidAssemblyError::ElementStiffness {
163 element_id: element.element_id.clone(),
164 message: err.to_string(),
165 })?;
166 scatter_tetrahedron4_csr_rows(&mut rows, dof_offsets, &element_stiffness);
167 }
168 Ok(rows_to_csr(rows))
169}
170
171fn scatter_tetrahedron4(
172 dense: &mut [f64],
173 dof_count: usize,
174 dof_offsets: [usize; 4],
175 element_stiffness: &[[f64; TETRAHEDRON4_ELEMENT_DOF_COUNT]; TETRAHEDRON4_ELEMENT_DOF_COUNT],
176) {
177 for local_row_node in 0..4 {
178 for local_row_axis in 0..TETRAHEDRON4_NODE_DOF_COUNT {
179 let local_row = local_row_node * TETRAHEDRON4_NODE_DOF_COUNT + local_row_axis;
180 let global_row = dof_offsets[local_row_node] + local_row_axis;
181 for (local_col_node, global_col_offset) in dof_offsets.iter().enumerate() {
182 for local_col_axis in 0..TETRAHEDRON4_NODE_DOF_COUNT {
183 let local_col = local_col_node * TETRAHEDRON4_NODE_DOF_COUNT + local_col_axis;
184 let global_col = global_col_offset + local_col_axis;
185 dense[global_row * dof_count + global_col] +=
186 element_stiffness[local_row][local_col];
187 }
188 }
189 }
190 }
191}
192
193fn scatter_tetrahedron4_csr_rows(
194 rows: &mut [BTreeMap<usize, f64>],
195 dof_offsets: [usize; 4],
196 element_stiffness: &[[f64; TETRAHEDRON4_ELEMENT_DOF_COUNT]; TETRAHEDRON4_ELEMENT_DOF_COUNT],
197) {
198 for local_row_node in 0..4 {
199 for local_row_axis in 0..TETRAHEDRON4_NODE_DOF_COUNT {
200 let local_row = local_row_node * TETRAHEDRON4_NODE_DOF_COUNT + local_row_axis;
201 let global_row = dof_offsets[local_row_node] + local_row_axis;
202 for (local_col_node, global_col_offset) in dof_offsets.iter().enumerate() {
203 for local_col_axis in 0..TETRAHEDRON4_NODE_DOF_COUNT {
204 let local_col = local_col_node * TETRAHEDRON4_NODE_DOF_COUNT + local_col_axis;
205 let global_col = global_col_offset + local_col_axis;
206 *rows[global_row].entry(global_col).or_insert(0.0) +=
207 element_stiffness[local_row][local_col];
208 }
209 }
210 }
211 }
212}
213
214fn rows_to_csr(rows: Vec<BTreeMap<usize, f64>>) -> CsrMatrix {
215 let mut row_offsets = Vec::with_capacity(rows.len() + 1);
216 let mut column_indices = Vec::new();
217 let mut values = Vec::new();
218 row_offsets.push(0);
219 for row in rows {
220 for (column, value) in row {
221 if value.abs() > 0.0 {
222 column_indices.push(column);
223 values.push(value);
224 }
225 }
226 row_offsets.push(values.len());
227 }
228 CsrMatrix {
229 row_offsets,
230 column_indices,
231 values,
232 }
233}
234
235#[cfg(test)]
236mod tests {
237 use super::*;
238 use runmat_meshing_core::{
239 AnalysisMeshNode, AnalysisMeshProvenance, AnalysisMeshQualityReport, AnalysisVolumeElement,
240 MeshSizingField,
241 };
242
243 fn mesh(kind: VolumeElementKind) -> AnalysisMeshArtifact {
244 AnalysisMeshArtifact {
245 schema_version: "analysis-mesh/v1".to_string(),
246 mesh_id: "mesh".to_string(),
247 nodes: vec![
248 AnalysisMeshNode {
249 node_id: 1,
250 coordinates_m: [0.0, 0.0, 0.0],
251 provenance: Vec::new(),
252 },
253 AnalysisMeshNode {
254 node_id: 2,
255 coordinates_m: [1.0, 0.0, 0.0],
256 provenance: Vec::new(),
257 },
258 AnalysisMeshNode {
259 node_id: 3,
260 coordinates_m: [0.0, 1.0, 0.0],
261 provenance: Vec::new(),
262 },
263 AnalysisMeshNode {
264 node_id: 4,
265 coordinates_m: [0.0, 0.0, 1.0],
266 provenance: Vec::new(),
267 },
268 ],
269 volume_elements: vec![AnalysisVolumeElement {
270 element_id: "tetrahedron_1".to_string(),
271 kind,
272 node_ids: vec![1, 2, 3, 4],
273 material_region_id: "region".to_string(),
274 provenance: Vec::new(),
275 }],
276 boundary_faces: Vec::new(),
277 boundary_edges: Vec::new(),
278 quality: AnalysisMeshQualityReport::default(),
279 sizing: MeshSizingField::default(),
280 field_topology: Vec::new(),
281 backend: Default::default(),
282 adaptive_iterations: Vec::new(),
283 provenance: AnalysisMeshProvenance {
284 algorithm: "test".to_string(),
285 source_geometry_id: "geo".to_string(),
286 source_geometry_revision: 1,
287 source_geometry_sha256: None,
288 },
289 }
290 }
291
292 #[test]
293 fn solid_topology_uses_analysis_mesh_nodes_and_tetrahedron4_elements() {
294 let topology =
295 solid_topology_from_analysis_mesh(&mesh(VolumeElementKind::Tetrahedron4), 3).unwrap();
296 assert_eq!(topology.dof_count, 12);
297 assert_eq!(topology.node_count, 4);
298 assert_eq!(topology.volume_element_count, 1);
299 }
300
301 #[test]
302 fn solid_topology_rejects_unsupported_volume_elements() {
303 let err = solid_topology_from_analysis_mesh(&mesh(VolumeElementKind::Hex8), 3)
304 .expect_err("hex solid assembly is not supported yet");
305 assert_eq!(
306 err,
307 SolidAssemblyError::UnsupportedVolumeElementKind {
308 element_id: "tetrahedron_1".to_string()
309 }
310 );
311 }
312
313 #[test]
314 fn solid_stiffness_scatter_assembles_tetrahedron4_dense_matrix() {
315 let mesh = mesh(VolumeElementKind::Tetrahedron4);
316 let dense = assemble_solid_stiffness_dense(
317 &mesh,
318 SolidMaterial {
319 youngs_modulus_pa: 200.0e9,
320 poisson_ratio: 0.3,
321 },
322 3,
323 )
324 .expect("Tetrahedron4 stiffness should assemble");
325 let dof_count = 12;
326 assert_eq!(dense.len(), dof_count * dof_count);
327 for row in 0..dof_count {
328 assert!(dense[row * dof_count + row] > 0.0);
329 for col in 0..dof_count {
330 assert!(
331 (dense[row * dof_count + col] - dense[col * dof_count + row]).abs() < 1.0e-5
332 );
333 }
334 }
335 }
336
337 #[test]
338 fn solid_stiffness_scatter_assembles_tetrahedron4_csr_matrix() {
339 let mesh = mesh(VolumeElementKind::Tetrahedron4);
340 let csr = assemble_solid_stiffness_csr(
341 &mesh,
342 SolidMaterial {
343 youngs_modulus_pa: 200.0e9,
344 poisson_ratio: 0.3,
345 },
346 3,
347 )
348 .expect("Tetrahedron4 stiffness should assemble");
349 let dof_count = 12;
350 assert_eq!(csr.row_offsets.len(), dof_count + 1);
351 assert_eq!(csr.row_offsets.last().copied(), Some(csr.values.len()));
352 assert_eq!(csr.column_indices.len(), csr.values.len());
353 assert!(csr.values.len() <= dof_count * dof_count);
354 for row in 0..dof_count {
355 let start = csr.row_offsets[row];
356 let end = csr.row_offsets[row + 1];
357 assert!(csr.column_indices[start..end].binary_search(&row).is_ok());
358 }
359 }
360}