fix(query): validate_orientation catches faces whose material normal points into the solid

main/opengeometry/src/query/measure.rs

@@ -1,3 +1,4 @@
+use super::orientation::validate_orientation_validated;
 use crate::brep::{BodyType, BrepEnvelope, GeometryError};
 
 #[derive(Clone, Debug, PartialEq)]
@@ -9,7 +10,10 @@ pub struct Validity {
 }
 
 pub fn validity(body: &BrepEnvelope) -> Validity {
-    match body.body_type() {
+    let checked = body
+        .body_type()
+        .and_then(|body_type| validate_orientation_validated(body).map(|()| body_type));
+    match checked {
         Ok(body_type) => Validity {
             valid: true,
             closed: closed(body, body_type),

main/opengeometry/src/query/mod.rs

@@ -2,6 +2,7 @@ mod classification;
 mod classify;
 mod face_trim;
 mod measure;
+mod orientation;
 mod ray_roots;
 #[cfg(test)]
 mod tests;
@@ -11,3 +12,4 @@ pub use classify::classify_point;
 pub(super) use classify::classify_point_in_shell;
 pub(super) use face_trim::face_contains_uv;
 pub use measure::{validity, Validity};
+pub use orientation::validate_orientation;

main/opengeometry/src/query/orientation.rs

@@ -0,0 +1,115 @@
+use super::classification::PointClassification;
+use super::classify::classify_point_validated;
+use super::face_trim::face_contains_uv;
+use crate::brep::{unit, BrepEnvelope, Face, GeometryError, Surface, UV};
+use crate::math::{add, norm, scale, sub};
+
+const GRID: usize = 6;
+const PROBE_OFFSET: f64 = 100.0;
+const RELIABLE_TOLERANCE: f64 = 1e-12;
+
+pub fn validate_orientation(brep: &BrepEnvelope) -> Result<(), GeometryError> {
+    brep.validate()?;
+    validate_orientation_validated(brep)
+}
+
+pub(super) fn validate_orientation_validated(brep: &BrepEnvelope) -> Result<(), GeometryError> {
+    if brep.solids.is_empty() {
+        return Ok(());
+    }
+    match check_faces(brep) {
+        Err(GeometryError::Math(_)) => Ok(()),
+        outcome => outcome,
+    }
+}
+
+fn check_faces(brep: &BrepEnvelope) -> Result<(), GeometryError> {
+    let Some(bounds) = brep.bounds_unchecked()? else {
+        return Ok(());
+    };
+    let diagonal = norm(std::array::from_fn(|axis| bounds.axes[axis].width()));
+    let coordinates = bounds
+        .axes
+        .iter()
+        .flat_map(|axis| [axis.lo, axis.hi])
+        .chain(brep.geometry.surfaces.iter().flat_map(|s| s.frame().origin))
+        .fold(diagonal, |scale, value| scale.max(value.abs()));
+    if !coordinates.is_finite() || brep.accuracy.geometric < RELIABLE_TOLERANCE * coordinates {
+        return Ok(());
+    }
+    for face in &brep.topology.faces {
+        if face_is_inverted(brep, face)? {
+            return Err(GeometryError::InvalidTopology(format!(
+                "face {} ('{}') has its material normal pointing into the solid",
+                face.id, face.key
+            )));
+        }
+    }
+    Ok(())
+}
+
+fn face_is_inverted(brep: &BrepEnvelope, face: &Face) -> Result<bool, GeometryError> {
+    let [first, second] = probe_points(brep, face)?;
+    let Some(first) = first else {
+        return Ok(false);
+    };
+    if probe(brep, face, first)? != Some(true) {
+        return Ok(false);
+    }
+    let Some(second) = second else {
+        return Ok(true);
+    };
+    Ok(probe(brep, face, second)? == Some(true))
+}
+
+fn probe(brep: &BrepEnvelope, face: &Face, uv: UV) -> Result<Option<bool>, GeometryError> {
+    let surface = brep.geometry.surface(face.surface)?;
+    let Ok(normal) = surface.normal_at(uv) else {
+        return Ok(None);
+    };
+    let step = scale(
+        unit(normal)?,
+        face.sense.multiplier() * PROBE_OFFSET * brep.accuracy.geometric,
+    );
+    let on_face = surface.point_at(uv)?;
+    let front = classify_point_validated(brep, add(on_face, step))?;
+    let back = classify_point_validated(brep, sub(on_face, step))?;
+    Ok(match (front, back) {
+        (PointClassification::Inside, PointClassification::Outside) => Some(true),
+        (PointClassification::Outside, PointClassification::Inside) => Some(false),
+        _ => None,
+    })
+}
+
+fn probe_points(brep: &BrepEnvelope, face: &Face) -> Result<[Option<UV>; 2], GeometryError> {
+    let cell =
+        |i: usize, j: usize| -> UV { [i, j].map(|index| (index as f64 + 0.5) / GRID as f64) };
+    let at = |fraction: UV| -> UV {
+        std::array::from_fn(|axis| {
+            let bounds = face.trim.uv_bounds[axis];
+            bounds.lo + bounds.width() * fraction[axis]
+        })
+    };
+    let mut inside = [[false; GRID]; GRID];
+    for (i, row) in inside.iter_mut().enumerate() {
+        for (j, value) in row.iter_mut().enumerate() {
+            *value = face_contains_uv(brep, face, at(cell(i, j)))? == Some(true);
+        }
+    }
+    let mut ranked: Vec<(usize, usize, usize)> = (0..GRID * GRID)
+        .map(|index| (index / GRID, index % GRID))
+        .filter(|&(i, j)| inside[i][j])
+        .map(|(i, j)| (depth(&inside, i, j), i, j))
+        .collect();
+    ranked.sort_by_key(|&(depth, i, j)| (std::cmp::Reverse(depth), i, j));
+    Ok([0, 1].map(|rank| ranked.get(rank).map(|&(_, i, j)| at(cell(i, j)))))
+}
+
+fn depth(inside: &[[bool; GRID]; GRID], i: usize, j: usize) -> usize {
+    let border = [i + 1, j + 1, GRID - i, GRID - j].into_iter().min();
+    let outside = (0..GRID * GRID)
+        .map(|index| (index / GRID, index % GRID))
+        .filter(|&(k, l)| !inside[k][l])
+        .map(|(k, l)| i.abs_diff(k).max(j.abs_diff(l)));
+    border.into_iter().chain(outside).min().unwrap_or(0)
+}

main/opengeometry/src/query/tests.rs

@@ -1,6 +1,9 @@
 use super::classification::PointClassification;
 use super::classify::{classify_point, classify_point_validated};
-use crate::brep::{coarse_accuracy, unit, Accuracy, BrepEnvelope, Frame3};
+use super::{validate_orientation, validity};
+use crate::brep::{
+    coarse_accuracy, fine_accuracy, unit, Accuracy, BrepEnvelope, Frame3, GeometryError,
+};
 use crate::math::{add, cross, dot, scale, Point3};
 use crate::primitives;
 
@@ -431,3 +434,75 @@ fn cuboid_corners_edges_and_face_points_are_boundary() {
             .collect(),
     );
 }
+
+fn assert_inverted(body: &BrepEnvelope, face: u32) {
+    body.validate().unwrap();
+    let key = &body.topology.faces[face as usize].key;
+    match validate_orientation(body) {
+        Err(GeometryError::InvalidTopology(message)) => {
+            assert!(message.contains(&format!("'{key}'")), "{message}")
+        }
+        other => panic!("face {face} not reported: {other:?}"),
+    }
+    assert!(!validity(body).valid);
+}
+
+#[test]
+fn a_face_whose_material_normal_points_into_the_solid_is_rejected() {
+    let accuracy = coarse_accuracy(1e-5);
+    let mut cube = primitives::cuboid("cube".into(), Frame3::IDENTITY, [1.0; 3], accuracy).unwrap();
+    validate_orientation(&cube).unwrap();
+    cube.topology.faces[5].sense = cube.topology.faces[5].sense.reversed();
+    assert_inverted(&cube, 5);
+
+    let mut sphere = primitives::sphere("ball".into(), Frame3::IDENTITY, 1.0, accuracy).unwrap();
+    for face in &mut sphere.topology.faces {
+        face.sense = face.sense.reversed();
+    }
+    assert_inverted(&sphere, 0);
+}
+
+#[test]
+fn reversed_bores_thin_plates_and_curved_solids_pass_the_orientation_check() {
+    let accuracy = coarse_accuracy(1e-5);
+    let bodies = [
+        primitives::annular_cylinder("tube".into(), Frame3::IDENTITY, 0.5, 1.0, 2.0, accuracy),
+        primitives::cuboid(
+            "plate".into(),
+            Frame3::IDENTITY,
+            [1000.0, 500.0, 1.0],
+            Accuracy::STANDARD,
+        ),
+        primitives::torus("ring".into(), Frame3::IDENTITY, 2.0, 0.5, accuracy),
+        primitives::cone("cone".into(), Frame3::IDENTITY, 1.0, 2.0, accuracy),
+    ];
+    for body in bodies {
+        let body = body.unwrap();
+        validate_orientation(&body).unwrap();
+        assert!(validity(&body).valid, "{}", body.id);
+    }
+}
+
+#[test]
+fn the_orientation_check_stands_aside_where_classification_is_unreliable() {
+    let near_cylinder = primitives::frustum(
+        "near".into(),
+        Frame3::IDENTITY,
+        1.0,
+        1.0 + 2.0_f64.powi(-48),
+        1.0,
+        fine_accuracy(),
+    )
+    .unwrap();
+    validate_orientation(&near_cylinder).unwrap();
+    for (radius, geometric) in [(1e300, 1e290), (1e-200, 1e-210)] {
+        let accuracy = Accuracy {
+            geometric,
+            intersection: geometric / 4.0,
+            ..fine_accuracy()
+        };
+        let sphere =
+            primitives::sphere("extreme".into(), Frame3::IDENTITY, radius, accuracy).unwrap();
+        validate_orientation(&sphere).unwrap();
+    }
+}

main/opengeometry/tests/cases/builders.rs

@@ -1,5 +1,6 @@
 use opengeometry::brep::{Accuracy, BodyType, CurveGeometry, Frame3};
 use opengeometry::primitives;
+use opengeometry::query::validate_orientation;
 use opengeometry_test_support::stored;
 use serde_json::json;
 
@@ -91,6 +92,7 @@ fn primitive_builders_match_stored_json() {
         ),
     ];
     for (name, body) in bodies {
+        validate_orientation(&body).unwrap();
         assert_eq!(
             body.body_type().unwrap(),
             if name == "circle" {
@@ -184,6 +186,7 @@ fn extrusions_match_stored_json() {
         ),
     ];
     for (name, body) in bodies {
+        validate_orientation(&body).unwrap();
         assert_eq!(body.body_type().unwrap(), BodyType::Solid);
         if stored::compared() {
             let expected = match name {

main/opengeometry/tests/cases/matrix.rs

@@ -3,6 +3,7 @@ use opengeometry::operations::modifying::boolean::{
     boolean_brep_outcome_with_handlers, boolean_brep_with_handlers, BooleanOp,
 };
 use opengeometry::primitives;
+use opengeometry::query::validate_orientation;
 use opengeometry_test_support::{stored, volume};
 use serde_json::json;
 
@@ -56,6 +57,7 @@ fn boolean_fixtures_match_stored_brep_and_handlers() {
         ($name:literal, $a:expr, $b:expr, $operation:expr) => {{
             let (result, handlers) =
                 boolean_brep_with_handlers($a, $b, $operation, $name.into()).unwrap();
+            validate_orientation(&result.brep).unwrap();
             let json = result.brep.to_json().unwrap();
             if stored::compared() {
                 assert_eq!(
@@ -110,7 +112,7 @@ fn boolean_support_matrix_matches_stored_outcomes_and_handlers() {
             boolean_brep_outcome_with_handlers(&a, &b, operation, format!("matrix-{name}"));
         let actual = match result {
             Ok(result) => {
-                result.brep.validate().unwrap();
+                validate_orientation(&result.brep).unwrap();
                 let expected_volume = match name {
                     "box-union" => Some(27.0),
                     "box-intersection" => Some(8.0),