diff --git a/crates/mesh/src/octree/mod.rs b/crates/mesh/src/octree/mod.rs index 6497781..b294579 100644 --- a/crates/mesh/src/octree/mod.rs +++ b/crates/mesh/src/octree/mod.rs @@ -11,7 +11,8 @@ //! 2. **Grid** ([`grid`]). Root cubes over the part's box, refined as a 27-tree where a B-rep //! face curves (its normals turn by more than [`OctreeOptions::curvature_angle`] inside a //! cell), where the part is thin (fewer than [`OctreeOptions::thickness_cells`] cells across a -//! wall, measured along the inward normal of each boundary triangle to a face that faces it) +//! wall, measured along the inward normal of a boundary triangle to a face that faces it, at +//! the triangle's centre and, on a triangle larger than the cell, where it is nearest the cell) //! and where a patch is narrow, then strongly balanced and closed under the transition //! templates. The levels the thinnest wall and narrowest patch need are found before the grid //! is built, and the root size is chosen with them (see [`octree_hex_mesh`]'s plan): a part @@ -500,17 +501,29 @@ fn choose_frame(g: &Geometry) -> [[f64; 3]; 3] { /// ray that starts near a convex edge leaves through the face beyond the edge after a distance /// that only measures how near the edge it started; the faces of a wall face each other. fn triangle_thickness(g: &Geometry) -> Vec { + (0..g.surface.triangles().len() as u32).map(|t| thickness_at(g, t, g.triangle_centre(t))).collect() +} + +/// The distance through the part from point `p` of boundary triangle `t` (see +/// [`triangle_thickness`]). +fn thickness_at(g: &Geometry, t: u32, p: P) -> f64 { let far = 2.0 * norm(sub(g.bounds.hi, g.bounds.lo)) + 1.0; let opposed = -libm::cos(WALL_OPENING.to_radians()); - (0..g.surface.triangles().len() as u32) - .map(|t| { - let n = g.normals.get(t as usize).copied().unwrap_or([0.0; 3]); - match g.ray(g.triangle_centre(t), n.map(|x| -x), far, t) { - Some((hit, s)) if g.normals.get(hit as usize).is_some_and(|&m| dot(m, n) <= opposed) => s, - _ => f64::INFINITY, - } - }) - .collect() + let n = g.normals.get(t as usize).copied().unwrap_or([0.0; 3]); + match g.ray(p, n.map(|x| -x), far, t) { + Some((hit, s)) if g.normals.get(hit as usize).is_some_and(|&m| dot(m, n) <= opposed) => s, + _ => f64::INFINITY, + } +} + +fn longest_edge(a: P, b: P, c: P) -> f64 { + norm(sub(b, a)).max(norm(sub(c, b))).max(norm(sub(a, c))) +} + +/// `p` on a triangle moved a millionth of the way to its centroid, so a ray from it does not +/// start on an edge it shares with a neighbour. +fn toward_centre(p: P, [a, b, c]: [P; 3]) -> P { + [0, 1, 2].map(|k| p[k] + 1e-6 * ((a[k] + b[k] + c[k]) / 3.0 - p[k])) } fn centroid(points: &[P], hx: &[u32; 8]) -> P { @@ -629,7 +642,8 @@ impl Context<'_> { let mut first: BTreeMap = BTreeMap::new(); for t in self.g.overlapping(&bx) { let [a, b, c] = self.g.triangle(t); - if geom::dist2(centre, geom::closest_on_triangle(centre, a, b, c)) > reach * reach { + let near = geom::closest_on_triangle(centre, a, b, c); + if geom::dist2(centre, near) > reach * reach { continue; } let n = self.g.normals.get(t as usize).copied().unwrap_or([0.0; 3]); @@ -638,8 +652,18 @@ impl Context<'_> { if key.0 < self.options.curvature_levels && dot(f, n) < cos { return Ok(true); } - // (With a relative margin, so cells planned at exactly the size a wall needs stay.) - if self.thickness.get(t as usize).is_some_and(|&th| th * (1.0 + 1e-9) < self.options.thickness_cells * e) { + // The wall at the triangle's centre, and, on a triangle larger than the cell, where it + // is nearest the cell: a planar face beside a pocket is a few triangles whose centres + // may face a far wall while the wall here is thin. That reading is held to the + // thinnest wall the plan chose the levels for (towards a wedge's blade the wall + // thins to nothing). (With a relative margin, so cells planned at exactly the size a + // wall needs stay.) + let need = self.options.thickness_cells * e; + let thin = |th: f64| th * (1.0 + 1e-9) < need; + if self.thickness.get(t as usize).is_some_and(|&th| thin(th)) { + return Ok(true); + } + if longest_edge(a, b, c) > e && thin(thickness_at(self.g, t, toward_centre(near, [a, b, c])).max(self.needs.wall)) { return Ok(true); } // Narrow patches (a hole's wall, a thin rim) need cells across them. diff --git a/crates/mesh/src/octree/tests.rs b/crates/mesh/src/octree/tests.rs index c500a04..140e290 100644 --- a/crates/mesh/src/octree/tests.rs +++ b/crates/mesh/src/octree/tests.rs @@ -182,3 +182,44 @@ fn a_size_box_refines_the_grid_inside_it() { let e = octree_hex_mesh(&s, &deep, &shallow).unwrap_err().to_string(); assert!(e.contains("max_depth"), "{e}"); } + +/// A prism: the polygon `profile` in the x–z plane (counter-clockwise, star-shaped from its +/// first point) extruded from y = 0 to `depth`. Faces 0 and 1 are the caps, 2 + i the side on +/// edge i. +fn prism_surface(profile: &[[f64; 2]], depth: f64) -> BoundarySurface { + let n = profile.len(); + let mut p: Vec

= profile.iter().map(|q| [q[0], 0.0, q[1]]).collect(); + p.extend(profile.iter().map(|q| [q[0], depth, q[1]])); + let (f, b) = (|i: usize| (i % n) as u32, |i: usize| (n + i % n) as u32); + let (mut tris, mut face) = (Vec::new(), Vec::new()); + for i in 1..n - 1 { + tris.extend([[f(0), f(i), f(i + 1)], [b(0), b(i + 1), b(i)]]); + face.extend([0, 1]); + } + for i in 0..n { + tris.extend([[f(i), b(i), b(i + 1)], [f(i), b(i + 1), f(i + 1)]]); + face.extend([2 + i as u32; 2]); + } + BoundarySurface::new(&p, &tris, &face, &crate::surface::SurfaceOptions::default()).unwrap() +} + +#[test] +fn a_wall_is_measured_where_a_large_triangle_meets_the_cell() { + // A 30 mm block with a 3 mm flange: the bottom face is two triangles whose centres lie under + // the block, so their centre rays measure 30, but under the flange the wall is 3 thick. + let s = prism_surface(&[[0.0, 0.0], [60.0, 0.0], [60.0, 3.0], [50.0, 3.0], [50.0, 30.0], [0.0, 30.0]], 40.0); + let g = Geometry::new(&s, 30.0, 160.0).unwrap(); + let thickness = triangle_thickness(&g); + let bottom: Vec = (0..s.triangles().len() as u32).filter(|&t| s.tri_face()[t as usize] == 2).collect(); + assert!(bottom.iter().all(|&t| thickness[t as usize] == 30.0), "{thickness:?}"); + let o = OctreeOptions::default(); + let needs = Needs::of(&g, &thickness, &o); + let identity = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]; + let field = UniformSize(100.0); + let ctx = Context::new(&g, &field, &o, &thickness, identity, 3.0, ([0.0; 3], [60.0, 40.0, 30.0]), needs).unwrap(); + assert_eq!(ctx.origin, [-1.5, -2.5, -1.5]); + // Root cells of 3 on the bottom face, clear of every other face: under the flange + // (x 52.5..55.5, 3 cells across need cells of 1) and under the block (x 28.5..31.5). + assert!(ctx.needs_split((0, [18, 7, 0])).unwrap(), "under the flange"); + assert!(!ctx.needs_split((0, [10, 7, 0])).unwrap(), "under the block"); +}