Skip to main content

runmat_meshing/analysis_prep/
mod.rs

1use runmat_geometry_core::{GeometryAsset, MeshKind};
2use std::collections::BTreeMap;
3
4mod metrics;
5mod types;
6
7use metrics::{
8    coordinate_active_dimension_count, coordinate_characteristic_length_m, mesh_coordinate_span_m,
9    mesh_element_geometry_metrics,
10};
11pub use types::*;
12
13pub fn prepare_geometry_for_analysis(
14    geometry: &GeometryAsset,
15    options: MeshingOptions,
16) -> Result<MeshingPrepResult, String> {
17    if geometry.meshes.is_empty() {
18        return Err("geometry has no meshes to prepare".to_string());
19    }
20    if options.target_element_budget == 0 {
21        return Err("target_element_budget must be greater than zero".to_string());
22    }
23
24    let mut source_meshes = geometry.meshes.clone();
25    source_meshes.sort_by(|a, b| a.mesh_id.cmp(&b.mesh_id));
26    let per_mesh_budget =
27        (options.target_element_budget / source_meshes.len().max(1)).max(1) as u64;
28
29    let mut prepared_meshes = Vec::with_capacity(source_meshes.len());
30    for mesh in source_meshes {
31        let profile_scale = match options.profile {
32            MeshingProfile::SurfaceOnly => 1.0,
33            MeshingProfile::AnalysisReady => {
34                if mesh.kind == MeshKind::Surface {
35                    1.4
36                } else {
37                    1.1
38                }
39            }
40            MeshingProfile::AdaptiveRefine => {
41                if mesh.kind == MeshKind::Surface {
42                    1.8
43                } else {
44                    1.35
45                }
46            }
47        };
48        let proposed = ((mesh.element_count as f64) * profile_scale).round() as u64;
49        let min_refined_elements = match options.profile {
50            MeshingProfile::AdaptiveRefine => mesh.element_count.max(64),
51            MeshingProfile::AnalysisReady | MeshingProfile::SurfaceOnly => 1,
52        };
53        let element_count = proposed
54            .max(min_refined_elements)
55            .min(per_mesh_budget.max(mesh.element_count));
56        let node_count = (mesh.vertex_count.max(3)).max(element_count / 2 + 2);
57        let connectivity_class = match mesh.kind {
58            MeshKind::Surface => {
59                if element_count > 20_000 {
60                    MeshConnectivityClass::SparseBand
61                } else {
62                    MeshConnectivityClass::SurfacePatch
63                }
64            }
65            MeshKind::Volume => MeshConnectivityClass::VolumeCore,
66        };
67        let element_family_hint = match mesh.kind {
68            MeshKind::Surface => {
69                if element_count % 2 == 0 {
70                    ElementFamilyHint::Triangle
71                } else {
72                    ElementFamilyHint::Quad
73                }
74            }
75            MeshKind::Volume => {
76                if element_count % 2 == 0 {
77                    ElementFamilyHint::Tetrahedron
78                } else {
79                    ElementFamilyHint::Hex
80                }
81            }
82        };
83        let coordinate_span_m = mesh_coordinate_span_m(geometry, &mesh.mesh_id);
84        let coordinate_active_dimension_count =
85            coordinate_active_dimension_count(coordinate_span_m);
86        let coordinate_characteristic_length_m = coordinate_characteristic_length_m(
87            coordinate_span_m,
88            coordinate_active_dimension_count,
89            node_count,
90        );
91        let element_geometry = mesh_element_geometry_metrics(geometry, &mesh.mesh_id);
92        let region_span_hint = (geometry.regions.len().max(1) as u32)
93            .clamp(1, 64)
94            .saturating_sub((prepared_meshes.len() as u32) % 2);
95        prepared_meshes.push(PreparedMeshDescriptor {
96            prepared_mesh_id: format!("prep_{}_{}", geometry.revision, mesh.mesh_id),
97            source_mesh_id: mesh.mesh_id,
98            kind: mesh.kind,
99            node_count,
100            element_count,
101            connectivity_class,
102            element_family_hint,
103            region_span_hint,
104            coordinate_span_m,
105            coordinate_active_dimension_count,
106            coordinate_characteristic_length_m,
107            element_geometry_node_count: element_geometry.node_count,
108            element_geometry_edge_count: element_geometry.edge_count,
109            mean_element_edge_length_m: element_geometry.mean_edge_length_m,
110            mean_element_area_m2: element_geometry.mean_area_m2,
111            element_geometry_coverage_ratio: element_geometry.coverage_ratio,
112            reference_element_coordinates_m: element_geometry.reference_coordinates_m,
113            reference_element_area_m2: element_geometry.reference_area_m2,
114            control_volume_cell_count: element_geometry.control_volume_cell_count,
115            control_volume_face_count: element_geometry.control_volume_face_count,
116            control_volume_internal_face_count: element_geometry.control_volume_internal_face_count,
117            control_volume_boundary_face_count: element_geometry.control_volume_boundary_face_count,
118            control_volume_connectivity_coverage_ratio: element_geometry.coverage_ratio,
119            element_topology_sample_element_count: element_geometry
120                .element_topology_sample
121                .element_count,
122            element_topology_sample_edge_count: element_geometry.element_topology_sample.edge_count,
123            element_topology_sample_edge_nodes: element_geometry.element_topology_sample.edge_nodes,
124            element_topology_sample_node_coordinates_m: element_geometry
125                .element_topology_sample
126                .node_coordinates_m,
127            element_topology_sample_element_edges: element_geometry
128                .element_topology_sample
129                .element_edges,
130            element_topology_sample_element_orientations: element_geometry
131                .element_topology_sample
132                .element_orientations,
133            element_topology_sample_element_areas_m2: element_geometry
134                .element_topology_sample
135                .element_areas_m2,
136            element_topology_node_coordinates_m: element_geometry
137                .element_topology_node_coordinates_m,
138            element_topology_edge_nodes: element_geometry.element_topology_edge_nodes,
139            element_topology_element_edges: element_geometry.element_topology_element_edges,
140            element_topology_element_orientations: element_geometry
141                .element_topology_element_orientations,
142            element_topology_element_areas_m2: element_geometry.element_topology_element_areas_m2,
143        });
144    }
145
146    let mut prepared_by_source = BTreeMap::<String, String>::new();
147    for prepared in &prepared_meshes {
148        prepared_by_source.insert(
149            prepared.source_mesh_id.clone(),
150            prepared.prepared_mesh_id.clone(),
151        );
152    }
153
154    let mut source_mesh_ids_by_region = BTreeMap::<String, Vec<String>>::new();
155    for mapping in &geometry.region_entity_mappings {
156        let entry = source_mesh_ids_by_region
157            .entry(mapping.region_id.clone())
158            .or_default();
159        if !entry.iter().any(|mesh_id| mesh_id == &mapping.mesh_id) {
160            entry.push(mapping.mesh_id.clone());
161        }
162    }
163    for mesh_ids in source_mesh_ids_by_region.values_mut() {
164        mesh_ids.sort();
165    }
166
167    let mut region_mappings = Vec::<RegionMeshMapping>::new();
168    for region in &geometry.regions {
169        let source_mesh_ids = source_mesh_ids_by_region
170            .get(&region.region_id)
171            .cloned()
172            .filter(|mesh_ids| !mesh_ids.is_empty())
173            .ok_or_else(|| {
174                format!(
175                    "region {} has no mesh entity mapping for analysis prep",
176                    region.region_id
177                )
178            })?;
179        let prepared_mesh_ids = source_mesh_ids
180            .iter()
181            .filter_map(|mesh_id| prepared_by_source.get(mesh_id).cloned())
182            .collect::<Vec<_>>();
183        if prepared_mesh_ids.len() != source_mesh_ids.len() {
184            let unknown_mesh_ids = source_mesh_ids
185                .iter()
186                .filter(|mesh_id| !prepared_by_source.contains_key(*mesh_id))
187                .cloned()
188                .collect::<Vec<_>>();
189            return Err(format!(
190                "region {} references unknown mesh entity mapping(s): {}",
191                region.region_id,
192                unknown_mesh_ids.join(", ")
193            ));
194        }
195        region_mappings.push(RegionMeshMapping {
196            region_id: region.region_id.clone(),
197            source_mesh_ids,
198            prepared_mesh_ids,
199        });
200    }
201    if region_mappings.is_empty() {
202        let source_mesh_ids = prepared_meshes
203            .iter()
204            .map(|mesh| mesh.source_mesh_id.clone())
205            .collect::<Vec<_>>();
206        let prepared_mesh_ids = prepared_meshes
207            .iter()
208            .map(|mesh| mesh.prepared_mesh_id.clone())
209            .collect::<Vec<_>>();
210        region_mappings.push(RegionMeshMapping {
211            region_id: "region_default".to_string(),
212            source_mesh_ids: source_mesh_ids.clone(),
213            prepared_mesh_ids: prepared_mesh_ids.clone(),
214        });
215    }
216    region_mappings.sort_by(|a, b| a.region_id.cmp(&b.region_id));
217
218    let total_elements = prepared_meshes
219        .iter()
220        .map(|mesh| mesh.element_count)
221        .sum::<u64>()
222        .max(1);
223    let total_nodes = prepared_meshes
224        .iter()
225        .map(|mesh| mesh.node_count)
226        .sum::<u64>()
227        .max(1);
228    let element_density = total_elements as f64 / total_nodes as f64;
229    let normalized_density = element_density.min(1.0);
230    let (min_scaled_jacobian, mean_aspect_ratio) = match options.profile {
231        MeshingProfile::SurfaceOnly => (
232            (0.89 - 0.1 * normalized_density).max(0.5),
233            1.4 + normalized_density,
234        ),
235        MeshingProfile::AnalysisReady => (
236            (0.92 - 0.1 * normalized_density).max(0.5),
237            1.2 + normalized_density,
238        ),
239        MeshingProfile::AdaptiveRefine => (
240            (0.95 - 0.03 * normalized_density).clamp(0.5, 0.99),
241            1.05 + 0.15 * normalized_density,
242        ),
243    };
244
245    Ok(MeshingPrepResult {
246        schema_version: "geometry-prep-for-analysis/v1".to_string(),
247        prepared_meshes,
248        region_mappings,
249        quality: MeshingQualityReport {
250            min_scaled_jacobian,
251            mean_aspect_ratio,
252            inverted_element_count: 0,
253        },
254        provenance: MeshingProvenance {
255            algorithm: "deterministic_topology_seed/v1".to_string(),
256            profile: options.profile,
257            source_geometry_id: geometry.geometry_id.clone(),
258            source_geometry_revision: geometry.revision,
259        },
260    })
261}
262
263#[cfg(test)]
264mod tests;