1use runmat_geometry_core::{GeometryAsset, MeshKind};
2use serde::{Deserialize, Serialize};
3use std::collections::{BTreeMap, BTreeSet};
4
5#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
6#[serde(rename_all = "snake_case")]
7pub enum MeshingProfile {
8 SurfaceOnly,
9 AnalysisReady,
10 AdaptiveRefine,
11}
12
13#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
14pub struct MeshingOptions {
15 pub profile: MeshingProfile,
16 pub target_element_budget: usize,
17}
18
19impl Default for MeshingOptions {
20 fn default() -> Self {
21 Self {
22 profile: MeshingProfile::AnalysisReady,
23 target_element_budget: 250_000,
24 }
25 }
26}
27
28#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
29#[serde(rename_all = "snake_case")]
30pub enum MeshConnectivityClass {
31 SparseBand,
32 SurfacePatch,
33 VolumeCore,
34}
35
36#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
37#[serde(rename_all = "snake_case")]
38pub enum ElementFamilyHint {
39 Triangle,
40 Quad,
41 Tet,
42 Hex,
43 Mixed,
44}
45
46fn default_coordinate_span_m() -> [f64; 3] {
47 [1.0, 0.0, 0.0]
48}
49
50fn default_coordinate_active_dimension_count() -> u8 {
51 1
52}
53
54fn default_coordinate_characteristic_length_m() -> f64 {
55 1.0
56}
57
58fn default_zero_u64() -> u64 {
59 0
60}
61
62fn default_zero_f64() -> f64 {
63 0.0
64}
65
66fn default_reference_element_coordinates_m() -> [[f64; 3]; 3] {
67 [[0.0; 3]; 3]
68}
69
70fn default_element_topology_sample_edge_nodes() -> [[u32; 2]; 8] {
71 [[0; 2]; 8]
72}
73
74fn default_element_topology_sample_node_coordinates_m() -> [[f64; 3]; 8] {
75 [[0.0; 3]; 8]
76}
77
78fn default_element_topology_sample_element_edges() -> [[u32; 3]; 4] {
79 [[0; 3]; 4]
80}
81
82fn default_element_topology_sample_element_orientations() -> [[i8; 3]; 4] {
83 [[0; 3]; 4]
84}
85
86fn default_element_topology_sample_element_areas_m2() -> [f64; 4] {
87 [0.0; 4]
88}
89
90fn default_element_topology_node_coordinates_m() -> Vec<[f64; 3]> {
91 Vec::new()
92}
93
94fn default_element_topology_edge_nodes() -> Vec<[u32; 2]> {
95 Vec::new()
96}
97
98fn default_element_topology_element_edges() -> Vec<[u32; 3]> {
99 Vec::new()
100}
101
102fn default_element_topology_element_orientations() -> Vec<[i8; 3]> {
103 Vec::new()
104}
105
106fn default_element_topology_element_areas_m2() -> Vec<f64> {
107 Vec::new()
108}
109
110#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
111pub struct PreparedMeshDescriptor {
112 pub prepared_mesh_id: String,
113 pub source_mesh_id: String,
114 pub kind: MeshKind,
115 pub node_count: u64,
116 pub element_count: u64,
117 pub connectivity_class: MeshConnectivityClass,
118 pub element_family_hint: ElementFamilyHint,
119 pub region_span_hint: u32,
120 #[serde(default = "default_coordinate_span_m")]
121 pub coordinate_span_m: [f64; 3],
122 #[serde(default = "default_coordinate_active_dimension_count")]
123 pub coordinate_active_dimension_count: u8,
124 #[serde(default = "default_coordinate_characteristic_length_m")]
125 pub coordinate_characteristic_length_m: f64,
126 #[serde(default = "default_zero_u64")]
127 pub element_geometry_node_count: u64,
128 #[serde(default = "default_zero_u64")]
129 pub element_geometry_edge_count: u64,
130 #[serde(default = "default_zero_f64")]
131 pub mean_element_edge_length_m: f64,
132 #[serde(default = "default_zero_f64")]
133 pub mean_element_area_m2: f64,
134 #[serde(default = "default_zero_f64")]
135 pub element_geometry_coverage_ratio: f64,
136 #[serde(default = "default_reference_element_coordinates_m")]
137 pub reference_element_coordinates_m: [[f64; 3]; 3],
138 #[serde(default = "default_zero_f64")]
139 pub reference_element_area_m2: f64,
140 #[serde(default = "default_zero_u64")]
141 pub control_volume_cell_count: u64,
142 #[serde(default = "default_zero_u64")]
143 pub control_volume_face_count: u64,
144 #[serde(default = "default_zero_u64")]
145 pub control_volume_internal_face_count: u64,
146 #[serde(default = "default_zero_u64")]
147 pub control_volume_boundary_face_count: u64,
148 #[serde(default = "default_zero_f64")]
149 pub control_volume_connectivity_coverage_ratio: f64,
150 #[serde(default = "default_zero_u64")]
151 pub element_topology_sample_element_count: u64,
152 #[serde(default = "default_zero_u64")]
153 pub element_topology_sample_edge_count: u64,
154 #[serde(default = "default_element_topology_sample_edge_nodes")]
155 pub element_topology_sample_edge_nodes: [[u32; 2]; 8],
156 #[serde(default = "default_element_topology_sample_node_coordinates_m")]
157 pub element_topology_sample_node_coordinates_m: [[f64; 3]; 8],
158 #[serde(default = "default_element_topology_sample_element_edges")]
159 pub element_topology_sample_element_edges: [[u32; 3]; 4],
160 #[serde(default = "default_element_topology_sample_element_orientations")]
161 pub element_topology_sample_element_orientations: [[i8; 3]; 4],
162 #[serde(default = "default_element_topology_sample_element_areas_m2")]
163 pub element_topology_sample_element_areas_m2: [f64; 4],
164 #[serde(default = "default_element_topology_node_coordinates_m")]
165 pub element_topology_node_coordinates_m: Vec<[f64; 3]>,
166 #[serde(default = "default_element_topology_edge_nodes")]
167 pub element_topology_edge_nodes: Vec<[u32; 2]>,
168 #[serde(default = "default_element_topology_element_edges")]
169 pub element_topology_element_edges: Vec<[u32; 3]>,
170 #[serde(default = "default_element_topology_element_orientations")]
171 pub element_topology_element_orientations: Vec<[i8; 3]>,
172 #[serde(default = "default_element_topology_element_areas_m2")]
173 pub element_topology_element_areas_m2: Vec<f64>,
174}
175
176#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
177pub struct RegionMeshMapping {
178 pub region_id: String,
179 pub source_mesh_ids: Vec<String>,
180 pub prepared_mesh_ids: Vec<String>,
181}
182
183#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
184pub struct MeshingQualityReport {
185 pub min_scaled_jacobian: f64,
186 pub mean_aspect_ratio: f64,
187 pub inverted_element_count: u64,
188}
189
190#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
191pub struct MeshingProvenance {
192 pub algorithm: String,
193 pub profile: MeshingProfile,
194 pub source_geometry_id: String,
195 pub source_geometry_revision: u32,
196}
197
198#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
199pub struct MeshingPrepResult {
200 pub schema_version: String,
201 pub prepared_meshes: Vec<PreparedMeshDescriptor>,
202 pub region_mappings: Vec<RegionMeshMapping>,
203 pub quality: MeshingQualityReport,
204 pub provenance: MeshingProvenance,
205}
206
207pub fn prepare_geometry_for_analysis(
208 geometry: &GeometryAsset,
209 options: MeshingOptions,
210) -> Result<MeshingPrepResult, String> {
211 if geometry.meshes.is_empty() {
212 return Err("geometry has no meshes to prepare".to_string());
213 }
214 if options.target_element_budget == 0 {
215 return Err("target_element_budget must be greater than zero".to_string());
216 }
217
218 let mut source_meshes = geometry.meshes.clone();
219 source_meshes.sort_by(|a, b| a.mesh_id.cmp(&b.mesh_id));
220 let per_mesh_budget =
221 (options.target_element_budget / source_meshes.len().max(1)).max(1) as u64;
222
223 let mut prepared_meshes = Vec::with_capacity(source_meshes.len());
224 for mesh in source_meshes {
225 let profile_scale = match options.profile {
226 MeshingProfile::SurfaceOnly => 1.0,
227 MeshingProfile::AnalysisReady => {
228 if mesh.kind == MeshKind::Surface {
229 1.4
230 } else {
231 1.1
232 }
233 }
234 MeshingProfile::AdaptiveRefine => {
235 if mesh.kind == MeshKind::Surface {
236 1.8
237 } else {
238 1.35
239 }
240 }
241 };
242 let proposed = ((mesh.element_count as f64) * profile_scale).round() as u64;
243 let min_refined_elements = match options.profile {
244 MeshingProfile::AdaptiveRefine => mesh.element_count.max(64),
245 MeshingProfile::AnalysisReady | MeshingProfile::SurfaceOnly => 1,
246 };
247 let element_count = proposed
248 .max(min_refined_elements)
249 .min(per_mesh_budget.max(mesh.element_count));
250 let node_count = (mesh.vertex_count.max(3)).max(element_count / 2 + 2);
251 let connectivity_class = match mesh.kind {
252 MeshKind::Surface => {
253 if element_count > 20_000 {
254 MeshConnectivityClass::SparseBand
255 } else {
256 MeshConnectivityClass::SurfacePatch
257 }
258 }
259 MeshKind::Volume => MeshConnectivityClass::VolumeCore,
260 };
261 let element_family_hint = match mesh.kind {
262 MeshKind::Surface => {
263 if element_count % 2 == 0 {
264 ElementFamilyHint::Triangle
265 } else {
266 ElementFamilyHint::Quad
267 }
268 }
269 MeshKind::Volume => {
270 if element_count % 2 == 0 {
271 ElementFamilyHint::Tet
272 } else {
273 ElementFamilyHint::Hex
274 }
275 }
276 };
277 let coordinate_span_m = mesh_coordinate_span_m(geometry, &mesh.mesh_id);
278 let coordinate_active_dimension_count =
279 coordinate_active_dimension_count(coordinate_span_m);
280 let coordinate_characteristic_length_m = coordinate_characteristic_length_m(
281 coordinate_span_m,
282 coordinate_active_dimension_count,
283 node_count,
284 );
285 let element_geometry = mesh_element_geometry_metrics(geometry, &mesh.mesh_id);
286 let region_span_hint = (geometry.regions.len().max(1) as u32)
287 .clamp(1, 64)
288 .saturating_sub((prepared_meshes.len() as u32) % 2);
289 prepared_meshes.push(PreparedMeshDescriptor {
290 prepared_mesh_id: format!("prep_{}_{}", geometry.revision, mesh.mesh_id),
291 source_mesh_id: mesh.mesh_id,
292 kind: mesh.kind,
293 node_count,
294 element_count,
295 connectivity_class,
296 element_family_hint,
297 region_span_hint,
298 coordinate_span_m,
299 coordinate_active_dimension_count,
300 coordinate_characteristic_length_m,
301 element_geometry_node_count: element_geometry.node_count,
302 element_geometry_edge_count: element_geometry.edge_count,
303 mean_element_edge_length_m: element_geometry.mean_edge_length_m,
304 mean_element_area_m2: element_geometry.mean_area_m2,
305 element_geometry_coverage_ratio: element_geometry.coverage_ratio,
306 reference_element_coordinates_m: element_geometry.reference_coordinates_m,
307 reference_element_area_m2: element_geometry.reference_area_m2,
308 control_volume_cell_count: element_geometry.control_volume_cell_count,
309 control_volume_face_count: element_geometry.control_volume_face_count,
310 control_volume_internal_face_count: element_geometry.control_volume_internal_face_count,
311 control_volume_boundary_face_count: element_geometry.control_volume_boundary_face_count,
312 control_volume_connectivity_coverage_ratio: element_geometry.coverage_ratio,
313 element_topology_sample_element_count: element_geometry
314 .element_topology_sample
315 .element_count,
316 element_topology_sample_edge_count: element_geometry.element_topology_sample.edge_count,
317 element_topology_sample_edge_nodes: element_geometry.element_topology_sample.edge_nodes,
318 element_topology_sample_node_coordinates_m: element_geometry
319 .element_topology_sample
320 .node_coordinates_m,
321 element_topology_sample_element_edges: element_geometry
322 .element_topology_sample
323 .element_edges,
324 element_topology_sample_element_orientations: element_geometry
325 .element_topology_sample
326 .element_orientations,
327 element_topology_sample_element_areas_m2: element_geometry
328 .element_topology_sample
329 .element_areas_m2,
330 element_topology_node_coordinates_m: element_geometry
331 .element_topology_node_coordinates_m,
332 element_topology_edge_nodes: element_geometry.element_topology_edge_nodes,
333 element_topology_element_edges: element_geometry.element_topology_element_edges,
334 element_topology_element_orientations: element_geometry
335 .element_topology_element_orientations,
336 element_topology_element_areas_m2: element_geometry.element_topology_element_areas_m2,
337 });
338 }
339
340 let mut prepared_by_source = BTreeMap::<String, String>::new();
341 for prepared in &prepared_meshes {
342 prepared_by_source.insert(
343 prepared.source_mesh_id.clone(),
344 prepared.prepared_mesh_id.clone(),
345 );
346 }
347
348 let mut source_mesh_ids_by_region = BTreeMap::<String, Vec<String>>::new();
349 for mapping in &geometry.region_entity_mappings {
350 let entry = source_mesh_ids_by_region
351 .entry(mapping.region_id.clone())
352 .or_default();
353 if !entry.iter().any(|mesh_id| mesh_id == &mapping.mesh_id) {
354 entry.push(mapping.mesh_id.clone());
355 }
356 }
357 for mesh_ids in source_mesh_ids_by_region.values_mut() {
358 mesh_ids.sort();
359 }
360
361 let fallback_source_mesh_ids = prepared_meshes
362 .iter()
363 .map(|mesh| mesh.source_mesh_id.clone())
364 .collect::<Vec<_>>();
365 let fallback_prepared_mesh_ids = prepared_meshes
366 .iter()
367 .map(|mesh| mesh.prepared_mesh_id.clone())
368 .collect::<Vec<_>>();
369
370 let mut region_mappings = Vec::<RegionMeshMapping>::new();
371 for region in &geometry.regions {
372 let source_mesh_ids = source_mesh_ids_by_region
373 .get(®ion.region_id)
374 .cloned()
375 .filter(|mesh_ids| !mesh_ids.is_empty())
376 .unwrap_or_else(|| fallback_source_mesh_ids.clone());
377 let mut prepared_mesh_ids = source_mesh_ids
378 .iter()
379 .filter_map(|mesh_id| prepared_by_source.get(mesh_id).cloned())
380 .collect::<Vec<_>>();
381 if prepared_mesh_ids.is_empty() {
382 prepared_mesh_ids = fallback_prepared_mesh_ids.clone();
383 }
384 region_mappings.push(RegionMeshMapping {
385 region_id: region.region_id.clone(),
386 source_mesh_ids,
387 prepared_mesh_ids,
388 });
389 }
390 if region_mappings.is_empty() {
391 let source_mesh_ids = prepared_meshes
392 .iter()
393 .map(|mesh| mesh.source_mesh_id.clone())
394 .collect::<Vec<_>>();
395 let prepared_mesh_ids = prepared_meshes
396 .iter()
397 .map(|mesh| mesh.prepared_mesh_id.clone())
398 .collect::<Vec<_>>();
399 region_mappings.push(RegionMeshMapping {
400 region_id: "region_default".to_string(),
401 source_mesh_ids: source_mesh_ids.clone(),
402 prepared_mesh_ids: prepared_mesh_ids.clone(),
403 });
404 }
405 region_mappings.sort_by(|a, b| a.region_id.cmp(&b.region_id));
406
407 let total_elements = prepared_meshes
408 .iter()
409 .map(|mesh| mesh.element_count)
410 .sum::<u64>()
411 .max(1);
412 let total_nodes = prepared_meshes
413 .iter()
414 .map(|mesh| mesh.node_count)
415 .sum::<u64>()
416 .max(1);
417 let element_density = total_elements as f64 / total_nodes as f64;
418 let normalized_density = element_density.min(1.0);
419 let (min_scaled_jacobian, mean_aspect_ratio) = match options.profile {
420 MeshingProfile::SurfaceOnly => (
421 (0.89 - 0.1 * normalized_density).max(0.5),
422 1.4 + normalized_density,
423 ),
424 MeshingProfile::AnalysisReady => (
425 (0.92 - 0.1 * normalized_density).max(0.5),
426 1.2 + normalized_density,
427 ),
428 MeshingProfile::AdaptiveRefine => (
429 (0.95 - 0.03 * normalized_density).clamp(0.5, 0.99),
430 1.05 + 0.15 * normalized_density,
431 ),
432 };
433
434 Ok(MeshingPrepResult {
435 schema_version: "geometry-prep-for-analysis/v1".to_string(),
436 prepared_meshes,
437 region_mappings,
438 quality: MeshingQualityReport {
439 min_scaled_jacobian,
440 mean_aspect_ratio,
441 inverted_element_count: 0,
442 },
443 provenance: MeshingProvenance {
444 algorithm: "deterministic_topology_seed/v1".to_string(),
445 profile: options.profile,
446 source_geometry_id: geometry.geometry_id.clone(),
447 source_geometry_revision: geometry.revision,
448 },
449 })
450}
451
452#[derive(Debug, Clone, Default)]
453struct ElementGeometryMetrics {
454 node_count: u64,
455 edge_count: u64,
456 mean_edge_length_m: f64,
457 mean_area_m2: f64,
458 coverage_ratio: f64,
459 reference_coordinates_m: [[f64; 3]; 3],
460 reference_area_m2: f64,
461 control_volume_cell_count: u64,
462 control_volume_face_count: u64,
463 control_volume_internal_face_count: u64,
464 control_volume_boundary_face_count: u64,
465 element_topology_sample: ElementTopologySample,
466 element_topology_node_coordinates_m: Vec<[f64; 3]>,
467 element_topology_edge_nodes: Vec<[u32; 2]>,
468 element_topology_element_edges: Vec<[u32; 3]>,
469 element_topology_element_orientations: Vec<[i8; 3]>,
470 element_topology_element_areas_m2: Vec<f64>,
471}
472
473#[derive(Debug, Clone, Copy, Default)]
474struct ElementTopologySample {
475 element_count: u64,
476 edge_count: u64,
477 edge_nodes: [[u32; 2]; 8],
478 node_coordinates_m: [[f64; 3]; 8],
479 element_edges: [[u32; 3]; 4],
480 element_orientations: [[i8; 3]; 4],
481 element_areas_m2: [f64; 4],
482}
483
484fn mesh_element_geometry_metrics(
485 geometry: &GeometryAsset,
486 mesh_id: &str,
487) -> ElementGeometryMetrics {
488 let Some(surface) = geometry
489 .surface_meshes
490 .iter()
491 .find(|surface| surface.mesh_id == mesh_id)
492 else {
493 return ElementGeometryMetrics::default();
494 };
495 let Some(descriptor) = geometry.meshes.iter().find(|mesh| mesh.mesh_id == mesh_id) else {
496 return ElementGeometryMetrics::default();
497 };
498
499 let mut referenced_nodes = BTreeSet::<u32>::new();
500 let mut unique_edges = BTreeSet::<(u32, u32)>::new();
501 let mut edge_length_sum = 0.0_f64;
502 let mut edge_length_count = 0_u64;
503 let mut area_sum = 0.0_f64;
504 let mut valid_triangle_count = 0_u64;
505 let mut edge_incidence = BTreeMap::<(u32, u32), u64>::new();
506 let mut edge_indices = BTreeMap::<(u32, u32), u32>::new();
507 let mut element_topology_sample = ElementTopologySample::default();
508 let element_topology_node_coordinates_m = surface.vertices.clone();
509 let mut element_topology_edge_nodes = Vec::<[u32; 2]>::new();
510 let mut element_topology_element_edges = Vec::<[u32; 3]>::new();
511 let mut element_topology_element_orientations = Vec::<[i8; 3]>::new();
512 let mut element_topology_element_areas_m2 = Vec::<f64>::new();
513 let mut reference_coordinates_m = [[0.0_f64; 3]; 3];
514 let mut reference_area_m2 = 0.0_f64;
515 for triangle in &surface.triangles {
516 let indices = [triangle[0], triangle[1], triangle[2]];
517 let Some(vertices) = triangle_vertices(&surface.vertices, indices) else {
518 continue;
519 };
520 valid_triangle_count += 1;
521 let triangle_area = triangle_area_m2(vertices);
522 if reference_area_m2 == 0.0 && triangle_area.is_finite() && triangle_area > 0.0 {
523 reference_coordinates_m = vertices;
524 reference_area_m2 = triangle_area;
525 }
526 for index in indices {
527 referenced_nodes.insert(index);
528 }
529 for (index, vertex) in indices.into_iter().zip(vertices) {
530 if (index as usize) < element_topology_sample.node_coordinates_m.len() {
531 element_topology_sample.node_coordinates_m[index as usize] = vertex;
532 }
533 }
534 for (left, right) in [
535 (indices[0], indices[1]),
536 (indices[1], indices[2]),
537 (indices[2], indices[0]),
538 ] {
539 let edge = (left.min(right), left.max(right));
540 unique_edges.insert(edge);
541 *edge_incidence.entry(edge).or_insert(0) += 1;
542 if !edge_indices.contains_key(&edge) {
543 let edge_index = edge_indices.len() as u32;
544 edge_indices.insert(edge, edge_index);
545 element_topology_edge_nodes.push([edge.0, edge.1]);
546 if (edge_index as usize) < element_topology_sample.edge_nodes.len() {
547 element_topology_sample.edge_nodes[edge_index as usize] = [edge.0, edge.1];
548 element_topology_sample.edge_count = (edge_index as u64 + 1)
549 .min(element_topology_sample.edge_nodes.len() as u64);
550 }
551 }
552 }
553 let mut full_element_edges = [0_u32; 3];
554 let mut full_element_orientations = [0_i8; 3];
555 for (local_index, (left, right)) in [
556 (indices[0], indices[1]),
557 (indices[1], indices[2]),
558 (indices[2], indices[0]),
559 ]
560 .into_iter()
561 .enumerate()
562 {
563 let edge = (left.min(right), left.max(right));
564 full_element_edges[local_index] = *edge_indices.get(&edge).unwrap_or(&0);
565 full_element_orientations[local_index] = if left <= right { 1 } else { -1 };
566 }
567 element_topology_element_edges.push(full_element_edges);
568 element_topology_element_orientations.push(full_element_orientations);
569 element_topology_element_areas_m2.push(triangle_area);
570 if (element_topology_sample.element_count as usize) < 4 {
571 let element_index = element_topology_sample.element_count as usize;
572 for (local_index, (left, right)) in [
573 (indices[0], indices[1]),
574 (indices[1], indices[2]),
575 (indices[2], indices[0]),
576 ]
577 .into_iter()
578 .enumerate()
579 {
580 let edge = (left.min(right), left.max(right));
581 element_topology_sample.element_edges[element_index][local_index] =
582 *edge_indices.get(&edge).unwrap_or(&0);
583 element_topology_sample.element_orientations[element_index][local_index] =
584 if left <= right { 1 } else { -1 };
585 }
586 element_topology_sample.element_areas_m2[element_index] = triangle_area;
587 element_topology_sample.element_count += 1;
588 }
589 for (left, right) in [
590 (vertices[0], vertices[1]),
591 (vertices[1], vertices[2]),
592 (vertices[2], vertices[0]),
593 ] {
594 edge_length_sum += distance_m(left, right);
595 edge_length_count += 1;
596 }
597 area_sum += triangle_area;
598 }
599
600 let control_volume_internal_face_count =
601 edge_incidence.values().filter(|count| **count > 1).count() as u64;
602 let control_volume_boundary_face_count =
603 edge_incidence.values().filter(|count| **count == 1).count() as u64;
604
605 ElementGeometryMetrics {
606 node_count: referenced_nodes.len() as u64,
607 edge_count: unique_edges.len() as u64,
608 mean_edge_length_m: if edge_length_count == 0 {
609 0.0
610 } else {
611 edge_length_sum / edge_length_count as f64
612 },
613 mean_area_m2: if valid_triangle_count == 0 {
614 0.0
615 } else {
616 area_sum / valid_triangle_count as f64
617 },
618 coverage_ratio: if descriptor.element_count == 0 {
619 0.0
620 } else {
621 (valid_triangle_count as f64 / descriptor.element_count as f64).clamp(0.0, 1.0)
622 },
623 reference_coordinates_m,
624 reference_area_m2,
625 control_volume_cell_count: valid_triangle_count,
626 control_volume_face_count: unique_edges.len() as u64,
627 control_volume_internal_face_count,
628 control_volume_boundary_face_count,
629 element_topology_sample,
630 element_topology_node_coordinates_m,
631 element_topology_edge_nodes,
632 element_topology_element_edges,
633 element_topology_element_orientations,
634 element_topology_element_areas_m2,
635 }
636}
637
638fn triangle_vertices(vertices: &[[f64; 3]], indices: [u32; 3]) -> Option<[[f64; 3]; 3]> {
639 let a = *vertices.get(indices[0] as usize)?;
640 let b = *vertices.get(indices[1] as usize)?;
641 let c = *vertices.get(indices[2] as usize)?;
642 Some([a, b, c])
643}
644
645fn distance_m(left: [f64; 3], right: [f64; 3]) -> f64 {
646 ((right[0] - left[0]).powi(2) + (right[1] - left[1]).powi(2) + (right[2] - left[2]).powi(2))
647 .sqrt()
648}
649
650fn triangle_area_m2(vertices: [[f64; 3]; 3]) -> f64 {
651 let ab = [
652 vertices[1][0] - vertices[0][0],
653 vertices[1][1] - vertices[0][1],
654 vertices[1][2] - vertices[0][2],
655 ];
656 let ac = [
657 vertices[2][0] - vertices[0][0],
658 vertices[2][1] - vertices[0][1],
659 vertices[2][2] - vertices[0][2],
660 ];
661 let cross = [
662 ab[1] * ac[2] - ab[2] * ac[1],
663 ab[2] * ac[0] - ab[0] * ac[2],
664 ab[0] * ac[1] - ab[1] * ac[0],
665 ];
666 0.5 * (cross[0] * cross[0] + cross[1] * cross[1] + cross[2] * cross[2]).sqrt()
667}
668
669fn mesh_coordinate_span_m(geometry: &GeometryAsset, mesh_id: &str) -> [f64; 3] {
670 let Some(surface) = geometry
671 .surface_meshes
672 .iter()
673 .find(|surface| surface.mesh_id == mesh_id)
674 else {
675 return default_coordinate_span_m();
676 };
677 let Some(first) = surface.vertices.first().copied() else {
678 return default_coordinate_span_m();
679 };
680 let mut min = first;
681 let mut max = first;
682 for vertex in &surface.vertices {
683 for axis in 0..3 {
684 min[axis] = min[axis].min(vertex[axis]);
685 max[axis] = max[axis].max(vertex[axis]);
686 }
687 }
688 [
689 finite_positive_or_default(max[0] - min[0], 0.0),
690 finite_positive_or_default(max[1] - min[1], 0.0),
691 finite_positive_or_default(max[2] - min[2], 0.0),
692 ]
693}
694
695fn coordinate_active_dimension_count(span_m: [f64; 3]) -> u8 {
696 span_m
697 .iter()
698 .filter(|span| span.is_finite() && **span > 1.0e-12)
699 .count()
700 .max(1) as u8
701}
702
703fn coordinate_characteristic_length_m(
704 span_m: [f64; 3],
705 active_dimension_count: u8,
706 node_count: u64,
707) -> f64 {
708 let active_spans = span_m
709 .into_iter()
710 .filter(|span| span.is_finite() && *span > 1.0e-12)
711 .collect::<Vec<_>>();
712 if active_spans.is_empty() {
713 return default_coordinate_characteristic_length_m();
714 }
715 let domain_measure = active_spans.iter().product::<f64>();
716 let node_scale = (node_count.max(2) as f64).powf(1.0 / active_dimension_count.max(1) as f64);
717 finite_positive_or_default(
718 domain_measure.powf(1.0 / active_dimension_count as f64) / node_scale,
719 1.0,
720 )
721}
722
723fn finite_positive_or_default(value: f64, default: f64) -> f64 {
724 if value.is_finite() && value > 0.0 {
725 value
726 } else {
727 default
728 }
729}
730
731#[cfg(test)]
732mod tests {
733 use super::*;
734 use runmat_geometry_core::{
735 GeometrySource, MaterialEvidence, MeshDescriptor, Region, RegionEntityMapping,
736 SourceGeometry, SourceGeometryKind, SurfaceMesh, TessellationProfile, UnitSystem,
737 };
738
739 fn sample_geometry() -> GeometryAsset {
740 GeometryAsset {
741 geometry_id: "geo_meshing_test".to_string(),
742 source: GeometrySource {
743 path: "/fixtures/cube.stl".to_string(),
744 sha256: "mesh-test".to_string(),
745 importer_version: "stl/v1".to_string(),
746 },
747 source_geometry: SourceGeometry {
748 kind: SourceGeometryKind::Mesh,
749 assembly: None,
750 material_evidence: vec![MaterialEvidence {
751 source_key: "fixture".to_string(),
752 normalized_key: "fixture".to_string(),
753 value: "steel".to_string(),
754 confidence: runmat_geometry_core::MaterialEvidenceConfidence::High,
755 unit_basis: None,
756 assumptions: Vec::new(),
757 }],
758 },
759 tessellation_profile: TessellationProfile::default(),
760 units: UnitSystem::Meter,
761 revision: 7,
762 meshes: vec![MeshDescriptor {
763 mesh_id: "mesh_a".to_string(),
764 kind: MeshKind::Surface,
765 vertex_count: 120,
766 element_count: 200,
767 }],
768 surface_meshes: Vec::new(),
769 regions: vec![Region {
770 region_id: "region_main".to_string(),
771 name: "main".to_string(),
772 tag: None,
773 cad_ownership: None,
774 }],
775 region_entity_mappings: vec![RegionEntityMapping::all_faces(
776 "region_main",
777 "mesh_a",
778 200,
779 )],
780 diagnostics: Vec::new(),
781 }
782 }
783
784 #[test]
785 fn meshing_prep_is_deterministic() {
786 let geometry = sample_geometry();
787 let first = prepare_geometry_for_analysis(&geometry, MeshingOptions::default())
788 .expect("first meshing prep should work");
789 let second = prepare_geometry_for_analysis(&geometry, MeshingOptions::default())
790 .expect("second meshing prep should work");
791 assert_eq!(first, second);
792 let descriptor = first
793 .prepared_meshes
794 .first()
795 .expect("prepared mesh descriptor should exist");
796 assert!(descriptor.region_span_hint >= 1);
797 }
798
799 #[test]
800 fn meshing_prep_carries_element_geometry_metrics() {
801 let mut geometry = sample_geometry();
802 geometry.meshes[0].vertex_count = 4;
803 geometry.meshes[0].element_count = 2;
804 geometry.surface_meshes = vec![SurfaceMesh::new(
805 "mesh_a",
806 vec![
807 [0.0, 0.0, 0.0],
808 [1.0, 0.0, 0.0],
809 [1.0, 1.0, 0.0],
810 [0.0, 1.0, 0.0],
811 ],
812 vec![[0, 1, 2], [0, 2, 3]],
813 )];
814
815 let prep = prepare_geometry_for_analysis(&geometry, MeshingOptions::default())
816 .expect("meshing prep should work");
817 let descriptor = prep
818 .prepared_meshes
819 .first()
820 .expect("prepared mesh descriptor should exist");
821 assert_eq!(descriptor.element_geometry_node_count, 4);
822 assert_eq!(descriptor.element_geometry_edge_count, 5);
823 assert!(descriptor.mean_element_edge_length_m > 1.0);
824 assert!((descriptor.mean_element_area_m2 - 0.5).abs() < 1.0e-12);
825 assert_eq!(descriptor.element_geometry_coverage_ratio, 1.0);
826 assert_eq!(
827 descriptor.reference_element_coordinates_m,
828 [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0]]
829 );
830 assert!((descriptor.reference_element_area_m2 - 0.5).abs() < 1.0e-12);
831 assert_eq!(descriptor.control_volume_cell_count, 2);
832 assert_eq!(descriptor.control_volume_face_count, 5);
833 assert_eq!(descriptor.control_volume_internal_face_count, 1);
834 assert_eq!(descriptor.control_volume_boundary_face_count, 4);
835 assert_eq!(descriptor.control_volume_connectivity_coverage_ratio, 1.0);
836 assert_eq!(descriptor.element_topology_sample_element_count, 2);
837 assert_eq!(descriptor.element_topology_sample_edge_count, 5);
838 assert_eq!(descriptor.element_topology_sample_edge_nodes[0], [0, 1]);
839 assert_eq!(
840 descriptor.element_topology_sample_node_coordinates_m[0],
841 [0.0, 0.0, 0.0]
842 );
843 assert_eq!(
844 descriptor.element_topology_sample_element_edges[0],
845 [0, 1, 2]
846 );
847 assert_eq!(
848 descriptor.element_topology_sample_element_orientations[0],
849 [1, 1, -1]
850 );
851 assert!((descriptor.element_topology_sample_element_areas_m2[0] - 0.5).abs() < 1.0e-12);
852 }
853
854 #[test]
855 fn meshing_prep_validates_options() {
856 let geometry = sample_geometry();
857 let error = prepare_geometry_for_analysis(
858 &geometry,
859 MeshingOptions {
860 profile: MeshingProfile::AnalysisReady,
861 target_element_budget: 0,
862 },
863 )
864 .expect_err("zero budget should fail");
865 assert!(error.contains("target_element_budget"));
866 }
867
868 #[test]
869 fn metadata_only_regions_map_to_prepared_meshes() {
870 let mut geometry = sample_geometry();
871 geometry.region_entity_mappings.clear();
872
873 let prep = prepare_geometry_for_analysis(&geometry, MeshingOptions::default())
874 .expect("meshing prep should work");
875 let mapping = prep
876 .region_mappings
877 .iter()
878 .find(|mapping| mapping.region_id == "region_main")
879 .expect("region mapping should exist");
880 assert_eq!(mapping.source_mesh_ids, vec!["mesh_a"]);
881 assert_eq!(mapping.prepared_mesh_ids, vec!["prep_7_mesh_a"]);
882 }
883
884 #[test]
885 fn adaptive_refine_profile_improves_quality_within_budget() {
886 let geometry = sample_geometry();
887 let analysis_ready = prepare_geometry_for_analysis(
888 &geometry,
889 MeshingOptions {
890 profile: MeshingProfile::AnalysisReady,
891 target_element_budget: 4_000,
892 },
893 )
894 .expect("analysis-ready prep should work");
895 let adaptive = prepare_geometry_for_analysis(
896 &geometry,
897 MeshingOptions {
898 profile: MeshingProfile::AdaptiveRefine,
899 target_element_budget: 4_000,
900 },
901 )
902 .expect("adaptive-refine prep should work");
903
904 let adaptive_total_elements = adaptive
905 .prepared_meshes
906 .iter()
907 .map(|mesh| mesh.element_count)
908 .sum::<u64>();
909 assert!(adaptive_total_elements <= 4_000);
910 assert!(adaptive.quality.min_scaled_jacobian >= analysis_ready.quality.min_scaled_jacobian);
911 assert!(adaptive.quality.mean_aspect_ratio <= analysis_ready.quality.mean_aspect_ratio);
912 }
913}