Skip to main content

runmat_meshing_surface/validate/
checks.rs

1use std::collections::{BTreeMap, BTreeSet};
2
3use crate::{
4    math::{dot, norm, point_triangle_distance, triangle_area, MeshingTolerance},
5    SurfaceDiscretization, SurfaceElement, INTERNAL_SOURCE_EDGE_ID,
6};
7use runmat_meshing_cad::{CadFace, CadTopologyModel, SourceTopologyEdge, SourceTopologyModel};
8
9use super::{
10    geometry::{surface_element_points, topology_face_points, unit_normal},
11    source_edges::{
12        count_closed_source_edge_loops, source_edge_is_recovered_by_chain, surface_edge_source_ids,
13    },
14    types::{SurfaceValidationError, SurfaceValidationOptions, SurfaceValidationReport},
15};
16
17pub fn validate_surface_discretization(
18    topology: &SourceTopologyModel,
19    surface: &SurfaceDiscretization,
20    options: SurfaceValidationOptions,
21) -> Result<SurfaceValidationReport, SurfaceValidationError> {
22    validate_options(options)?;
23    if surface.elements.is_empty() {
24        return Err(SurfaceValidationError::EmptySurface);
25    }
26
27    let tolerance = MeshingTolerance::from_bounds(topology.bounds_min_m, topology.bounds_max_m);
28    let source_faces = topology
29        .faces
30        .iter()
31        .map(|face| (face.face_id, face))
32        .collect::<BTreeMap<_, _>>();
33    let source_edges = topology
34        .edges
35        .iter()
36        .map(|edge| (edge.edge_id, edge))
37        .collect::<BTreeMap<_, _>>();
38    let surface_edges = surface_edge_source_ids(surface);
39
40    let mut covered_source_faces = BTreeSet::<u32>::new();
41    let mut conforming_source_edges = BTreeSet::<u32>::new();
42    let mut max_projection_error_m = 0.0_f64;
43    let mut min_orientation_alignment = f64::INFINITY;
44
45    for element in &surface.elements {
46        let source_face = source_faces.get(&element.source_face_id).ok_or(
47            SurfaceValidationError::MissingSourceFace {
48                source_face_id: element.source_face_id,
49            },
50        )?;
51        let points = surface_element_points(surface, element)?;
52        let computed_area_m2 = triangle_area(points);
53        if computed_area_m2 <= tolerance.length_epsilon(1.0).powi(2) {
54            return Err(SurfaceValidationError::DegenerateElement {
55                element_id: element.element_id,
56            });
57        }
58        validate_element_parametric_evidence(element)?;
59        validate_element_area_evidence(element, computed_area_m2, tolerance)?;
60        validate_element_projection_evidence(element)?;
61        let source_points = topology_face_points(topology, source_face.node_ids)?;
62        let projection_error_m = points
63            .into_iter()
64            .map(|point| point_triangle_distance(point, source_points))
65            .fold(0.0_f64, f64::max);
66        max_projection_error_m = max_projection_error_m.max(projection_error_m);
67        if projection_error_m > options.max_projection_error_m.max(tolerance.absolute_m) {
68            return Err(SurfaceValidationError::ProjectionError {
69                element_id: element.element_id,
70                error_m: projection_error_m,
71                max_error_m: options.max_projection_error_m,
72            });
73        }
74
75        let surface_normal =
76            unit_normal(points).ok_or(SurfaceValidationError::DegenerateElement {
77                element_id: element.element_id,
78            })?;
79        let element_normal = validate_element_unit_normal(element)?;
80        let alignment = dot(surface_normal, element_normal);
81        min_orientation_alignment = min_orientation_alignment.min(alignment);
82        if alignment < options.min_orientation_alignment {
83            return Err(SurfaceValidationError::OrientationMismatch {
84                element_id: element.element_id,
85                source_face_id: element.source_face_id,
86                alignment,
87                min_alignment: options.min_orientation_alignment,
88            });
89        }
90
91        covered_source_faces.insert(element.source_face_id);
92        for source_edge_id in element.source_edge_ids {
93            if source_edge_id == INTERNAL_SOURCE_EDGE_ID {
94                continue;
95            }
96            source_edges
97                .get(&source_edge_id)
98                .ok_or(SurfaceValidationError::MissingSourceEdge { source_edge_id })?;
99        }
100    }
101
102    for (source_edge_id, source_edge) in &source_edges {
103        if source_edge_is_recovered_by_chain(
104            surface_edges
105                .get(source_edge_id)
106                .map(Vec::as_slice)
107                .unwrap_or(&[]),
108            source_edge,
109        ) {
110            conforming_source_edges.insert(*source_edge_id);
111        }
112    }
113
114    if options.require_source_edge_conformity {
115        for edge in &topology.edges {
116            if !conforming_source_edges.contains(&edge.edge_id) {
117                return Err(SurfaceValidationError::EdgeConformityFailed {
118                    source_edge_id: edge.edge_id,
119                    source_edge_node_ids: edge.node_ids,
120                    recovered_segment_count: surface_edges
121                        .get(&edge.edge_id)
122                        .map(Vec::len)
123                        .unwrap_or_default(),
124                });
125            }
126        }
127    }
128
129    for source_face in source_faces.keys() {
130        if !covered_source_faces.contains(source_face) {
131            return Err(SurfaceValidationError::UncoveredSourceFace {
132                source_face_id: *source_face,
133            });
134        }
135    }
136
137    let (source_edge_loop_count, closed_source_edge_loop_count) =
138        count_closed_source_edge_loops(&topology.edges)?;
139
140    Ok(SurfaceValidationReport {
141        source_face_count: topology.faces.len(),
142        surface_element_count: surface.elements.len(),
143        source_edge_loop_count,
144        closed_source_edge_loop_count,
145        conforming_source_edge_count: conforming_source_edges.len(),
146        missing_source_edge_count: topology
147            .edges
148            .len()
149            .saturating_sub(conforming_source_edges.len()),
150        max_projection_error_m,
151        min_orientation_alignment: if min_orientation_alignment.is_finite() {
152            min_orientation_alignment
153        } else {
154            1.0
155        },
156        face_coverage_ratio: covered_source_faces.len() as f64 / topology.faces.len() as f64,
157    })
158}
159
160pub fn validate_cad_topology_surface_discretization(
161    cad_topology: &CadTopologyModel,
162    topology: &SourceTopologyModel,
163    surface: &SurfaceDiscretization,
164    options: SurfaceValidationOptions,
165) -> Result<SurfaceValidationReport, SurfaceValidationError> {
166    validate_options(options)?;
167    if surface.elements.is_empty() {
168        return Err(SurfaceValidationError::EmptySurface);
169    }
170
171    let tolerance = MeshingTolerance::from_bounds(topology.bounds_min_m, topology.bounds_max_m);
172    let source_edges = topology
173        .edges
174        .iter()
175        .map(|edge| (edge.edge_id, edge))
176        .collect::<BTreeMap<_, _>>();
177    let cad_faces = cad_topology
178        .faces
179        .iter()
180        .map(|face| (face.entity_id.id.clone(), face))
181        .collect::<BTreeMap<_, _>>();
182    let surface_edges = surface_edge_source_ids(surface);
183
184    let mut cad_face_loop_edges = BTreeMap::<String, BTreeSet<u32>>::new();
185    let mut all_loop_source_edges = BTreeSet::<u32>::new();
186    for cad_face in &cad_topology.faces {
187        let loop_edge_ids = cad_face_loop_source_edge_ids(cad_face)?;
188        all_loop_source_edges.extend(loop_edge_ids.iter().copied());
189        cad_face_loop_edges.insert(cad_face.entity_id.id.clone(), loop_edge_ids);
190    }
191
192    let mut covered_cad_faces = BTreeSet::<String>::new();
193    let mut conforming_source_edges = BTreeSet::<u32>::new();
194    let mut max_projection_error_m = 0.0_f64;
195    let mut min_orientation_alignment = f64::INFINITY;
196
197    for element in &surface.elements {
198        let cad_face_id =
199            element
200                .cad_face_id
201                .as_ref()
202                .ok_or(SurfaceValidationError::MissingCadFaceId {
203                    element_id: element.element_id,
204                })?;
205        let _cad_face =
206            cad_faces
207                .get(cad_face_id)
208                .ok_or_else(|| SurfaceValidationError::MissingCadFace {
209                    cad_face_id: cad_face_id.clone(),
210                })?;
211        let loop_edges = cad_face_loop_edges.get(cad_face_id).ok_or_else(|| {
212            SurfaceValidationError::MissingCadFace {
213                cad_face_id: cad_face_id.clone(),
214            }
215        })?;
216
217        let points = surface_element_points(surface, element)?;
218        let computed_area_m2 = triangle_area(points);
219        if computed_area_m2 <= tolerance.length_epsilon(1.0).powi(2) {
220            return Err(SurfaceValidationError::DegenerateElement {
221                element_id: element.element_id,
222            });
223        }
224        validate_element_parametric_evidence(element)?;
225        validate_element_area_evidence(element, computed_area_m2, tolerance)?;
226        validate_element_projection_evidence(element)?;
227        max_projection_error_m = max_projection_error_m.max(element.max_projection_error_m);
228        if element.max_projection_error_m > options.max_projection_error_m.max(tolerance.absolute_m)
229        {
230            return Err(SurfaceValidationError::ProjectionError {
231                element_id: element.element_id,
232                error_m: element.max_projection_error_m,
233                max_error_m: options.max_projection_error_m,
234            });
235        }
236
237        let surface_normal =
238            unit_normal(points).ok_or(SurfaceValidationError::DegenerateElement {
239                element_id: element.element_id,
240            })?;
241        let element_normal = validate_element_unit_normal(element)?;
242        let alignment = dot(surface_normal, element_normal);
243        min_orientation_alignment = min_orientation_alignment.min(alignment);
244        if alignment < options.min_orientation_alignment {
245            return Err(SurfaceValidationError::OrientationMismatch {
246                element_id: element.element_id,
247                source_face_id: element.source_face_id,
248                alignment,
249                min_alignment: options.min_orientation_alignment,
250            });
251        }
252
253        covered_cad_faces.insert(cad_face_id.clone());
254        for source_edge_id in element.source_edge_ids {
255            if source_edge_id == INTERNAL_SOURCE_EDGE_ID {
256                continue;
257            }
258            source_edges
259                .get(&source_edge_id)
260                .ok_or(SurfaceValidationError::MissingSourceEdge { source_edge_id })?;
261            if !loop_edges.contains(&source_edge_id) {
262                return Err(SurfaceValidationError::SourceEdgeOutsideCadFace {
263                    cad_face_id: cad_face_id.clone(),
264                    source_edge_id,
265                });
266            }
267        }
268    }
269
270    for source_edge_id in &all_loop_source_edges {
271        let source_edge =
272            source_edges
273                .get(source_edge_id)
274                .ok_or(SurfaceValidationError::MissingSourceEdge {
275                    source_edge_id: *source_edge_id,
276                })?;
277        if source_edge_is_recovered_by_chain(
278            surface_edges
279                .get(source_edge_id)
280                .map(Vec::as_slice)
281                .unwrap_or(&[]),
282            source_edge,
283        ) {
284            conforming_source_edges.insert(*source_edge_id);
285        }
286    }
287
288    if options.require_source_edge_conformity {
289        for source_edge_id in &all_loop_source_edges {
290            if !conforming_source_edges.contains(source_edge_id) {
291                let source_edge = source_edges.get(source_edge_id).ok_or(
292                    SurfaceValidationError::MissingSourceEdge {
293                        source_edge_id: *source_edge_id,
294                    },
295                )?;
296                return Err(SurfaceValidationError::EdgeConformityFailed {
297                    source_edge_id: *source_edge_id,
298                    source_edge_node_ids: source_edge.node_ids,
299                    recovered_segment_count: surface_edges
300                        .get(source_edge_id)
301                        .map(Vec::len)
302                        .unwrap_or_default(),
303                });
304            }
305        }
306    }
307
308    for cad_face_id in cad_faces.keys() {
309        if !covered_cad_faces.contains(cad_face_id) {
310            return Err(SurfaceValidationError::UncoveredCadFace {
311                cad_face_id: cad_face_id.clone(),
312            });
313        }
314    }
315
316    let loop_source_edges = all_loop_source_edges
317        .iter()
318        .map(|source_edge_id| {
319            source_edges
320                .get(source_edge_id)
321                .map(|edge| (*edge).clone())
322                .ok_or(SurfaceValidationError::MissingSourceEdge {
323                    source_edge_id: *source_edge_id,
324                })
325        })
326        .collect::<Result<Vec<SourceTopologyEdge>, SurfaceValidationError>>()?;
327    let (source_edge_loop_count, closed_source_edge_loop_count) =
328        count_closed_source_edge_loops(&loop_source_edges)?;
329
330    Ok(SurfaceValidationReport {
331        source_face_count: cad_topology.faces.len(),
332        surface_element_count: surface.elements.len(),
333        source_edge_loop_count,
334        closed_source_edge_loop_count,
335        conforming_source_edge_count: conforming_source_edges.len(),
336        missing_source_edge_count: all_loop_source_edges
337            .len()
338            .saturating_sub(conforming_source_edges.len()),
339        max_projection_error_m,
340        min_orientation_alignment: if min_orientation_alignment.is_finite() {
341            min_orientation_alignment
342        } else {
343            1.0
344        },
345        face_coverage_ratio: covered_cad_faces.len() as f64 / cad_topology.faces.len() as f64,
346    })
347}
348
349fn cad_face_loop_source_edge_ids(
350    cad_face: &CadFace,
351) -> Result<BTreeSet<u32>, SurfaceValidationError> {
352    cad_face
353        .loop_edge_ids
354        .iter()
355        .map(|loop_edge_id| {
356            loop_edge_id
357                .strip_prefix("cad_edge_")
358                .and_then(|edge_id| edge_id.parse::<u32>().ok())
359                .ok_or_else(|| SurfaceValidationError::InvalidCadLoopEdgeId {
360                    cad_face_id: cad_face.entity_id.id.clone(),
361                    loop_edge_id: loop_edge_id.clone(),
362                })
363        })
364        .collect()
365}
366
367fn validate_options(options: SurfaceValidationOptions) -> Result<(), SurfaceValidationError> {
368    if !options.max_projection_error_m.is_finite()
369        || options.max_projection_error_m < 0.0
370        || !options.min_orientation_alignment.is_finite()
371        || !(0.0..=1.0).contains(&options.min_orientation_alignment)
372    {
373        return Err(SurfaceValidationError::InvalidOptions);
374    }
375    Ok(())
376}
377
378fn validate_element_unit_normal(
379    element: &SurfaceElement,
380) -> Result<[f64; 3], SurfaceValidationError> {
381    if !element
382        .unit_normal
383        .iter()
384        .all(|component| component.is_finite())
385    {
386        return Err(SurfaceValidationError::InvalidElementNormal {
387            element_id: element.element_id,
388        });
389    }
390    let normal_length = norm(element.unit_normal);
391    if (normal_length - 1.0).abs() > 1.0e-8 {
392        return Err(SurfaceValidationError::InvalidElementNormal {
393            element_id: element.element_id,
394        });
395    }
396    Ok(element.unit_normal)
397}
398
399fn validate_element_area_evidence(
400    element: &SurfaceElement,
401    computed_area_m2: f64,
402    tolerance: MeshingTolerance,
403) -> Result<(), SurfaceValidationError> {
404    if !element.area_m2.is_finite() || element.area_m2 <= tolerance.length_epsilon(1.0).powi(2) {
405        return Err(SurfaceValidationError::InvalidElementArea {
406            element_id: element.element_id,
407        });
408    }
409    let area_tolerance_m2 = tolerance
410        .length_epsilon(1.0)
411        .powi(2)
412        .max(computed_area_m2.abs() * 1.0e-8);
413    if (element.area_m2 - computed_area_m2).abs() > area_tolerance_m2 {
414        return Err(SurfaceValidationError::InvalidElementArea {
415            element_id: element.element_id,
416        });
417    }
418    Ok(())
419}
420
421fn validate_element_parametric_evidence(
422    element: &SurfaceElement,
423) -> Result<(), SurfaceValidationError> {
424    if !element
425        .parametric_node_uv
426        .iter()
427        .flatten()
428        .all(|coordinate| coordinate.is_finite())
429    {
430        return Err(SurfaceValidationError::InvalidParametricEvidence {
431            element_id: element.element_id,
432        });
433    }
434    Ok(())
435}
436
437fn validate_element_projection_evidence(
438    element: &SurfaceElement,
439) -> Result<(), SurfaceValidationError> {
440    if !element.max_projection_error_m.is_finite() || element.max_projection_error_m < 0.0 {
441        return Err(SurfaceValidationError::InvalidProjectionEvidence {
442            element_id: element.element_id,
443        });
444    }
445    Ok(())
446}