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