diff --git a/main/opengeometry/src/query/measure.rs b/main/opengeometry/src/query/measure.rs index a5101c23..9ba83795 100644 --- a/main/opengeometry/src/query/measure.rs +++ b/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), diff --git a/main/opengeometry/src/query/mod.rs b/main/opengeometry/src/query/mod.rs index 88746abc..af82e880 100644 --- a/main/opengeometry/src/query/mod.rs +++ b/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; diff --git a/main/opengeometry/src/query/orientation.rs b/main/opengeometry/src/query/orientation.rs new file mode 100644 index 00000000..e3ae788c --- /dev/null +++ b/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 { + 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, 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; 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) +} diff --git a/main/opengeometry/src/query/tests.rs b/main/opengeometry/src/query/tests.rs index a1330868..58d255dd 100644 --- a/main/opengeometry/src/query/tests.rs +++ b/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(); + } +} diff --git a/main/opengeometry/tests/cases/builders.rs b/main/opengeometry/tests/cases/builders.rs index 9ea68654..19d5fe78 100644 --- a/main/opengeometry/tests/cases/builders.rs +++ b/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 { diff --git a/main/opengeometry/tests/cases/matrix.rs b/main/opengeometry/tests/cases/matrix.rs index d9c3eeec..01a5c433 100644 --- a/main/opengeometry/tests/cases/matrix.rs +++ b/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),