diff --git a/ATTRIBUTION.md b/ATTRIBUTION.md index ed41d30..f383ab0 100644 --- a/ATTRIBUTION.md +++ b/ATTRIBUTION.md @@ -47,8 +47,8 @@ is committed and listed here with its exports. | assets/icons/mesh.element.add.svg, assets/icons/mesh.element.add.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.element.add.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/mesh.hex.decompose.svg, assets/icons/mesh.hex.decompose.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.hex.decompose.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/mesh.hex.dominant.svg, assets/icons/mesh.hex.dominant.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.hex.dominant.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | -| assets/icons/mesh.hex.generate.svg, assets/icons/mesh.hex.generate.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.hex.generate.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/mesh.hex.mapped.svg, assets/icons/mesh.hex.mapped.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.hex.mapped.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | +| assets/icons/mesh.hex.octree.svg, assets/icons/mesh.hex.octree.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.hex.octree.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/mesh.hex.polycube.svg, assets/icons/mesh.hex.polycube.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.hex.polycube.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/mesh.hex.sweep.svg, assets/icons/mesh.hex.sweep.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.hex.sweep.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/mesh.import.orphan.svg, assets/icons/mesh.import.orphan.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.import.orphan.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | diff --git a/ROADMAP.md b/ROADMAP.md index 869756b..6bd860a 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -23,7 +23,7 @@ hex meshes today and it gives the best elements where stresses are read. | M2 Tet meshing and size field | Size field (global, per face and edge seeds, curvature and proximity; the `modelcraft_mesh::size::SizeField` trait and edge-division helpers are in); Bowyer–Watson Delaunay (exact, symbolic perturbation, verified; `modelcraft-mesh::delaunay`, done); robust CDT boundary recovery (Diazzi et al. 2023: implicit LNC Steiner points on segments, gift-wrapping facet recovery, inside/outside labelling, B-rep face tags on boundary tet faces, few-ulp rounding repair checked by the validity gate; `modelcraft-mesh::cdt` and `mesh.tet.generate`, done: all 20 valid corpus parts); Delaunay refinement, ODT smoothing, sliver removal; tet4 and tet10 with mid-side nodes on the B-rep; mesh verify per set | in progress | | M3 First deck | Sets and surfaces bound to SolidCraft face, edge, vertex and body names, resolved after meshing; materials (isotropic and orthotropic elastic, density) with a unit-system label and consistency check; sections (solid, shell, beam profiles) assigned by set; element formulation per region mapped to each solver; BCs (fixed, displacement, symmetry); loads (force and moment on reference points, pressure, gravity); reference points with kinematic and distributing coupling; steps (linear static, frequency) with explicit per-step load and BC tables and output presets; `model.check` naming every missing or bad item; `model.*` and `mesh.tag` commands (model definition: done); writers for Abaqus input (done: `deck.write`), Nastran SOL 101/103 (done: `deck.write`) and LS-DYNA implicit (done); the `.modelcraft` file (text definition, binary mesh, geometry references reopened on load, journal; atomic saves, autosave and recovery; undo/redo on the journal: done) | in progress | | M4 Own solver and first results | Element library per docs/elements.md (1D and 2D first); Linear static and modal solver (tet4/10, wedge6, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions (solver: done — `modelcraft-fem` kernels, `modelcraft-solver` with nested dissection, sparse Cholesky and block Lanczos, `solve.static`/`solve.frequency`/`solve.probe`, patch, cantilever, Lamé, Kirsch and beam-frequency tests; reference points with kinematic and distributing couplings and moments by master–slave constraint elimination, verified by rigid motion, clamping, Saint-Venant tension, bending and tube torsion; 1D elements I: rod2/rod3 trusses, point masses, rotary inertia, springs and dashpots in both solvers and all three decks, verified by the method of joints, the three-bar truss, spring–mass, two-mass chain and torsional pendulum frequencies, and damped decay; 2D elements: plane stress, plane strain and axisymmetric tri3, tri6, quad4 (incompatible modes QM6), quad8 with two degrees of freedom per node, centrifugal loads, `mesh.planar` 2D meshes of polygons and planar faces, verified by patch tests, Kirsch, Lamé (plane strain and axisymmetric), rotating discs and cylinders, and a plane stress cantilever; deck mapping to Abaqus CPS/CPE/CAX, Nastran membranes, LS-DYNA ELFORM 12/13); shells and beams (done: quad4 MITC4 and tri3 MITC3 flat shells with a Hughes–Brezzi drilling rotation, two-node Timoshenko beams with profile properties and orientations, six degrees of freedom at shell and beam nodes shared with solids and trusses, section forces and surface stresses, verified by shell patch tests, Navier's plate, cantilever and portal-frame beams, and the Scordelis–Lo and pinched-cylinder benchmarks); results display (done: `results.*` commands, deformed shape to scale, banded contours of displacement, von Mises, principal and component stresses with a legend, mode shapes with their frequencies and animation, probing a node by clicking it, reaction force and moment sums per set, stale results marked after an edit; explicit steps: field frames by slider or animation with contours per frame, energy histories plotted) | started | -| M5 Structured hex | 2D quad meshing of faces (`modelcraft-mesh::quad`: transfinite mapped meshes, free all-quad meshes of polygons with holes on prescribed boundary nodes, exact planar quad validity, shared `SizeField` and even edge divisions; done for planar faces); mapped hex of six-sided regions; sweep meshing; conforming interfaces and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | started | +| M5 Structured hex | 2D quad meshing of faces (`modelcraft-mesh::quad`: transfinite mapped meshes, free all-quad meshes of polygons with holes on prescribed boundary nodes, exact planar quad validity, shared `SizeField` and even edge divisions; done for planar faces); mapped hex of six-sided regions (trivariate TFI, curved faces) and sweep meshing (multi-level extrusions with imprinted cross-sections, revolutions, one-to-one sweeps through linking faces; interval assignment, hex8/hex20 on the geometry, B-rep face tags, `mesh.hex.generate`/`sweep`/`mapped`: done, 16 of 21 corpus parts all-hex); conforming interfaces and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | started | | M6 Automatic hex | Octree/grid-dual all-hex for any part; hex-dominant by tet recombination; polycube for blocky parts; smoothing, untangling, boundary padding, layers through thin walls; mesh editing and orphan-mesh import | planned | | M7 Assemblies and contact | Part instances from SolidCraft occurrences; tie between non-matching meshes; general and pairwise contact with automatic pair detection, friction; bolt pretension, connectors; amplitudes, multiple steps; nonlinear static steps in decks | planned | | M8 Materials, explicit, shells, composites | Plasticity, hyperelasticity, damage data per solver; explicit dynamics decks (Abaqus/Explicit, LS-DYNA); midsurface shells; composite layups and material orientations, continuum shells on swept hex; hybrid (mixed u–p) hex8/hex20/tet10 for nearly incompressible materials (done); heat transfer (solids done: steady and θ-method transient conduction with flux, convection and sources, thermal strain in static steps (solids and 2D) and explicit dynamics (with plasticity), coupled temperature–displacement, `solve.heat`; deck output, 2D/shell/beam conduction planned) | started | diff --git a/assets/icons/mesh.hex.generate.vectorcraft b/assets/icons/mesh.hex.generate.vectorcraft deleted file mode 100644 index 8625d88..0000000 Binary files a/assets/icons/mesh.hex.generate.vectorcraft and /dev/null differ diff --git a/assets/icons/mesh.hex.generate.svg b/assets/icons/mesh.hex.octree.svg similarity index 94% rename from assets/icons/mesh.hex.generate.svg rename to assets/icons/mesh.hex.octree.svg index d83bdbe..feb35c6 100644 --- a/assets/icons/mesh.hex.generate.svg +++ b/assets/icons/mesh.hex.octree.svg @@ -1,6 +1,6 @@ - mesh.hex.generate + mesh.hex.octree diff --git a/assets/icons/mesh.hex.octree.vectorcraft b/assets/icons/mesh.hex.octree.vectorcraft new file mode 100644 index 0000000..f74d552 Binary files /dev/null and b/assets/icons/mesh.hex.octree.vectorcraft differ diff --git a/assets/icons/src/mesh.hex.generate.json b/assets/icons/src/mesh.hex.octree.json similarity index 100% rename from assets/icons/src/mesh.hex.generate.json rename to assets/icons/src/mesh.hex.octree.json diff --git a/crates/engine/src/catalogue.rs b/crates/engine/src/catalogue.rs index c3dadb1..e02b5e0 100644 --- a/crates/engine/src/catalogue.rs +++ b/crates/engine/src/catalogue.rs @@ -256,11 +256,25 @@ pub static CATALOGUE: &[Capability] = &[ "Quad or triangle mesh of a planar face or polygon in the XY plane for plane stress, plane strain and axisymmetric models; edges tagged", ), soon("mesh.quad.faces", "Quad Faces", Mesh, "Structured Hex", "M5", "All-quad meshes of faces on shared edge nodes"), - soon("mesh.hex.mapped", "Mapped Hex", Mesh, "Structured Hex", "M5", "Mapped hex mesh of six-sided regions"), - soon("mesh.hex.sweep", "Sweep Hex", Mesh, "Structured Hex", "M5", "Hex mesh swept along extrusions and revolutions"), + ready_more( + "mesh.hex.mapped", + "Mapped Hex", + Mesh, + "Structured Hex", + "Mapped hex mesh of a six-sided block (curved faces allowed)", + &["mesh.hex.generate"], + ), + ready_more( + "mesh.hex.sweep", + "Sweep Hex", + Mesh, + "Structured Hex", + "Hex mesh swept along extrusions (several levels), around revolutions, or from a source face to a target", + &["mesh.hex.generate"], + ), soon("mesh.hex.decompose", "Decompose", Mesh, "Structured Hex", "M5", "Find sweepable and mappable regions and propose the splits"), soon("mesh.transition", "Transition", Mesh, "Structured Hex", "M5", "Conforming tet–hex interfaces with pyramids"), - soon("mesh.hex.generate", "All-Hex", Mesh, "Automatic Hex", "M6", "Octree dual all-hex mesh of any part"), + soon("mesh.hex.octree", "All-Hex", Mesh, "Automatic Hex", "M6", "Octree dual all-hex mesh of any part"), soon("mesh.hex.dominant", "Hex-Dominant", Mesh, "Automatic Hex", "M6", "Hex-dominant mesh by recombining tets"), soon("mesh.hex.polycube", "Polycube", Mesh, "Automatic Hex", "M6", "Polycube hex mesh of blocky parts"), soon("mesh.smooth", "Smooth", Mesh, "Improve", "M6", "Smooth node positions, keeping every element valid"), diff --git a/crates/engine/src/cmd/hex.rs b/crates/engine/src/cmd/hex.rs new file mode 100644 index 0000000..565706c --- /dev/null +++ b/crates/engine/src/cmd/hex.rs @@ -0,0 +1,338 @@ +//! Structured hexahedral meshing of the open design's bodies or of imported surfaces. + +use std::borrow::Cow; + +use modelcraft_geometry::{FaceSurface, Shape, SurfaceMesh}; +use modelcraft_mesh::hex::{HexMesh, HexMethod, HexOptions, Surface, hex_mesh}; +use modelcraft_mesh::size::UniformSize; +use modelcraft_mesh::surface::BoundarySurface; +use modelcraft_model::{GeometryRef, Model}; +use serde_json::{Value, json}; + +use crate::{CommandKind, CommandSpec, EngineError, Result, Session}; + +pub const HEX_GENERATE: CommandSpec = CommandSpec { + id: "mesh.hex.generate", + label: "Generate Hex Mesh", + params: r#"{"body": "Body1" (default: every body of the open design) or "surface": 0 (an imported surface), "size": 2.0 (element size, mm; default 1/25 of the bounding-box diagonal), "order": 1 (hex8) or 2 (hex20, mid-edge nodes on the geometry), "method": "auto" (default) | "swept" (any sweep) | "extrude" | "sweep" | "revolve" | "mapped", "tolerance": 0.05 (chordal tolerance of the tessellation, mm; default 1e-3 of the diagonal)} — replaces the mesh with structured hexahedra: multi-level sweeps along a direction or around an axis, one-to-one sweeps through linking faces, or mapped blocks. Boundary element faces are tagged with their B-rep faces (named as in geometry.bodies; by index when a face has no name; `weld`, `weld_distance` and `feature_angle` as in mesh.tet.generate), so sets by geometry resolve. Reports per body the method, counts and quality; a body no structured method fits is an error that says why"#, + kind: CommandKind::Model, + run: hex_generate, +}; + +/// Bodies at most a command meshes. +const MAX_BODIES: usize = 1000; +/// Smallest element size and chordal tolerance, relative to the bounding-box diagonal. +const MIN_RELATIVE_SIZE: f64 = 1e-4; +const MIN_RELATIVE_TOLERANCE: f64 = 1e-5; +/// Default element size and chordal tolerance, relative to the bounding-box diagonal. +const DEFAULT_RELATIVE_SIZE: f64 = 1.0 / 25.0; +const DEFAULT_RELATIVE_TOLERANCE: f64 = 1e-3; + +/// One body to mesh. +struct Input<'a> { + name: String, + mesh: Cow<'a, SurfaceMesh>, + surfaces: Vec>, + face_names: Vec, + diagonal: f64, +} + +fn analytic(s: FaceSurface) -> Surface { + match s { + FaceSurface::Plane { normal, point } => Surface::Plane { point, normal }, + FaceSurface::Cylinder { origin, axis, radius } => Surface::Cylinder { origin, axis, radius }, + FaceSurface::Cone { apex, axis, half_angle } => Surface::Cone { apex, axis, half_angle }, + FaceSurface::Sphere { center, radius } => Surface::Sphere { center, radius }, + FaceSurface::Torus { center, axis, major, minor } => Surface::Torus { center, axis, major, minor }, + } +} + +fn number(p: &Value, key: &str) -> Result> { + match p.get(key) { + None => Ok(None), + Some(v) => v + .as_f64() + .filter(|x| x.is_finite() && *x > 0.0) + .map(Some) + .ok_or_else(|| EngineError::Params(format!("`{key}` must be a positive number (mm)"))), + } +} + +fn diagonal(lo: [f64; 3], hi: [f64; 3]) -> f64 { + (0..3).map(|i| (hi[i] - lo[i]).powi(2)).sum::().sqrt() +} + +fn shape_input(shape: &Shape, tolerance: Option) -> Result> { + let diag = shape.size()?; + let tol = tolerance.unwrap_or(DEFAULT_RELATIVE_TOLERANCE * diag); + if tol < MIN_RELATIVE_TOLERANCE * diag { + return Err(EngineError::Params(format!("`tolerance` must be at least {} mm for body `{}`", MIN_RELATIVE_TOLERANCE * diag, shape.name))); + } + let surfaces = shape.face_surfaces()?.into_iter().map(|f| f.map(analytic)).collect(); + let face_names = shape.faces().iter().map(|f| f.name.clone()).collect(); + Ok(Input { name: shape.name.clone(), mesh: Cow::Owned(shape.tessellate(tol)?), surfaces, face_names, diagonal: diag }) +} + +fn inputs<'a>(s: &'a Session, p: &Value) -> Result>> { + let tolerance = number(p, "tolerance")?; + match (p.get("body"), p.get("surface")) { + (Some(b), None) => { + let name = b.as_str().ok_or_else(|| EngineError::Params("`body` must be a body name".into()))?; + let shape = s + .design + .as_deref() + .and_then(|d| d.body(name)) + .ok_or_else(|| EngineError::Params(format!("no body `{name}` (open a design first)")))?; + Ok(vec![shape_input(shape, tolerance)?]) + } + (None, Some(i)) => { + let n = s.surfaces.len(); + let index = i.as_u64().and_then(|i| usize::try_from(i).ok()).filter(|&i| i < n); + let index = index.ok_or_else(|| EngineError::Params(format!("`surface` must be the index of an imported surface (there are {n})")))?; + let t = s.surfaces.get(index).ok_or_else(|| EngineError::Params("no such surface".into()))?; + let (mut lo, mut hi) = ([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]); + for q in &t.mesh.positions { + for k in 0..3 { + lo[k] = lo[k].min(q[k]); + hi[k] = hi[k].max(q[k]); + } + } + let name = if t.name.trim().is_empty() { format!("Surface{index}") } else { t.name.clone() }; + // Faces named by index, as mesh.tet.generate names them. + let face_names = (0..t.faces).map(|f| f.to_string()).collect(); + Ok(vec![Input { name, mesh: Cow::Borrowed(&t.mesh), surfaces: Vec::new(), face_names, diagonal: diagonal(lo, hi) }]) + } + (None, None) => { + let design = s.design.as_deref().ok_or_else(|| EngineError::Params("open a design first, or give `surface`".into()))?; + if design.bodies.is_empty() || design.bodies.len() > MAX_BODIES { + return Err(EngineError::Params(format!("the design has {} bodies; meshing takes 1 to {MAX_BODIES}", design.bodies.len()))); + } + design.bodies.iter().map(|b| shape_input(b, tolerance)).collect() + } + _ => Err(EngineError::Params("give `body` or `surface`, not both".into())), + } +} + +pub const HEX_SWEEP: CommandSpec = CommandSpec { + id: "mesh.hex.sweep", + label: "Sweep Hex Mesh", + params: r#"the parameters of mesh.hex.generate, with "method": "swept" (default: extrude, sweep or revolve, whichever fits first) | "extrude" | "sweep" | "revolve" — mesh.hex.generate limited to the sweeping methods"#, + kind: CommandKind::Model, + run: hex_sweep, +}; + +pub const HEX_MAPPED: CommandSpec = CommandSpec { + id: "mesh.hex.mapped", + label: "Mapped Hex Mesh", + params: r#"the parameters of mesh.hex.generate without "method" — mesh.hex.generate with the mapped method: a six-sided block (six faces of four edges, three edges at each of eight corners)"#, + kind: CommandKind::Model, + run: hex_mapped, +}; + +fn hex_generate(s: &mut Session, p: &Value) -> Result { + let all = [HexMethod::Auto, HexMethod::Swept, HexMethod::Extrude, HexMethod::Sweep, HexMethod::Revolve, HexMethod::Mapped]; + generate(s, p, HexMethod::Auto, &all) +} + +fn hex_sweep(s: &mut Session, p: &Value) -> Result { + generate(s, p, HexMethod::Swept, &[HexMethod::Swept, HexMethod::Extrude, HexMethod::Sweep, HexMethod::Revolve]) +} + +fn hex_mapped(s: &mut Session, p: &Value) -> Result { + generate(s, p, HexMethod::Mapped, &[HexMethod::Mapped]) +} + +/// Mesh with `default` or the `method` parameter, which must be one of `allowed`. +fn generate(s: &mut Session, p: &Value, default: HexMethod, allowed: &[HexMethod]) -> Result { + let size = number(p, "size")?; + let quadratic = match p.get("order") { + None => false, + Some(v) => match v.as_u64() { + Some(1) => false, + Some(2) => true, + _ => return Err(EngineError::Params("`order` must be 1 (hex8) or 2 (hex20)".into())), + }, + }; + let method = match p.get("method") { + None => default, + Some(v) => v.as_str().and_then(HexMethod::from_name).filter(|m| allowed.contains(m)).ok_or_else(|| { + let names: Vec<&str> = allowed.iter().map(|m| m.name()).collect(); + EngineError::Params(format!("`method` must be one of {}", names.join(", "))) + })?, + }; + let bodies = inputs(s, p)?; + let mut meshes: Vec<(Input, HexMesh)> = Vec::with_capacity(bodies.len()); + for b in bodies { + if !(b.diagonal > 0.0 && b.diagonal.is_finite()) { + return Err(EngineError::Params(format!("body `{}` has no extent", b.name))); + } + let h = size.unwrap_or(DEFAULT_RELATIVE_SIZE * b.diagonal); + if h < MIN_RELATIVE_SIZE * b.diagonal { + return Err(EngineError::Params(format!("`size` must be at least {} mm for body `{}`", MIN_RELATIVE_SIZE * b.diagonal, b.name))); + } + let surface_options = super::mesh::surface_options(p, &b.mesh)?; + let surface = BoundarySurface::new(&b.mesh.positions, &b.mesh.triangles, &b.mesh.tri_face, &surface_options) + .map_err(|e| EngineError::Params(format!("body `{}`: {e}", b.name)))?; + let chord = if b.surfaces.is_empty() { 0.0 } else { number(p, "tolerance")?.unwrap_or(DEFAULT_RELATIVE_TOLERANCE * b.diagonal) }; + let options = + HexOptions { method, quadratic, chordal_tolerance: chord, max_elements: modelcraft_model::MAX_ELEMENTS / 8, ..HexOptions::default() }; + let mesh = hex_mesh(&surface, &b.surfaces, &UniformSize(h), &options).map_err(|e| EngineError::Params(format!("body `{}`: {e}", b.name)))?; + meshes.push((b, mesh)); + } + // Replace the mesh: definitions stay, nodes, elements and geometry tags are new. + let mut model = s.model.without_mesh(); + let mut report = Vec::with_capacity(meshes.len()); + let (mut node_label, mut element_label) = (1u32, 1u32); + for (b, m) in &meshes { + let first_node = node_label; + for xyz in &m.nodes { + model.set_node(node_label, *xyz)?; + node_label = node_label.checked_add(1).ok_or_else(|| EngineError::Params("too many nodes".into()))?; + } + let body_tag = model.geometry_tag(&GeometryRef::body(b.name.as_str())?)?; + let first_element = element_label; + for e in 0..m.element_count() { + let nodes: Vec = m.element(e).unwrap_or(&[]).iter().map(|&i| first_node + i).collect(); + model.set_element(element_label, m.kind, nodes)?; + model.tag_element(element_label, body_tag)?; + element_label = element_label.checked_add(1).ok_or_else(|| EngineError::Params("too many elements".into()))?; + } + tag_faces(&mut model, b, m, first_element)?; + report.push(json!({ + "body": b.name, + "method": m.report.method.name(), + "description": m.report.description, + "kind": m.kind.name(), + "elements": m.element_count(), + "nodes": m.nodes.len(), + "boundary_faces": m.boundary.len(), + "quality": {"min_scaled_jacobian": m.report.min_scaled_jacobian, "mean_scaled_jacobian": m.report.mean_scaled_jacobian}, + "volume": m.report.volume, + "surface_volume": m.report.surface_volume, + })); + } + s.model = model; + Ok(json!({"bodies": report, "nodes": node_label - 1, "elements": element_label - 1})) +} + +/// Tag every boundary element face with its B-rep face. +fn tag_faces(model: &mut Model, b: &Input<'_>, m: &HexMesh, first_element: u32) -> Result<()> { + let mut tags: Vec> = Vec::new(); + for f in &m.boundary { + let i = f.brep_face as usize; + if tags.len() <= i { + tags.resize(i + 1, None); + } + let tag = match tags.get(i).copied().flatten() { + Some(t) => t, + None => { + let name = b.face_names.get(i).cloned().unwrap_or_else(|| i.to_string()); + let t = model.geometry_tag(&GeometryRef::face(b.name.as_str(), name)?)?; + if let Some(slot) = tags.get_mut(i) { + *slot = Some(t); + } + t + } + }; + model.tag_element_face(first_element + f.element, f.face, tag)?; + } + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::*; + use modelcraft_geometry::sample::Plate; + + fn opened() -> (crate::tests::TempDir, Session) { + let (dir, path) = crate::tests::sample_design_file(&Plate::default()); + let mut s = Session::default(); + s.execute("geometry.open", &json!({"path": path})).unwrap(); + (dir, s) + } + + #[test] + fn plate_with_a_hole_meshes_with_hexahedra_and_face_sets_resolve() { + let (_dir, mut s) = opened(); + let r = s.execute("mesh.hex.generate", &json!({"body": "Body1", "size": 2.5})).unwrap(); + let body = &r["bodies"][0]; + assert_eq!(body["method"], "extrude", "{r}"); + assert_eq!(body["kind"], "hex8"); + let n = body["elements"].as_u64().unwrap() as usize; + assert!(n > 100, "{r}"); + assert!(body["quality"]["min_scaled_jacobian"].as_f64().unwrap() > 0.3, "{r}"); + assert_eq!(s.model.element_counts(), vec![(modelcraft_model::ElementKind::Hex8, n)]); + // Volume: the plate minus the hole, to the tessellation's accuracy. + let v = 40.0 * 30.0 * 10.0 - std::f64::consts::PI * 25.0 * 10.0; + assert!((body["volume"].as_f64().unwrap() - v).abs() < 0.01 * v, "{r}"); + // A set on the top face resolves to element faces. + let bodies = s.execute("geometry.bodies", &json!({})).unwrap(); + let faces: Vec = bodies[0]["faces"].as_array().unwrap().iter().map(|f| f["name"].as_str().unwrap().to_string()).collect(); + let top = faces.iter().find(|f| f.ends_with("top") || f.contains("end")).unwrap_or(&faces[0]).clone(); + s.execute("model.set.create", &json!({"name": "top", "geometry": [{"body": "Body1", "face": top}]})).unwrap(); + let set = s.execute("model.set.resolve", &json!({"name": "top"})).unwrap(); + assert!(set["faces"]["count"].as_u64().unwrap() > 0, "{set}"); + // Every boundary face is tagged: the skin of the body set. + s.execute("model.set.create", &json!({"name": "skin", "geometry": [{"body": "Body1"}]})).unwrap(); + let skin = s.execute("model.set.resolve", &json!({"name": "skin"})).unwrap(); + assert_eq!(skin["faces"]["count"], body["boundary_faces"], "{skin}"); + assert_eq!(skin["elements"]["count"].as_u64().unwrap() as usize, n); + // Hex20 with mid-edge nodes on the hole. + let r2 = s.execute("mesh.hex.generate", &json!({"order": 2, "size": 2.5})).unwrap(); + assert_eq!(r2["bodies"][0]["kind"], "hex20"); + assert_eq!(s.model.element_counts(), vec![(modelcraft_model::ElementKind::Hex20, n)]); + let v2 = r2["bodies"][0]["volume"].as_f64().unwrap(); + assert!((v2 - v).abs() < (body["volume"].as_f64().unwrap() - v).abs() + 1e-6 * v, "hex20 {v2} vs hex8 {} vs {v}", body["volume"]); + // The set still resolves after the remesh. + let set = s.execute("model.set.resolve", &json!({"name": "top"})).unwrap(); + assert!(set["faces"]["count"].as_u64().unwrap() > 0); + // Undo returns the hex8 mesh. + s.execute("edit.undo", &json!({})).unwrap(); + assert_eq!(s.model.element_counts(), vec![(modelcraft_model::ElementKind::Hex8, n)]); + } + + #[test] + fn hostile_params_are_errors_and_change_nothing() { + let (_dir, mut s) = opened(); + let before = s.model.clone(); + for p in [ + json!({"body": "Body1", "size": 0}), + json!({"body": "Body1", "size": -1}), + json!({"body": "Body1", "size": 1e-300}), + json!({"body": "Body1", "size": "big"}), + json!({"body": "Body1", "size": 1e300}), + json!({"body": "Body1", "order": 3}), + json!({"body": "Body1", "order": "2"}), + json!({"body": "Body1", "method": "octree"}), + json!({"body": "Body1", "method": null}), + json!({"body": "Body1", "tolerance": 1e-300}), + json!({"body": "Body1", "tolerance": -1}), + json!({"body": "Nope"}), + json!({"body": 7}), + json!({"surface": 0}), + json!({"surface": -1}), + json!({"body": "Body1", "surface": 0}), + json!({"body": "Body1", "method": "revolve"}), + json!({"body": "Body1", "method": "mapped"}), + ] { + let r = s.execute("mesh.hex.generate", &p); + if p.get("size") == Some(&json!(1e300)) { + // One element through the whole body is still a valid request. + if r.is_ok() { + s.model = before.clone(); + } + continue; + } + assert!(r.is_err(), "{p}: {r:?}"); + assert_eq!(s.model, before, "{p}"); + } + assert!(Session::default().execute("mesh.hex.generate", &json!({})).is_err()); + assert!(s.execute("mesh.hex.sweep", &json!({"body": "Body1", "method": "mapped"})).is_err()); + assert!(s.execute("mesh.hex.mapped", &json!({"body": "Body1", "method": "sweep"})).is_err()); + // The plate is not a six-sided block; it is a sweep. + assert!(s.execute("mesh.hex.mapped", &json!({"body": "Body1"})).unwrap_err().to_string().contains("logical cube")); + assert_eq!(s.model, before); + let r = s.execute("mesh.hex.sweep", &json!({"body": "Body1", "size": 4})).unwrap(); + assert_eq!(r["bodies"][0]["method"], "extrude"); + } +} diff --git a/crates/engine/src/cmd/mesh.rs b/crates/engine/src/cmd/mesh.rs index 7cc78d0..e11678c 100644 --- a/crates/engine/src/cmd/mesh.rs +++ b/crates/engine/src/cmd/mesh.rs @@ -104,8 +104,8 @@ fn optional_count(p: &Value, key: &str, default: u64, range: std::ops::RangeIncl usize::try_from(v).map_err(|_| EngineError::Params(format!("`{key}` is too large"))) } -/// The boundary options of `mesh.surface.check` and `mesh.tet.generate`. -fn surface_options(p: &Value, mesh: &SurfaceMesh) -> Result { +/// The boundary options of `mesh.surface.check`, `mesh.tet.generate` and the hex meshers. +pub(super) fn surface_options(p: &Value, mesh: &SurfaceMesh) -> Result { let defaults = SurfaceOptions::default(); Ok(SurfaceOptions { weld_tolerance: optional_number(p, "weld", defaults.weld_tolerance)?, diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index ec6bcbd..2d73866 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -6,7 +6,8 @@ pub(crate) mod define; mod edit; pub(crate) mod file; mod geometry; -mod mesh; +mod hex; +pub(crate) mod mesh; mod model; mod results; pub(crate) mod solve; @@ -36,6 +37,9 @@ pub static ALL: &[CommandSpec] = &[ mesh::PLANAR, mesh::TET_GENERATE, block::BLOCK, + hex::HEX_GENERATE, + hex::HEX_SWEEP, + hex::HEX_MAPPED, tags::TAG, tags::TAG_CLEAR, model::UNITS, diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs index 83071b8..2d24615 100644 --- a/crates/engine/src/lib.rs +++ b/crates/engine/src/lib.rs @@ -375,6 +375,11 @@ mod tests { json!({"polygon": [["a", 0], [1, null], {}], "size": "big", "kind": 7}), json!({"polygon": [[0, 0], [1, 0], [0, 1]], "size": 1e300, "kind": "tri6", "max_elements": 1e300}), json!({"body": "Body1", "face": "", "tolerance": 0.05, "size": 1}), + // Structured hex meshing. + json!({"size": -1, "order": 9, "method": "x"}), + json!({"size": 1e-300}), + json!({"order": 2, "method": "mapped", "tolerance": 1e300}), + json!({"surface": 0, "size": 1e300, "order": 2}), // Model definitions. json!({"name": "x".repeat(100_000)}), json!({"name": "A\u{0}"}), @@ -506,6 +511,10 @@ mod tests { json!({"body": "Body1", "tolerance": 0.5, "weld_distance": 0.05}), json!({"body": "Body1", "tolerance": 0.5, "feature_angle": 1e-300}), json!({"body": "Body1", "tolerance": 0.5, "feature_angle": null}), + json!({"body": "Body1", "size": 1e-300}), + json!({"body": "Body1", "size": 1e300, "order": 2}), + json!({"body": "Body1", "size": 4, "method": "sweep", "order": 2}), + json!({"body": "Body1", "method": "revolve"}), ]; for p in tess.iter().chain(&hostile) { for spec in commands() { @@ -516,14 +525,14 @@ mod tests { } /// A design file in a fresh temporary directory (removed when the guard drops). - fn sample_design_file(plate: &Plate) -> (TempDir, String) { + pub(crate) fn sample_design_file(plate: &Plate) -> (TempDir, String) { let dir = TempDir::new(); let path = dir.0.join("plate.solidcraft"); std::fs::write(&path, plate_design(plate).unwrap()).unwrap(); (dir, path.to_string_lossy().into_owned()) } - struct TempDir(std::path::PathBuf); + pub(crate) struct TempDir(std::path::PathBuf); impl TempDir { fn new() -> Self { static N: std::sync::atomic::AtomicUsize = std::sync::atomic::AtomicUsize::new(0); diff --git a/crates/geometry/src/lib.rs b/crates/geometry/src/lib.rs index 84b5a4b..6cc6b8b 100644 --- a/crates/geometry/src/lib.rs +++ b/crates/geometry/src/lib.rs @@ -78,6 +78,39 @@ pub struct Face { pub exact: bool, } +/// The analytic surface a B-rep face lies on, when it is one of the elementary kinds. Directions +/// are unit vectors; lengths are millimetres and angles radians. +#[derive(Debug, Clone, Copy, PartialEq)] +pub enum FaceSurface { + /// The outward unit normal and a point on the plane. + Plane { + normal: [f64; 3], + point: [f64; 3], + }, + Cylinder { + origin: [f64; 3], + axis: [f64; 3], + radius: f64, + }, + /// The apex, the axis pointing into the cone and the half-angle. + Cone { + apex: [f64; 3], + axis: [f64; 3], + half_angle: f64, + }, + Sphere { + center: [f64; 3], + radius: f64, + }, + /// The centre of the ring, its axis, and the ring (major) and tube (minor) radii. + Torus { + center: [f64; 3], + axis: [f64; 3], + major: f64, + minor: f64, + }, +} + /// Mass properties and extent of a [`Shape`]. #[derive(Debug, Clone, Copy, PartialEq)] pub struct Measure { @@ -127,6 +160,37 @@ impl Shape { self.faces.iter().position(|f| f.name == name).and_then(|i| u32::try_from(i).ok()) } + /// The analytic surface of each B-rep face, by face index (the [`SurfaceMesh::tri_face`] + /// tags): `None` for a free-form face. Meshers project nodes onto these surfaces. + pub fn face_surfaces(&self) -> Result>> { + guard("face surfaces", || { + Ok((0..self.body.face_count()) + .map(|f| { + solidcraft_kernel::face_surface(&self.body, f).map(|s| match s { + solidcraft_kernel::FaceSurface::Plane { normal, point } => FaceSurface::Plane { normal: xyz(normal), point: xyz(point) }, + solidcraft_kernel::FaceSurface::Cylinder { origin, axis, radius } => { + FaceSurface::Cylinder { origin: xyz(origin), axis: xyz(axis), radius } + } + solidcraft_kernel::FaceSurface::Cone { apex, axis, half_angle } => { + FaceSurface::Cone { apex: xyz(apex), axis: xyz(axis), half_angle } + } + solidcraft_kernel::FaceSurface::Sphere { center, radius } => FaceSurface::Sphere { center: xyz(center), radius }, + solidcraft_kernel::FaceSurface::Torus { center, axis, major, minor } => { + FaceSurface::Torus { center: xyz(center), axis: xyz(axis), major, minor } + } + }) + }) + .collect()) + }) + .map_err(Error::from) + } + + /// A characteristic size: the diagonal of the bounding box of the body's B-rep vertices + /// (cached by the kernel; curved faces may bulge past it). For default tolerances and sizes. + pub fn size(&self) -> Result { + Ok(guard("size", || Ok(self.body.size()))?) + } + pub fn measure(&self) -> Result { let m = guard("measure", || solidcraft_kernel::measure(&self.body))?; Ok(Measure { @@ -184,6 +248,19 @@ mod tests { assert_eq!(m.bbox_max, [10., 20., 30.]); } + #[test] + fn block_faces_are_planes_through_the_box() { + let surfaces = block([0., 0., 0.], [10., 20., 30.]).unwrap().face_surfaces().unwrap(); + assert_eq!(surfaces.len(), 6); + for s in surfaces { + let Some(FaceSurface::Plane { normal, point }) = s else { panic!("{s:?}") }; + // Each plane holds one of the box's faces: the point's coordinate along the normal is + // 0 or the box's extent there. + let along: f64 = normal.iter().zip(point).map(|(n, p)| n * p).sum(); + assert!([0.0, 10.0, 20.0, 30.0].iter().any(|v| (along.abs() - v).abs() < 1e-9), "{normal:?} {point:?}"); + } + } + #[test] fn tessellation_covers_every_face_and_encloses_the_volume() { let mesh = block([0., 0., 0.], [2., 3., 4.]).unwrap().tessellate(0.01).unwrap(); diff --git a/crates/mesh/src/hex/arrangement.rs b/crates/mesh/src/hex/arrangement.rs new file mode 100644 index 0000000..3966f72 --- /dev/null +++ b/crates/mesh/src/hex/arrangement.rs @@ -0,0 +1,623 @@ +//! The planar arrangement of the curves a multi-level sweep projects into its cross-section +//! plane: points merged within a tolerance, segments split where a point lies on them or where +//! two cross, then the faces of the planar graph (regions with their holes) and the curves +//! between curve ends. +//! +//! This is the imprinting step of many-to-many sweeping (Lai, Benzley, White, *Automated +//! hexahedral mesh generation by generalized multiple source to multiple target sweeping*, +//! IJNME 49(1–2):261–275, 2000): projecting every level's boundary into one plane splits the +//! cross-section into regions whose columns are each a one-to-one sweep. Crossing and side tests +//! use the exact `orient2d` predicate on the merged coordinates; the merge itself is a tolerance +//! (the projected copies of one point differ by rounding). + +use std::collections::{BTreeMap, HashMap}; + +use crate::predicates::{Point2, Sign, orient2d, snap_to_domain}; +use crate::{Error, Result}; + +pub(crate) type Q = [f64; 2]; + +/// The most segments one arrangement may hold. +pub(crate) const MAX_SEGMENTS: usize = 2_000_000; +/// The most rounds of splitting (each round settles every T-junction and crossing found). +const MAX_ROUNDS: usize = 16; + +/// A curve of the arrangement: a chain of points between curve ends (or a closed chain). +#[derive(Debug, Clone)] +pub(crate) struct Curve { + /// Point ids in order; closed curves repeat the first at the end. + pub points: Vec, + /// The tags of the input polylines it came from (sorted, unique). + pub tags: Vec, +} + +/// A bounded face of the arrangement. +#[derive(Debug, Clone)] +pub(crate) struct Region { + /// Loops of oriented curves `(curve, forward)`: the outer loop first (counter-clockwise), then + /// the holes (clockwise). + pub loops: Vec>, +} + +#[derive(Debug, Clone)] +pub(crate) struct Arrangement { + pub points: Vec, + pub curves: Vec, + pub regions: Vec, +} + +pub(super) fn d2(a: Q, b: Q) -> f64 { + (a[0] - b[0]).powi(2) + (a[1] - b[1]).powi(2) +} + +/// Squared distance from `p` to segment `ab`, and the parameter of the closest point. +pub(super) fn seg_dist2(p: Q, a: Q, b: Q) -> (f64, f64) { + let (dx, dy) = (b[0] - a[0], b[1] - a[1]); + let len2 = dx * dx + dy * dy; + if !(len2 > 0.0) { + return (d2(p, a), 0.0); + } + let t = (((p[0] - a[0]) * dx + (p[1] - a[1]) * dy) / len2).clamp(0.0, 1.0); + (d2(p, [a[0] + t * dx, a[1] + t * dy]), t) +} + +pub(super) fn sign(a: Q, b: Q, c: Q) -> Result { + Ok(orient2d(&Point2::new(a[0], a[1])?, &Point2::new(b[0], b[1])?, &Point2::new(c[0], c[1])?)) +} + +/// Points merged within a tolerance through a uniform grid; the first point of a cluster is its +/// representative, so the result follows insertion order. +struct Snapper { + tol: f64, + cell: f64, + grid: HashMap<(i64, i64), Vec>, + points: Vec, +} + +impl Snapper { + fn new(tol: f64) -> Self { + Snapper { tol, cell: tol * 4.0, grid: HashMap::new(), points: Vec::new() } + } + + fn key(&self, p: Q) -> (i64, i64) { + ((p[0] / self.cell).floor() as i64, (p[1] / self.cell).floor() as i64) + } + + fn add(&mut self, p: Q) -> Result { + let p = [snap_to_domain(p[0]), snap_to_domain(p[1])]; + Point2::new(p[0], p[1])?; + let (i, j) = self.key(p); + let mut best: Option<(f64, usize)> = None; + for di in -1..=1 { + for dj in -1..=1 { + for &k in self.grid.get(&(i + di, j + dj)).map(Vec::as_slice).unwrap_or(&[]) { + let d = d2(p, self.points.get(k).copied().unwrap_or([f64::INFINITY; 2])); + if d <= self.tol * self.tol && best.is_none_or(|(bd, bk)| d < bd || (d == bd && k < bk)) { + best = Some((d, k)); + } + } + } + } + if let Some((_, k)) = best { + return Ok(k); + } + self.points.push(p); + let id = self.points.len() - 1; + self.grid.entry((i, j)).or_default().push(id); + Ok(id) + } +} + +/// A grid over boxes, for candidate queries. +struct BoxGrid { + cell: f64, + cells: HashMap<(i64, i64), Vec>, +} + +/// The most grid cells one box may cover before the grid is judged too fine. +const MAX_CELLS_PER_BOX: i64 = 1 << 16; + +impl BoxGrid { + fn range(&self, lo: Q, hi: Q) -> Option<(i64, i64, i64, i64)> { + let (i0, j0) = ((lo[0] / self.cell).floor() as i64, (lo[1] / self.cell).floor() as i64); + let (i1, j1) = ((hi[0] / self.cell).floor() as i64, (hi[1] / self.cell).floor() as i64); + ((i1 - i0 + 1).saturating_mul(j1 - j0 + 1) <= MAX_CELLS_PER_BOX).then_some((i0, j0, i1, j1)) + } + + fn insert(&mut self, id: usize, lo: Q, hi: Q) -> Result<()> { + let (i0, j0, i1, j1) = self.range(lo, hi).ok_or_else(|| Error::Limit("a segment covers too many grid cells".into()))?; + for i in i0..=i1 { + for j in j0..=j1 { + self.cells.entry((i, j)).or_default().push(id); + } + } + Ok(()) + } + + /// Ids in the cells the box covers (sorted, unique). + fn query(&self, lo: Q, hi: Q) -> Result> { + let (i0, j0, i1, j1) = self.range(lo, hi).ok_or_else(|| Error::Limit("a query covers too many grid cells".into()))?; + let mut out = Vec::new(); + for i in i0..=i1 { + for j in j0..=j1 { + out.extend(self.cells.get(&(i, j)).map(Vec::as_slice).unwrap_or(&[])); + } + } + out.sort_unstable(); + out.dedup(); + Ok(out) + } +} + +fn bbox(a: Q, b: Q, pad: f64) -> (Q, Q) { + ([a[0].min(b[0]) - pad, a[1].min(b[1]) - pad], [a[0].max(b[0]) + pad, a[1].max(b[1]) + pad]) +} + +/// Build the arrangement of `polylines` (each with a tag), with `vertices` (curve ends that must +/// stay ends) merged into it. Points closer than `tol` are one point. Two input polylines with +/// the same end points that stay within `coincide` of each other are one curve (the copies of a +/// wall's outline at different levels). +pub(crate) fn build(vertices: &[Q], polylines: &[(Vec, usize)], tol: f64, coincide: f64) -> Result { + if !(tol > 0.0 && tol.is_finite() && coincide >= 0.0 && coincide.is_finite()) { + return Err(Error::Invalid("arrangement tolerances must be positive and finite".into())); + } + let mut snap = Snapper::new(tol); + let mut is_end = Vec::new(); + for &v in vertices { + let id = snap.add(v)?; + if is_end.len() <= id { + is_end.resize(id + 1, false); + } + if let Some(x) = is_end.get_mut(id) { + *x = true; + } + } + // Polylines as point ids, dropping repeats; then coincident duplicates. + let mut lines: Vec<(Vec, Vec)> = Vec::new(); + let mut by_ends: BTreeMap<(usize, usize), Vec> = BTreeMap::new(); + for (pl, tag) in polylines { + let mut ids: Vec = Vec::with_capacity(pl.len()); + for &p in pl { + let id = snap.add(p)?; + if ids.last() != Some(&id) { + ids.push(id); + } + } + if ids.len() < 2 { + continue; + } + let (a, b) = (ids.first().copied().unwrap_or(0), ids.last().copied().unwrap_or(0)); + let key = (a.min(b), a.max(b)); + let pts = &snap.points; + let at = |i: usize| pts.get(i).copied().unwrap_or([0.0; 2]); + let near = |p: Q, other: &[usize]| other.windows(2).any(|w| seg_dist2(p, at(w[0]), at(w[1])).0 <= coincide * coincide); + let same = by_ends.get(&key).and_then(|cands| { + cands + .iter() + .copied() + .find(|&k| lines.get(k).is_some_and(|(other, _)| ids.iter().all(|&i| near(at(i), other)) && other.iter().all(|&i| near(at(i), &ids)))) + }); + match same { + Some(k) => { + if let Some(l) = lines.get_mut(k) { + l.1.push(*tag); + } + } + None => { + by_ends.entry(key).or_default().push(lines.len()); + lines.push((ids, vec![*tag])); + } + } + } + // Segments, unique by end points. + let mut segs: BTreeMap<(usize, usize), Vec> = BTreeMap::new(); + for (ids, tags) in &lines { + for w in ids.windows(2) { + let key = (w[0].min(w[1]), w[0].max(w[1])); + segs.entry(key).or_default().extend(tags); + if segs.len() > MAX_SEGMENTS { + return Err(Error::Limit(format!("the arrangement exceeds {MAX_SEGMENTS} segments"))); + } + } + } + let mut segs: Vec<((usize, usize), Vec)> = segs.into_iter().collect(); + // Split at points on segments and at crossings until nothing changes. + let mut settled = false; + for _ in 0..MAX_ROUNDS { + let pts = snap.points.clone(); + let at = |i: usize| pts.get(i).copied().unwrap_or([0.0; 2]); + let total: f64 = segs.iter().map(|((a, b), _)| d2(at(*a), at(*b)).sqrt()).sum(); + let cell = (total / segs.len().max(1) as f64).max(tol * 16.0); + let mut seg_grid = BoxGrid { cell, cells: HashMap::new() }; + for (k, ((a, b), _)) in segs.iter().enumerate() { + let (lo, hi) = bbox(at(*a), at(*b), tol); + seg_grid.insert(k, lo, hi)?; + } + // Crossings add points. + let mut added = false; + for (k, ((a, b), _)) in segs.iter().enumerate() { + let (lo, hi) = bbox(at(*a), at(*b), tol); + for other in seg_grid.query(lo, hi)? { + if other <= k { + continue; + } + let Some(((c, d), _)) = segs.get(other) else { continue }; + if a == c || a == d || b == c || b == d { + continue; + } + let (pa, pb, pc, pd) = (at(*a), at(*b), at(*c), at(*d)); + let (o1, o2) = (sign(pa, pb, pc)?, sign(pa, pb, pd)?); + let (o3, o4) = (sign(pc, pd, pa)?, sign(pc, pd, pb)?); + let opposite = |x: Sign, y: Sign| matches!((x, y), (Sign::Positive, Sign::Negative) | (Sign::Negative, Sign::Positive)); + if opposite(o1, o2) && opposite(o3, o4) { + let den = (pb[0] - pa[0]) * (pd[1] - pc[1]) - (pb[1] - pa[1]) * (pd[0] - pc[0]); + let t = ((pc[0] - pa[0]) * (pd[1] - pc[1]) - (pc[1] - pa[1]) * (pd[0] - pc[0])) / den; + if t.is_finite() { + let before = snap.points.len(); + snap.add([pa[0] + t * (pb[0] - pa[0]), pa[1] + t * (pb[1] - pa[1])])?; + added |= snap.points.len() > before; + } + } + } + } + // Points on segments split them. + let pts = snap.points.clone(); + let at = |i: usize| pts.get(i).copied().unwrap_or([0.0; 2]); + let mut pt_grid = BoxGrid { cell, cells: HashMap::new() }; + for (i, p) in pts.iter().enumerate() { + pt_grid.insert(i, *p, *p)?; + } + let mut split = false; + let mut next: BTreeMap<(usize, usize), Vec> = BTreeMap::new(); + for ((a, b), tags) in &segs { + let (pa, pb) = (at(*a), at(*b)); + let (lo, hi) = bbox(pa, pb, tol); + let mut on: Vec<(f64, usize)> = Vec::new(); + for p in pt_grid.query(lo, hi)? { + if p == *a || p == *b { + continue; + } + let (dd, t) = seg_dist2(at(p), pa, pb); + if dd <= tol * tol && t > 0.0 && t < 1.0 { + on.push((t, p)); + } + } + on.sort_by(|x, y| x.0.total_cmp(&y.0).then(x.1.cmp(&y.1))); + split |= !on.is_empty(); + let chain: Vec = std::iter::once(*a).chain(on.into_iter().map(|x| x.1)).chain(std::iter::once(*b)).collect(); + for w in chain.windows(2) { + if w[0] != w[1] { + next.entry((w[0].min(w[1]), w[0].max(w[1]))).or_default().extend(tags); + } + } + if next.len() > MAX_SEGMENTS { + return Err(Error::Limit(format!("the arrangement exceeds {MAX_SEGMENTS} segments"))); + } + } + segs = next.into_iter().collect(); + if !split && !added { + settled = true; + break; + } + } + if !settled { + return Err(Error::Failed("the projected curves did not settle into a planar arrangement".into())); + } + for (_, tags) in segs.iter_mut() { + tags.sort_unstable(); + tags.dedup(); + } + faces(snap.points, is_end, segs) +} + +/// The planar graph's curves and bounded faces. +fn faces(points: Vec, mut is_end: Vec, segs: Vec<((usize, usize), Vec)>) -> Result { + let n = points.len(); + is_end.resize(n, false); + let at = |i: usize| points.get(i).copied().unwrap_or([0.0; 2]); + // Outgoing half-edges per point, counter-clockwise. Half-edge 2s goes a → b, 2s+1 b → a. + let head = |h: usize| segs.get(h / 2).map_or(0, |((a, b), _)| if h.is_multiple_of(2) { *b } else { *a }); + let tail = |h: usize| segs.get(h / 2).map_or(0, |((a, b), _)| if h.is_multiple_of(2) { *a } else { *b }); + let mut out: Vec> = vec![Vec::new(); n]; + for h in 0..segs.len() * 2 { + out.get_mut(tail(h)).ok_or_else(|| Error::Failed("arrangement point".into()))?.push(h); + } + for (p, list) in out.iter_mut().enumerate() { + let o = at(p); + let half = |h: usize| { + let q = at(head(h)); + let (dx, dy) = (q[0] - o[0], q[1] - o[1]); + if dy > 0.0 || (dy == 0.0 && dx > 0.0) { 0 } else { 1 } + }; + let mut err = None; + list.sort_by(|&x, &y| { + half(x).cmp(&half(y)).then_with(|| match sign(o, at(head(x)), at(head(y))) { + Ok(Sign::Positive) => std::cmp::Ordering::Less, + Ok(Sign::Negative) => std::cmp::Ordering::Greater, + Ok(Sign::Zero) => x.cmp(&y), + Err(e) => { + err = Some(e); + x.cmp(&y) + } + }) + }); + if let Some(e) = err { + return Err(e); + } + } + let pos: Vec = { + let mut pos = vec![0; segs.len() * 2]; + for list in &out { + for (k, &h) in list.iter().enumerate() { + if let Some(x) = pos.get_mut(h) { + *x = k; + } + } + } + pos + }; + let next = |h: usize| -> usize { + let b = head(h); + let twin = h ^ 1; + let list = out.get(b).map(Vec::as_slice).unwrap_or(&[]); + let k = pos.get(twin).copied().unwrap_or(0); + list.get((k + list.len() - 1) % list.len().max(1)).copied().unwrap_or(twin) + }; + // Curves: chains between ends (degree ≠ 2 or marked). + let degree = |p: usize| out.get(p).map_or(0, Vec::len); + let end = |p: usize| degree(p) != 2 || is_end.get(p).copied().unwrap_or(false); + let mut curve_of: Vec<(usize, usize, bool)> = vec![(usize::MAX, 0, true); segs.len() * 2]; + let mut curves: Vec = Vec::new(); + let mut seen = vec![false; segs.len()]; + let walk = |start: usize, curves: &mut Vec, curve_of: &mut Vec<(usize, usize, bool)>, seen: &mut Vec| { + let id = curves.len(); + let mut pts = vec![tail(start)]; + let mut tags: Vec = Vec::new(); + let mut h = start; + for k in 0..=segs.len() { + if let Some(s) = seen.get_mut(h / 2) { + *s = true; + } + if let Some(x) = curve_of.get_mut(h) { + *x = (id, k, true); + } + if let Some(x) = curve_of.get_mut(h ^ 1) { + *x = (id, k, false); + } + tags.extend(segs.get(h / 2).map(|s| s.1.as_slice()).unwrap_or(&[])); + let b = head(h); + pts.push(b); + if end(b) || b == tail(start) { + break; + } + // Continue through the other half-edge leaving b. + match out.get(b).and_then(|l| l.iter().copied().find(|&o| o / 2 != h / 2)) { + Some(o) if !seen.get(o / 2).copied().unwrap_or(true) => h = o, + _ => break, + } + } + tags.sort_unstable(); + tags.dedup(); + curves.push(Curve { points: pts, tags }); + }; + for p in 0..n { + if !end(p) { + continue; + } + for &h in out.get(p).map(Vec::as_slice).unwrap_or(&[]) { + if !seen.get(h / 2).copied().unwrap_or(true) { + walk(h, &mut curves, &mut curve_of, &mut seen); + } + } + } + for s in 0..segs.len() { + if !seen.get(s).copied().unwrap_or(true) { + walk(2 * s, &mut curves, &mut curve_of, &mut seen); + } + } + // Face cycles. + let mut cycle_of = vec![usize::MAX; segs.len() * 2]; + let mut cycles: Vec> = Vec::new(); + for h0 in 0..segs.len() * 2 { + if cycle_of.get(h0).copied() != Some(usize::MAX) { + continue; + } + let id = cycles.len(); + let mut cyc = Vec::new(); + let mut h = h0; + for _ in 0..=segs.len() * 2 { + if let Some(c) = cycle_of.get_mut(h) { + *c = id; + } + cyc.push(h); + h = next(h); + if h == h0 { + break; + } + } + if h != h0 { + return Err(Error::Failed("an arrangement face cycle does not close".into())); + } + cycles.push(cyc); + } + let area = |cyc: &[usize]| -> f64 { + 0.5 * cyc + .iter() + .map(|&h| { + let (a, b) = (at(tail(h)), at(head(h))); + a[0] * b[1] - a[1] * b[0] + }) + .sum::() + }; + // Connected components (union-find over points). + let mut parent: Vec = (0..n).collect(); + fn find(parent: &mut [usize], mut x: usize) -> usize { + for _ in 0..parent.len() { + let p = parent.get(x).copied().unwrap_or(x); + if p == x { + break; + } + let gp = parent.get(p).copied().unwrap_or(p); + if let Some(slot) = parent.get_mut(x) { + *slot = gp; + } + x = gp; + } + x + } + for ((a, b), _) in &segs { + let (ra, rb) = (find(&mut parent, *a), find(&mut parent, *b)); + if ra != rb + && let Some(slot) = parent.get_mut(ra.max(rb)) + { + *slot = ra.min(rb); + } + } + let comp: Vec = (0..n).map(|p| find(&mut parent, p)).collect(); + let areas: Vec = cycles.iter().map(|c| area(c)).collect(); + // Is point p inside the polygon of cycle c (crossing number with exact side tests)? + let inside = |p: Q, cyc: &[usize]| -> Result { + let mut wn = 0i64; + for &h in cyc { + let (a, b) = (at(tail(h)), at(head(h))); + if a[1] <= p[1] { + if b[1] > p[1] && sign(a, b, p)? == Sign::Positive { + wn += 1; + } + } else if b[1] <= p[1] && sign(a, b, p)? == Sign::Negative { + wn -= 1; + } + } + Ok(wn != 0) + }; + let mut holes: Vec> = vec![Vec::new(); cycles.len()]; + for (c, cyc) in cycles.iter().enumerate() { + if areas.get(c).copied().unwrap_or(0.0) >= 0.0 { + continue; + } + let Some(&h) = cyc.first() else { continue }; + let p = tail(h); + let mut best: Option<(f64, usize)> = None; + for (o, other) in cycles.iter().enumerate() { + let a = areas.get(o).copied().unwrap_or(0.0); + let same = other.first().is_some_and(|&x| comp.get(tail(x)) == comp.get(p)); + if a <= 0.0 || same || best.is_some_and(|(ba, _)| a >= ba) { + continue; + } + if inside(at(p), other)? { + best = Some((a, o)); + } + } + if let Some((_, o)) = best + && let Some(hs) = holes.get_mut(o) + { + hs.push(c); + } + } + // Cycles as loops of curves: one entry per curve use, in cycle order. A use starts at the + // curve's first segment going forwards, or at its last going backwards. + let starts_use = |h: usize| match curve_of.get(h) { + Some(&(c, k, true)) => k == 0 && c != usize::MAX, + Some(&(c, k, false)) => curves.get(c).is_some_and(|cv| k + 2 == cv.points.len()), + None => false, + }; + let as_loop = |cyc: &[usize]| -> Vec<(usize, bool)> { + let start = cyc.iter().position(|&h| starts_use(h)).unwrap_or(0); + (0..cyc.len()) + .filter_map(|k| cyc.get((start + k) % cyc.len()).copied()) + .filter(|&h| starts_use(h)) + .filter_map(|h| curve_of.get(h).map(|&(c, _, fwd)| (c, fwd))) + .collect() + }; + let mut regions = Vec::new(); + for (c, cyc) in cycles.iter().enumerate() { + let a = areas.get(c).copied().unwrap_or(0.0); + if a <= 0.0 { + continue; + } + let mut loops = vec![as_loop(cyc)]; + for &h in holes.get(c).map(Vec::as_slice).unwrap_or(&[]) { + loops.push(as_loop(cycles.get(h).map(Vec::as_slice).unwrap_or(&[]))); + } + regions.push(Region { loops }); + } + Ok(Arrangement { points, curves, regions }) +} + +#[cfg(test)] +mod tests { + use super::*; + + /// A region's area: its loops' signed areas summed (holes run clockwise). + fn area(a: &Arrangement, r: &Region) -> f64 { + r.loops + .iter() + .map(|lp| { + let pts: Vec = lp + .iter() + .flat_map(|&(c, fwd)| { + let mut v: Vec = a.curves[c].points.clone(); + if !fwd { + v.reverse(); + } + v.pop(); + v + }) + .map(|p| a.points[p]) + .collect(); + 0.5 * (0..pts.len()).map(|k| pts[k][0] * pts[(k + 1) % pts.len()][1] - pts[(k + 1) % pts.len()][0] * pts[k][1]).sum::() + }) + .sum() + } + + fn square(x0: f64, y0: f64, s: f64) -> Vec { + vec![[x0, y0], [x0 + s, y0], [x0 + s, y0 + s], [x0, y0 + s], [x0, y0]] + } + + #[test] + fn square_with_a_hole() { + let polys = vec![(square(0.0, 0.0, 10.0), 0), (square(4.0, 4.0, 2.0), 1)]; + let verts: Vec = square(0.0, 0.0, 10.0).into_iter().chain(square(4.0, 4.0, 2.0)).collect(); + let a = build(&verts, &polys, 1e-9, 0.0).unwrap(); + assert_eq!(a.points.len(), 8); + assert_eq!(a.curves.len(), 8); + assert_eq!(a.regions.len(), 2); + let ring = a.regions.iter().find(|r| r.loops.len() == 2).unwrap(); + assert!((area(&a, ring) - 96.0).abs() < 1e-9, "{}", area(&a, ring)); + assert_eq!(ring.loops[0].len(), 4); + assert_eq!(ring.loops[1].len(), 4); + let inner = a.regions.iter().find(|r| r.loops.len() == 1).unwrap(); + assert!((area(&a, inner) - 4.0).abs() < 1e-9); + } + + #[test] + fn t_junctions_and_crossings_split() { + // A square cut by a line from the middle of the bottom to the middle of the top, and a + // line crossing it. + let polys = vec![(square(0.0, 0.0, 10.0), 0), (vec![[5.0, 0.0], [5.0, 10.0]], 1), (vec![[0.0, 5.0], [10.0, 5.0]], 2)]; + let a = build(&[], &polys, 1e-9, 0.0).unwrap(); + assert_eq!(a.regions.len(), 4); + for r in &a.regions { + assert!((area(&a, r) - 25.0).abs() < 1e-9); + assert_eq!(r.loops.len(), 1); + } + } + + #[test] + fn coincident_copies_merge_and_tags_join() { + let mut copy = square(0.0, 0.0, 10.0); + copy[1][1] += 1e-12; + let polys = vec![(square(0.0, 0.0, 10.0), 0), (copy, 1)]; + let a = build(&[[0.0, 0.0]], &polys, 1e-9, 1e-6).unwrap(); + assert_eq!(a.regions.len(), 1); + assert_eq!(a.curves.len(), 1); + assert_eq!(a.curves[0].tags, vec![0, 1]); + } + + #[test] + fn hostile_input_is_an_error() { + assert!(build(&[[f64::NAN, 0.0]], &[], 1e-9, 0.0).is_err()); + assert!(build(&[], &[], 0.0, 0.0).is_err()); + assert!(build(&[], &[(vec![[0.0, 0.0], [1e300, 0.0]], 0)], 1e-9, 0.0).is_err()); + } +} diff --git a/crates/mesh/src/hex/brep.rs b/crates/mesh/src/hex/brep.rs new file mode 100644 index 0000000..04bcc04 --- /dev/null +++ b/crates/mesh/src/hex/brep.rs @@ -0,0 +1,302 @@ +//! The B-rep topology recovered from a face-tagged boundary surface: faces (their triangles, +//! geometry and boundary loops), edges (chains of face-boundary feature edges between the same +//! two faces) and vertices (where three or more faces meet, or where the pair of faces along a +//! chain changes). +//! +//! The hex meshers reason on this topology the way B-rep sweepers do (Blacker, *The Cooper +//! tool*, Proc. 5th International Meshing Roundtable, 1996): which faces are sources, targets and +//! linking faces, and which edges bound them. + +use std::collections::BTreeMap; + +use super::geom::{Surface, add, dist, scale, unit}; +use crate::surface::{BoundarySurface, FeatureKind}; +use crate::vec3::{P, cross, dot, norm, sub}; +use crate::{Error, Result}; + +/// A B-rep face. +#[derive(Debug, Clone)] +pub(crate) struct Face { + /// The B-rep face tag (`tri_face`). + pub id: u32, + pub triangles: Vec, + /// Surface vertices of the face (sorted). + pub vertices: Vec, + pub area: f64, + /// Its analytic surface, when the geometry gave one. + pub surface: Option, + /// When the face is planar: a point and the outward unit normal. + pub plane: Option<(P, P)>, + /// Boundary loops: oriented edges `(edge, forward)` with the face on the left. + pub loops: Vec>, +} + +/// A B-rep edge: a chain of surface vertices along the boundary between two faces. +#[derive(Debug, Clone)] +pub(crate) struct Edge { + /// In order; a closed edge repeats its first point at the end. + pub points: Vec, + /// Face indices: `faces[0]` traverses the edge forwards (has it on its counter-clockwise + /// boundary seen from outside), `faces[1]` backwards. + pub faces: [usize; 2], +} + +impl Edge { + pub fn first(&self) -> u32 { + self.points.first().copied().unwrap_or(0) + } + pub fn last(&self) -> u32 { + self.points.last().copied().unwrap_or(0) + } + pub fn closed(&self) -> bool { + self.first() == self.last() + } + /// The other face of the edge. + pub fn other(&self, face: usize) -> usize { + if self.faces[0] == face { self.faces[1] } else { self.faces[0] } + } +} + +pub(crate) struct Brep<'a> { + pub surface: &'a BoundarySurface, + pub faces: Vec, + pub edges: Vec, + /// Whether each surface vertex is a B-rep vertex. + pub is_vertex: Vec, + /// The two triangles of each face-boundary segment, by its vertices (low, high). + edge_triangles: BTreeMap<(u32, u32), [u32; 2]>, + /// Length of the bounding-box diagonal. + pub diagonal: f64, +} + +/// Relative geometric tolerance (of the bounding-box diagonal). +pub(crate) const RELATIVE_EPS: f64 = 1e-7; + +impl<'a> Brep<'a> { + /// Recover the topology. `surfaces[tag]` is the analytic surface of B-rep face `tag`, if any. + pub(crate) fn new(surface: &'a BoundarySurface, surfaces: &[Option]) -> Result { + let pts = surface.points(); + let mut lo = [f64::INFINITY; 3]; + let mut hi = [f64::NEG_INFINITY; 3]; + for p in pts { + for ((l, h), x) in lo.iter_mut().zip(hi.iter_mut()).zip(p) { + *l = l.min(*x); + *h = h.max(*x); + } + } + let diagonal = dist(lo, hi); + if !(diagonal > 0.0 && diagonal.is_finite()) { + return Err(Error::Invalid("the boundary surface is empty or degenerate".into())); + } + let eps = RELATIVE_EPS * diagonal; + // Faces in tag order. + let mut by_tag: BTreeMap> = BTreeMap::new(); + for (t, &f) in surface.tri_face().iter().enumerate() { + by_tag.entry(f).or_default().push(t as u32); + } + let index: BTreeMap = by_tag.keys().enumerate().map(|(i, &f)| (f, i)).collect(); + let mut faces = Vec::with_capacity(by_tag.len()); + for (&id, tris) in &by_tag { + let mut verts = Vec::with_capacity(tris.len() * 3); + let (mut area, mut nsum) = (0.0, [0.0; 3]); + for &t in tris { + let tri = surface.triangles().get(t as usize).copied().unwrap_or([0; 3]); + verts.extend(tri); + let [a, b, c] = tri.map(|v| pts.get(v as usize).copied().unwrap_or([0.0; 3])); + let n = cross(sub(b, a), sub(c, a)); + area += 0.5 * norm(n); + nsum = add(nsum, scale(n, 0.5)); + } + verts.sort_unstable(); + verts.dedup(); + let analytic = surfaces.get(id as usize).copied().flatten().and_then(|s| s.checked()); + let normal = unit(nsum).unwrap_or([0.0; 3]); + let plane = planarity(pts, &verts, normal, analytic, eps); + faces.push(Face { id, triangles: tris.clone(), vertices: verts, area, surface: analytic, plane, loops: Vec::new() }); + } + let face_of = |t: u32| surface.tri_face().get(t as usize).and_then(|f| index.get(f)).copied(); + // Oriented boundary segments: `from → to` with face `left` traversing it. + struct Seg { + from: u32, + to: u32, + left: usize, + right: usize, + } + let mut segs = Vec::new(); + let mut edge_triangles = BTreeMap::new(); + for fe in surface.feature_edges().iter().filter(|e| e.kind == FeatureKind::FaceBoundary) { + let (Some(left), Some(right)) = (face_of(fe.triangles[0]), face_of(fe.triangles[1])) else { + return Err(Error::Failed("feature edge without faces".into())); + }; + // Orient each segment as the lower-indexed face traverses it, so that the segments of + // one edge chain head to tail. + let [lo, hi] = fe.vertices; + edge_triangles.insert((lo, hi), fe.triangles); + segs.push(if left < right { Seg { from: lo, to: hi, left, right } } else { Seg { from: hi, to: lo, left: right, right: left } }); + } + let mut incident: Vec> = vec![Vec::new(); pts.len()]; + for (i, s) in segs.iter().enumerate() { + for v in [s.from, s.to] { + incident.get_mut(v as usize).ok_or_else(|| Error::Failed("vertex index".into()))?.push(i); + } + } + let pair = |s: &Seg| (s.left, s.right); + let mut is_vertex = vec![false; pts.len()]; + for (v, inc) in incident.iter().enumerate() { + let vertex = match inc.as_slice() { + [] => false, + [a, b] => match (segs.get(*a), segs.get(*b)) { + (Some(a), Some(b)) => { + pair(a) != pair(b) || !((a.to == v as u32 && b.from == v as u32) || (b.to == v as u32 && a.from == v as u32)) + } + _ => true, + }, + _ => true, + }; + if let Some(x) = is_vertex.get_mut(v) { + *x = vertex; + } + } + // The segment leaving `v` along the chain that arrived through segment `s`. + let next = |s: usize, v: u32| -> Option { + incident.get(v as usize)?.iter().copied().find(|&o| o != s && segs.get(o).is_some_and(|x| x.from == v)) + }; + let mut visited = vec![false; segs.len()]; + let mut edges = Vec::new(); + for pass in 0..2 { + for start in 0..segs.len() { + let Some(s0) = segs.get(start) else { continue }; + if visited.get(start).copied().unwrap_or(true) || (pass == 0 && !is_vertex.get(s0.from as usize).copied().unwrap_or(false)) { + continue; + } + let mut points = vec![s0.from]; + let mut s = start; + // A chain has at most one step per segment. + for _ in 0..=segs.len() { + if let Some(x) = visited.get_mut(s) { + *x = true; + } + let Some(seg) = segs.get(s) else { break }; + points.push(seg.to); + if is_vertex.get(seg.to as usize).copied().unwrap_or(true) || seg.to == s0.from { + break; + } + match next(s, seg.to) { + Some(n) if !visited.get(n).copied().unwrap_or(true) => s = n, + _ => break, + } + } + edges.push(Edge { points, faces: [s0.left, s0.right] }); + } + } + // Loops of each face from its oriented edges. + let mut uses: Vec> = vec![Vec::new(); faces.len()]; + for (e, edge) in edges.iter().enumerate() { + if let Some(u) = uses.get_mut(edge.faces[0]) { + u.push((e, true)); + } + if let Some(u) = uses.get_mut(edge.faces[1]) { + u.push((e, false)); + } + } + for (f, u) in uses.iter().enumerate() { + let ends = |(e, fwd): (usize, bool)| -> (u32, u32) { + let edge = edges.get(e); + let (a, b) = (edge.map_or(0, Edge::first), edge.map_or(0, Edge::last)); + if fwd { (a, b) } else { (b, a) } + }; + let mut used = vec![false; u.len()]; + let mut loops = Vec::new(); + for start in 0..u.len() { + if used.get(start).copied().unwrap_or(true) { + continue; + } + let Some(&first) = u.get(start) else { continue }; + let origin = ends(first).0; + let mut lp = Vec::new(); + let mut cur = start; + for _ in 0..=u.len() { + if let Some(x) = used.get_mut(cur) { + *x = true; + } + let Some(&oe) = u.get(cur) else { break }; + lp.push(oe); + let end = ends(oe).1; + if end == origin { + break; + } + match (0..u.len()).find(|&k| !used.get(k).copied().unwrap_or(true) && u.get(k).is_some_and(|&x| ends(x).0 == end)) { + Some(k) => cur = k, + None => return Err(Error::Failed(format!("face {} has an open boundary loop", faces.get(f).map_or(0, |x| x.id)))), + } + } + loops.push(lp); + } + if let Some(face) = faces.get_mut(f) { + face.loops = loops; + } + } + Ok(Brep { surface, faces, edges, is_vertex, edge_triangles, diagonal }) + } + + pub(crate) fn point(&self, v: u32) -> P { + self.surface.points().get(v as usize).copied().unwrap_or([0.0; 3]) + } + + /// The points of an oriented edge. + pub(crate) fn edge_points(&self, e: usize, forward: bool) -> Vec

{ + let mut pts: Vec

= self.edges.get(e).map(|x| x.points.iter().map(|&v| self.point(v)).collect()).unwrap_or_default(); + if !forward { + pts.reverse(); + } + pts + } + + /// The oriented edge's surface vertices. + pub(crate) fn edge_vertices(&self, e: usize, forward: bool) -> Vec { + let mut v = self.edges.get(e).map(|x| x.points.clone()).unwrap_or_default(); + if !forward { + v.reverse(); + } + v + } + + /// The vertex of face `f`'s triangle on edge `e`'s first segment that is not on the edge. + pub(crate) fn third_vertex(&self, e: usize, f: usize) -> Option { + let edge = self.edges.get(e)?; + let (a, b) = (*edge.points.first()?, *edge.points.get(1)?); + let tri = self.edge_triangles.get(&(a.min(b), a.max(b)))?.iter().copied().find(|&t| self.face_of_triangle(t) == Some(f))?; + self.surface.triangles().get(tri as usize)?.iter().copied().find(|&v| v != a && v != b) + } + + /// The face's name for messages: its tag. + pub(crate) fn face_tag(&self, f: usize) -> u32 { + self.faces.get(f).map_or(u32::MAX, |x| x.id) + } + + /// The faces that use the triangle's tag (index into `faces`). + pub(crate) fn face_of_triangle(&self, t: u32) -> Option { + let tag = *self.surface.tri_face().get(t as usize)?; + self.faces.binary_search_by_key(&tag, |f| f.id).ok() + } +} + +/// A point and outward normal of a planar face: the analytic plane, or the plane fitted through +/// the face's vertices when they all lie within `eps` of it. +fn planarity(pts: &[P], verts: &[u32], normal: P, analytic: Option, eps: f64) -> Option<(P, P)> { + let first = pts.get(*verts.first()? as usize).copied()?; + let (point, n) = match analytic { + Some(Surface::Plane { point, normal: an }) => (point, if dot(an, normal) < 0.0 { scale(an, -1.0) } else { an }), + Some(_) => return None, + None => { + if norm(normal) == 0.0 { + return None; + } + (first, normal) + } + }; + let flat = verts.iter().all(|&v| pts.get(v as usize).is_some_and(|p| dot(sub(*p, point), n).abs() <= eps.max(1e-9 * norm(sub(*p, point))))); + // An analytic plane is trusted within a looser bound: its vertices lie on it to the kernel's + // tolerance, not to ours. + (flat || matches!(analytic, Some(Surface::Plane { .. }))).then_some((point, n)) +} diff --git a/crates/mesh/src/hex/geom.rs b/crates/mesh/src/hex/geom.rs new file mode 100644 index 0000000..e17652f --- /dev/null +++ b/crates/mesh/src/hex/geom.rs @@ -0,0 +1,408 @@ +//! Geometry the hex meshers place nodes on: the analytic surface of a B-rep face (exact closest +//! point), the triangles of the boundary surface (closest point through a bounding-volume tree, +//! for free-form faces), and the inside test by the generalized winding number. +//! +//! Closest points on the elementary surfaces are closed form. The point-to-triangle distance is +//! the Voronoi-region method of C. Ericson, *Real-Time Collision Detection*, Morgan Kaufmann 2005, +//! §5.1.5. The winding number is the sum of the solid angles the triangles subtend (A. Van +//! Oosterom, J. Strackee, *The solid angle of a plane triangle*, IEEE Trans. Biomed. Eng. +//! 30(2):125–126, 1983; A. Jacobson, L. Kavan, O. Sorkine-Hornung, *Robust inside-outside +//! segmentation using generalized winding numbers*, ACM TOG 32(4):33, 2013). +//! +//! These are floating-point measures for placing nodes and choosing tags; no topological decision +//! of the meshers rests on a bare sign here except the inside test, which is only asked at points +//! far from the boundary (region centres), where the winding number is close to 0 or 1. + +use crate::surface::BoundarySurface; +use crate::vec3::{P, cross, dot, norm, sub}; +use crate::{Error, Result}; + +pub(crate) fn add(a: P, b: P) -> P { + [a[0] + b[0], a[1] + b[1], a[2] + b[2]] +} +pub(crate) fn scale(a: P, s: f64) -> P { + [a[0] * s, a[1] * s, a[2] * s] +} +pub(crate) fn dist(a: P, b: P) -> f64 { + norm(sub(a, b)) +} +/// The unit vector along `a`, or `None` when it has no direction. +pub(crate) fn unit(a: P) -> Option

{ + let n = norm(a); + (n > 0.0 && n.is_finite()).then(|| scale(a, 1.0 / n)) +} + +/// The analytic surface a B-rep face lies on. Directions need not be unit length; they are +/// normalized when the surface is checked ([`Surface::checked`]). Lengths are model units and +/// angles radians. +#[derive(Debug, Clone, Copy, PartialEq)] +pub enum Surface { + Plane { + point: [f64; 3], + normal: [f64; 3], + }, + Cylinder { + origin: [f64; 3], + axis: [f64; 3], + radius: f64, + }, + /// The apex, the axis pointing into the cone, and the half-angle (0, π/2). + Cone { + apex: [f64; 3], + axis: [f64; 3], + half_angle: f64, + }, + Sphere { + center: [f64; 3], + radius: f64, + }, + /// The centre of the ring, its axis, the ring (major) and tube (minor) radii. + Torus { + center: [f64; 3], + axis: [f64; 3], + major: f64, + minor: f64, + }, +} + +impl Surface { + /// The surface with unit directions, or `None` when a number is not finite or a direction or + /// radius is degenerate (such a face is treated as free-form). + pub(crate) fn checked(&self) -> Option { + let fin = |p: &[f64; 3]| p.iter().all(|x| x.is_finite()); + let pos = |r: f64| r > 0.0 && r.is_finite(); + Some(match *self { + Surface::Plane { point, normal } if fin(&point) => Surface::Plane { point, normal: unit(normal)? }, + Surface::Cylinder { origin, axis, radius } if fin(&origin) && pos(radius) => Surface::Cylinder { origin, axis: unit(axis)?, radius }, + Surface::Cone { apex, axis, half_angle } if fin(&apex) && half_angle > 0.0 && half_angle < std::f64::consts::FRAC_PI_2 => { + Surface::Cone { apex, axis: unit(axis)?, half_angle } + } + Surface::Sphere { center, radius } if fin(¢er) && pos(radius) => Surface::Sphere { center, radius }, + Surface::Torus { center, axis, major, minor } if fin(¢er) && pos(major) && pos(minor) => { + Surface::Torus { center, axis: unit(axis)?, major, minor } + } + _ => return None, + }) + } + + /// The closest point on the (untrimmed) surface, or `None` where it is not unique (on a + /// cylinder's axis, at a sphere's centre, …). + pub(crate) fn project(&self, p: P) -> Option

{ + let q = match *self { + Surface::Plane { point, normal } => sub(p, scale(normal, dot(sub(p, point), normal))), + Surface::Cylinder { origin, axis, radius } => { + let w = sub(p, origin); + let h = dot(w, axis); + let radial = unit(sub(w, scale(axis, h)))?; + add(add(origin, scale(axis, h)), scale(radial, radius)) + } + Surface::Cone { apex, axis, half_angle } => { + let w = sub(p, apex); + let h = dot(w, axis); + let radial = unit(sub(w, scale(axis, h)))?; + let rho = norm(sub(w, scale(axis, h))); + let (c, s) = (half_angle.cos(), half_angle.sin()); + let t = (h * c + rho * s).max(0.0); + add(apex, add(scale(axis, t * c), scale(radial, t * s))) + } + Surface::Sphere { center, radius } => add(center, scale(unit(sub(p, center))?, radius)), + Surface::Torus { center, axis, major, minor } => { + let w = sub(p, center); + let h = dot(w, axis); + let ring = add(center, scale(unit(sub(w, scale(axis, h)))?, major)); + add(ring, scale(unit(sub(p, ring))?, minor)) + } + }; + q.iter().all(|x| x.is_finite()).then_some(q) + } +} + +/// The closest point to `p` on triangle `abc`. +pub(crate) fn closest_on_triangle(p: P, a: P, b: P, c: P) -> P { + let (ab, ac, ap) = (sub(b, a), sub(c, a), sub(p, a)); + let (d1, d2) = (dot(ab, ap), dot(ac, ap)); + if d1 <= 0.0 && d2 <= 0.0 { + return a; + } + let bp = sub(p, b); + let (d3, d4) = (dot(ab, bp), dot(ac, bp)); + if d3 >= 0.0 && d4 <= d3 { + return b; + } + let vc = d1 * d4 - d3 * d2; + if vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0 { + return add(a, scale(ab, d1 / (d1 - d3))); + } + let cp = sub(p, c); + let (d5, d6) = (dot(ab, cp), dot(ac, cp)); + if d6 >= 0.0 && d5 <= d6 { + return c; + } + let vb = d5 * d2 - d1 * d6; + if vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0 { + return add(a, scale(ac, d2 / (d2 - d6))); + } + let va = d3 * d6 - d5 * d4; + if va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0 { + return add(b, scale(sub(c, b), (d4 - d3) / ((d4 - d3) + (d5 - d6)))); + } + let den = va + vb + vc; + if !(den.abs() > 0.0) { + return a; + } + add(a, add(scale(ab, vb / den), scale(ac, vc / den))) +} + +#[derive(Debug, Clone, Copy)] +struct Bounds { + lo: P, + hi: P, +} + +impl Bounds { + const EMPTY: Bounds = Bounds { lo: [f64::INFINITY; 3], hi: [f64::NEG_INFINITY; 3] }; + fn add(&mut self, p: P) { + for ((l, h), x) in self.lo.iter_mut().zip(self.hi.iter_mut()).zip(p) { + *l = l.min(x); + *h = h.max(x); + } + } + /// Squared distance from `p` to the box (0 inside). + fn dist2(&self, p: P) -> f64 { + (0..3).map(|i| (self.lo[i] - p[i]).max(0.0).max(p[i] - self.hi[i]).powi(2)).sum() + } +} + +enum Node { + Leaf { bounds: Bounds, first: usize, count: usize }, + Inner { bounds: Bounds, left: usize, right: usize }, +} + +impl Node { + fn bounds(&self) -> &Bounds { + match self { + Node::Leaf { bounds, .. } | Node::Inner { bounds, .. } => bounds, + } + } +} + +const LEAF_SIZE: usize = 8; +/// Nodes a closest-point query may visit before it settles for the best triangle found. +const MAX_QUERY_VISITS: usize = 1 << 20; + +/// A bounding-volume tree over the boundary surface's triangles (median splits along the longest +/// axis), for closest-point queries restricted to a set of B-rep faces. +pub(crate) struct TriangleTree<'a> { + surface: &'a BoundarySurface, + order: Vec, + nodes: Vec, + normals: Vec

, +} + +/// The closest triangle found by [`TriangleTree::closest`]. +#[derive(Debug, Clone, Copy)] +pub(crate) struct Hit { + pub triangle: u32, + pub point: P, + pub dist: f64, +} + +impl<'a> TriangleTree<'a> { + pub(crate) fn new(surface: &'a BoundarySurface) -> Result { + let n = surface.triangles().len(); + let mut tree = TriangleTree { surface, order: (0..n as u32).collect(), nodes: Vec::new(), normals: Vec::with_capacity(n) }; + let mut centres = Vec::with_capacity(n); + for t in 0..n as u32 { + let [a, b, c] = tree.corners(t)?; + centres.push(scale(add(add(a, b), c), 1.0 / 3.0)); + tree.normals.push(unit(cross(sub(b, a), sub(c, a))).unwrap_or([0.0; 3])); + } + if n > 0 { + tree.build(0, n, ¢res, 0)?; + } + Ok(tree) + } + + pub(crate) fn corners(&self, t: u32) -> Result<[P; 3]> { + let pts = self.surface.points(); + let tri = self.surface.triangles().get(t as usize).ok_or_else(|| Error::Failed(format!("no triangle {t}")))?; + let p = |v: u32| pts.get(v as usize).copied().ok_or_else(|| Error::Failed(format!("no vertex {v}"))); + Ok([p(tri[0])?, p(tri[1])?, p(tri[2])?]) + } + + /// The unit normal of triangle `t` (outward), zero when degenerate. + pub(crate) fn normal(&self, t: u32) -> P { + self.normals.get(t as usize).copied().unwrap_or([0.0; 3]) + } + + fn build(&mut self, first: usize, count: usize, centres: &[P], depth: usize) -> Result { + let mut bounds = Bounds::EMPTY; + let mut cb = Bounds::EMPTY; + for &t in self.order.get(first..first + count).unwrap_or(&[]) { + for p in self.corners(t)? { + bounds.add(p); + } + cb.add(centres.get(t as usize).copied().unwrap_or([0.0; 3])); + } + let id = self.nodes.len(); + // Depth is bounded by log2 of the count (median splits); the cap only guards the stack. + if count <= LEAF_SIZE || depth > 64 { + self.nodes.push(Node::Leaf { bounds, first, count }); + return Ok(id); + } + let ext = [cb.hi[0] - cb.lo[0], cb.hi[1] - cb.lo[1], cb.hi[2] - cb.lo[2]]; + let axis = if ext[0] >= ext[1] && ext[0] >= ext[2] { + 0 + } else if ext[1] >= ext[2] { + 1 + } else { + 2 + }; + let key = |t: &u32| centres.get(*t as usize).map_or(0.0, |c| c[axis]); + if let Some(slice) = self.order.get_mut(first..first + count) { + slice.sort_by(|a, b| key(a).total_cmp(&key(b)).then(a.cmp(b))); + } + self.nodes.push(Node::Leaf { bounds, first, count }); + let half = count / 2; + let left = self.build(first, half, centres, depth + 1)?; + let right = self.build(first + half, count - half, centres, depth + 1)?; + if let Some(n) = self.nodes.get_mut(id) { + *n = Node::Inner { bounds, left, right }; + } + Ok(id) + } + + /// The closest triangle to `p` among those `accept` admits (by triangle index). + pub(crate) fn closest(&self, p: P, accept: &dyn Fn(u32) -> bool) -> Option { + let mut best: Option = None; + let mut best2 = f64::INFINITY; + let mut stack = vec![0usize]; + let mut visits = 0usize; + while let Some(id) = stack.pop() { + visits += 1; + if visits > MAX_QUERY_VISITS { + break; + } + let Some(node) = self.nodes.get(id) else { continue }; + if node.bounds().dist2(p) > best2 { + continue; + } + match node { + Node::Leaf { first, count, .. } => { + for &t in self.order.get(*first..*first + *count).unwrap_or(&[]) { + if !accept(t) { + continue; + } + let Ok([a, b, c]) = self.corners(t) else { continue }; + let q = closest_on_triangle(p, a, b, c); + let d2 = dot(sub(p, q), sub(p, q)); + if d2 < best2 || (d2 == best2 && best.is_some_and(|h| t < h.triangle)) { + best2 = d2; + best = Some(Hit { triangle: t, point: q, dist: d2.sqrt() }); + } + } + } + Node::Inner { left, right, .. } => { + let (dl, dr) = ( + self.nodes.get(*left).map_or(f64::INFINITY, |n| n.bounds().dist2(p)), + self.nodes.get(*right).map_or(f64::INFINITY, |n| n.bounds().dist2(p)), + ); + // Visit the nearer child first (pushed last). + if dl <= dr { + stack.push(*right); + stack.push(*left); + } else { + stack.push(*left); + stack.push(*right); + } + } + } + } + best + } +} + +/// The most triangle evaluations the inside tests of one mesh may spend. +pub(crate) const MAX_INSIDE_WORK: u64 = 2_000_000_000; + +/// Inside tests by the generalized winding number, with a shared work budget. +pub(crate) struct Inside<'a> { + surface: &'a BoundarySurface, + work: u64, +} + +impl<'a> Inside<'a> { + pub(crate) fn new(surface: &'a BoundarySurface) -> Self { + Inside { surface, work: 0 } + } + + /// The generalized winding number of the boundary at `x` (1 inside, 0 outside). + pub(crate) fn winding(&mut self, x: P) -> Result { + let n = self.surface.triangles().len() as u64; + self.work = self.work.saturating_add(n); + if self.work > MAX_INSIDE_WORK { + return Err(Error::Limit(format!("inside tests exceed {MAX_INSIDE_WORK} triangle evaluations"))); + } + let pts = self.surface.points(); + let mut sum = 0.0; + for t in self.surface.triangles() { + let (Some(a), Some(b), Some(c)) = (pts.get(t[0] as usize), pts.get(t[1] as usize), pts.get(t[2] as usize)) else { continue }; + let (a, b, c) = (sub(*a, x), sub(*b, x), sub(*c, x)); + let (la, lb, lc) = (norm(a), norm(b), norm(c)); + let den = la * lb * lc + dot(a, b) * lc + dot(a, c) * lb + dot(b, c) * la; + sum += 2.0 * dot(a, cross(b, c)).atan2(den); + } + Ok(sum / (4.0 * std::f64::consts::PI)) + } + + /// Whether `x` is inside the solid. A winding number away from both 0 and 1 means `x` is on + /// or next to the boundary, which callers avoid; it is an error rather than a guess. + pub(crate) fn contains(&mut self, x: P) -> Result { + let w = self.winding(x)?; + if (w - 1.0).abs() < 0.25 { + Ok(true) + } else if w.abs() < 0.25 { + Ok(false) + } else { + Err(Error::Failed(format!("cannot tell whether {x:?} is inside the solid (winding number {w:.3})"))) + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn close(a: P, b: P) -> bool { + dist(a, b) < 1e-12 + } + + #[test] + fn analytic_projections() { + let cyl = Surface::Cylinder { origin: [1.0, 0.0, 0.0], axis: [0.0, 0.0, 2.0], radius: 2.0 }.checked().unwrap(); + assert!(close(cyl.project([4.0, 0.0, 5.0]).unwrap(), [3.0, 0.0, 5.0])); + assert!(cyl.project([1.0, 0.0, 7.0]).is_none()); + let sph = Surface::Sphere { center: [0.0; 3], radius: 3.0 }.checked().unwrap(); + assert!(close(sph.project([0.0, 0.0, 1.0]).unwrap(), [0.0, 0.0, 3.0])); + let tor = Surface::Torus { center: [0.0; 3], axis: [0.0, 0.0, 1.0], major: 10.0, minor: 2.0 }.checked().unwrap(); + assert!(close(tor.project([13.0, 0.0, 0.0]).unwrap(), [12.0, 0.0, 0.0])); + assert!(close(tor.project([10.0, 0.0, 5.0]).unwrap(), [10.0, 0.0, 2.0])); + let cone = Surface::Cone { apex: [0.0; 3], axis: [0.0, 0.0, 1.0], half_angle: std::f64::consts::FRAC_PI_4 }.checked().unwrap(); + let q = cone.project([2.0, 0.0, 0.0]).unwrap(); + assert!(close(q, [1.0, 0.0, 1.0]), "{q:?}"); + let plane = Surface::Plane { point: [0.0, 0.0, 1.0], normal: [0.0, 0.0, 5.0] }.checked().unwrap(); + assert!(close(plane.project([3.0, 4.0, 9.0]).unwrap(), [3.0, 4.0, 1.0])); + assert!(Surface::Sphere { center: [0.0; 3], radius: f64::NAN }.checked().is_none()); + assert!(Surface::Plane { point: [0.0; 3], normal: [0.0; 3] }.checked().is_none()); + } + + #[test] + fn triangle_closest_point_regions() { + let (a, b, c) = ([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]); + assert!(close(closest_on_triangle([0.2, 0.2, 3.0], a, b, c), [0.2, 0.2, 0.0])); + assert!(close(closest_on_triangle([-1.0, -1.0, 0.0], a, b, c), a)); + assert!(close(closest_on_triangle([2.0, -1.0, 0.0], a, b, c), b)); + assert!(close(closest_on_triangle([0.5, -1.0, 0.0], a, b, c), [0.5, 0.0, 0.0])); + assert!(close(closest_on_triangle([1.0, 1.0, 0.0], a, b, c), [0.5, 0.5, 0.0])); + } +} diff --git a/crates/mesh/src/hex/interval.rs b/crates/mesh/src/hex/interval.rs new file mode 100644 index 0000000..cb604b7 --- /dev/null +++ b/crates/mesh/src/hex/interval.rs @@ -0,0 +1,229 @@ +//! Interval assignment: the number of mesh divisions of every curve, chosen together so that +//! every face mesher downstream can succeed with conforming meshes on shared curves. +//! +//! The problem is the integer program of Tam and Armstrong (*Finite element mesh control by +//! integer programming*, IJNME 36(15):2581–2605, 1993) and Mitchell (*High fidelity interval +//! assignment*, Proc. 6th International Meshing Roundtable, 1997; IJCGA 10(4):399–415, 2000): +//! integer counts `x_i ≥ 1` close to goals `g_i` (the curve length over the target size), subject to +//! - **equalities** `Σ_A x = Σ_B x`: opposite sides of a mapped face or of a sweep's linking +//! face; +//! - **parity**: `Σ_L x` even on every loop of a face meshed by the free all-quad mesher (a +//! multiple of 4 around a disk meshed as an O-grid); +//! - **minimum sums** `Σ_L x ≥ m`: a loop needs enough points to be meshed. +//! +//! Mitchell's *incremental interval assignment* (Mitchell, *Incremental interval assignment by +//! integer linear algebra*, International Meshing Roundtable 2021; Computer-Aided Design 158:103485, 2023) starts from a +//! rounded solution and repairs violated constraints by unit steps chosen by a priority on the +//! change in relative error. This module uses that scheme with increments only: start from the +//! rounded goals, then repeatedly take the first violated constraint and raise the variable whose +//! increase costs least in squared relative error `((x − g)/g)²`, preferring, for parity, variables that +//! no equality couples to other faces. Increments only make the counts never coarser than asked +//! and make termination provable for feasible systems up to the step cap; an infeasible system +//! (one that only grows) hits the cap and is reported, never looped on. + +use crate::{Error, Result}; + +/// The constraint system. +#[derive(Debug, Clone, Default)] +pub(crate) struct Intervals { + /// Goal count of each variable (curve length over target size). + pub goals: Vec, + /// Lower bound of each variable (at least 1). + pub mins: Vec, + /// `Σ lhs = Σ rhs` (variables may repeat). + pub equal: Vec<(Vec, Vec)>, + /// `Σ vars` a multiple of the modulus (2: even loops for the free mesher; 4: a disk's + /// O-grid). + pub multiple: Vec<(Vec, usize)>, + /// `Σ vars ≥ min`. + pub at_least: Vec<(Vec, usize)>, +} + +/// The largest count a variable may reach. +const MAX_COUNT: usize = 1 << 20; + +impl Intervals { + /// Add a variable with its goal and minimum; returns its index. + pub(crate) fn var(&mut self, goal: f64, min: usize) -> usize { + self.goals.push(goal); + self.mins.push(min.max(1)); + self.goals.len() - 1 + } + + pub(crate) fn solve(&self) -> Result> { + let n = self.goals.len(); + if self.mins.len() != n { + return Err(Error::Failed("interval assignment: goals and minimums differ in length".into())); + } + let mut x = Vec::with_capacity(n); + for (i, &g) in self.goals.iter().enumerate() { + if !(g >= 0.0 && g.is_finite()) { + return Err(Error::Invalid(format!("interval goal {g} is not a finite, non-negative number"))); + } + if g > MAX_COUNT as f64 { + return Err(Error::Limit(format!("a curve would get {g:.0} divisions; the limit is {MAX_COUNT}"))); + } + x.push((g.round() as usize).max(self.mins.get(i).copied().unwrap_or(1))); + } + let check = |v: &[usize]| v.iter().all(|&i| i < n); + if !self.equal.iter().all(|(a, b)| check(a) && check(b)) + || !self.multiple.iter().all(|(l, m)| check(l) && *m >= 2) + || !self.at_least.iter().all(|(l, _)| check(l)) + { + return Err(Error::Failed("interval assignment: a constraint names an unknown curve".into())); + } + let mut coupled = vec![false; n]; + for (a, b) in &self.equal { + for &i in a.iter().chain(b) { + if let Some(c) = coupled.get_mut(i) { + *c = true; + } + } + } + let error = |i: usize, c: usize| -> f64 { + let g = self.goals.get(i).copied().unwrap_or(1.0).max(1e-3); + // Squared relative deviation: strictly convex, so repeated increments spread over + // the candidates instead of piling onto one. + ((c as f64 - g) / g).powi(2) + }; + let step_cost = |x: &[usize], i: usize| -> f64 { + let c = x.get(i).copied().unwrap_or(0); + error(i, c + 1) - error(i, c) + }; + let total: usize = x.iter().sum(); + let constraints = self.equal.len() + self.multiple.len() + self.at_least.len(); + let max_steps = total.saturating_mul(8).saturating_add(64 * (n + constraints)).saturating_add(1000); + for _ in 0..max_steps { + let sum = |v: &[usize], x: &[usize]| v.iter().map(|&i| x.get(i).copied().unwrap_or(0)).sum::(); + // The candidates to raise and whether raising them is what the constraint needs. + let mut raise: Option> = None; + let mut parity = false; + for (a, b) in &self.equal { + let (sa, sb) = (sum(a, &x), sum(b, &x)); + if sa == sb { + continue; + } + // Raise a variable whose net coefficient is on the smaller side. + let (small, large) = if sa < sb { (a, b) } else { (b, a) }; + let net = |i: usize| small.iter().filter(|&&k| k == i).count() as i64 - large.iter().filter(|&&k| k == i).count() as i64; + let cands: Vec = small.iter().copied().filter(|&i| net(i) > 0).collect(); + if cands.is_empty() { + return Err(Error::Failed("interval assignment: a side constraint cannot be met (the same curves on both sides)".into())); + } + raise = Some(cands); + break; + } + if raise.is_none() { + for (l, m) in &self.at_least { + if sum(l, &x) < *m { + raise = Some(l.clone()); + break; + } + } + } + if raise.is_none() { + for (l, m) in &self.multiple { + if !sum(l, &x).is_multiple_of(*m) { + // Only a variable listed a number of times that is not a multiple of the + // modulus moves the remainder. + let moving: Vec = l.iter().copied().filter(|&i| !l.iter().filter(|&&k| k == i).count().is_multiple_of(*m)).collect(); + if moving.is_empty() { + return Err(Error::Failed("interval assignment: a loop's count cannot reach the multiple it needs".into())); + } + raise = Some(moving); + parity = true; + break; + } + } + } + let Some(cands) = raise else { return Ok(x) }; + let pick = cands + .iter() + .copied() + .min_by(|&a, &b| { + let key = |i: usize| (parity && coupled.get(i).copied().unwrap_or(false), step_cost(&x, i)); + let (ka, kb) = (key(a), key(b)); + ka.0.cmp(&kb.0).then(ka.1.total_cmp(&kb.1)).then(a.cmp(&b)) + }) + .ok_or_else(|| Error::Failed("interval assignment: empty constraint".into()))?; + let v = x.get_mut(pick).ok_or_else(|| Error::Failed("interval index".into()))?; + *v += 1; + if *v > MAX_COUNT { + return Err(Error::Limit(format!("interval assignment drove a curve past {MAX_COUNT} divisions"))); + } + } + Err(Error::Failed(format!( + "interval assignment did not settle in {max_steps} steps: the side and parity constraints of the faces conflict (the region cannot be mapped or swept as decomposed)" + ))) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn rounds_goals_without_constraints() { + let mut p = Intervals::default(); + p.var(2.4, 1); + p.var(0.2, 1); + p.var(7.6, 3); + assert_eq!(p.solve().unwrap(), vec![2, 1, 8]); + } + + #[test] + fn mapped_sides_become_equal_and_loops_even() { + // A rectangle: sides 0 and 2 opposite, 1 and 3 opposite; side 0 shared with a free face + // whose loop is 0 + 4 + 5. + let mut p = Intervals::default(); + for g in [4.2, 1.1, 3.6, 1.4, 2.0, 2.0] { + p.var(g, 1); + } + p.equal.push((vec![0], vec![2])); + p.equal.push((vec![1], vec![3])); + p.multiple.push((vec![0, 4, 5], 2)); + let x = p.solve().unwrap(); + assert_eq!(x[0], x[2]); + assert_eq!(x[1], x[3]); + assert_eq!((x[0] + x[4] + x[5]) % 2, 0); + // The parity fix went to an uncoupled curve. + assert_eq!(x[0], 4); + } + + #[test] + fn sums_and_minimums() { + let mut p = Intervals::default(); + let a = p.var(1.0, 1); + let b = p.var(1.0, 1); + let c = p.var(5.0, 1); + p.equal.push((vec![a, b], vec![c])); + p.at_least.push((vec![a, b, c], 12)); + p.multiple.push((vec![a, b, c], 2)); + let x = p.solve().unwrap(); + assert_eq!(x[a] + x[b], x[c]); + assert!(x[a] + x[b] + x[c] >= 12); + assert_eq!((x[a] + x[b] + x[c]) % 2, 0); + } + + #[test] + fn infeasible_systems_are_errors_not_loops() { + let mut p = Intervals::default(); + let a = p.var(3.0, 1); + let b = p.var(3.0, 1); + // a = a + b has no solution with b ≥ 1. + p.equal.push((vec![a], vec![a, b])); + assert!(p.solve().is_err()); + let mut q = Intervals::default(); + let a = q.var(3.0, 1); + q.multiple.push((vec![a, a], 2)); + q.multiple.push((vec![a, a, a], 2)); + q.multiple.push((vec![a, a, a], 4)); + assert!(q.solve().is_ok()); + let mut r = Intervals::default(); + r.var(f64::NAN, 1); + assert!(r.solve().is_err()); + let mut s = Intervals::default(); + s.var(1e30, 1); + assert!(s.solve().is_err()); + } +} diff --git a/crates/mesh/src/hex/layered.rs b/crates/mesh/src/hex/layered.rs new file mode 100644 index 0000000..7f81b40 --- /dev/null +++ b/crates/mesh/src/hex/layered.rs @@ -0,0 +1,793 @@ +//! Multi-level sweeping along a straight direction (extrusions) or around an axis (revolutions). +//! +//! **Recognition.** In a sweep frame every point has cross-section coordinates `q` and a sweep +//! coordinate `h` (height along the direction, or angle about the axis). A body is sweepable in +//! the frame when every face is a *cap* (a plane across the sweep: normal along the direction, or +//! a half-plane through the axis) or a *wall* (along the sweep: planes and cylinders parallel to +//! the direction; planes across the axis and coaxial cylinders, cones, spheres and tori), and +//! every edge runs either across the sweep (constant `h`) or along it (constant `q`). The caps +//! and the across edges sit at a finite set of *levels*. +//! +//! **Decomposition.** Projecting into the cross-section every edge that starts a wall upward +//! (the lower outline of each wall) and imprinting them on each other gives a planar arrangement +//! ([`super::arrangement`]) whose regions are the columns of the body: over each region the body +//! is a union of whole level intervals. This is many-to-many sweeping by imprinting all sources +//! and targets into one projection plane (Lai, Benzley, White, IJNME 49(1–2):261–275, 2000; +//! Knupp, *Next-generation sweep tool: a method for generating all-hex meshes on two-and-one-half +//! dimensional geometries*, Proc. 7th International Meshing Roundtable, 1998). Which intervals a +//! region fills is decided by the inside test at its centre ([`super::geom::Inside`]). +//! +//! **Meshing.** Interval assignment ([`super::interval`]) gives every curve its count so that +//! four-sided regions are mapped (opposite sides equal) and the others have even loops for the +//! free all-quad mesher; each level interval gets its own layer count. Each region's quads are +//! stacked through its intervals; shared curves and levels share nodes, so the hexahedra +//! conform. Nodes of curves on analytic walls are projected onto them; in a revolution, the +//! rotated nodes then lie exactly on every surface of revolution. + +use std::collections::HashMap; +use std::f64::consts::{FRAC_PI_4, PI, TAU}; + +use super::arrangement::{self, Arrangement, Q}; +use super::brep::{Brep, Face}; +use super::geom::{Inside, Surface, add, dist, scale, unit}; +use super::interval::Intervals; +use super::region::{self, Scheme}; +use super::{Ctx, Raw}; +use crate::predicates::Sign; +use crate::size::{SizeField, divide_polyline, edge_divisions, ideal_divisions}; +use crate::vec3::{P, cross, dot, norm, sub}; +use crate::{Error, Result}; + +/// Angular tolerance (radians) of "parallel" and "perpendicular" for the recognition. +const ANGLE_TOL: f64 = 1e-6; +/// Distance (relative to the diagonal) within which two axes of revolution are one. +const RELATIVE_AXIS: f64 = 1e-6; +/// The most candidate frames tried per method. +pub(crate) const MAX_CANDIDATES: usize = 6; + +/// A sweep frame. +#[derive(Debug, Clone, Copy, PartialEq)] +pub(crate) enum Frame { + /// Cross-section `(p·u, p·v)` relative to `o`, height `p·d`; `u × v = d`. + Translate { o: P, u: P, v: P, d: P }, + /// Cross-section `(z, r)` (along the axis, from the axis), angle about `a` from `e1` + /// counter-clockwise; `e2 = a × e1`. + Rotate { o: P, a: P, e1: P, e2: P }, +} + +fn perpendicular(d: P) -> P { + let k = if d[0].abs() <= d[1].abs() && d[0].abs() <= d[2].abs() { + [1.0, 0.0, 0.0] + } else if d[1].abs() <= d[2].abs() { + [0.0, 1.0, 0.0] + } else { + [0.0, 0.0, 1.0] + }; + unit(cross(d, k)).unwrap_or([1.0, 0.0, 0.0]) +} + +impl Frame { + pub(crate) fn translate(o: P, d: P) -> Option { + let d = unit(d)?; + let u = perpendicular(d); + Some(Frame::Translate { o, u, v: cross(d, u), d }) + } + + pub(crate) fn rotate(o: P, a: P) -> Option { + let a = unit(a)?; + let e1 = perpendicular(a); + Some(Frame::Rotate { o, a, e1, e2: cross(a, e1) }) + } + + pub(crate) fn coords(&self, p: P) -> (Q, f64) { + match *self { + Frame::Translate { o, u, v, d } => { + let w = sub(p, o); + ([dot(w, u), dot(w, v)], dot(w, d)) + } + Frame::Rotate { o, a, e1, e2 } => { + let w = sub(p, o); + let z = dot(w, a); + let (x, y) = (dot(w, e1), dot(w, e2)); + let t = y.atan2(x); + ([z, x.hypot(y)], if t < 0.0 { t + TAU } else { t }) + } + } + } + + pub(crate) fn point(&self, q: Q, h: f64) -> P { + match *self { + Frame::Translate { o, u, v, d } => add(o, add(add(scale(u, q[0]), scale(v, q[1])), scale(d, h))), + Frame::Rotate { o, a, e1, e2 } => add(o, add(scale(a, q[0]), scale(add(scale(e1, h.cos()), scale(e2, h.sin())), q[1]))), + } + } + + fn periodic(&self) -> bool { + matches!(self, Frame::Rotate { .. }) + } + + /// The signed sweep distance from `h0` to `h1` (angles wrapped to (−π, π]). + fn dh(&self, h1: f64, h0: f64) -> f64 { + let d = h1 - h0; + match self { + Frame::Translate { .. } => d, + Frame::Rotate { .. } => { + let w = d.rem_euclid(TAU); + if w > PI { w - TAU } else { w } + } + } + } + + /// The length a sweep distance `dh` spans at cross-section point `q`. + fn span(&self, q: Q, dh: f64) -> f64 { + match self { + Frame::Translate { .. } => dh.abs(), + Frame::Rotate { .. } => dh.abs() * q[1].abs(), + } + } + + pub(crate) fn describe(&self) -> String { + let f = |p: P| format!("({:.6}, {:.6}, {:.6})", p[0], p[1], p[2]); + match *self { + Frame::Translate { d, .. } => format!("extruded along {}", f(d)), + Frame::Rotate { o, a, .. } => format!("revolved about the axis through {} along {}", f(o), f(a)), + } + } +} + +#[derive(Debug, Clone, Copy, PartialEq)] +enum Role { + Cap, + Wall, +} + +fn parallel(a: P, b: P) -> bool { + norm(cross(a, b)) <= ANGLE_TOL +} + +fn perpendicular_to(a: P, b: P) -> bool { + dot(a, b).abs() <= ANGLE_TOL +} + +/// Distance from `p` to the line through `o` along unit `a`. +fn line_distance(p: P, o: P, a: P) -> f64 { + let w = sub(p, o); + norm(sub(w, scale(a, dot(w, a)))) +} + +/// The face's role in the frame, or why it has none. +fn role(brep: &Brep, face: &Face, frame: &Frame, tol: f64) -> std::result::Result { + let name = face.id; + match *frame { + Frame::Translate { d, .. } => { + if let Some((_, n)) = face.plane { + if parallel(n, d) { + return Ok(Role::Cap); + } + if perpendicular_to(n, d) { + return Ok(Role::Wall); + } + return Err(format!("planar face {name} is slanted to the direction")); + } + match face.surface { + Some(Surface::Cylinder { axis, .. }) if parallel(axis, d) => Ok(Role::Wall), + Some(Surface::Cylinder { .. }) => Err(format!("cylindrical face {name} is not parallel to the direction")), + Some(_) => Err(format!("face {name} is curved across the direction")), + None => { + // A free-form face is a wall when every triangle is parallel to the direction. + let ok = face.triangles.iter().all(|&t| { + let Some(tri) = brep.surface.triangles().get(t as usize) else { return false }; + let [a, b, c] = tri.map(|v| brep.point(v)); + unit(cross(sub(b, a), sub(c, a))).is_some_and(|n| perpendicular_to(n, d)) + }); + if ok { Ok(Role::Wall) } else { Err(format!("free-form face {name} is not parallel to the direction")) } + } + } + } + Frame::Rotate { o, a, .. } => { + let coaxial = |origin: P, axis: P| parallel(axis, a) && line_distance(origin, o, a) <= tol; + if let Some((p, n)) = face.plane { + if perpendicular_to(n, a) && dot(sub(p, o), n).abs() <= tol { + return Ok(Role::Cap); + } + if parallel(n, a) { + return Ok(Role::Wall); + } + return Err(format!("planar face {name} neither contains nor crosses the axis")); + } + match face.surface { + Some(Surface::Cylinder { origin, axis, .. }) + | Some(Surface::Cone { apex: origin, axis, .. }) + | Some(Surface::Torus { center: origin, axis, .. }) + if coaxial(origin, axis) => + { + Ok(Role::Wall) + } + Some(Surface::Sphere { center, .. }) if line_distance(center, o, a) <= tol => Ok(Role::Wall), + Some(_) => Err(format!("face {name} is not a surface of revolution about the axis")), + None => Err(format!("free-form face {name}: a revolution needs analytic faces")), + } + } + } +} + +/// Candidate extrusion directions: the normals of planar faces and the axes of cylinders, most +/// common first. +pub(crate) fn translations(brep: &Brep) -> Vec { + let mut dirs: Vec<(P, f64)> = Vec::new(); + let mut push = |d: P, weight: f64| { + let Some(mut d) = unit(d) else { return }; + // A canonical sign: the largest component positive. + let k = (0..3).max_by(|&i, &j| d[i].abs().total_cmp(&d[j].abs()).then(j.cmp(&i))).unwrap_or(0); + if d.get(k).copied().unwrap_or(0.0) < 0.0 { + d = scale(d, -1.0); + } + match dirs.iter_mut().find(|(e, _)| parallel(*e, d)) { + Some(slot) => slot.1 += weight, + None => dirs.push((d, weight)), + } + }; + for f in &brep.faces { + if let Some((_, n)) = f.plane { + push(n, f.area); + } else if let Some(Surface::Cylinder { axis, .. }) = f.surface { + push(axis, 0.0); + } + } + dirs.sort_by(|a, b| b.1.total_cmp(&a.1)); + let o = brep.point(0); + dirs.into_iter().filter_map(|(d, _)| Frame::translate(o, d)).collect() +} + +/// Candidate axes of revolution: the axes of the analytic surfaces of revolution. +pub(crate) fn rotations(brep: &Brep) -> Vec { + let mut axes: Vec<(P, P)> = Vec::new(); + for f in &brep.faces { + let (o, a) = match f.surface { + Some(Surface::Cylinder { origin, axis, .. }) + | Some(Surface::Cone { apex: origin, axis, .. }) + | Some(Surface::Torus { center: origin, axis, .. }) => (origin, axis), + _ => continue, + }; + if !axes.iter().any(|(p, b)| parallel(*b, a) && line_distance(o, *p, *b) <= RELATIVE_AXIS * brep.diagonal) { + axes.push((o, a)); + } + } + axes.into_iter().filter_map(|(o, a)| Frame::rotate(o, a)).collect() +} + +/// The frame's classification of the body: face roles and levels. +pub(crate) struct Setup { + frame: Frame, + roles: Vec, + /// Sorted levels (heights or angles in `[0, 2π)`). + levels: Vec, +} + +/// Classify faces and edges and find the levels; an `Unsupported` error says why the body does +/// not sweep in this frame. +pub(crate) fn setup(ctx: &Ctx, frame: &Frame) -> Result { + let brep = ctx.brep; + let tol = ctx.tol; + let mut roles = Vec::with_capacity(brep.faces.len()); + for f in &brep.faces { + roles.push(role(brep, f, frame, tol).map_err(Error::Unsupported)?); + } + if frame.periodic() { + // A revolved region touching the axis would collapse its elements there. + if let Some(p) = brep.surface.points().iter().find(|&&p| frame.coords(p).0[1] <= tol) { + return Err(Error::Unsupported(format!("the body reaches the axis of revolution (at {p:?}); revolved hexahedra would collapse there"))); + } + } + let r_max = brep.surface.points().iter().map(|&p| frame.coords(p).0[1].abs()).fold(0.0, f64::max); + let h_tol = match frame { + Frame::Translate { .. } => tol, + Frame::Rotate { .. } => tol / r_max.max(tol), + }; + let mut heights = Vec::new(); + for (e, edge) in brep.edges.iter().enumerate() { + let pts: Vec<(Q, f64)> = edge.points.iter().map(|&v| frame.coords(brep.point(v))).collect(); + let Some(&(q0, h0)) = pts.first() else { continue }; + let across = pts.iter().all(|&(q, h)| frame.span(q, frame.dh(h, h0)) <= tol); + let along = pts.iter().all(|&(q, _)| dist([q[0], q[1], 0.0], [q0[0], q0[1], 0.0]) <= tol); + if across { + heights.push(h0); + } else if !along { + let [f0, f1] = edge.faces.map(|f| brep.face_tag(f)); + return Err(Error::Unsupported(format!("edge {e} between faces {f0} and {f1} runs neither along nor across the sweep"))); + } + } + for (f, r) in brep.faces.iter().zip(&roles) { + if *r == Role::Cap + && let Some(&v) = f.vertices.first() + { + heights.push(frame.coords(brep.point(v)).1); + } + } + heights.sort_by(f64::total_cmp); + let mut levels: Vec = Vec::new(); + for h in heights { + if levels.last().is_none_or(|&l| h - l > h_tol) { + levels.push(h); + } + } + if frame.periodic() + && levels.len() > 1 + && let (Some(&first), Some(&last)) = (levels.first(), levels.last()) + && first + TAU - last <= h_tol + { + levels.pop(); + } + let need = if frame.periodic() { 1 } else { 2 }; + if levels.len() < need { + return Err(Error::Unsupported("no source and target faces across the sweep".into())); + } + Ok(Setup { frame: *frame, roles, levels }) +} + +impl Setup { + /// The number of levels the frame sweeps through (for ordering candidates). + pub(crate) fn level_count(&self) -> usize { + self.levels.len() + } +} + +/// A size field on the cross-section plane at sweep coordinate `h`. +pub(crate) struct PlaneField<'a> { + pub field: &'a dyn SizeField, + pub frame: Frame, + pub h: f64, +} + +impl SizeField for PlaneField<'_> { + fn size(&self, p: [f64; 3]) -> f64 { + self.field.size(self.frame.point([p[0], p[1]], self.h)) + } +} + +/// The cross-section of a frame the body sweeps in: the arrangement (cut at reflex corners), +/// the level intervals and which of them each region fills. It does not depend on the element +/// size, so refinement retries reuse it. +pub(crate) struct Prepared { + setup: Setup, + arr: Arrangement, + intervals: Vec<(f64, f64)>, + filled: Vec>, + used: Vec, +} + +/// Lay out the cross-section of `setup`'s frame. +pub(crate) fn prepare(ctx: &Ctx, setup: Setup) -> Result { + let brep = ctx.brep; + let tol = ctx.tol; + let frame = &setup.frame; + let (roles, levels) = (&setup.roles, &setup.levels); + // Lower outlines of the walls, and every vertex. + let mut polylines: Vec<(Vec, usize)> = Vec::new(); + for (e, edge) in brep.edges.iter().enumerate() { + let pts: Vec<(Q, f64)> = edge.points.iter().map(|&v| frame.coords(brep.point(v))).collect(); + let Some(&(_, h0)) = pts.first() else { continue }; + if !pts.iter().all(|&(q, h)| frame.span(q, frame.dh(h, h0)) <= tol) { + continue; + } + for &f in &edge.faces { + if roles.get(f) != Some(&Role::Wall) { + continue; + } + // The wall's triangle on the edge's first segment: is its third vertex above? + let Some(c) = brep.third_vertex(e, f) else { continue }; + let (qc, hc) = frame.coords(brep.point(c)); + if frame.dh(hc, h0) > 0.0 && frame.span(qc, frame.dh(hc, h0)) > tol { + polylines.push((pts.iter().map(|x| x.0).collect(), f)); + } + } + } + let vertices: Vec = brep + .edges + .iter() + .flat_map(|e| [e.first(), e.last()]) + .chain(brep.is_vertex.iter().enumerate().filter(|(_, v)| **v).map(|(i, _)| i as u32)) + .map(|v| frame.coords(brep.point(v)).0) + .collect(); + let coincide = 2.0 * ctx.opts.chordal_tolerance.max(0.0) + 4.0 * tol; + // Level intervals; a revolution's last one wraps round to the first level. + let (Some(&first_level), Some(rest)) = (levels.first(), levels.get(1..)) else { + return Err(Error::Unsupported("no levels".into())); + }; + let ends: Vec = if frame.periodic() { rest.iter().copied().chain(std::iter::once(first_level + TAU)).collect() } else { rest.to_vec() }; + let intervals: Vec<(f64, f64)> = levels.iter().copied().zip(ends).collect(); + let mut inside = Inside::new(brep.surface); + let first = arrangement::build(&vertices, &polylines, tol, coincide)?; + let (_, used_first) = occupancy(&first, frame, &intervals, &mut inside)?; + // Cut the filled regions at their reflex corners, then lay the arrangement out again. + let cut_lines = cuts(&first, &used_first, tol); + let arr = if cut_lines.is_empty() { + first + } else { + polylines.extend(cut_lines); + arrangement::build(&vertices, &polylines, tol, coincide)? + }; + let (filled, used) = occupancy(&arr, frame, &intervals, &mut inside)?; + Ok(Prepared { setup, arr, intervals, filled, used }) +} + +/// Mesh the prepared cross-section, with its curves `2^refine` times finer than the size asks. +pub(crate) fn mesh(ctx: &Ctx, prepared: &Prepared, refine: u32) -> Result { + let brep = ctx.brep; + let Prepared { setup, arr, intervals, filled, used } = prepared; + let (frame, levels) = (&setup.frame, &setup.levels); + let h_ref = levels.first().copied().unwrap_or(0.0); + // Interval assignment over the curves of filled regions, and the layer counts. + let mut ia = Intervals::default(); + let plane_field = PlaneField { field: ctx.field, frame: *frame, h: h_ref }; + let curve3d = |c: usize| -> Vec

{ + arr.curves + .get(c) + .map(|cv| cv.points.iter().map(|&p| frame.point(arr.points.get(p).copied().unwrap_or([0.0; 2]), h_ref)).collect()) + .unwrap_or_default() + }; + let fine = f64::from(1u32 << refine.min(8)); + for c in 0..arr.curves.len() { + let goal = ideal_divisions(ctx.field, &curve3d(c))?; + ia.var(goal * fine, 1); + } + let oriented = |c: usize, fwd: bool| -> Vec { + let mut v: Vec = + arr.curves.get(c).map(|x| x.points.iter().map(|&p| arr.points.get(p).copied().unwrap_or([0.0; 2])).collect()).unwrap_or_default(); + if !fwd { + v.reverse(); + } + v + }; + let mut schemes: Vec = vec![Scheme::Free; arr.regions.len()]; + for (r, region) in arr.regions.iter().enumerate() { + if !used.get(r).copied().unwrap_or(false) { + continue; + } + let polylines: Vec>> = region.loops.iter().map(|lp| lp.iter().map(|&(c, f)| oriented(c, f)).collect()).collect(); + let scheme = region::classify(&polylines); + let vars: Vec> = region.loops.iter().map(|lp| lp.iter().map(|x| x.0).collect()).collect(); + region::constrain(&mut ia, &vars, scheme); + if let Some(slot) = schemes.get_mut(r) { + *slot = scheme; + } + } + let counts = ia.solve()?; + let r_max = arr.points.iter().map(|q| q[1].abs()).fold(0.0, f64::max); + let mut layers = Vec::with_capacity(intervals.len()); + for &(lo, hi) in intervals.iter() { + let path: Vec

= match frame { + Frame::Translate { .. } => vec![frame.point([0.0, 0.0], lo), frame.point([0.0, 0.0], hi)], + Frame::Rotate { .. } => (0..=16).map(|k| frame.point([0.0, r_max], lo + (hi - lo) * k as f64 / 16.0)).collect(), + }; + let mut n = edge_divisions(ctx.field, &path)?; + if frame.periodic() { + n = n.max(((hi - lo) / FRAC_PI_4).ceil() as usize); + } + layers.push(n); + } + // 2D nodes: curve ends, then curve interiors (projected onto analytic walls). + let mut nodes2: Vec = Vec::new(); + let mut point_node: HashMap = HashMap::new(); + let mut curve_nodes: Vec> = Vec::with_capacity(arr.curves.len()); + for (c, curve) in arr.curves.iter().enumerate() { + let n = counts.get(c).copied().unwrap_or(1); + let pts3 = divide_polyline(ctx.field, &curve3d(c), n)?; + let wall = curve.tags.iter().find_map(|&f| brep.faces.get(f).and_then(|x| x.surface)).filter(|s| !matches!(s, Surface::Plane { .. })); + let mut ids = Vec::with_capacity(pts3.len()); + let last = pts3.len().saturating_sub(1); + for (k, p) in pts3.iter().enumerate() { + let end = if k == 0 { + curve.points.first() + } else if k == last { + curve.points.last() + } else { + None + }; + let id = match end { + Some(&pid) => *point_node.entry(pid).or_insert_with(|| { + nodes2.push(arr.points.get(pid).copied().unwrap_or([0.0; 2])); + nodes2.len() as u32 - 1 + }), + None => { + let p = wall.and_then(|s| s.project(*p)).unwrap_or(*p); + nodes2.push(frame.coords(p).0); + nodes2.len() as u32 - 1 + } + }; + ids.push(id); + } + curve_nodes.push(ids); + } + // Quad meshes of the filled regions. + let mut region_quads: Vec> = vec![Vec::new(); arr.regions.len()]; + let quad_opts = &ctx.opts.quad; + for (r, region) in arr.regions.iter().enumerate() { + if !used.get(r).copied().unwrap_or(false) { + continue; + } + let mut loops_ids: Vec> = Vec::with_capacity(region.loops.len()); + let mut starts: Vec> = Vec::new(); + for lp in ®ion.loops { + let mut ids = Vec::new(); + let mut st = Vec::new(); + for &(c, fwd) in lp { + st.push(ids.len()); + let cn = curve_nodes.get(c).map(Vec::as_slice).unwrap_or(&[]); + let mut seq: Vec = cn.to_vec(); + if !fwd { + seq.reverse(); + } + seq.pop(); + ids.extend(seq); + } + loops_ids.push(ids); + starts.push(st); + } + let at = |i: u32| nodes2.get(i as usize).copied().unwrap_or([0.0; 2]); + let loops: Vec> = loops_ids.iter().map(|l| l.iter().map(|&i| at(i)).collect()).collect(); + let scheme = schemes.get(r).copied().unwrap_or(Scheme::Free); + let (mesh, used_scheme) = region::mesh(&loops, &starts, scheme, &plane_field, quad_opts).map_err(|e| region_error(r, e))?; + if let Some(slot) = schemes.get_mut(r) { + *slot = used_scheme; + } + let boundary: Vec = loops_ids.concat(); + let mut local: Vec = boundary.clone(); + for p in mesh.nodes.iter().skip(mesh.boundary_nodes) { + nodes2.push(*p); + local.push(nodes2.len() as u32 - 1); + } + let quads = region_quads.get_mut(r).ok_or_else(|| Error::Failed("region index".into()))?; + for q in &mesh.quads { + quads.push(q.map(|i| local.get(i as usize).copied().unwrap_or(0))); + } + } + // Layer coordinates. + let mut heights: Vec = Vec::new(); + let mut first_layer: Vec = Vec::with_capacity(intervals.len()); + for (k, &(lo, hi)) in intervals.iter().enumerate() { + first_layer.push(heights.len()); + let n = layers.get(k).copied().unwrap_or(1); + for j in 0..n { + heights.push(lo + (hi - lo) * j as f64 / n as f64); + } + } + if !frame.periodic() + && let Some(&(_, hi)) = intervals.last() + { + heights.push(hi); + } + let total_layers = heights.len(); + // Element budget. + let mut count = 0usize; + for (r, quads) in region_quads.iter().enumerate() { + for (k, f) in filled.get(r).map(Vec::as_slice).unwrap_or(&[]).iter().enumerate() { + if *f { + count = count.saturating_add(quads.len().saturating_mul(layers.get(k).copied().unwrap_or(0))); + } + } + } + if count > ctx.opts.max_elements { + return Err(Error::Limit(format!("the mesh would have {count} hexahedra; the limit is {}", ctx.opts.max_elements))); + } + // Hexahedra, with nodes made as first used. + let mut nodes: Vec

= Vec::new(); + let mut node_at: HashMap<(u32, usize), u32> = HashMap::new(); + let mut node = |n2: u32, layer: usize, nodes: &mut Vec

| -> u32 { + let layer = if frame.periodic() { layer % total_layers.max(1) } else { layer }; + *node_at.entry((n2, layer)).or_insert_with(|| { + let q = nodes2.get(n2 as usize).copied().unwrap_or([0.0; 2]); + nodes.push(frame.point(q, heights.get(layer).copied().unwrap_or(0.0))); + nodes.len() as u32 - 1 + }) + }; + let mut hexes: Vec<[u32; 8]> = Vec::with_capacity(count); + for (r, quads) in region_quads.iter().enumerate() { + for (k, f) in filled.get(r).map(Vec::as_slice).unwrap_or(&[]).iter().enumerate() { + if !*f { + continue; + } + let (s, n) = (first_layer.get(k).copied().unwrap_or(0), layers.get(k).copied().unwrap_or(0)); + for q in quads { + for j in s..s + n { + let b = q.map(|v| node(v, j, &mut nodes)); + let t = q.map(|v| node(v, j + 1, &mut nodes)); + hexes.push([b[0], b[1], b[2], b[3], t[0], t[1], t[2], t[3]]); + } + } + } + } + let filled_regions = used.iter().filter(|&&u| u).count(); + let structured = schemes.iter().zip(used.iter()).filter(|(s, u)| **u && **s != Scheme::Free).count(); + let description = format!( + "{}: {} level{}, {} region{} ({} structured, {} free), {} layer{}", + frame.describe(), + levels.len(), + if levels.len() == 1 { "" } else { "s" }, + filled_regions, + if filled_regions == 1 { "" } else { "s" }, + structured, + filled_regions - structured, + total_layers.saturating_sub(usize::from(!frame.periodic())), + if total_layers == 2 { "" } else { "s" }, + ); + Ok(Raw { nodes, hexes, description }) +} + +/// Which level intervals each region of the arrangement fills (inside test at the region's +/// centre, mid-interval), and which regions fill any. +fn occupancy(arr: &Arrangement, frame: &Frame, intervals: &[(f64, f64)], inside: &mut Inside) -> Result<(Vec>, Vec)> { + if arr.regions.is_empty() { + return Err(Error::Unsupported("the projected outline encloses no region".into())); + } + let mut filled: Vec> = Vec::with_capacity(arr.regions.len()); + for region in &arr.regions { + let c = interior_point(arr, region)?; + let mut row = Vec::with_capacity(intervals.len()); + for &(lo, hi) in intervals { + row.push(inside.contains(frame.point(c, 0.5 * (lo + hi)))?); + } + filled.push(row); + } + let used: Vec = filled.iter().map(|r| r.iter().any(|&x| x)).collect(); + if !used.iter().any(|&u| u) { + return Err(Error::Unsupported("no region of the cross-section lies inside the body".into())); + } + Ok((filled, used)) +} + +/// The tag of cut lines (no wall behind them). +const CUT: usize = usize::MAX; + +/// Straight cuts from the reflex corners of the filled regions into them, each along one of the +/// corner's edges (the shorter way) to the first boundary it meets. A rectilinear region falls +/// into rectangles this way; others lose their reflex corners. This is the explicit form of the +/// virtual decomposition behind submapping (Whiteley, White, Benzley, Blacker 1996), done in the +/// arrangement so that interval assignment sees every piece. +fn cuts(arr: &Arrangement, used: &[bool], tol: f64) -> Vec<(Vec, usize)> { + let at = |p: usize| arr.points.get(p).copied().unwrap_or([0.0; 2]); + let segments: Vec<(usize, usize)> = arr.curves.iter().flat_map(|c| c.points.windows(2).map(|w| (w[0], w[1]))).collect(); + // The nearest boundary hit of the ray from point v along unit d (ignoring segments at v). + let hit = |v: usize, d: Q| -> Option<(f64, Q, Option)> { + let o = at(v); + let mut best: Option<(f64, Q, Option)> = None; + for &(a, b) in &segments { + if a == v || b == v { + continue; + } + let (pa, pb) = (at(a), at(b)); + let e = [pb[0] - pa[0], pb[1] - pa[1]]; + let den = d[0] * e[1] - d[1] * e[0]; + if den.abs() <= f64::EPSILON * (e[0].abs() + e[1].abs()) { + continue; + } + let w = [pa[0] - o[0], pa[1] - o[1]]; + let t = (w[0] * e[1] - w[1] * e[0]) / den; + let u = (w[0] * d[1] - w[1] * d[0]) / den; + if t > tol && (-1e-12..=1.0 + 1e-12).contains(&u) && best.is_none_or(|(bt, _, _)| t < bt) { + let x = [o[0] + t * d[0], o[1] + t * d[1]]; + // End at a segment end that is this close to the hit. + let near = [(a, pa), (b, pb)].into_iter().find(|(_, p)| ((p[0] - x[0]).powi(2) + (p[1] - x[1]).powi(2)).sqrt() <= 0.1 * t); + best = Some((t, near.map_or(x, |n| n.1), near.map(|n| n.0))); + } + } + best + }; + let mut out: Vec<(Vec, usize)> = Vec::new(); + let mut seen: Vec<(usize, usize)> = Vec::new(); + for (r, region) in arr.regions.iter().enumerate() { + if !used.get(r).copied().unwrap_or(false) { + continue; + } + for lp in ®ion.loops { + let n = lp.len(); + for k in 0..n { + let (Some(&(pc, pf)), Some(&(nc, nf))) = (lp.get((k + n - 1) % n), lp.get(k)) else { continue }; + let pts = |c: usize, f: bool| { + let mut v = arr.curves.get(c).map(|x| x.points.clone()).unwrap_or_default(); + if !f { + v.reverse(); + } + v + }; + let (prev, next) = (pts(pc, pf), pts(nc, nf)); + let (Some(&a), Some(&v), Some(&b)) = (prev.get(prev.len().wrapping_sub(2)), next.first(), next.get(1)) else { continue }; + let (pa, pv, pb) = (at(a), at(v), at(b)); + let (u, w) = ([pv[0] - pa[0], pv[1] - pa[1]], [pb[0] - pv[0], pb[1] - pv[1]]); + let turn = (u[0] * w[1] - u[1] * w[0]).atan2(u[0] * w[0] + u[1] * w[1]); + // The region is on the left: a right turn of more than 45° is a reflex corner. + if PI - turn < 1.25 * PI { + continue; + } + let norm2 = |x: Q| (x[0] * x[0] + x[1] * x[1]).sqrt(); + let (lu, lw) = (norm2(u), norm2(w)); + if !(lu > 0.0 && lw > 0.0) { + continue; + } + let rays = [[u[0] / lu, u[1] / lu], [-w[0] / lw, -w[1] / lw]]; + let best = rays.iter().filter_map(|&d| hit(v, d)).min_by(|x, y| x.0.total_cmp(&y.0)); + if let Some((_, end, end_point)) = best { + if let Some(e) = end_point { + let key = (v.min(e), v.max(e)); + if seen.contains(&key) { + continue; + } + seen.push(key); + } + out.push((vec![pv, end], CUT)); + } + } + } + } + out +} + +fn region_error(r: usize, e: Error) -> Error { + match e { + Error::Limit(m) => Error::Limit(m), + other => Error::Failed(format!("cross-section region {r}: {other}")), + } +} + +/// A point well inside a region: off the middle of one of its outer edges, inside the outer +/// loop, outside the holes and clear of every boundary segment. +fn interior_point(arr: &Arrangement, region: &arrangement::Region) -> Result { + let at = |p: usize| arr.points.get(p).copied().unwrap_or([0.0; 2]); + let polygon = |lp: &[(usize, bool)]| -> Vec { + let mut out = Vec::new(); + for &(c, fwd) in lp { + let mut v: Vec = arr.curves.get(c).map(|x| x.points.clone()).unwrap_or_default(); + if !fwd { + v.reverse(); + } + v.pop(); + out.extend(v.into_iter().map(at)); + } + out + }; + let loops: Vec> = region.loops.iter().map(|l| polygon(l)).collect(); + // Crossing-number winding with exact side tests (a probe point off the plane's domain + // counts as outside). + let winding = |p: Q, poly: &[Q]| -> i64 { + let mut wn = 0; + for k in 0..poly.len() { + let (a, b) = (poly[k], poly[(k + 1) % poly.len()]); + let side = arrangement::sign(a, b, p).unwrap_or(Sign::Zero); + if a[1] <= p[1] { + if b[1] > p[1] && side == Sign::Positive { + wn += 1; + } + } else if b[1] <= p[1] && side == Sign::Negative { + wn -= 1; + } + } + wn + }; + let clearance = |p: Q| -> f64 { + loops + .iter() + .flat_map(|poly| (0..poly.len()).map(move |k| (poly[k], poly[(k + 1) % poly.len()]))) + .map(|(a, b)| arrangement::seg_dist2(p, a, b).0.sqrt()) + .fold(f64::INFINITY, f64::min) + }; + let outer = loops.first().ok_or_else(|| Error::Failed("region without a loop".into()))?; + let mut edges: Vec<(f64, usize)> = (0..outer.len()).map(|k| (arrangement::d2(outer[k], outer[(k + 1) % outer.len()]), k)).collect(); + edges.sort_by(|a, b| b.0.total_cmp(&a.0).then(a.1.cmp(&b.1))); + let mut best: Option<(f64, Q)> = None; + for &(len2, k) in edges.iter().take(64) { + let (a, b) = (outer[k], outer[(k + 1) % outer.len()]); + let len = len2.sqrt(); + let (m, n) = ([(a[0] + b[0]) * 0.5, (a[1] + b[1]) * 0.5], [-(b[1] - a[1]) / len, (b[0] - a[0]) / len]); + for f in [0.5, 0.25, 0.1, 0.03, 0.01, 1e-3] { + let p = [m[0] + n[0] * f * len, m[1] + n[1] * f * len]; + if winding(p, outer) == 0 || loops.iter().skip(1).any(|h| winding(p, h) != 0) { + continue; + } + let c = clearance(p); + if best.is_none_or(|(bc, _)| c > bc) { + best = Some((c, p)); + } + } + } + best.filter(|(c, _)| *c > 0.0).map(|(_, p)| p).ok_or_else(|| Error::Failed("no point found inside a cross-section region".into())) +} diff --git a/crates/mesh/src/hex/mapped.rs b/crates/mesh/src/hex/mapped.rs new file mode 100644 index 0000000..a31458d --- /dev/null +++ b/crates/mesh/src/hex/mapped.rs @@ -0,0 +1,246 @@ +//! Mapped meshing of a logical cube: six faces with one loop of four edges each, three edges at +//! each of eight corners. +//! +//! The corners get logical coordinates `(i, j, k) ∈ {0, 1}³` from the edge graph; the four edges +//! along each logical direction share one division count (interval assignment). Edges are divided +//! on the geometry (nodes on both faces of each edge), each face gets the transfinite grid of its +//! four edges projected onto its surface (Gordon and Hall 1973), and the volume is the trivariate +//! transfinite interpolation of the six face grids: faces blended linearly, minus edges blended +//! bilinearly, plus corners blended trilinearly (W. J. Gordon, L. C. Thiel, *Transfinite +//! mappings and their application to grid generation*, Applied Mathematics and Computation +//! 10–11:171–233, 1982). + +use std::collections::BTreeMap; + +use super::geom::{add, scale}; +use super::interval::Intervals; +use super::{Ctx, Raw, project_to_face, project_to_faces}; +use crate::quad::transfinite_grid; +use crate::size::{divide_polyline, ideal_divisions}; +use crate::vec3::{P, det, sub}; +use crate::{Error, Result}; + +/// The most nodes of a mapped block. +const MAX_NODES: usize = 1 << 26; + +pub(crate) fn mesh(ctx: &Ctx) -> Result { + let brep = ctx.brep; + if brep.faces.len() != 6 || brep.edges.len() != 12 { + return Err(Error::Unsupported(format!("not a logical cube: {} faces and {} edges (6 and 12 needed)", brep.faces.len(), brep.edges.len()))); + } + if let Some(f) = brep.faces.iter().find(|f| f.loops.len() != 1 || f.loops.first().map(Vec::len) != Some(4)) { + return Err(Error::Unsupported(format!("not a logical cube: face {} is not bounded by one loop of four edges", f.id))); + } + // Corner graph. + let mut adj: BTreeMap> = BTreeMap::new(); + for (e, edge) in brep.edges.iter().enumerate() { + if edge.closed() { + return Err(Error::Unsupported("not a logical cube: a closed edge".into())); + } + adj.entry(edge.first()).or_default().push((edge.last(), e)); + adj.entry(edge.last()).or_default().push((edge.first(), e)); + } + if adj.len() != 8 || adj.values().any(|v| v.len() != 3) { + return Err(Error::Unsupported("not a logical cube: the corners do not each have three edges".into())); + } + let neighbours = |v: u32| -> Vec { adj.get(&v).map(|l| l.iter().map(|x| x.0).collect()).unwrap_or_default() }; + let common = |a: u32, b: u32, not: u32| -> Option { + let nb = neighbours(b); + neighbours(a).into_iter().find(|x| *x != not && nb.contains(x)) + }; + let c000 = *adj.keys().next().ok_or_else(|| Error::Failed("no corner".into()))?; + let n = neighbours(c000); + let (Some(&a), Some(&b), Some(&c)) = (n.first(), n.get(1), n.get(2)) else { return Err(Error::Failed("corner".into())) }; + // Right-handed: (a, b, c) along (i, j, k) when the triple product is positive. + let p = |v: u32| brep.point(v); + let (a, b, c) = if det(sub(p(a), p(c000)), sub(p(b), p(c000)), sub(p(c), p(c000))) > 0.0 { (a, b, c) } else { (a, c, b) }; + let c110 = common(a, b, c000); + let c101 = common(a, c, c000); + let c011 = common(b, c, c000); + let (Some(c110), Some(c101), Some(c011)) = (c110, c101, c011) else { + return Err(Error::Unsupported("not a logical cube: the corner graph is not a cube".into())); + }; + let c111 = common(c110, c101, a).filter(|&x| neighbours(c011).contains(&x)); + let Some(c111) = c111 else { return Err(Error::Unsupported("not a logical cube: the corner graph is not a cube".into())) }; + // corner[i][j][k] + let corner = [[[c000, c], [b, c011]], [[a, c101], [c110, c111]]]; + let at = |i: usize, j: usize, k: usize| -> u32 { corner[i][j][k] }; + let mut distinct: Vec = corner.iter().flatten().flatten().copied().collect(); + distinct.sort_unstable(); + distinct.dedup(); + if distinct.len() != 8 { + return Err(Error::Unsupported("not a logical cube: corners repeat".into())); + } + // The edge between two corners, oriented from the first. + let edge_between = |u: u32, v: u32| -> Option<(usize, bool)> { + adj.get(&u)?.iter().find(|x| x.0 == v).map(|x| (x.1, brep.edges.get(x.1).is_some_and(|e| e.first() == u))) + }; + // The four edges along each direction, from the low corner to the high one. + let mut groups: [Vec<(usize, bool)>; 3] = [Vec::new(), Vec::new(), Vec::new()]; + for (d, group) in groups.iter_mut().enumerate() { + for s in 0..2 { + for t in 0..2 { + let (lo, hi) = match d { + 0 => (at(0, s, t), at(1, s, t)), + 1 => (at(s, 0, t), at(s, 1, t)), + _ => (at(s, t, 0), at(s, t, 1)), + }; + group.push(edge_between(lo, hi).ok_or_else(|| Error::Unsupported("not a logical cube: a corner pair has no edge".into()))?); + } + } + } + let mut ia = Intervals::default(); + for e in 0..brep.edges.len() { + ia.var(ideal_divisions(ctx.field, &brep.edge_points(e, true))?, 1); + } + for g in &groups { + let first = g.first().map(|x| x.0).unwrap_or(0); + for &(e, _) in g.iter().skip(1) { + ia.equal.push((vec![first], vec![e])); + } + } + let counts = ia.solve()?; + let count = |g: &Vec<(usize, bool)>| g.first().and_then(|x| counts.get(x.0)).copied().unwrap_or(1); + let (ni, nj, nk) = (count(&groups[0]), count(&groups[1]), count(&groups[2])); + let total = (ni + 1).checked_mul(nj + 1).and_then(|x| x.checked_mul(nk + 1)).filter(|&x| x <= MAX_NODES); + let total = total.ok_or_else(|| Error::Limit(format!("a {ni} × {nj} × {nk} block exceeds {MAX_NODES} nodes")))?; + if ni * nj * nk > ctx.opts.max_elements { + return Err(Error::Limit(format!("the mesh would have {} hexahedra; the limit is {}", ni * nj * nk, ctx.opts.max_elements))); + } + // Edge node positions, low corner to high. + let edge_pts = |(e, fwd): (usize, bool), n: usize| -> Result> { + let pts = divide_polyline(ctx.field, &brep.edge_points(e, fwd), n)?; + let faces = brep.edges.get(e).map(|x| x.faces.to_vec()).unwrap_or_default(); + let last = pts.len().saturating_sub(1); + Ok(pts.iter().enumerate().map(|(k, q)| if k == 0 || k == last { *q } else { project_to_faces(ctx, &faces, *q) }).collect()) + }; + let idx = |i: usize, j: usize, k: usize| (k * (nj + 1) + j) * (ni + 1) + i; + let mut x: Vec> = vec![None; total]; + // Edges: E[d][s][t] along direction d at the other two coordinates (s, t) ∈ {0, 1}². + for (d, g) in groups.iter().enumerate() { + let n = [ni, nj, nk][d]; + for (m, &use_) in g.iter().enumerate() { + let (s, t) = (m / 2, m % 2); + for (q, pt) in edge_pts(use_, n)?.into_iter().enumerate() { + let (i, j, k) = match d { + 0 => (q, s * nj, t * nk), + 1 => (s * ni, q, t * nk), + _ => (s * ni, t * nj, q), + }; + if let Some(slot) = x.get_mut(idx(i, j, k)) { + slot.get_or_insert(pt); + } + } + } + } + let get = |x: &[Option

], i: usize, j: usize, k: usize| -> Result

{ + x.get(idx(i, j, k)).copied().flatten().ok_or_else(|| Error::Failed("mapped grid node".into())) + }; + // Faces: grids from their edges, projected onto the face found at the face's centre. + for d in 0..3 { + for side in 0..2 { + // Face d = side·n_d; its parameters (a, b) run along the other two directions. + let (da, db) = match d { + 0 => (1, 2), + 1 => (0, 2), + _ => (0, 1), + }; + let (na, nb) = ([ni, nj, nk][da], [ni, nj, nk][db]); + let ijk = |a: usize, b: usize| -> (usize, usize, usize) { + let mut c = [0usize; 3]; + c[d] = side * [ni, nj, nk][d]; + c[da] = a; + c[db] = b; + (c[0], c[1], c[2]) + }; + let line = |f: &dyn Fn(usize) -> (usize, usize, usize), n: usize| -> Result> { + (0..=n) + .map(|q| { + let (i, j, k) = f(q); + get(&x, i, j, k) + }) + .collect() + }; + let bottom = line(&|a| ijk(a, 0), na)?; + let top = line(&|a| ijk(a, nb), na)?; + let left = line(&|b| ijk(0, b), nb)?; + let right = line(&|b| ijk(na, b), nb)?; + let grid = transfinite_grid(&bottom, &right, &top, &left)?; + // The B-rep face bounded by this side's four edges. + let logical = |a: usize, b: usize| { + let (i, j, k) = ijk(a, b); + at(usize::from(i > 0), usize::from(j > 0), usize::from(k > 0)) + }; + let ring = [logical(0, 0), logical(na, 0), logical(na, nb), logical(0, nb)]; + let mut bounding: Vec = (0..4).filter_map(|m| edge_between(ring[m], ring[(m + 1) % 4]).map(|x| x.0)).collect(); + bounding.sort_unstable(); + let face = brep.faces.iter().position(|f| { + let mut es: Vec = f.loops.iter().flatten().map(|x| x.0).collect(); + es.sort_unstable(); + es == bounding + }); + for b in 1..nb { + for a in 1..na { + let q = grid.get(b * (na + 1) + a).copied().unwrap_or([0.0; 3]); + let q = face.and_then(|f| project_to_face(ctx, f, q)).unwrap_or(q); + let (i, j, k) = ijk(a, b); + if let Some(slot) = x.get_mut(idx(i, j, k)) { + *slot = Some(q); + } + } + } + } + } + // Volume: trivariate transfinite interpolation. + let (fi, fj, fk) = (ni as f64, nj as f64, nk as f64); + for k in 1..nk { + for j in 1..nj { + for i in 1..ni { + let (u, v, w) = (i as f64 / fi, j as f64 / fj, k as f64 / fk); + let g = |i: usize, j: usize, k: usize| get(&x, i, j, k); + let lin = |t: f64, a: P, b: P| add(scale(a, 1.0 - t), scale(b, t)); + // Faces. + let faces = add(add(lin(u, g(0, j, k)?, g(ni, j, k)?), lin(v, g(i, 0, k)?, g(i, nj, k)?)), lin(w, g(i, j, 0)?, g(i, j, nk)?)); + // Edges. + let bil = |s: f64, t: f64, p00: P, p10: P, p01: P, p11: P| lin(t, lin(s, p00, p10), lin(s, p01, p11)); + let edges = add( + add( + bil(v, w, g(i, 0, 0)?, g(i, nj, 0)?, g(i, 0, nk)?, g(i, nj, nk)?), + bil(u, w, g(0, j, 0)?, g(ni, j, 0)?, g(0, j, nk)?, g(ni, j, nk)?), + ), + bil(u, v, g(0, 0, k)?, g(ni, 0, k)?, g(0, nj, k)?, g(ni, nj, k)?), + ); + // Corners. + let corners = lin( + w, + bil(u, v, g(0, 0, 0)?, g(ni, 0, 0)?, g(0, nj, 0)?, g(ni, nj, 0)?), + bil(u, v, g(0, 0, nk)?, g(ni, 0, nk)?, g(0, nj, nk)?, g(ni, nj, nk)?), + ); + let p = add(sub(faces, edges), corners); + if let Some(slot) = x.get_mut(idx(i, j, k)) { + *slot = Some(p); + } + } + } + } + let nodes: Vec

= x.into_iter().map(|p| p.ok_or_else(|| Error::Failed("mapped grid node missing".into()))).collect::>()?; + let id = |i: usize, j: usize, k: usize| idx(i, j, k) as u32; + let mut hexes = Vec::with_capacity(ni * nj * nk); + for k in 0..nk { + for j in 0..nj { + for i in 0..ni { + hexes.push([ + id(i, j, k), + id(i + 1, j, k), + id(i + 1, j + 1, k), + id(i, j + 1, k), + id(i, j, k + 1), + id(i + 1, j, k + 1), + id(i + 1, j + 1, k + 1), + id(i, j + 1, k + 1), + ]); + } + } + } + Ok(Raw { nodes, hexes, description: format!("mapped block of {ni} × {nj} × {nk}") }) +} diff --git a/crates/mesh/src/hex/mod.rs b/crates/mesh/src/hex/mod.rs new file mode 100644 index 0000000..431bf34 --- /dev/null +++ b/crates/mesh/src/hex/mod.rs @@ -0,0 +1,553 @@ +//! Structured all-hexahedral meshing of sweepable and mappable bodies. +//! +//! [`hex_mesh`] meshes one body's boundary surface (face-tagged, from surface preparation) with +//! hexahedra only, by the first method that fits: +//! 1. **Extrude** ([`layered`]): a straight multi-level sweep: every face is across the direction +//! (a source or target level) or along it (a linking wall), possibly at several levels with +//! regions of different heights (stepped blocks, ribs on plates, shelled boxes). The cross- +//! section is imprinted, quad-meshed region by region and stacked layer by layer. +//! 2. **Sweep** ([`sweep`]): a one-to-one sweep between a planar source and a planar target face +//! through rings of logically rectangular linking faces (bent tubes, lofts, partial +//! revolutions): linking faces get transfinite grids projected onto their surfaces, and every +//! layer's interior nodes come from the source mesh by the least-squares affine map of the +//! boundary rings (Roca and Sarrate) with boundary residual correction (Knupp). +//! 3. **Revolve** ([`layered`] in a rotating frame): bodies of revolution whose cross-section does +//! not reach the axis (rings, tori, grooved pulleys). +//! 4. **Mapped** ([`mapped`]): a logical cube (six four-sided faces, three edges at each of eight +//! corners), meshed by transfinite interpolation of its face grids (Gordon and Thiel). +//! +//! A body none of these fits is [`Error::Unsupported`] with each method's reason: the input for +//! decomposition into sweepable and mappable regions. +//! +//! **Common steps.** Interval assignment ([`interval`]) picks the division counts of shared curves +//! together. Boundary nodes on analytic faces are placed on the exact surface; quadratic hexahedra +//! (hex20) get their mid-edge nodes projected onto the faces (and onto both faces of a B-rep +//! edge). Every boundary quadrilateral is tagged with the B-rep face it lies on, so model sets +//! named by face resolve. The output passes the certified validity gate (every element's +//! Jacobian determinant proven positive, [`crate::quality`]), a tiling check (the elements' +//! volume equals the volume their boundary encloses, so nothing overlaps), and a fidelity check +//! (the volume differs from the boundary surface's by no more than the gap between the mesh +//! boundary and the surface accounts for). + +mod arrangement; +mod brep; +mod geom; +mod interval; +mod layered; +mod mapped; +mod region; +mod sweep; + +#[cfg(test)] +mod tests; + +use std::collections::HashMap; + +pub use geom::Surface; + +use brep::Brep; +use geom::{TriangleTree, add, dist, scale, unit}; + +use crate::quad::QuadOptions; +use crate::quality::{GateOptions, GateReport, Metric, element_quality, validity_gate_checks}; +use crate::size::SizeField; +use crate::surface::BoundarySurface; +use crate::vec3::{P, cross, dot, sub}; +use crate::{Error, Result}; +use modelcraft_model::ElementKind; + +/// Which structured method to use. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)] +pub enum HexMethod { + /// The first method that fits, in the order extrude, sweep, revolve, mapped. + Auto, + /// The first sweeping method that fits: extrude, sweep, revolve. + Swept, + Extrude, + Sweep, + Revolve, + Mapped, +} + +impl HexMethod { + pub fn name(self) -> &'static str { + match self { + HexMethod::Auto => "auto", + HexMethod::Swept => "swept", + HexMethod::Extrude => "extrude", + HexMethod::Sweep => "sweep", + HexMethod::Revolve => "revolve", + HexMethod::Mapped => "mapped", + } + } + + pub fn from_name(name: &str) -> Option { + [HexMethod::Auto, HexMethod::Swept, HexMethod::Extrude, HexMethod::Sweep, HexMethod::Revolve, HexMethod::Mapped] + .into_iter() + .find(|m| m.name() == name) + } +} + +/// Options of [`hex_mesh`]. +#[derive(Debug, Clone, PartialEq)] +pub struct HexOptions { + pub method: HexMethod, + /// Hex20 (mid-edge nodes on the geometry) instead of hex8. + pub quadratic: bool, + /// The chordal tolerance the boundary surface was tessellated with (0 when exact): copies of a + /// curve at different levels closer than this are one curve, and the fidelity check allows + /// it. + pub chordal_tolerance: f64, + /// The most hexahedra ([`Error::Limit`] beyond). + pub max_elements: usize, + /// Options of the face meshers. + pub quad: QuadOptions, +} + +impl Default for HexOptions { + fn default() -> Self { + HexOptions { method: HexMethod::Auto, quadratic: false, chordal_tolerance: 0.0, max_elements: 5_000_000, quad: QuadOptions::default() } + } +} + +/// A boundary face of the mesh and the B-rep face it lies on. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct BoundaryFace { + /// Element index (0-based, into [`HexMesh::element`]). + pub element: u32, + /// Face number (1-based, as in [`ElementKind::faces`]). + pub face: u8, + /// The B-rep face tag (`tri_face`). + pub brep_face: u32, +} + +/// How the mesh was made and how good it is. +#[derive(Debug, Clone, PartialEq)] +pub struct HexReport { + pub method: HexMethod, + /// What the method found (direction, levels, regions, layers, …). + pub description: String, + /// Smallest and mean scaled Jacobian (Knupp) over the elements. + pub min_scaled_jacobian: f64, + pub mean_scaled_jacobian: f64, + /// The elements' total volume, and the volume the boundary surface encloses. + pub volume: f64, + pub surface_volume: f64, + /// The validity gate's report. + pub gate: GateReport, +} + +/// An all-hexahedral mesh. +#[derive(Debug, Clone, PartialEq)] +pub struct HexMesh { + pub nodes: Vec<[f64; 3]>, + /// [`ElementKind::Hex8`] or [`ElementKind::Hex20`]. + pub kind: ElementKind, + /// Element node indices, `kind.node_count()` per element, in the model's node order. + pub connectivity: Vec, + pub boundary: Vec, + pub report: HexReport, +} + +impl HexMesh { + pub fn element_count(&self) -> usize { + self.connectivity.len() / self.kind.node_count() + } + + /// The node indices of element `e`. + pub fn element(&self, e: usize) -> Option<&[u32]> { + let n = self.kind.node_count(); + self.connectivity.get(e * n..(e + 1) * n) + } +} + +/// What every method gets. +pub(crate) struct Ctx<'a> { + pub brep: &'a Brep<'a>, + pub tree: &'a TriangleTree<'a>, + pub field: &'a dyn SizeField, + pub opts: &'a HexOptions, + /// Classification tolerance (model units). + pub tol: f64, +} + +/// What a method produces: linear hexahedra in the model's node order. +pub(crate) struct Raw { + pub nodes: Vec

, + pub hexes: Vec<[u32; 8]>, + pub description: String, +} + +/// Refinements of the cross-section curves tried when a face mesher fails (each halves them). +const MAX_REFINE: u32 = 2; + +/// Classification tolerance relative to the bounding-box diagonal. +const RELATIVE_TOLERANCE: f64 = 1e-6; + +/// Mesh the solid bounded by `surface` with hexahedra. `faces[tag]` is the analytic surface of +/// B-rep face `tag` when known (exact node placement; free-form faces use the triangles); +/// `field` sets the element size. A body no structured method fits is [`Error::Unsupported`] +/// with the reasons. +pub fn hex_mesh(surface: &BoundarySurface, faces: &[Option], field: &dyn SizeField, opts: &HexOptions) -> Result { + if !(opts.chordal_tolerance >= 0.0 && opts.chordal_tolerance.is_finite()) { + return Err(Error::Invalid(format!("the chordal tolerance must be a non-negative number, got {}", opts.chordal_tolerance))); + } + if surface.shells().len() != 1 { + return Err(Error::Unsupported(format!( + "the boundary has {} shells; structured meshing takes one solid without internal voids", + surface.shells().len() + ))); + } + // A mesh far beyond the element limit is refused before any work: the solid's volume over + // the cube of the smallest size the field asks for at the boundary. + let mut h_min = f64::INFINITY; + for p in surface.points() { + h_min = h_min.min(crate::size::checked_size(field, *p)?); + } + let estimate = surface.volume().abs() / h_min.powi(3); + if !(estimate <= 8.0 * opts.max_elements as f64) { + return Err(Error::Limit(format!("about {estimate:.3e} hexahedra at this size; the limit is {}", opts.max_elements))); + } + let brep = Brep::new(surface, faces)?; + let tree = TriangleTree::new(surface)?; + let ctx = Ctx { brep: &brep, tree: &tree, field, opts, tol: RELATIVE_TOLERANCE * brep.diagonal }; + let methods: &[HexMethod] = match opts.method { + HexMethod::Auto => &[HexMethod::Extrude, HexMethod::Sweep, HexMethod::Revolve, HexMethod::Mapped], + HexMethod::Swept => &[HexMethod::Extrude, HexMethod::Sweep, HexMethod::Revolve], + HexMethod::Extrude => &[HexMethod::Extrude], + HexMethod::Sweep => &[HexMethod::Sweep], + HexMethod::Revolve => &[HexMethod::Revolve], + HexMethod::Mapped => &[HexMethod::Mapped], + }; + let mut reasons: Vec = Vec::new(); + for &method in methods { + let mut last = None; + match method { + HexMethod::Extrude | HexMethod::Revolve => { + let setups = match candidate_setups(&ctx, method) { + Ok(s) => s, + Err(why) => { + reasons.push(format!("{}: {why}", method.name())); + continue; + } + }; + for setup in setups { + // The cross-section is laid out once per frame; refinement retries reuse it. + let prepared = match layered::prepare(&ctx, setup) { + Ok(p) => p, + Err(e @ Error::Limit(_)) => return Err(e), + Err(e) => { + last = Some(e); + continue; + } + }; + match with_refinement(&ctx, method, |refine| layered::mesh(&ctx, &prepared, refine)) { + Ok(mesh) => return Ok(mesh), + Err(e @ Error::Limit(_)) => return Err(e), + Err(e) => last = Some(e), + } + } + } + HexMethod::Sweep => match with_refinement(&ctx, method, |refine| sweep::mesh(&ctx, refine)) { + Ok(mesh) => return Ok(mesh), + Err(e @ Error::Limit(_)) => return Err(e), + Err(e) => last = Some(e), + }, + _ => match mapped::mesh(&ctx).and_then(|raw| finish(&ctx, raw, method)) { + Ok(mesh) => return Ok(mesh), + Err(e @ Error::Limit(_)) => return Err(e), + Err(e) => last = Some(e), + }, + } + if let Some(e) = last { + reasons.push(format!("{}: {}", method.name(), plain(&e))); + } + } + Err(Error::Unsupported(format!("no structured hex method fits this body ({})", reasons.join("; ")))) +} + +/// The frames a layered method can try, fewest levels first, or the first reason none fits. +fn candidate_setups(ctx: &Ctx, method: HexMethod) -> std::result::Result, String> { + let frames = if method == HexMethod::Extrude { layered::translations(ctx.brep) } else { layered::rotations(ctx.brep) }; + let mut fitting: Vec<(usize, usize, layered::Setup)> = Vec::new(); + let mut first_reason = None; + for (k, f) in frames.into_iter().enumerate() { + match layered::setup(ctx, &f) { + Ok(s) => fitting.push((s.level_count(), k, s)), + Err(e) => { + first_reason.get_or_insert_with(|| plain(&e)); + } + } + } + if fitting.is_empty() { + return Err(first_reason.unwrap_or_else(|| if method == HexMethod::Extrude { "no direction to try".into() } else { "no axis to try".into() })); + } + fitting.sort_by_key(|x| (x.0, x.1)); + Ok(fitting.into_iter().take(layered::MAX_CANDIDATES).map(|x| x.2).collect()) +} + +/// Run a method's mesher and [`finish`], and when that fails (a face mesher on a boundary too +/// coarse for it, or curved elements the gate rejects) again with the cross-section's curves +/// halved, up to [`MAX_REFINE`] times. Limits and unsupported shapes end it at once. +fn with_refinement(ctx: &Ctx, method: HexMethod, run: impl Fn(u32) -> Result) -> Result { + let mut last = Error::Failed("no attempt".into()); + for refine in 0..=MAX_REFINE { + match run(refine).and_then(|raw| finish(ctx, raw, method)) { + Ok(mesh) => return Ok(mesh), + Err(e @ Error::Failed(_)) => last = e, + Err(e) => return Err(e), + } + } + Err(last) +} + +/// An error's message without the "not supported yet" prefix (for the combined reasons). +fn plain(e: &Error) -> String { + match e { + Error::Unsupported(m) => m.clone(), + other => other.to_string(), + } +} + +/// Hex face `k` (0-based) corner positions in the element. +fn hex_face(k: usize) -> [usize; 4] { + ElementKind::Hex8.faces().get(k).and_then(|f| <[usize; 4]>::try_from(*f).ok()).unwrap_or([0, 1, 2, 3]) +} + +/// Tag, check, upgrade to hex20 when asked, gate and measure. +fn finish(ctx: &Ctx, raw: Raw, method: HexMethod) -> Result { + let Raw { mut nodes, hexes, description } = raw; + if hexes.is_empty() { + return Err(Error::Failed("the method produced no elements".into())); + } + if hexes.len() > ctx.opts.max_elements { + return Err(Error::Limit(format!("the mesh has {} hexahedra; the limit is {}", hexes.len(), ctx.opts.max_elements))); + } + modelcraft_model::check_coordinates(nodes.as_flattened())?; + let at = |nodes: &[P], i: u32| nodes.get(i as usize).copied().ok_or_else(|| Error::Failed(format!("node {i} out of range"))); + // Boundary faces: hex faces used once. + let mut uses: HashMap<[u32; 4], u32> = HashMap::with_capacity(hexes.len() * 3); + for h in &hexes { + for k in 0..6 { + let mut key = hex_face(k).map(|i| h.get(i).copied().unwrap_or(0)); + key.sort_unstable(); + *uses.entry(key).or_insert(0) += 1; + } + } + if let Some(n) = uses.values().copied().max().filter(|&n| n > 2) { + return Err(Error::Failed(format!("a hexahedron face is shared by {n} elements (the mesh is not a manifold)"))); + } + let mut boundary = Vec::new(); + let mut fidelity = 0.0; + for (e, h) in hexes.iter().enumerate() { + for k in 0..6 { + let corners = hex_face(k).map(|i| h.get(i).copied().unwrap_or(0)); + let mut key = corners; + key.sort_unstable(); + if uses.get(&key) != Some(&1) { + continue; + } + let [a, b, c, d] = [at(&nodes, corners[0])?, at(&nodes, corners[1])?, at(&nodes, corners[2])?, at(&nodes, corners[3])?]; + let centre = scale(add(add(a, b), add(c, d)), 0.25); + let n = cross(sub(c, a), sub(d, b)); + let area = 0.5 * dot(n, n).sqrt(); + let Some(n) = unit(n) else { return Err(Error::Failed(format!("element {e} has a degenerate boundary face"))) }; + let hit = ctx + .tree + .closest(centre, &|t| dot(ctx.tree.normal(t), n) > 0.0) + .ok_or_else(|| Error::Failed(format!("no boundary face faces the same way as element {e}'s face {}", k + 1)))?; + let f = ctx.brep.face_of_triangle(hit.triangle).ok_or_else(|| Error::Failed("triangle without a face".into()))?; + fidelity += area * (hit.dist + ctx.opts.chordal_tolerance); + boundary.push(BoundaryFace { element: e as u32, face: k as u8 + 1, brep_face: ctx.brep.face_tag(f) }); + } + } + let kind = if ctx.opts.quadratic { ElementKind::Hex20 } else { ElementKind::Hex8 }; + let mut connectivity: Vec = Vec::with_capacity(hexes.len() * kind.node_count()); + if ctx.opts.quadratic { + let mids = mid_edge_nodes(ctx, &mut nodes, &hexes, &boundary)?; + for (h, m) in hexes.iter().zip(&mids) { + connectivity.extend_from_slice(h); + connectivity.extend_from_slice(m); + } + } else { + for h in &hexes { + connectivity.extend_from_slice(h); + } + } + let npe = kind.node_count(); + let coords = |e: usize| -> Result> { connectivity.get(e * npe..(e + 1) * npe).unwrap_or(&[]).iter().map(|&i| at(&nodes, i)).collect() }; + // The volume the mesh's own boundary encloses: equal to the elements' total exactly when they + // tile it without overlapping. + let mut enclosed = 0.0; + for b in &boundary { + let x = coords(b.element as usize)?; + let local = kind.faces().get(b.face as usize - 1).copied().unwrap_or(&[]); + let face: Vec

= local.iter().filter_map(|&i| x.get(i).copied()).collect(); + enclosed += face_flux(&face); + } + // One certified pass: every element's quality, with its Jacobian check for the gate. + let (mut min_sj, mut sum_sj) = (f64::INFINITY, 0.0); + let mut checks = Vec::with_capacity(hexes.len()); + for e in 0..hexes.len() { + let q = element_quality(kind, &coords(e)?)?; + let sj = q.get(Metric::ScaledJacobian).unwrap_or(-1.0); + min_sj = min_sj.min(sj); + sum_sj += sj; + checks.push((e as u32 + 1, q.jacobian.ok_or_else(|| Error::Failed("an element was not certified".into()))?)); + } + let gate = validity_gate_checks(checks.iter().map(|(l, c)| (*l, c)), &GateOptions { expected_volume: Some(enclosed), volume_tolerance: 1e-9 })? + .into_result() + .map_err(|e| Error::Failed(format!("{} ({description})", plain(&e))))?; + let surface_volume = ctx.brep.surface.volume(); + let slack = 2.0 * fidelity + 1e-9 * surface_volume.abs(); + if (gate.volume - surface_volume).abs() > slack { + return Err(Error::Failed(format!( + "the mesh volume {} differs from the boundary's {surface_volume} by more than the gap between them allows ({slack}): the mesh does not fill the body ({description})", + gate.volume + ))); + } + let report = HexReport { + method, + description, + min_scaled_jacobian: min_sj, + mean_scaled_jacobian: sum_sj / hexes.len().max(1) as f64, + volume: gate.volume, + surface_volume, + gate, + }; + Ok(HexMesh { nodes, kind, connectivity, boundary, report }) +} + +/// `∫ x n_x dA` over a quad4 or quad8 face (outward normal): its contribution to the enclosed +/// volume. The 3 × 3 Gauss rule integrates it exactly (degree at most 5 per direction). +fn face_flux(face: &[P]) -> f64 { + const G: [(f64, f64); 3] = [(-0.774_596_669_241_483_4, 5.0 / 9.0), (0.0, 8.0 / 9.0), (0.774_596_669_241_483_4, 5.0 / 9.0)]; + let mut sum = 0.0; + for (r, wr) in G { + for (s, ws) in G { + let [n, dr, ds] = quad_shape(face.len() == 8, r, s); + let (mut x, mut xr, mut xs) = ([0.0; 3], [0.0; 3], [0.0; 3]); + for (k, p) in face.iter().enumerate().take(8) { + x = add(x, scale(*p, n[k])); + xr = add(xr, scale(*p, dr[k])); + xs = add(xs, scale(*p, ds[k])); + } + sum += wr * ws * x[0] * cross(xr, xs)[0]; + } + } + sum +} + +/// Shape functions and their `r`, `s` derivatives of quad4 (the first four entries) or quad8 at +/// `(r, s)`. +fn quad_shape(quadratic: bool, r: f64, s: f64) -> [[f64; 8]; 3] { + const C: [(f64, f64); 4] = [(-1.0, -1.0), (1.0, -1.0), (1.0, 1.0), (-1.0, 1.0)]; + let mut out = [[0.0; 8]; 3]; + for (k, (ri, si)) in C.iter().enumerate() { + if quadratic { + out[0][k] = 0.25 * (1.0 + ri * r) * (1.0 + si * s) * (ri * r + si * s - 1.0); + out[1][k] = 0.25 * ri * (1.0 + si * s) * (2.0 * ri * r + si * s); + out[2][k] = 0.25 * si * (1.0 + ri * r) * (ri * r + 2.0 * si * s); + } else { + out[0][k] = 0.25 * (1.0 + ri * r) * (1.0 + si * s); + out[1][k] = 0.25 * ri * (1.0 + si * s); + out[2][k] = 0.25 * si * (1.0 + ri * r); + } + } + if quadratic { + // Mid-side nodes on sides 0-1 (s = −1), 1-2 (r = 1), 2-3 (s = 1), 3-0 (r = −1). + out[0][4..].copy_from_slice(&[ + 0.5 * (1.0 - r * r) * (1.0 - s), + 0.5 * (1.0 + r) * (1.0 - s * s), + 0.5 * (1.0 - r * r) * (1.0 + s), + 0.5 * (1.0 - r) * (1.0 - s * s), + ]); + out[1][4..].copy_from_slice(&[-r * (1.0 - s), 0.5 * (1.0 - s * s), -r * (1.0 + s), -0.5 * (1.0 - s * s)]); + out[2][4..].copy_from_slice(&[-0.5 * (1.0 - r * r), -s * (1.0 + r), 0.5 * (1.0 - r * r), -s * (1.0 - r)]); + } + out +} + +/// Projection rounds for a point on several faces (alternating projections onto each). +const PROJECTION_ROUNDS: usize = 24; + +/// The closest point to `p` on B-rep face index `f`: its analytic surface, or its triangles. +pub(crate) fn project_to_face(ctx: &Ctx, f: usize, p: P) -> Option

{ + let face = ctx.brep.faces.get(f)?; + if let Some(q) = face.surface.and_then(|s| s.project(p)) { + return Some(q); + } + let tag = face.id; + let tri_face = ctx.brep.surface.tri_face(); + ctx.tree.closest(p, &|t| tri_face.get(t as usize) == Some(&tag)).map(|h| h.point) +} + +/// The point on every listed face (their common curve or corner) nearest `p`, by alternating +/// projections. +pub(crate) fn project_to_faces(ctx: &Ctx, faces: &[usize], p: P) -> P { + let mut x = p; + for _ in 0..PROJECTION_ROUNDS { + let before = x; + for &f in faces { + if let Some(q) = project_to_face(ctx, f, x) { + x = q; + } + } + if faces.len() < 2 || dist(before, x) <= ctx.tol * 1e-3 { + break; + } + } + x +} + +/// Mid-edge nodes of every hexahedron (Hex20 order: positions 8..20), shared between elements; +/// boundary edges get theirs on the B-rep faces of the boundary quads around them. +fn mid_edge_nodes(ctx: &Ctx, nodes: &mut Vec

, hexes: &[[u32; 8]], boundary: &[BoundaryFace]) -> Result> { + // B-rep faces at each boundary edge. + let mut edge_faces: HashMap<(u32, u32), Vec> = HashMap::new(); + for b in boundary { + let Some(h) = hexes.get(b.element as usize) else { continue }; + let c = hex_face(b.face as usize - 1).map(|i| h.get(i).copied().unwrap_or(0)); + let f = ctx.brep.faces.binary_search_by_key(&b.brep_face, |x| x.id).map_err(|_| Error::Failed("boundary tag".into()))?; + for k in 0..4 { + let (a, z) = (c[k], c[(k + 1) % 4]); + let list = edge_faces.entry((a.min(z), a.max(z))).or_default(); + if !list.contains(&f) { + list.push(f); + } + } + } + let edges = ElementKind::Hex20.edges(); + let mut mid: HashMap<(u32, u32), u32> = HashMap::new(); + let mut out = Vec::with_capacity(hexes.len()); + for h in hexes { + let mut m = [0u32; 12]; + for (k, e) in edges.iter().enumerate() { + let (Some(&i), Some(&j)) = (e.first().and_then(|&x| h.get(x)), e.get(1).and_then(|&x| h.get(x))) else { + return Err(Error::Failed("hex20 edge table".into())); + }; + let key = (i.min(j), i.max(j)); + let id = match mid.get(&key) { + Some(&id) => id, + None => { + let (a, b) = (nodes.get(i as usize).copied().unwrap_or([0.0; 3]), nodes.get(j as usize).copied().unwrap_or([0.0; 3])); + let p = scale(add(a, b), 0.5); + let p = match edge_faces.get(&key) { + Some(fs) => project_to_faces(ctx, fs, p), + None => p, + }; + nodes.push(p); + let id = nodes.len() as u32 - 1; + mid.insert(key, id); + id + } + }; + if let Some(slot) = m.get_mut(k) { + *slot = id; + } + } + out.push(m); + } + modelcraft_model::check_coordinates(nodes.as_flattened())?; + Ok(out) +} diff --git a/crates/mesh/src/hex/region.rs b/crates/mesh/src/hex/region.rs new file mode 100644 index 0000000..f5e5ed3 --- /dev/null +++ b/crates/mesh/src/hex/region.rs @@ -0,0 +1,387 @@ +//! Planar face meshing for the sweepers: the scheme a region gets, the interval constraints the +//! scheme needs, and the quadrilateral mesh on the assigned boundary nodes. +//! +//! - **Mapped**: one loop with four corners (vertex classification by interior angle, as in +//! submapping: Whiteley, White, Benzley, Blacker, *Two and three-quarter dimensional meshing +//! facilitators*, Engineering with Computers 12:144–154, 1996): opposite sides equal, +//! transfinite interpolation ([`crate::quad::mapped_quad_mesh`]). +//! - **Annulus**: two smooth loops (no corners), one inside the other: the two loops get equal +//! counts and the ring between them is a periodic structured grid, each radial line joining a +//! node of the outer loop to its partner on the inner loop (the rotation that pairs them most +//! closely) — the O-grid of a tube's cross-section. +//! - **Disk**: one smooth loop (no corners): an O-grid ("butterfly"): a mapped square in the +//! middle, its corners halfway from the centroid to the loop's quarter points, and rings of +//! quadrilaterals from the square out to the loop. The loop's count is a multiple of four. +//! - **Free**: anything else: even loops for the free all-quad mesher +//! ([`crate::quad::free_quad_mesh`]). +//! +//! A structured scheme whose grid would fold falls back to the free mesher (its constraints imply +//! even loops). + +use std::f64::consts::PI; + +use super::interval::Intervals; +use crate::quad::{QuadMesh, QuadOptions, QuadStats, free_quad_mesh, mapped_quad_mesh, transfinite_grid}; +use crate::quality::{Validity, planar_quad_validity}; +use crate::size::SizeField; +use crate::{Error, Result}; + +type Q = [f64; 2]; + +/// The most nodes the annulus pairing compares (its cost is quadratic in the loop length). +const MAX_ANNULUS_NODES: usize = 4096; + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub(crate) enum Scheme { + /// Corner positions (curve uses) in the single loop. + Mapped([usize; 4]), + Annulus, + Disk, + Free, +} + +/// The interior angle at every curve end of a loop, given each curve use's oriented polyline +/// (with the loop counter-clockwise when `ccw`). +fn angles(lp: &[Vec], ccw: bool) -> Option> { + let n = lp.len(); + let mut out = Vec::with_capacity(n); + for k in 0..n { + let prev = lp.get((k + n - 1) % n)?; + let next = lp.get(k)?; + let (a, o, b) = (*prev.get(prev.len().checked_sub(2)?)?, *next.first()?, *next.get(1)?); + let (u, v) = ([o[0] - a[0], o[1] - a[1]], [b[0] - o[0], b[1] - o[1]]); + let turn = (u[0] * v[1] - u[1] * v[0]).atan2(u[0] * v[0] + u[1] * v[1]); + out.push(PI - if ccw { turn } else { -turn }); + } + Some(out) +} + +fn signed_area(pts: &[Q]) -> f64 { + let n = pts.len(); + 0.5 * (0..n).map(|k| pts[k][0] * pts[(k + 1) % n][1] - pts[(k + 1) % n][0] * pts[k][1]).sum::() +} + +/// Choose the scheme of a region from its loops: per loop, the oriented polyline of each curve +/// use (loop 0 outer). +pub(crate) fn classify(loops: &[Vec>]) -> Scheme { + let ccw = |lp: &Vec>| { + let pts: Vec = lp.iter().flat_map(|c| c.iter().take(c.len().saturating_sub(1)).copied()).collect(); + signed_area(&pts) > 0.0 + }; + match loops { + [lp] => { + let Some(a) = angles(lp, ccw(lp)) else { return Scheme::Free }; + if a.iter().any(|&x| x >= 1.25 * PI) { + return Scheme::Free; + } + let corners: Vec = (0..a.len()).filter(|&k| a.get(k).is_some_and(|&x| x <= 0.75 * PI)).collect(); + if corners.is_empty() { + return Scheme::Disk; + } + corners.try_into().map(Scheme::Mapped).unwrap_or(Scheme::Free) + } + [outer, inner] => { + let smooth = |lp: &Vec>| angles(lp, ccw(lp)).is_some_and(|a| a.iter().all(|&x| x > 0.75 * PI && x < 1.25 * PI)); + if smooth(outer) && smooth(inner) { Scheme::Annulus } else { Scheme::Free } + } + _ => Scheme::Free, + } +} + +/// Add the scheme's constraints on the loops' curve variables. +pub(crate) fn constrain(ia: &mut Intervals, loops: &[Vec], scheme: Scheme) { + match scheme { + Scheme::Mapped(c) => { + let Some(lp) = loops.first() else { return }; + let side = |a: usize, b: usize| (a..b).filter_map(|k| lp.get(k).copied()).collect::>(); + ia.equal.push((side(c[0], c[1]), side(c[2], c[3]))); + ia.equal.push((side(c[1], c[2]), side(c[3], lp.len()).into_iter().chain(side(0, c[0])).collect())); + } + Scheme::Annulus => { + if let (Some(o), Some(i)) = (loops.first(), loops.get(1)) { + ia.equal.push((o.clone(), i.clone())); + ia.multiple.push((o.clone(), 2)); + ia.at_least.push((o.clone(), 6)); + } + } + Scheme::Disk => { + if let Some(lp) = loops.first() { + ia.multiple.push((lp.clone(), 4)); + ia.at_least.push((lp.clone(), 8)); + } + } + Scheme::Free => { + let min = if loops.len() > 1 { 6 } else { 4 }; + for lp in loops { + ia.multiple.push((lp.clone(), 2)); + ia.at_least.push((lp.clone(), min)); + } + } + } +} + +/// Mesh a region on its boundary nodes: `loops[k]` are loop `k`'s nodes in order, `starts[k]` the +/// positions where its curve uses start. Returns the mesh (boundary nodes first, in the order of +/// the loops) and the scheme used. +pub(crate) fn mesh(loops: &[Vec], starts: &[Vec], scheme: Scheme, size: &dyn SizeField, opts: &QuadOptions) -> Result<(QuadMesh, Scheme)> { + let structured = match scheme { + Scheme::Mapped(c) => { + let s0 = starts.first(); + let pos = c.map(|k| s0.and_then(|s| s.get(k)).copied().unwrap_or(0)); + match loops.first().map(|b| mapped_quad_mesh(b, pos, opts)) { + Some(Ok(m)) => Some(m.mesh), + Some(Err(e @ Error::Limit(_))) => return Err(e), + _ => None, + } + } + Scheme::Annulus => annulus(loops, size, opts)?, + Scheme::Disk => disk(loops, size, opts)?, + Scheme::Free => None, + }; + match structured { + Some(m) => Ok((m, scheme)), + None => Ok((free_quad_mesh(loops, size, opts)?, Scheme::Free)), + } +} + +/// The O-grid between two loops of equal length, or `None` when it would fold. +fn annulus(loops: &[Vec], size: &dyn SizeField, opts: &QuadOptions) -> Result> { + let [outer, inner] = loops else { return Ok(None) }; + let n = outer.len(); + if n != inner.len() || !(3..=MAX_ANNULUS_NODES).contains(&n) { + return Ok(None); + } + // Index maps that make both loops counter-clockwise. + let ccw_order = |pts: &[Q]| -> Vec { if signed_area(pts) > 0.0 { (0..pts.len()).collect() } else { (0..pts.len()).rev().collect() } }; + let (oi, ii) = (ccw_order(outer), ccw_order(inner)); + let o = |k: usize| oi.get(k % n).and_then(|&i| outer.get(i)).copied().unwrap_or([0.0; 2]); + let inn = |k: usize| ii.get(k % n).and_then(|&i| inner.get(i)).copied().unwrap_or([0.0; 2]); + let d2 = |a: Q, b: Q| (a[0] - b[0]).powi(2) + (a[1] - b[1]).powi(2); + let shift = (0..n) + .min_by(|&s, &t| { + let cost = |s: usize| (0..n).map(|k| d2(o(k), inn(k + s))).sum::(); + cost(s).total_cmp(&cost(t)).then(s.cmp(&t)) + }) + .unwrap_or(0); + let mean_gap = (0..n).map(|k| d2(o(k), inn(k + shift)).sqrt()).sum::() / n as f64; + let mid = |k: usize| { + let (a, b) = (o(k), inn(k + shift)); + [(a[0] + b[0]) * 0.5, (a[1] + b[1]) * 0.5, 0.0] + }; + let h = crate::size::checked_size(size, mid(0))?; + let m = (mean_gap / h).round().max(1.0) as usize; + let quads_count = n.saturating_mul(m); + if quads_count > opts.max_elements { + return Err(Error::Limit(format!("the mesh would have {quads_count} quadrilaterals; the limit is {}", opts.max_elements))); + } + // Node ids: outer loop 0..n, inner loop n..2n (input order), then the rings between. + let outer_ids: Vec = (0..n).map(|k| oi.get(k % n).copied().unwrap_or(0) as u32).collect(); + let inner_ids: Vec = (0..n).map(|k| (n + ii.get((k + shift) % n).copied().unwrap_or(0)) as u32).collect(); + let mut nodes: Vec = outer.iter().chain(inner.iter()).copied().collect(); + let mut quads = Vec::with_capacity(quads_count); + if !rings(&mut nodes, &outer_ids, &inner_ids, m, &mut quads)? { + return Ok(None); + } + Ok(Some(QuadMesh { nodes, quads, boundary_nodes: 2 * n, stats: QuadStats::default() })) +} + +/// Rings of quadrilaterals between a closed outer node sequence (counter-clockwise) and an inner +/// one of the same length, node `k` of each joined by a straight radial line divided into `m`. +/// `false` when a quadrilateral would not be strictly convex. +fn rings(nodes: &mut Vec, outer: &[u32], inner: &[u32], m: usize, quads: &mut Vec<[u32; 4]>) -> Result { + let n = outer.len(); + let at = |nodes: &[Q], i: u32| nodes.get(i as usize).copied().unwrap_or([0.0; 2]); + let mut ring_ids: Vec> = vec![outer.to_vec()]; + for r in 1..m { + let t = r as f64 / m as f64; + let mut ids = Vec::with_capacity(n); + for k in 0..n { + let (a, b) = (at(nodes, outer.get(k).copied().unwrap_or(0)), at(nodes, inner.get(k).copied().unwrap_or(0))); + nodes.push([a[0] + t * (b[0] - a[0]), a[1] + t * (b[1] - a[1])]); + ids.push(nodes.len() as u32 - 1); + } + ring_ids.push(ids); + } + ring_ids.push(inner.to_vec()); + for w in ring_ids.windows(2) { + let (a, b) = (&w[0], &w[1]); + for k in 0..n { + let (Some(&p0), Some(&p1), Some(&q1), Some(&q0)) = (a.get(k), a.get((k + 1) % n), b.get((k + 1) % n), b.get(k)) else { + return Ok(false); + }; + // Outer ring counter-clockwise, inner ring inside it: (p0, p1, q1, q0) turns left. + let q = [p0, p1, q1, q0]; + if planar_quad_validity(&q.map(|i| at(nodes, i)))? != Validity::Valid { + return Ok(false); + } + quads.push(q); + } + } + Ok(true) +} + +/// The O-grid of one smooth loop whose length is a multiple of four (at least 8), or `None` when +/// it would fold. +fn disk(loops: &[Vec], size: &dyn SizeField, opts: &QuadOptions) -> Result> { + let [outer] = loops else { return Ok(None) }; + let n = outer.len(); + if n < 8 || !n.is_multiple_of(4) || n > MAX_ANNULUS_NODES { + return Ok(None); + } + let q = n / 4; + let ccw = signed_area(outer) > 0.0; + let order: Vec = if ccw { (0..n as u32).collect() } else { (0..n as u32).rev().collect() }; + let at = |nodes: &[Q], i: u32| nodes.get(i as usize).copied().unwrap_or([0.0; 2]); + // The area centroid. + let (mut cx, mut cy, mut a2) = (0.0, 0.0, 0.0); + for k in 0..n { + let (p, r) = (outer[k], outer[(k + 1) % n]); + let c = p[0] * r[1] - r[0] * p[1]; + a2 += c; + cx += (p[0] + r[0]) * c; + cy += (p[1] + r[1]) * c; + } + if !(a2.abs() > 0.0) { + return Ok(None); + } + let centre = [cx / (3.0 * a2), cy / (3.0 * a2)]; + let mut nodes: Vec = outer.clone(); + // The middle square's corners, halfway to the loop's quarter points, and its perimeter. + let corner = |k: usize| { + let b = at(outer, order.get((k * q) % n).copied().unwrap_or(0)); + [centre[0] + 0.5 * (b[0] - centre[0]), centre[1] + 0.5 * (b[1] - centre[1])] + }; + let mut square: Vec = Vec::with_capacity(n); + for k in 0..4 { + let (c0, c1) = (corner(k), corner(k + 1)); + for j in 0..q { + let t = j as f64 / q as f64; + nodes.push([c0[0] + t * (c1[0] - c0[0]), c0[1] + t * (c1[1] - c0[1])]); + square.push(nodes.len() as u32 - 1); + } + } + let gap = (0..n).map(|k| { + let (a, b) = (at(&nodes, order.get(k).copied().unwrap_or(0)), at(&nodes, square.get(k).copied().unwrap_or(0))); + ((a[0] - b[0]).powi(2) + (a[1] - b[1]).powi(2)).sqrt() + }); + let mean_gap = gap.sum::() / n as f64; + let h = crate::size::checked_size(size, [centre[0], centre[1], 0.0])?; + let m = (mean_gap / h).round().max(1.0) as usize; + let count = n.saturating_mul(m).saturating_add(q * q); + if count > opts.max_elements { + return Err(Error::Limit(format!("the mesh would have {count} quadrilaterals; the limit is {}", opts.max_elements))); + } + let mut quads = Vec::with_capacity(count); + if !rings(&mut nodes, &order, &square, m, &mut quads)? { + return Ok(None); + } + // The square: a q × q transfinite grid on its perimeter (counter-clockwise from corner 0). + let side = |k: usize, rev: bool| -> Vec { + let mut v: Vec = (0..=q).map(|j| square.get((k * q + j) % n).copied().unwrap_or(0)).collect(); + if rev { + v.reverse(); + } + v + }; + let (bottom, right, top, left) = (side(0, false), side(1, false), side(2, true), side(3, true)); + let pts = |ids: &[u32], nodes: &[Q]| ids.iter().map(|&i| at(nodes, i)).collect::>(); + let grid = transfinite_grid(&pts(&bottom, &nodes), &pts(&right, &nodes), &pts(&top, &nodes), &pts(&left, &nodes))?; + let mut id = vec![0u32; (q + 1) * (q + 1)]; + for j in 0..=q { + for i in 0..=q { + let v = if j == 0 { + bottom.get(i).copied() + } else if j == q { + top.get(i).copied() + } else if i == 0 { + left.get(j).copied() + } else if i == q { + right.get(j).copied() + } else { + nodes.push(grid.get(j * (q + 1) + i).copied().unwrap_or([0.0; 2])); + Some(nodes.len() as u32 - 1) + }; + if let (Some(slot), Some(v)) = (id.get_mut(j * (q + 1) + i), v) { + *slot = v; + } + } + } + let g = |i: usize, j: usize| id.get(j * (q + 1) + i).copied().unwrap_or(0); + for j in 0..q { + for i in 0..q { + let quad = [g(i, j), g(i + 1, j), g(i + 1, j + 1), g(i, j + 1)]; + if planar_quad_validity(&quad.map(|v| at(&nodes, v)))? != Validity::Valid { + return Ok(None); + } + quads.push(quad); + } + } + Ok(Some(QuadMesh { nodes, quads, boundary_nodes: n, stats: QuadStats::default() })) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::size::UniformSize; + + fn circle(r: f64, n: usize, ccw: bool) -> Vec { + (0..n) + .map(|k| { + let t = 2.0 * PI * k as f64 / n as f64 * if ccw { 1.0 } else { -1.0 }; + [r * t.cos(), r * t.sin()] + }) + .collect() + } + + #[test] + fn rings_get_an_o_grid() { + let (outer, inner) = (circle(5.0, 16, true), circle(4.0, 16, false)); + let polyline = |pts: &Vec| { + let mut p = pts.clone(); + p.push(pts[0]); + vec![p] + }; + assert_eq!(classify(&[polyline(&outer), polyline(&inner)]), Scheme::Annulus); + let (m, s) = mesh(&[outer.clone(), inner.clone()], &[vec![0], vec![0]], Scheme::Annulus, &UniformSize(0.5), &QuadOptions::default()).unwrap(); + assert_eq!(s, Scheme::Annulus); + assert_eq!(m.quads.len(), 16 * 2); + assert_eq!(m.boundary_nodes, 32); + assert_eq!(&m.nodes[..16], &outer[..]); + // Unequal loops fall back to the free mesher. + let (_, s) = + mesh(&[circle(5.0, 16, true), circle(2.0, 8, true)], &[vec![0], vec![0]], Scheme::Annulus, &UniformSize(1.0), &QuadOptions::default()) + .unwrap(); + assert_eq!(s, Scheme::Free); + } + + #[test] + fn disks_get_an_o_grid() { + let c = circle(5.0, 16, false); + let mut closed = c.clone(); + closed.push(c[0]); + assert_eq!(classify(&[vec![closed]]), Scheme::Disk); + let (m, s) = mesh(std::slice::from_ref(&c), &[vec![0]], Scheme::Disk, &UniformSize(1.5), &QuadOptions::default()).unwrap(); + assert_eq!(s, Scheme::Disk); + assert_eq!(&m.nodes[..16], &c[..]); + // 16 · 2 ring quadrilaterals and the 4 × 4 square. + assert_eq!(m.quads.len(), 16 * 2 + 16); + let area: f64 = m.quads.iter().map(|q| signed_area(&q.map(|i| m.nodes[i as usize]))).sum(); + assert!((area - signed_area(&circle(5.0, 16, true))).abs() < 1e-9, "{area}"); + } + + #[test] + fn squares_are_mapped() { + let sq: Vec> = + vec![vec![[0.0, 0.0], [1.0, 0.0]], vec![[1.0, 0.0], [1.0, 1.0]], vec![[1.0, 1.0], [0.0, 1.0]], vec![[0.0, 1.0], [0.0, 0.0]]]; + assert_eq!(classify(std::slice::from_ref(&sq)), Scheme::Mapped([0, 1, 2, 3])); + let l: Vec> = vec![ + vec![[0.0, 0.0], [2.0, 0.0]], + vec![[2.0, 0.0], [2.0, 1.0]], + vec![[2.0, 1.0], [1.0, 1.0]], + vec![[1.0, 1.0], [1.0, 2.0]], + vec![[1.0, 2.0], [0.0, 2.0]], + vec![[0.0, 2.0], [0.0, 0.0]], + ]; + assert_eq!(classify(&[l]), Scheme::Free); + } +} diff --git a/crates/mesh/src/hex/sweep.rs b/crates/mesh/src/hex/sweep.rs new file mode 100644 index 0000000..1eaa38a --- /dev/null +++ b/crates/mesh/src/hex/sweep.rs @@ -0,0 +1,506 @@ +//! One-to-one sweeping from a planar source face to a planar target face through rings of +//! linking faces. +//! +//! **Recognition** (Blacker, *The Cooper tool*, Proc. 5th International Meshing Roundtable, +//! 1996). Starting from the source's loops, each loop edge leads to the linking face across it. +//! A linking face has one loop made of a *bottom* chain (the edges it shares with the previous +//! ring), one *side* edge, a *top* chain and another side edge; its top chains form the next ring. +//! Rings are followed until every top edge of a loop borders one planar face, the target. Every +//! face but the source and the target must be used exactly once. +//! +//! **Intervals** ([`super::interval`]): on each linking face the bottom and top chains have equal +//! sums and the two sides are equal; the source is mapped when it is a four-cornered single loop, +//! otherwise its loops are even for the free all-quad mesher. Every ring has the same number of +//! layers per stage. +//! +//! **Placement.** Edge nodes lie on both faces of their edge; each linking face gets the +//! transfinite grid of its four sides (Gordon and Hall 1973) projected onto its surface. The +//! interior nodes of each layer follow from the source mesh by the least-squares affine map that +//! carries the source boundary ring onto the layer's ring (Roca and Sarrate, *An automatic and +//! general least-squares projection procedure for sweep meshing*, Engineering with Computers +//! 26(4):391–406, 2010), plus the map's residuals at the ring interpolated by inverse-distance +//! weights (the boundary-residual correction of Knupp, *Next-generation sweep tool: a method for +//! generating all-hex meshes on two-and-one-half dimensional geometries*, Proc. 7th International +//! Meshing Roundtable, 1998). The last layer is projected onto the target face. + +use std::collections::{BTreeMap, BTreeSet}; + +use super::brep::Brep; +use super::geom::{add, scale}; +use super::interval::Intervals; +use super::layered::{Frame, PlaneField}; +use super::region; +use super::{Ctx, Raw, project_to_face, project_to_faces}; +use crate::quad::transfinite_grid; +use crate::size::{divide_polyline, ideal_divisions}; +use crate::vec3::{P, sub}; +use crate::{Error, Result}; + +/// The most source faces tried. +const MAX_SOURCES: usize = 8; +/// The most interior-node × ring-node × layer products the residual correction may spend. +const MAX_CORRECTION_WORK: u64 = 2_000_000_000; + +/// A linking face of one stage: its bottom chain, sides and top chain (oriented edges as the face +/// traverses them). +#[derive(Debug, Clone)] +struct Link { + face: usize, + bottom: Vec<(usize, bool)>, + right: (usize, bool), + top: Vec<(usize, bool)>, + left: (usize, bool), +} + +/// The sweep's topology: per source loop, its stages, each a ring of linking faces in loop order. +struct Topology { + source: usize, + target: usize, + rings: Vec>>, +} + +/// Sweep the body, with the source's curves `2^refine` times finer than the size asks. +pub(crate) fn mesh(ctx: &Ctx, refine: u32) -> Result { + let brep = ctx.brep; + let mut planar: Vec = (0..brep.faces.len()).filter(|&f| brep.faces.get(f).is_some_and(|x| x.plane.is_some())).collect(); + if planar.is_empty() { + return Err(Error::Unsupported("no planar face to sweep from".into())); + } + let area = |f: usize| brep.faces.get(f).map_or(0.0, |x| x.area); + planar.sort_by(|&a, &b| area(b).total_cmp(&area(a)).then(a.cmp(&b))); + let mut reasons = Vec::new(); + for &s in planar.iter().take(MAX_SOURCES) { + match topology(brep, s) { + Ok(topo) => return build(ctx, &topo, refine), + Err(e) => reasons.push(format!("from face {}: {e}", brep.face_tag(s))), + } + } + Err(Error::Unsupported(reasons.into_iter().take(2).collect::>().join("; "))) +} + +/// Follow the rings of linking faces from source face `s`. +fn topology(brep: &Brep, s: usize) -> std::result::Result { + let source = brep.faces.get(s).ok_or("no such face")?; + let mut used: BTreeSet = BTreeSet::from([s]); + let mut rings = Vec::new(); + let mut target: Option = None; + for front0 in &source.loops { + // Front edges as the face below them traverses them, with that face. + let mut front: Vec<(usize, bool, usize)> = front0.iter().map(|&(e, f)| (e, f, s)).collect(); + let mut stages: Vec> = Vec::new(); + loop { + if stages.len() > brep.faces.len() { + return Err("the rings do not end".into()); + } + let across: Vec = front.iter().map(|&(e, _, below)| brep.edges.get(e).map_or(usize::MAX, |x| x.other(below))).collect(); + let first = across.first().copied().ok_or("an empty loop")?; + if across.iter().all(|&f| f == first) && brep.faces.get(first).is_some_and(|x| x.plane.is_some()) && is_loop_of(brep, first, &front) { + match target { + None => target = Some(first), + Some(t) if t == first => {} + Some(_) => return Err("the loops end on different faces (not one-to-one)".into()), + } + break; + } + // Rotate the front to start where the face across changes, then take runs. + let n = front.len(); + let start = (0..n).find(|&k| across.get((k + n - 1) % n) != across.get(k)).unwrap_or(0); + let mut ring: Vec = Vec::new(); + let mut k = 0; + while k < n { + let f = across.get((start + k) % n).copied().ok_or("front index")?; + if used.contains(&f) || ring.iter().any(|l| l.face == f) { + return Err(format!("face {} is reached twice (not a one-to-one sweep)", brep.face_tag(f))); + } + let mut run = Vec::new(); + while k < n && across.get((start + k) % n) == Some(&f) { + let &(e, fwd, _) = front.get((start + k) % n).ok_or("front index")?; + run.push((e, fwd)); + k += 1; + } + ring.push(link(brep, f, &run)?); + } + used.extend(ring.iter().map(|l| l.face)); + front = ring.iter().flat_map(|l| l.top.iter().map(move |&(e, fwd)| (e, fwd, l.face))).collect(); + stages.push(ring); + } + if stages.is_empty() { + return Err("the source touches the target".into()); + } + rings.push(stages); + } + let target = target.ok_or("no target face")?; + used.insert(target); + if used.len() != brep.faces.len() { + return Err(format!("{} faces are not on the sweep path", brep.faces.len() - used.len())); + } + if brep.faces.get(target).map(|f| f.loops.len()) != Some(source.loops.len()) { + return Err("the source and target have different numbers of loops".into()); + } + Ok(Topology { source: s, target, rings }) +} + +/// Whether the front's edges are exactly one loop of face `t`. +fn is_loop_of(brep: &Brep, t: usize, front: &[(usize, bool, usize)]) -> bool { + let edges: BTreeSet = front.iter().map(|x| x.0).collect(); + brep.faces.get(t).is_some_and(|f| f.loops.iter().any(|l| l.iter().map(|x| x.0).collect::>() == edges)) +} + +/// Split linking face `f`'s loop into bottom (the front `run`), sides and top. +fn link(brep: &Brep, f: usize, run: &[(usize, bool)]) -> std::result::Result { + let face = brep.faces.get(f).ok_or("face index")?; + let tag = face.id; + let [lp] = face.loops.as_slice() else { + return Err(format!("linking face {tag} has {} boundary loops (periodic or holed linking faces are not swept)", face.loops.len())); + }; + let n = lp.len(); + let bottom_edges: BTreeSet = run.iter().map(|x| x.0).collect(); + let in_bottom = |k: usize| lp.get(k % n).is_some_and(|x| bottom_edges.contains(&x.0)); + let start = (0..n) + .find(|&k| in_bottom(k) && !in_bottom(k + n - 1)) + .ok_or_else(|| format!("linking face {tag}: the shared edges do not form one chain"))?; + let len = (0..n).take_while(|&k| in_bottom(start + k)).count(); + if len != run.len() || (0..n).filter(|&k| in_bottom(k)).count() != len { + return Err(format!("linking face {tag}: the shared edges do not form one chain")); + } + if n < len + 3 { + return Err(format!("linking face {tag} has no top side")); + } + let at = |k: usize| lp.get(k % n).copied().ok_or_else(|| "loop index".to_string()); + let bottom: Vec<(usize, bool)> = (0..len).map(|k| at(start + k)).collect::>()?; + let right = at(start + len)?; + let left = at(start + n - 1)?; + let top: Vec<(usize, bool)> = (len + 1..n - 1).map(|k| at(start + k)).collect::>()?; + Ok(Link { face: f, bottom, right, top, left }) +} + +/// Node positions while building. +struct Nodes { + xyz: Vec

, +} + +impl Nodes { + fn add(&mut self, p: P) -> u32 { + self.xyz.push(p); + self.xyz.len() as u32 - 1 + } + fn get(&self, i: u32) -> P { + self.xyz.get(i as usize).copied().unwrap_or([0.0; 3]) + } +} + +/// Node ids along oriented edge uses, end to end (each use's last node is the next one's first). +fn chain(edge_nodes: &[Vec], uses: &[(usize, bool)]) -> Vec { + let mut out: Vec = Vec::new(); + for &(e, fwd) in uses { + let mut ids = edge_nodes.get(e).cloned().unwrap_or_default(); + if !fwd { + ids.reverse(); + } + if !out.is_empty() && !ids.is_empty() { + ids.remove(0); + } + out.extend(ids); + } + out +} + +fn build(ctx: &Ctx, topo: &Topology, refine: u32) -> Result { + let brep = ctx.brep; + let source = brep.faces.get(topo.source).ok_or_else(|| Error::Failed("source face".into()))?; + let (s_point, s_normal) = source.plane.ok_or_else(|| Error::Failed("the source is not planar".into()))?; + // Interval assignment over the B-rep edges. + let mut ia = Intervals::default(); + let sides: BTreeSet = topo.rings.iter().flatten().flatten().flat_map(|l| [l.right.0, l.left.0]).collect(); + let fine = f64::from(1u32 << refine.min(8)); + for e in 0..brep.edges.len() { + let goal = ideal_divisions(ctx.field, &brep.edge_points(e, true))?; + ia.var(if sides.contains(&e) { goal } else { goal * fine }, 1); + } + let vars = |uses: &[(usize, bool)]| uses.iter().map(|x| x.0).collect::>(); + for stages in &topo.rings { + for ring in stages { + for l in ring { + ia.equal.push((vars(&l.bottom), vars(&l.top))); + ia.equal.push((vec![l.right.0], vec![l.left.0])); + } + } + } + // The same layers on every ring: stage by stage when the rings have as many stages. + let side_of = |stages: &Vec>, k: usize| stages.get(k).and_then(|r| r.first()).map(|l| l.right.0); + if let Some(first) = topo.rings.first() { + for other in topo.rings.iter().skip(1) { + if other.len() == first.len() { + for k in 0..first.len() { + if let (Some(a), Some(b)) = (side_of(first, k), side_of(other, k)) { + ia.equal.push((vec![a], vec![b])); + } + } + } else { + let all = |st: &Vec>| (0..st.len()).filter_map(|k| side_of(st, k)).collect::>(); + ia.equal.push((all(first), all(other))); + } + } + } + // The source's frame: u × v along the sweep (into the solid). + let frame = Frame::translate(s_point, scale(s_normal, -1.0)).ok_or_else(|| Error::Failed("source normal".into()))?; + let q = |p: P| frame.coords(p).0; + let polylines: Vec>> = + source.loops.iter().map(|lp| lp.iter().map(|&(e, f)| brep.edge_points(e, f).into_iter().map(q).collect()).collect()).collect(); + let scheme = region::classify(&polylines); + let loop_vars: Vec> = source.loops.iter().map(|lp| vars(lp)).collect(); + region::constrain(&mut ia, &loop_vars, scheme); + let counts = ia.solve()?; + // Edge nodes: B-rep vertices once, interiors on both faces of the edge. + let mut nodes = Nodes { xyz: Vec::new() }; + let mut vertex_node: BTreeMap = BTreeMap::new(); + let mut edge_nodes: Vec> = Vec::with_capacity(brep.edges.len()); + for (e, edge) in brep.edges.iter().enumerate() { + let n = counts.get(e).copied().unwrap_or(1); + let pts = divide_polyline(ctx.field, &brep.edge_points(e, true), n)?; + let last = pts.len().saturating_sub(1); + let mut ids = Vec::with_capacity(pts.len()); + for (k, p) in pts.iter().enumerate() { + let end = if k == 0 { + Some(edge.first()) + } else if k == last { + Some(edge.last()) + } else { + None + }; + ids.push(match end { + Some(v) => *vertex_node.entry(v).or_insert_with(|| nodes.add(brep.point(v))), + None => nodes.add(project_to_faces(ctx, &edge.faces, *p)), + }); + } + edge_nodes.push(ids); + } + // Linking-face grids and the ring of every layer: rings_by_layer[loop][layer]. + let mut rings_by_layer: Vec>> = Vec::with_capacity(topo.rings.len()); + for stages in &topo.rings { + let mut layers: Vec> = Vec::new(); + for (k, ring) in stages.iter().enumerate() { + let n_side = ring.first().map_or(0, |l| chain(&edge_nodes, &[l.right]).len().saturating_sub(1)); + let mut rows: Vec> = vec![Vec::new(); n_side + 1]; + for l in ring { + let grid = face_grid(ctx, &mut nodes, &edge_nodes, l)?; + if grid.len() != n_side + 1 { + return Err(Error::Failed(format!("linking face {} has sides of different lengths", brep.face_tag(l.face)))); + } + // Rows run along the face's bottom traversal, which is the front's direction. + for (row, slot) in grid.iter().zip(rows.iter_mut()) { + let mut r = row.clone(); + r.pop(); + slot.extend(r); + } + } + // A later stage's rings may start elsewhere around the loop: turn them to continue the + // previous stage's columns. + if let (Some(prev), Some(first)) = (layers.last(), rows.first()) { + let n = first.len(); + let turn = (0..n) + .find(|&o| (0..n).all(|i| first.get((i + o) % n) == prev.get(i))) + .ok_or_else(|| Error::Failed(format!("stage {k} of the sweep does not continue the previous stage's ring")))?; + for r in rows.iter_mut() { + r.rotate_left(turn); + } + } + layers.extend(rows.into_iter().enumerate().filter(|(j, _)| k == 0 || *j > 0).map(|(_, r)| r)); + } + rings_by_layer.push(layers); + } + let n_layers = rings_by_layer.first().map_or(0, Vec::len).saturating_sub(1); + if n_layers == 0 || rings_by_layer.iter().any(|l| l.len() != n_layers + 1) { + return Err(Error::Failed("the rings have different numbers of layers".into())); + } + // Source mesh on the layer-0 rings, which run as the source's loops do. + let ring0: Vec> = rings_by_layer.iter().map(|l| l.first().cloned().unwrap_or_default()).collect(); + let loops2: Vec> = ring0.iter().map(|r| r.iter().map(|&i| q(nodes.get(i))).collect()).collect(); + let plane_field = PlaneField { field: ctx.field, frame, h: 0.0 }; + // Where each loop's edge uses start in its ring: at the edges' first vertices. + let starts: Vec> = source + .loops + .iter() + .zip(&ring0) + .map(|(lp, ring)| { + lp.iter() + .filter_map(|&(e, fwd)| { + brep.edge_vertices(e, fwd).first().and_then(|v| vertex_node.get(v)).and_then(|id| ring.iter().position(|x| x == id)) + }) + .collect() + }) + .collect(); + // The rings may start mid-loop: give the mapped scheme its corners as ring positions. + let (scheme, starts) = match scheme { + region::Scheme::Mapped(c) => { + let s0 = starts.first().cloned().unwrap_or_default(); + let mut pos: Vec = c.iter().filter_map(|&k| s0.get(k).copied()).collect(); + pos.sort_unstable(); + if pos.len() == 4 { (region::Scheme::Mapped([0, 1, 2, 3]), vec![pos]) } else { (region::Scheme::Free, starts) } + } + other => (other, starts), + }; + let (mesh, used_scheme) = region::mesh(&loops2, &starts, scheme, &plane_field, &ctx.opts.quad)?; + let boundary_ids: Vec = ring0.concat(); + let n_boundary = boundary_ids.len(); + if mesh.boundary_nodes != n_boundary { + return Err(Error::Failed("the source mesh changed its boundary".into())); + } + let interior2: Vec<[f64; 2]> = mesh.nodes.iter().skip(n_boundary).copied().collect(); + let count = mesh.quads.len().saturating_mul(n_layers); + if count > ctx.opts.max_elements { + return Err(Error::Limit(format!("the mesh would have {count} hexahedra; the limit is {}", ctx.opts.max_elements))); + } + let work = (interior2.len() as u64).saturating_mul(n_boundary as u64).saturating_mul(n_layers as u64); + if work > MAX_CORRECTION_WORK { + return Err(Error::Limit(format!("placing the sweep's interior nodes would take {work} steps; the limit is {MAX_CORRECTION_WORK}"))); + } + let s_ring: Vec<[f64; 2]> = boundary_ids.iter().map(|&i| q(nodes.get(i))).collect(); + // Interior nodes, layer by layer. The residual weights depend only on the source plane. + let weights: Vec> = interior2.iter().map(|p| shepard_weights(*p, &s_ring)).collect(); + let mut layer_nodes: Vec> = Vec::with_capacity(n_layers + 1); + for j in 0..=n_layers { + let ring: Vec = rings_by_layer.iter().flat_map(|l| l.get(j).cloned().unwrap_or_default()).collect(); + if ring.len() != n_boundary { + return Err(Error::Failed(format!("layer {j}'s ring has {} nodes, the source's {n_boundary}", ring.len()))); + } + let mut ids = ring.clone(); + if j == 0 { + ids.extend(interior2.iter().map(|p| nodes.add(frame.point(*p, 0.0)))); + } else { + let x: Vec

= ring.iter().map(|&i| nodes.get(i)).collect(); + let map = Affine::fit(&s_ring, &x)?; + let residual: Vec

= s_ring.iter().zip(&x).map(|(s, xi)| sub(*xi, map.at(*s))).collect(); + for (p, w) in interior2.iter().zip(&weights) { + let correction = w.iter().fold([0.0; 3], |acc, &(k, wk)| add(acc, scale(residual.get(k).copied().unwrap_or([0.0; 3]), wk))); + let mut y = add(map.at(*p), correction); + if j == n_layers { + y = project_to_face(ctx, topo.target, y).unwrap_or(y); + } + ids.push(nodes.add(y)); + } + } + layer_nodes.push(ids); + } + let mut hexes = Vec::with_capacity(count); + for quad in &mesh.quads { + for w in layer_nodes.windows(2) { + let (lo, hi) = (&w[0], &w[1]); + let b = quad.map(|v| lo.get(v as usize).copied().unwrap_or(0)); + let t = quad.map(|v| hi.get(v as usize).copied().unwrap_or(0)); + hexes.push([b[0], b[1], b[2], b[3], t[0], t[1], t[2], t[3]]); + } + } + let stages: usize = topo.rings.first().map_or(0, Vec::len); + let description = format!( + "swept from face {} to face {} through {stages} stage{} of linking faces: {n_layers} layers of {} {} source quadrilaterals", + brep.face_tag(topo.source), + brep.face_tag(topo.target), + if stages == 1 { "" } else { "s" }, + mesh.quads.len(), + match used_scheme { + region::Scheme::Mapped(_) => "mapped", + region::Scheme::Annulus => "annulus", + region::Scheme::Disk => "disk", + region::Scheme::Free => "free", + } + ); + Ok(Raw { nodes: nodes.xyz, hexes, description }) +} + +/// The transfinite grid of a linking face projected onto it: `grid[j][i]`, `i` along the bottom +/// as the front runs, `j` from bottom to top. +fn face_grid(ctx: &Ctx, nodes: &mut Nodes, edge_nodes: &[Vec], l: &Link) -> Result>> { + // The face traverses bottom → right → top → left; the front runs against its bottom. + let mut bottom = chain(edge_nodes, &l.bottom); + bottom.reverse(); + let mut right = chain(edge_nodes, &[l.left]); + right.reverse(); + let top = chain(edge_nodes, &l.top); + let left = chain(edge_nodes, &[l.right]); + let pts = |ids: &[u32]| ids.iter().map(|&i| nodes.get(i)).collect::>(); + let grid = transfinite_grid(&pts(&bottom), &pts(&right), &pts(&top), &pts(&left)) + .map_err(|e| Error::Failed(format!("linking face {}: {e}", ctx.brep.face_tag(l.face))))?; + let (nu, nv) = (bottom.len().saturating_sub(1), right.len().saturating_sub(1)); + let mut out = Vec::with_capacity(nv + 1); + for j in 0..=nv { + let mut row = Vec::with_capacity(nu + 1); + for i in 0..=nu { + let id = if j == 0 { + bottom.get(i).copied() + } else if j == nv { + top.get(i).copied() + } else if i == 0 { + left.get(j).copied() + } else if i == nu { + right.get(j).copied() + } else { + let p = grid.get(j * (nu + 1) + i).copied().unwrap_or([0.0; 3]); + Some(nodes.add(project_to_face(ctx, l.face, p).unwrap_or(p))) + }; + row.push(id.ok_or_else(|| Error::Failed("grid side".into()))?); + } + out.push(row); + } + Ok(out) +} + +/// `x ≈ c + M s`: the least-squares affine map from the source plane into space. +struct Affine { + c: P, + m: [P; 2], +} + +impl Affine { + fn fit(s: &[[f64; 2]], x: &[P]) -> Result { + if s.len() != x.len() || s.len() < 3 { + return Err(Error::Failed("a sweep layer needs at least three ring nodes".into())); + } + // Normal equations of [1, s_u, s_v], centred for conditioning. + let n = s.len() as f64; + let (mu, mv) = (s.iter().map(|p| p[0]).sum::() / n, s.iter().map(|p| p[1]).sum::() / n); + let (mut suu, mut suv, mut svv) = (0.0, 0.0, 0.0); + let (mut xu, mut xv, mut xm) = ([0.0; 3], [0.0; 3], [0.0; 3]); + for (p, xi) in s.iter().zip(x) { + let (u, v) = (p[0] - mu, p[1] - mv); + suu += u * u; + suv += u * v; + svv += v * v; + xu = add(xu, scale(*xi, u)); + xv = add(xv, scale(*xi, v)); + xm = add(xm, scale(*xi, 1.0 / n)); + } + let det = suu * svv - suv * suv; + if !(det > 1e-12 * suu * svv) { + return Err(Error::Failed("the source ring is collinear; the sweep's layers cannot be mapped".into())); + } + let mu_vec = scale(sub(scale(xu, svv), scale(xv, suv)), 1.0 / det); + let mv_vec = scale(sub(scale(xv, suu), scale(xu, suv)), 1.0 / det); + let c = sub(xm, add(scale(mu_vec, mu), scale(mv_vec, mv))); + Ok(Affine { c, m: [mu_vec, mv_vec] }) + } + + fn at(&self, s: [f64; 2]) -> P { + add(self.c, add(scale(self.m[0], s[0]), scale(self.m[1], s[1]))) + } +} + +/// Inverse-distance (Shepard, power 4) weights of the ring nodes at `p`, normalized: the +/// residual correction at `p` is their weighted sum of the ring's residuals. +fn shepard_weights(p: [f64; 2], ring: &[[f64; 2]]) -> Vec<(usize, f64)> { + let mut w: Vec<(usize, f64)> = Vec::with_capacity(ring.len()); + for (k, s) in ring.iter().enumerate() { + let d2 = (p[0] - s[0]).powi(2) + (p[1] - s[1]).powi(2); + if d2 == 0.0 { + return vec![(k, 1.0)]; + } + w.push((k, 1.0 / (d2 * d2))); + } + let sum: f64 = w.iter().map(|x| x.1).sum(); + if sum > 0.0 { + for x in &mut w { + x.1 /= sum; + } + } + w +} diff --git a/crates/mesh/src/hex/tests.rs b/crates/mesh/src/hex/tests.rs new file mode 100644 index 0000000..c1c00f2 --- /dev/null +++ b/crates/mesh/src/hex/tests.rs @@ -0,0 +1,294 @@ +use super::*; +use crate::size::UniformSize; +use crate::surface::{BoundarySurface, SurfaceOptions}; +use crate::vec3::{cross, dot, sub}; + +/// Ear clipping of a simple planar polygon (counter-clockwise about `normal`). +fn ears(points: &[[f64; 3]], poly: &[u32]) -> Vec<[u32; 3]> { + let p = |i: u32| points[i as usize]; + let normal = (0..poly.len()).fold([0.0; 3], |n, k| { + let (a, b) = (p(poly[k]), p(poly[(k + 1) % poly.len()])); + [n[0] + (a[1] - b[1]) * (a[2] + b[2]), n[1] + (a[2] - b[2]) * (a[0] + b[0]), n[2] + (a[0] - b[0]) * (a[1] + b[1])] + }); + let turn = |a: [f64; 3], b: [f64; 3], c: [f64; 3]| dot(cross(sub(b, a), sub(c, a)), normal); + let mut v: Vec = poly.to_vec(); + let mut out = Vec::new(); + while v.len() > 3 { + let n = v.len(); + let k = (0..n) + .find(|&k| { + let (a, b, c) = (p(v[(k + n - 1) % n]), p(v[k]), p(v[(k + 1) % n])); + turn(a, b, c) > 0.0 + && v.iter().all(|&o| { + let q = p(o); + o == v[(k + n - 1) % n] + || o == v[k] + || o == v[(k + 1) % n] + || !(turn(a, b, q) >= 0.0 && turn(b, c, q) >= 0.0 && turn(c, a, q) >= 0.0) + }) + }) + .unwrap(); + out.push([v[(k + n - 1) % n], v[k], v[(k + 1) % n]]); + v.remove(k); + } + out.push([v[0], v[1], v[2]]); + out +} + +/// A closed surface from planar polygonal faces (counter-clockwise from outside), tagged by face. +pub(crate) fn polyhedron(points: &[[f64; 3]], faces: &[Vec]) -> BoundarySurface { + let mut tris = Vec::new(); + let mut tags = Vec::new(); + for (f, poly) in faces.iter().enumerate() { + for t in ears(points, poly) { + tris.push(t); + tags.push(f as u32); + } + } + BoundarySurface::new(points, &tris, &tags, &SurfaceOptions::default()).unwrap() +} + +/// A prism over a counter-clockwise polygon from z0 to z1. +pub(crate) fn prism(poly: &[[f64; 2]], z0: f64, z1: f64) -> BoundarySurface { + let n = poly.len() as u32; + let mut pts: Vec<[f64; 3]> = poly.iter().map(|p| [p[0], p[1], z0]).collect(); + pts.extend(poly.iter().map(|p| [p[0], p[1], z1])); + let mut faces = vec![(0..n).rev().collect::>(), (n..2 * n).collect()]; + for k in 0..n { + let j = (k + 1) % n; + faces.push(vec![k, j, n + j, n + k]); + } + polyhedron(&pts, &faces) +} + +fn check(m: &HexMesh, volume: f64) { + assert!(m.report.gate.passed()); + assert!((m.report.volume - volume).abs() < 1e-9 * volume, "{} vs {volume}", m.report.volume); + assert!(m.report.min_scaled_jacobian > 0.0); + // Every boundary face is tagged and the boundary faces are exactly the hex faces used once. + assert!(!m.boundary.is_empty()); +} + +#[test] +fn box_extrudes_into_a_grid() { + let s = prism(&[[0.0, 0.0], [4.0, 0.0], [4.0, 2.0], [0.0, 2.0]], 0.0, 3.0); + let m = hex_mesh(&s, &[], &UniformSize(1.0), &HexOptions::default()).unwrap(); + check(&m, 24.0); + assert_eq!(m.element_count(), 24, "{}", m.report.description); + assert!((m.report.min_scaled_jacobian - 1.0).abs() < 1e-12); + // Boundary: 2·(4·2 + 4·3 + 2·3) = 52 quads. + assert_eq!(m.boundary.len(), 52); +} + +#[test] +fn l_shape_extrudes_with_a_free_or_mapped_section() { + let l = [[0.0, 0.0], [6.0, 0.0], [6.0, 2.0], [2.0, 2.0], [2.0, 6.0], [0.0, 6.0]]; + let s = prism(&l, 0.0, 2.0); + let m = hex_mesh(&s, &[], &UniformSize(1.0), &HexOptions::default()).unwrap(); + check(&m, 40.0); +} + +#[test] +fn quadratic_box() { + let s = prism(&[[0.0, 0.0], [2.0, 0.0], [2.0, 2.0], [0.0, 2.0]], 0.0, 2.0); + let m = hex_mesh(&s, &[], &UniformSize(1.0), &HexOptions { quadratic: true, ..HexOptions::default() }).unwrap(); + check(&m, 8.0); + assert_eq!(m.kind, ElementKind::Hex20); + // 3³ corner nodes + 3·(2·3·3) = 54 edge nodes. + assert_eq!(m.nodes.len(), 27 + 54); +} + +/// A block whose top is the curved surface z = 2 + 0.3 sin(πx/4), tessellated finely. +fn curved_block() -> BoundarySurface { + let n = 16usize; + let top = |x: f64| 2.0 + 0.3 * (std::f64::consts::PI * x / 4.0).sin(); + let mut pts: Vec<[f64; 3]> = vec![[0.0, 0.0, 0.0], [4.0, 0.0, 0.0], [4.0, 2.0, 0.0], [0.0, 2.0, 0.0]]; + // Top grid rows at y = 0 and y = 2: ids 4..4+n+1 and 5+n.. + for y in [0.0, 2.0] { + for k in 0..=n { + let x = 4.0 * k as f64 / n as f64; + pts.push([x, y, top(x)]); + } + } + let t0 = |k: usize| (4 + k) as u32; + let t2 = |k: usize| (5 + n + k) as u32; + let mut faces: Vec> = vec![vec![0, 3, 2, 1]]; + // y = 0 side: 0, 1, then the top edge from x = 4 back to 0. + faces.push([0u32, 1].into_iter().chain((0..=n).rev().map(t0)).collect()); + faces.push([2u32, 3].into_iter().chain((0..=n).map(t2)).collect()); + faces.push(vec![1, 2, t2(n), t0(n)]); + faces.push(vec![3, 0, t0(0), t2(0)]); + // The top: strips of two triangles, one face. + let mut tris: Vec<[u32; 3]> = Vec::new(); + let mut tags: Vec = Vec::new(); + for (f, poly) in faces.iter().enumerate() { + for t in ears(&pts, poly) { + tris.push(t); + tags.push(f as u32); + } + } + for k in 0..n { + tris.push([t0(k), t0(k + 1), t2(k + 1)]); + tris.push([t0(k), t2(k + 1), t2(k)]); + tags.extend([5, 5]); + } + BoundarySurface::new(&pts, &tris, &tags, &SurfaceOptions::default()).unwrap() +} + +#[test] +fn mapped_block_with_a_curved_top() { + let s = curved_block(); + let opts = HexOptions { method: HexMethod::Mapped, ..HexOptions::default() }; + let m = hex_mesh(&s, &[], &UniformSize(0.5), &opts).unwrap(); + assert!(m.report.gate.passed()); + assert!(m.report.min_scaled_jacobian > 0.8, "{}", m.report.min_scaled_jacobian); + assert!(m.report.description.contains("mapped block of 8 × 4 × 4"), "{}", m.report.description); + // Top nodes lie on the tessellated top. + let quadratic = hex_mesh(&s, &[], &UniformSize(0.5), &HexOptions { quadratic: true, ..opts }).unwrap(); + assert!(quadratic.report.gate.passed()); + assert!((quadratic.report.volume - s.volume()).abs() < 1e-3 * s.volume(), "{} vs {}", quadratic.report.volume, s.volume()); +} + +/// A thick-walled pipe along z (radii `ri` < `ro`, height `h`, `n` segments around, `n` even): +/// eight faces, each cylinder and annulus in two halves, with their analytic surfaces. +fn pipe(ri: f64, ro: f64, h: f64, n: usize) -> (BoundarySurface, Vec>) { + let ring = |r: f64, z: f64| -> Vec<[f64; 3]> { + (0..n) + .map(|k| { + let t = std::f64::consts::TAU * k as f64 / n as f64; + [r * t.cos(), r * t.sin(), z] + }) + .collect() + }; + let mut pts = Vec::new(); + for (r, z) in [(ro, 0.0), (ri, 0.0), (ro, h), (ri, h)] { + pts.extend(ring(r, z)); + } + let id = |ring: usize, k: usize| (ring * n + k % n) as u32; + let (bo, bi, to, ti) = (0, 1, 2, 3); + let mut tris: Vec<[u32; 3]> = Vec::new(); + let mut tags: Vec = Vec::new(); + // A quad as two triangles, turned to face `outward`. + let mut quad = |q: [u32; 4], outward: [f64; 3], tag: u32| { + let p = |i: u32| pts[i as usize]; + let n0 = cross(sub(p(q[1]), p(q[0])), sub(p(q[2]), p(q[0]))); + let q = if dot(n0, outward) < 0.0 { [q[0], q[3], q[2], q[1]] } else { q }; + tris.push([q[0], q[1], q[2]]); + tris.push([q[0], q[2], q[3]]); + tags.extend([tag, tag]); + }; + for k in 0..n { + let half = u32::from(k >= n / 2); + let t = std::f64::consts::TAU * (k as f64 + 0.5) / n as f64; + let radial = [t.cos(), t.sin(), 0.0]; + quad([id(bo, k), id(bo, k + 1), id(to, k + 1), id(to, k)], radial, half); + quad([id(bi, k), id(bi, k + 1), id(ti, k + 1), id(ti, k)], [-radial[0], -radial[1], 0.0], 2 + half); + quad([id(bo, k), id(bo, k + 1), id(bi, k + 1), id(bi, k)], [0.0, 0.0, -1.0], 4 + half); + quad([id(to, k), id(to, k + 1), id(ti, k + 1), id(ti, k)], [0.0, 0.0, 1.0], 6 + half); + } + let surface = BoundarySurface::new(&pts, &tris, &tags, &SurfaceOptions::default()).unwrap(); + let cyl = |r: f64| Some(Surface::Cylinder { origin: [0.0; 3], axis: [0.0, 0.0, 1.0], radius: r }); + let plane = |z: f64, nz: f64| Some(Surface::Plane { point: [0.0, 0.0, z], normal: [0.0, 0.0, nz] }); + (surface, vec![cyl(ro), cyl(ro), cyl(ri), cyl(ri), plane(0.0, -1.0), plane(0.0, -1.0), plane(h, 1.0), plane(h, 1.0)]) +} + +/// Every node of the boundary faces tagged with `faces` lies at `radius` from the z axis. +fn on_radius(m: &HexMesh, faces: &[u32], radius: f64) -> bool { + m.boundary.iter().filter(|b| faces.contains(&b.brep_face)).all(|b| { + let local = m.kind.faces()[b.face as usize - 1]; + let e = m.element(b.element as usize).unwrap(); + local.iter().all(|&i| { + let p = m.nodes[e[i] as usize]; + ((p[0] * p[0] + p[1] * p[1]).sqrt() - radius).abs() < 1e-9 * radius + }) + }) +} + +#[test] +fn pipe_extrudes_and_revolves_with_nodes_on_the_cylinders() { + let (s, faces) = pipe(2.0, 3.0, 4.0, 32); + for method in [HexMethod::Auto, HexMethod::Revolve] { + for quadratic in [false, true] { + let opts = + HexOptions { method, quadratic, chordal_tolerance: 3.0 * (1.0 - (std::f64::consts::PI / 32.0).cos()), ..HexOptions::default() }; + let m = hex_mesh(&s, &faces, &UniformSize(0.5), &opts).unwrap(); + let expected = if method == HexMethod::Auto { HexMethod::Extrude } else { HexMethod::Revolve }; + assert_eq!(m.report.method, expected, "{}", m.report.description); + assert!(m.report.gate.passed()); + assert!(m.report.min_scaled_jacobian > 0.5, "{method:?} {quadratic}: {} {}", m.report.min_scaled_jacobian, m.report.description); + assert!(on_radius(&m, &[0, 1], 3.0) && on_radius(&m, &[2, 3], 2.0), "{method:?} {quadratic}"); + let exact = std::f64::consts::PI * (9.0 - 4.0) * 4.0; + let tolerance = if quadratic { 1e-3 } else { 2e-2 }; + assert!((m.report.volume - exact).abs() < tolerance * exact, "{method:?} {quadratic}: {} vs {exact}", m.report.volume); + // Every boundary face is tagged with one of the pipe's faces. + assert!(m.boundary.iter().all(|b| b.brep_face < 8)); + } + } +} + +#[test] +fn box_faces_are_tagged_and_meshes_are_deterministic() { + let s = prism(&[[0.0, 0.0], [4.0, 0.0], [4.0, 2.0], [0.0, 2.0]], 0.0, 3.0); + let a = hex_mesh(&s, &[], &UniformSize(1.0), &HexOptions::default()).unwrap(); + let b = hex_mesh(&s, &[], &UniformSize(1.0), &HexOptions::default()).unwrap(); + assert_eq!(a, b); + // prism(): face 0 bottom, 1 top, then the sides along the polygon's edges. + let count = |f: u32| a.boundary.iter().filter(|x| x.brep_face == f).count(); + assert_eq!([count(0), count(1), count(2), count(3), count(4), count(5)], [8, 8, 12, 6, 12, 6]); + // Each boundary face lies in its B-rep face's plane. + let plane = |f: u32| -> (usize, f64) { + match f { + 0 => (2, 0.0), + 1 => (2, 3.0), + 2 => (1, 0.0), + 3 => (0, 4.0), + 4 => (1, 2.0), + _ => (0, 0.0), + } + }; + for b in &a.boundary { + let (axis, value) = plane(b.brep_face); + let e = a.element(b.element as usize).unwrap(); + for &i in ElementKind::Hex8.faces()[b.face as usize - 1] { + assert!((a.nodes[e[i] as usize][axis] - value).abs() < 1e-12, "{b:?}"); + } + } +} + +#[test] +fn an_l_prism_sweeps_one_to_one() { + let l = [[0.0, 0.0], [6.0, 0.0], [6.0, 2.0], [2.0, 2.0], [2.0, 6.0], [0.0, 6.0]]; + let s = prism(&l, 0.0, 2.0); + let m = hex_mesh(&s, &[], &UniformSize(1.0), &HexOptions { method: HexMethod::Sweep, ..HexOptions::default() }).unwrap(); + assert_eq!(m.report.method, HexMethod::Sweep); + check(&m, 40.0); +} + +#[test] +fn unsweepable_bodies_say_why_and_hostile_input_is_an_error() { + // A tetrahedron: every face is slanted to every other. + let tet = polyhedron( + &[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]], + &[vec![0, 2, 1], vec![0, 1, 3], vec![1, 2, 3], vec![0, 3, 2]], + ); + match hex_mesh(&tet, &[], &UniformSize(0.2), &HexOptions::default()) { + Err(Error::Unsupported(why)) => { + for m in ["extrude", "sweep", "revolve", "mapped"] { + assert!(why.contains(m), "{why}"); + } + } + other => panic!("{other:?}"), + } + let s = prism(&[[0.0, 0.0], [4.0, 0.0], [4.0, 2.0], [0.0, 2.0]], 0.0, 3.0); + let bad_sizes: [&dyn crate::size::SizeField; 3] = [&UniformSize(f64::NAN), &UniformSize(0.0), &UniformSize(-1.0)]; + for f in bad_sizes { + assert!(hex_mesh(&s, &[], f, &HexOptions::default()).is_err()); + } + assert!(matches!(hex_mesh(&s, &[], &UniformSize(1e-4), &HexOptions::default()), Err(Error::Limit(_)))); + assert!(matches!(hex_mesh(&s, &[], &UniformSize(0.1), &HexOptions { max_elements: 100, ..HexOptions::default() }), Err(Error::Limit(_)))); + assert!(hex_mesh(&s, &[], &UniformSize(1.0), &HexOptions { chordal_tolerance: f64::NAN, ..HexOptions::default() }).is_err()); + // Analytic surfaces with garbage numbers are treated as free-form, not trusted. + let garbage = vec![Some(Surface::Cylinder { origin: [f64::NAN; 3], axis: [0.0; 3], radius: -1.0 }); 6]; + let m = hex_mesh(&s, &garbage, &UniformSize(1.0), &HexOptions::default()).unwrap(); + assert_eq!(m.element_count(), 24); +} diff --git a/crates/mesh/src/lib.rs b/crates/mesh/src/lib.rs index 54ffb2f..a6eb594 100644 --- a/crates/mesh/src/lib.rs +++ b/crates/mesh/src/lib.rs @@ -10,6 +10,7 @@ pub mod block; pub mod cdt; pub mod delaunay; +pub mod hex; pub mod join; pub mod planar; pub mod quad; @@ -26,6 +27,10 @@ pub enum Error { Limit(String), #[error("meshing failed: {0}")] Failed(String), + /// The input is valid but outside what the mesher handles yet (a body that is not + /// sweepable, …); the message says why. + #[error("not supported yet: {0}")] + Unsupported(String), /// The boundary surface has defects ([`surface::BoundarySurface::new`]). #[error("invalid boundary surface: {0}")] Surface(surface::Defects), diff --git a/crates/mesh/src/quality/mod.rs b/crates/mesh/src/quality/mod.rs index efc99e3..6a7ba98 100644 --- a/crates/mesh/src/quality/mod.rs +++ b/crates/mesh/src/quality/mod.rs @@ -26,7 +26,9 @@ pub(crate) mod shape; use modelcraft_model::{ElementKind, check_coordinates}; pub use jacobian::{JacobianCheck, Validity}; -pub use report::{GateOptions, GateReport, MetricStats, QualityReport, enclosed_volume, quality_report, validity_gate, validity_gate_elements}; +pub use report::{ + GateOptions, GateReport, MetricStats, QualityReport, enclosed_volume, quality_report, validity_gate, validity_gate_checks, validity_gate_elements, +}; pub(crate) use shape::{planar_corner, quad_axes}; use crate::predicates::{Point2, Point3, Sign, orient2d, snap_to_domain, tet_volume_sign}; diff --git a/crates/mesh/src/quality/report.rs b/crates/mesh/src/quality/report.rs index 3835aa0..e02b2ce 100644 --- a/crates/mesh/src/quality/report.rs +++ b/crates/mesh/src/quality/report.rs @@ -4,7 +4,7 @@ use std::collections::{BTreeMap, BTreeSet}; use modelcraft_model::{Element, ElementKind, Model, check_coordinates}; -use super::{ElementQuality, Metric, Validity, element_quality, element_validity, is_measured}; +use super::{ElementQuality, JacobianCheck, Metric, Validity, element_quality, element_validity, is_measured}; use crate::{Error, Result}; /// The most histogram bins a report may ask for. @@ -218,7 +218,11 @@ impl Gate { if kind.dimension() != 3 { return Ok(()); } - let check = element_validity(kind, xyz)?; + self.add_check(label, &element_validity(kind, xyz)?); + Ok(()) + } + + fn add_check(&mut self, label: u32, check: &JacobianCheck) { let r = &mut self.0; r.checked += 1; r.volume += check.volume; @@ -227,10 +231,23 @@ impl Gate { Validity::Invalid => r.invalid.push(label), Validity::Undetermined => r.undetermined.push(label), } - Ok(()) } } +/// [`validity_gate_elements`] over certified checks already made, for a mesher that measured its +/// elements with [`super::element_quality`] (which certifies them too) and should not certify +/// them twice: `(label, check)` per solid element. +pub fn validity_gate_checks<'a, I>(checks: I, options: &GateOptions) -> Result +where + I: IntoIterator, +{ + let mut gate = Gate::new(options)?; + for (label, check) in checks { + gate.add_check(label, check); + } + Ok(gate.0) +} + /// The validity gate every mesher runs on its output before it reaches the model: each solid /// element's Jacobian determinant is proven positive everywhere, and the elements' total volume /// equals the volume the boundary encloses. Elements come as `(label, kind, node coordinates in diff --git a/crates/mesh/src/quality/tests.rs b/crates/mesh/src/quality/tests.rs index 34ab06a..cf34f4d 100644 --- a/crates/mesh/src/quality/tests.rs +++ b/crates/mesh/src/quality/tests.rs @@ -284,7 +284,13 @@ fn report_and_gate() { // The same gate straight from a mesher's arrays. let raw = validity_gate_elements([(7, ElementKind::Hex8, inside_out_hex()), (8, ElementKind::Hex8, UNIT_CUBE.to_vec())], &GateOptions::default()) .unwrap(); - assert_eq!((raw.checked, raw.invalid), (2, vec![7])); + assert_eq!((raw.checked, raw.invalid.clone()), (2, vec![7])); + // And from checks a mesher already made with element_quality. + let checks: Vec<(u32, JacobianCheck)> = [(7, inside_out_hex()), (8, UNIT_CUBE.to_vec())] + .into_iter() + .map(|(l, x)| (l, element_quality(ElementKind::Hex8, &x).unwrap().jacobian.unwrap())) + .collect(); + assert_eq!(validity_gate_checks(checks.iter().map(|(l, c)| (*l, c)), &GateOptions::default()).unwrap(), raw); // A unit cube split into six tetrahedra around the 1-7 diagonal passes against its boundary. let c = UNIT_CUBE; diff --git a/crates/testkit/src/runner.rs b/crates/testkit/src/runner.rs index 8c8f1f2..9e8e05c 100644 --- a/crates/testkit/src/runner.rs +++ b/crates/testkit/src/runner.rs @@ -13,8 +13,10 @@ use std::path::Path; use std::time::Instant; use crate::corpus::{self, CATEGORIES, Part}; -use modelcraft_geometry::{Design, Measure, SurfaceMesh}; +use modelcraft_geometry::{Design, FaceSurface, Measure, SurfaceMesh}; use modelcraft_mesh::cdt::{CdtOptions, tetrahedralize}; +use modelcraft_mesh::hex::{HexOptions, Surface, hex_mesh}; +use modelcraft_mesh::size::UniformSize; use modelcraft_mesh::surface::{BoundarySurface, SurfaceOptions}; pub const FILE: &str = "docs/meshing/corpus.md"; @@ -42,7 +44,7 @@ pub trait Stage { /// The stages, in column order. Meshers add theirs here. fn stages() -> Vec> { - vec![Box::new(SurfaceVolume), Box::new(BoundaryCheck), Box::new(TetCdt)] + vec![Box::new(SurfaceVolume), Box::new(BoundaryCheck), Box::new(TetCdt), Box::new(StructuredHex)] } /// The boundary surface's enclosed volume against SolidCraft's measured volume and the closed @@ -114,6 +116,81 @@ impl Stage for TetCdt { } } +/// Element size of the structured hex stage, relative to the part's bounding-box diagonal. +const HEX_RELATIVE_SIZE: f64 = 1.0 / 25.0; +/// Characters of a "not meshed" reason kept in the table. +const REASON_CHARS: usize = 240; + +/// Structured all-hex meshing (`modelcraft_mesh::hex`): the method that meshed each body, the +/// element count and the smallest scaled Jacobian of the hex8 mesh and of the hex20 mesh (mid-edge +/// nodes on the geometry), or why no structured method fits. +struct StructuredHex; + +impl Stage for StructuredHex { + fn headers(&self) -> &'static [&'static str] { + &["hex", "hex elements", "hex8 min SJ", "hex20 min SJ", "hex ms"] + } + + fn run(&self, l: &Loaded) -> Result, String> { + let start = Instant::now(); + let (lo, hi) = l.measures.iter().fold(([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]), |(mut lo, mut hi), m| { + for i in 0..3 { + lo[i] = lo[i].min(m.bbox_min[i]); + hi[i] = hi[i].max(m.bbox_max[i]); + } + (lo, hi) + }); + let diagonal = (0..3).map(|i| (hi[i] - lo[i]).powi(2)).sum::().sqrt(); + let size = UniformSize(HEX_RELATIVE_SIZE * diagonal); + let linear = HexOptions { chordal_tolerance: l.tolerance, ..HexOptions::default() }; + let quadratic = HexOptions { quadratic: true, ..linear.clone() }; + let (mut methods, mut elements, mut sj8, mut sj20) = (Vec::new(), 0usize, f64::INFINITY, Ok(f64::INFINITY)); + for (body, mesh) in l.design.bodies.iter().zip(&l.surfaces) { + let surface_options = SurfaceOptions { weld_distance: mesh.weld_distance, ..SurfaceOptions::default() }; + let s = BoundarySurface::new(&mesh.positions, &mesh.triangles, &mesh.tri_face, &surface_options) + .map_err(|e| format!("{}: {e}", body.name))?; + let faces: Vec> = + body.face_surfaces().map_err(|e| format!("{}: {e}", body.name))?.into_iter().map(|f| f.map(analytic)).collect(); + let m_method = match hex_mesh(&s, &faces, &size, &linear) { + Ok(m) => { + methods.push(m.report.method.name()); + elements += m.element_count(); + sj8 = sj8.min(m.report.min_scaled_jacobian); + m.report.method + } + Err(e) => { + let reason: String = e.to_string().chars().take(REASON_CHARS).collect(); + let ms = start.elapsed().as_secs_f64() * 1e3; + return Ok(vec![format!("no: {}: {reason}", body.name), String::new(), String::new(), String::new(), format!("{ms:.0}")]); + } + }; + // The hex20 run goes straight to the method that meshed the body. + let quadratic = HexOptions { method: m_method, ..quadratic.clone() }; + sj20 = match (sj20, hex_mesh(&s, &faces, &size, &quadratic)) { + (Ok(a), Ok(m)) => Ok(a.min(m.report.min_scaled_jacobian)), + (Err(e), _) => Err(e), + (_, Err(e)) => Err(format!("{}: {}", body.name, e.to_string().chars().take(REASON_CHARS).collect::())), + }; + } + methods.sort_unstable(); + methods.dedup(); + let ms = start.elapsed().as_secs_f64() * 1e3; + let sj20 = sj20.map_or_else(|e| format!("error: {e}"), |v| format!("{v:.3}")); + Ok(vec![methods.join(", "), elements.to_string(), format!("{sj8:.3}"), sj20, format!("{ms:.0}")]) + } +} + +/// The geometry's analytic face surface in the mesher's terms. +pub fn analytic(s: FaceSurface) -> Surface { + match s { + FaceSurface::Plane { normal, point } => Surface::Plane { point, normal }, + FaceSurface::Cylinder { origin, axis, radius } => Surface::Cylinder { origin, axis, radius }, + FaceSurface::Cone { apex, axis, half_angle } => Surface::Cone { apex, axis, half_angle }, + FaceSurface::Sphere { center, radius } => Surface::Sphere { center, radius }, + FaceSurface::Torus { center, axis, major, minor } => Surface::Torus { center, axis, major, minor }, + } +} + fn load(part: &Part) -> Result<(Loaded<'_>, f64), String> { let start = Instant::now(); let design = corpus::load(part).map_err(|e| e.to_string())?; diff --git a/docs/architecture.md b/docs/architecture.md index f568fcf..9aed28d 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -79,7 +79,17 @@ Rules beyond the layer order: recombination (Lee and Lo 1994) and one uniform split that lands on the boundary nodes, then smoothing and untangling; every quad is checked with exact predicates. One `SizeField` trait (`modelcraft-mesh::size`) sets interior sizes and edge divisions, with even counts per loop. - Mapped and swept hex build on these face meshes. Without geometry, `mesh.block` + Mapped and swept hex build on these face meshes (`modelcraft-mesh::hex`, `mesh.hex.generate`): + the B-rep topology is recovered from the face-tagged surface; a body is extruded (straight + multi-level sweep: the levels' outlines imprinted into one cross-section, cut at reflex corners, + regions mapped, O-gridded or free-meshed, then stacked), revolved (the same in a rotating + frame), swept one-to-one through rings of linking faces (least-squares affine layers, Roca and + Sarrate 2010), or mapped as a logical cube (Gordon and Thiel 1982); interval assignment + couples the counts of shared curves. Nodes on analytic faces (`Shape::face_surfaces`) are + exact, hex20 mid-edge nodes go onto the faces, boundary quads carry B-rep face tags, and the + result passes the certified validity gate plus tiling and fidelity checks. A body no method + fits gets each method's reason. + Without geometry, `mesh.block` (`modelcraft-mesh::block`) meshes one six-sided block from its eight corners and optional circular or elliptical edges by the same edge-based transfinite interpolation, as hex8, hex20, tet4/tet10 (Kuhn's six tetrahedra per cell) or wedge6, with grading; blocks join where their diff --git a/docs/features.md b/docs/features.md index 2fce278..bfa3880 100644 --- a/docs/features.md +++ b/docs/features.md @@ -6,7 +6,7 @@ Every entry is a toolbar button. **Ready** entries run their command; **coming s shown greyed with the roadmap milestone that brings them (ROADMAP.md). Scripts, the command line, the control channel and AI agents over MCP (docs/mcp.md) run the same commands. -52 ready, 51 coming soon. +54 ready, 49 coming soon. ## Geometry @@ -33,8 +33,8 @@ the control channel and AI agents over MCP (docs/mcp.md) run the same commands. | Tet | Delaunay Points — Delaunay tetrahedralization of a point list (debugging aid) | ready | `mesh.delaunay.points` | | 2D | 2D Mesh — Quad or triangle mesh of a planar face or polygon in the XY plane for plane stress, plane strain and axisymmetric models; edges tagged | ready | `mesh.planar` | | Structured Hex | Quad Faces — All-quad meshes of faces on shared edge nodes | coming soon (M5) | | -| Structured Hex | Mapped Hex — Mapped hex mesh of six-sided regions | coming soon (M5) | | -| Structured Hex | Sweep Hex — Hex mesh swept along extrusions and revolutions | coming soon (M5) | | +| Structured Hex | Mapped Hex — Mapped hex mesh of a six-sided block (curved faces allowed) | ready | `mesh.hex.mapped`, `mesh.hex.generate` | +| Structured Hex | Sweep Hex — Hex mesh swept along extrusions (several levels), around revolutions, or from a source face to a target | ready | `mesh.hex.sweep`, `mesh.hex.generate` | | Structured Hex | Decompose — Find sweepable and mappable regions and propose the splits | coming soon (M5) | | | Structured Hex | Transition — Conforming tet–hex interfaces with pyramids | coming soon (M5) | | | Automatic Hex | All-Hex — Octree dual all-hex mesh of any part | coming soon (M6) | | diff --git a/docs/meshing/corpus.md b/docs/meshing/corpus.md index 4a09133..86741d5 100644 --- a/docs/meshing/corpus.md +++ b/docs/meshing/corpus.md @@ -8,29 +8,29 @@ clock on the machine that ran it; every other column is deterministic. ## Results -| part | categories | bodies | faces | triangles | tol. mm | build ms | enclosed mm³ | SolidCraft mm³ | rel. diff | exact mm³ | boundary | genus | feature edges | check ms | CDT tets | Steiner | cavities | tet gate | CDT ms | -|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---| -| bent-tube | sweepable, thin-wall, genus | 1 | 18 | 4296 | 0.121 | 435 | 4213.1 | 4159.4 | 1.3e-2 | 4159.8 | ok, 2148 vertices | 1 | 268 | 10.9 | 15550 | 2338 | 0 | pass | 550 | -| bolted-assembly | assembly, blocky | 3 | 22 | 284 | 0.117 | 54 | 33004.7 | 32999.0 | 1.7e-4 | 32998.9 | ok, 144 vertices | 0, 1, 1 | 104 | 0.7 | 690 | 181 | 0 | pass | 31 | -| cross-cylinders | boolean, sharp | 1 | 10 | 1395 | 0.087 | 7682 | 32197.4 | 32357.8 | -5.0e-3 | 32365.8 | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … | | | | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … | | | | | -| filleted-bracket | fillet, genus, sweepable | 1 | 13 | 540 | 0.088 | 381 | 32162.3 | 32145.0 | 5.4e-4 | 32144.8 | ok, 268 vertices | 2 | 135 | 1.4 | 5582 | 1613 | 0 | pass | 187 | -| i-beam | sweepable, thin-wall | 1 | 14 | 44 | 0.208 | 4 | 61600.0 | 61600.0 | 0.0e0 | 61600.0 | ok, 24 vertices | 0 | 36 | 0.1 | 94 | 24 | 0 | pass | 3 | -| loft-transition | organic | 1 | 6 | 292 | 0.075 | 55 | 56390.2 | 56450.6 | -1.1e-3 | – | ok, 148 vertices | 0 | 40 | 0.7 | 1329 | 256 | 1 | pass | 33 | -| perforated-plate | genus, thin-wall, sweepable | 1 | 54 | 1644 | 0.144 | 265 | 26822.2 | 26764.6 | 2.2e-3 | 26764.2 | ok, 776 vertices | 24 | 828 | 4.6 | 33611 | 9798 | 0 | pass | 1195 | -| pinhole-block | small-feature, genus | 1 | 38 | 4044 | 0.073 | 202 | 49919.4 | 49909.6 | 2.0e-4 | 49909.5 | ok, 1992 vertices | 16 | 1548 | 12.4 | 96411 | 30221 | 0 | pass | 4408 | -| rounded-box | fillet | 1 | 26 | 2540 | 0.054 | 156 | 23325.7 | 23340.5 | -6.3e-4 | 23340.8 | ok, 1272 vertices | 0 | 216 | 6.3 | 11666 | 1913 | 0 | pass | 371 | -| shallow-groove | near-tangent, boolean | 1 | 8 | 62 | 0.075 | 11 | 47795.2 | 47782.3 | 2.7e-4 | 47781.5 | ok, 33 vertices | 0 | 24 | 0.2 | 267 | 87 | 0 | pass | 11 | -| sharp-wedge | sharp, thin-wall, sweepable | 1 | 5 | 8 | 0.083 | 2 | 3200.0 | 3200.0 | 0.0e0 | 3200.0 | ok, 6 vertices | 0 | 9 | 0.0 | 11 | 3 | 0 | pass | 0 | -| sphere-on-plate | near-tangent, boolean | 1 | 12 | 2826 | 0.098 | 703 | 69353.3 | 69509.0 | -2.2e-3 | 69509.7 | ok, 1415 vertices | 0 | 124 | 5.7 | 4898 | 309 | 9 | pass | 214 | -| spline-vase | organic | 1 | 11 | 3890 | 0.106 | 230 | 94891.2 | 94948.5 | -6.0e-4 | – | ok, 1947 vertices | 0 | 243 | 8.5 | 12461 | 1512 | 0 | pass | 419 | -| stepped-block | blocky, sharp | 1 | 14 | 96 | 0.078 | 30 | 47000.0 | 47000.0 | 0.0e0 | 47000.0 | ok, 48 vertices | 1 | 60 | 0.4 | 223 | 41 | 0 | pass | 6 | -| thin-ribs | thin-wall, blocky | 1 | 26 | 92 | 0.119 | 62 | 32000.0 | 32000.0 | -4.5e-16 | 32000.0 | ok, 48 vertices | 0 | 72 | 0.5 | 171 | 33 | 0 | pass | 6 | -| thin-wall-enclosure | thin-wall, blocky | 1 | 11 | 36 | 0.108 | 8 | 15564.0 | 15564.0 | 0.0e0 | 15564.0 | ok, 20 vertices | 0 | 28 | 0.1 | 191 | 47 | 2 | pass | 5 | -| analytic-sphere | analytic | 1 | 6 | 2784 | 0.069 | 263 | 33353.4 | 33507.9 | -4.6e-3 | 33510.3 | ok, 1394 vertices | 0 | 96 | 5.8 | 4211 | 58 | 0 | pass | 258 | -| analytic-torus | analytic, genus | 1 | 6 | 7230 | 0.109 | 345 | 37644.1 | 37893.7 | -6.6e-3 | 37899.3 | ok, 3615 vertices | 1 | 180 | 15.4 | 19790 | 466 | 0 | pass | 717 | -| analytic-cylinder | analytic, sweepable | 1 | 4 | 1004 | 0.049 | 25 | 12533.1 | 12564.7 | -2.5e-3 | 12566.4 | ok, 504 vertices | 0 | 106 | 2.7 | 7437 | 1568 | 0 | pass | 259 | -| analytic-pipe | analytic, genus, sweepable | 1 | 6 | 1448 | 0.104 | 43 | 94150.9 | 94238.0 | -9.3e-4 | 94247.8 | ok, 724 vertices | 1 | 220 | 3.7 | 6238 | 1340 | 0 | pass | 264 | -| analytic-plate-hole | analytic, genus, sweepable | 1 | 8 | 112 | 0.224 | 11 | 98449.1 | 98429.5 | 2.0e-4 | 98429.2 | ok, 56 vertices | 1 | 62 | 0.4 | 1868 | 543 | 0 | pass | 53 | +| part | categories | bodies | faces | triangles | tol. mm | build ms | enclosed mm³ | SolidCraft mm³ | rel. diff | exact mm³ | boundary | genus | feature edges | check ms | CDT tets | Steiner | cavities | tet gate | CDT ms | hex | hex elements | hex8 min SJ | hex20 min SJ | hex ms | +|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---| +| bent-tube | sweepable, thin-wall, genus | 1 | 18 | 4296 | 0.121 | 454 | 4213.1 | 4159.4 | 1.3e-2 | 4159.8 | ok, 2148 vertices | 1 | 268 | 11.5 | 15550 | 2338 | 0 | pass | 606 | sweep | 192 | 0.859 | 0.997 | 63 | +| bolted-assembly | assembly, blocky | 3 | 22 | 284 | 0.117 | 206 | 33004.7 | 32999.0 | 1.7e-4 | 32998.9 | ok, 144 vertices | 0, 1, 1 | 104 | 0.7 | 690 | 181 | 0 | pass | 34 | extrude | 357 | 0.398 | 0.222 | 40 | +| cross-cylinders | boolean, sharp | 1 | 10 | 1395 | 0.087 | 7920 | 32197.4 | 32357.8 | -5.0e-3 | 32365.8 | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … | | | | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … | | | | | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … | | | | | +| filleted-bracket | fillet, genus, sweepable | 1 | 13 | 540 | 0.088 | 388 | 32162.3 | 32145.0 | 5.4e-4 | 32144.8 | ok, 268 vertices | 2 | 135 | 1.4 | 5582 | 1613 | 0 | pass | 208 | no: Bracket: not supported yet: no structured hex method fits this body (extrude: cylindrical face 8 is not parallel to the direction; sweep: from face 1: linking face 6 has 2 boundary loops (periodic or holed linking faces are not swept); from face 6: | | | | 5 | +| i-beam | sweepable, thin-wall | 1 | 14 | 44 | 0.208 | 3 | 61600.0 | 61600.0 | 0.0e0 | 61600.0 | ok, 24 vertices | 0 | 36 | 0.1 | 94 | 24 | 0 | pass | 4 | extrude | 336 | 0.707 | 0.707 | 34 | +| loft-transition | organic | 1 | 6 | 292 | 0.075 | 57 | 56390.2 | 56450.6 | -1.1e-3 | – | ok, 148 vertices | 0 | 40 | 0.7 | 1329 | 256 | 1 | pass | 38 | sweep | 2873 | 0.189 | 0.189 | 328 | +| perforated-plate | genus, thin-wall, sweepable | 1 | 54 | 1644 | 0.144 | 273 | 26822.2 | 26764.6 | 2.2e-3 | 26764.2 | ok, 776 vertices | 24 | 828 | 4.8 | 33611 | 9798 | 0 | pass | 1345 | extrude | 467 | 0.304 | 0.033 | 98 | +| pinhole-block | small-feature, genus | 1 | 38 | 4044 | 0.073 | 201 | 49919.4 | 49909.6 | 2.0e-4 | 49909.5 | ok, 1992 vertices | 16 | 1548 | 12.9 | 96411 | 30221 | 0 | pass | 8142 | extrude | 5236 | 0.408 | 0.392 | 735 | +| rounded-box | fillet | 1 | 26 | 2540 | 0.054 | 235 | 23325.7 | 23340.5 | -6.3e-4 | 23340.8 | ok, 1272 vertices | 0 | 216 | 8.4 | 11663 | 1916 | 0 | pass | 469 | no: Box: not supported yet: no structured hex method fits this body (extrude: cylindrical face 8 is not parallel to the direction; sweep: from face 0: 12 faces are not on the sweep path; from face 5: 12 faces are not on the sweep path; revolve: plan | | | | 15 | +| shallow-groove | near-tangent, boolean | 1 | 8 | 62 | 0.075 | 11 | 47795.2 | 47782.3 | 2.7e-4 | 47781.5 | ok, 33 vertices | 0 | 24 | 0.2 | 267 | 87 | 0 | pass | 13 | extrude | 1960 | 0.991 | 0.986 | 201 | +| sharp-wedge | sharp, thin-wall, sweepable | 1 | 5 | 8 | 0.083 | 2 | 3200.0 | 3200.0 | 0.0e0 | 3200.0 | ok, 6 vertices | 0 | 9 | 0.0 | 11 | 3 | 0 | pass | 0 | extrude | 294 | 0.050 | 0.050 | 30 | +| sphere-on-plate | near-tangent, boolean | 1 | 12 | 2826 | 0.098 | 974 | 69353.3 | 69509.0 | -2.2e-3 | 69509.7 | ok, 1415 vertices | 0 | 124 | 7.2 | 4898 | 309 | 9 | pass | 260 | no: Plate: not supported yet: no structured hex method fits this body (extrude: face 6 is curved across the direction; sweep: from face 0: 6 faces are not on the sweep path; from face 5: linking face 11 has 2 boundary loops (periodic or holed linking | | | | 11 | +| spline-vase | organic | 1 | 11 | 3890 | 0.106 | 226 | 94891.2 | 94948.5 | -6.0e-4 | – | ok, 1947 vertices | 0 | 243 | 9.0 | 12461 | 1512 | 0 | pass | 456 | sweep | 4554 | 0.649 | 0.649 | 547 | +| stepped-block | blocky, sharp | 1 | 14 | 96 | 0.078 | 264 | 47000.0 | 47000.0 | 0.0e0 | 47000.0 | ok, 48 vertices | 1 | 60 | 0.4 | 223 | 41 | 0 | pass | 7 | extrude | 1416 | 0.999 | 0.999 | 142 | +| thin-ribs | thin-wall, blocky | 1 | 26 | 92 | 0.119 | 1004 | 32000.0 | 32000.0 | -4.5e-16 | 32000.0 | ok, 48 vertices | 0 | 72 | 0.6 | 171 | 33 | 0 | pass | 8 | extrude | 777 | 0.844 | 0.844 | 100 | +| thin-wall-enclosure | thin-wall, blocky | 1 | 11 | 36 | 0.108 | 340 | 15564.0 | 15564.0 | 0.0e0 | 15564.0 | ok, 20 vertices | 0 | 28 | 0.2 | 191 | 47 | 2 | pass | 11 | extrude | 894 | 0.297 | 0.297 | 112 | +| analytic-sphere | analytic | 1 | 6 | 2784 | 0.069 | 719 | 33353.4 | 33507.9 | -4.6e-3 | 33510.3 | ok, 1394 vertices | 0 | 96 | 12.6 | 4211 | 58 | 0 | pass | 440 | no: Sphere: not supported yet: no structured hex method fits this body (extrude: no direction to try; sweep: no planar face to sweep from; revolve: no axis to try; mapped: meshing failed: meshing failed: mesh validity gate: 8 inverted or degenerate ele | | | | 31 | +| analytic-torus | analytic, genus | 1 | 6 | 7230 | 0.109 | 944 | 37644.1 | 37893.7 | -6.6e-3 | 37899.3 | ok, 3615 vertices | 1 | 180 | 25.6 | 19790 | 466 | 0 | pass | 1412 | revolve | 1134 | 0.706 | 0.706 | 280 | +| analytic-cylinder | analytic, sweepable | 1 | 4 | 1004 | 0.049 | 190 | 12533.1 | 12564.7 | -2.5e-3 | 12566.4 | ok, 504 vertices | 0 | 106 | 3.4 | 7437 | 1568 | 0 | pass | 452 | extrude | 3200 | 0.707 | 0.707 | 585 | +| analytic-pipe | analytic, genus, sweepable | 1 | 6 | 1448 | 0.104 | 336 | 94150.9 | 94238.0 | -9.3e-4 | 94247.8 | ok, 724 vertices | 1 | 220 | 9.0 | 6238 | 1340 | 0 | pass | 479 | extrude | 1288 | 0.998 | 0.998 | 290 | +| analytic-plate-hole | analytic, genus, sweepable | 1 | 8 | 112 | 0.224 | 231 | 98449.1 | 98429.5 | 2.0e-4 | 98429.2 | ok, 56 vertices | 1 | 62 | 0.9 | 1868 | 543 | 0 | pass | 72 | extrude | 321 | 0.320 | 0.142 | 55 | ## Categories diff --git a/docs/meshing/literature.md b/docs/meshing/literature.md index f595e0a..0c29137 100644 --- a/docs/meshing/literature.md +++ b/docs/meshing/literature.md @@ -878,6 +878,65 @@ input points), Steiner points kept half the local spacing from the boundary. --- +### B8. Mapped and swept hexahedral meshing (implemented in `modelcraft-mesh::hex`) + +`mesh.hex.generate` (and `mesh.hex.sweep`, `mesh.hex.mapped`). Every method below is implemented from +its paper; no code was read. DOIs marked "checked" were confirmed on 2026-10-13 through a repository +record or the publisher; the others are cited by venue only. + +**B8.1 Blacker 1996, the Cooper tool.** **Verdict: BUILD ON (sweep recognition).** +- **Citation.** T. D. Blacker. *The Cooper tool.* Proc. 5th International Meshing Roundtable, 1996, pp. 13–29. +- **Use.** Source, target and linking faces from B-rep topology: rings of logically rectangular + linking faces followed from the source loops until they close on one target (`hex::sweep`). + +**B8.2 Lai, Benzley, White 2000, many-to-many sweeping.** **Verdict: BUILD ON (multi-level extrusions).** +- **Citation.** M. Lai, S. Benzley, D. White. *Automated hexahedral mesh generation by generalized multiple + source to multiple target sweeping.* IJNME 49(1–2):261–275, 2000. +- **Use.** Imprinting every level's outline into one projection plane so each region of the arrangement + is a one-to-one column (`hex::layered`, `hex::arrangement`); here done in a straight or rotating + frame, with the regions' fill decided by the generalized winding number (Jacobson et al. 2013). + +**B8.3 Knupp 1998, next-generation sweep tool.** **Verdict: BUILD ON (residual correction).** +- **Citation.** P. M. Knupp. *Next-generation sweep tool: a method for generating all-hex meshes on + two-and-one-half dimensional geometries.* Proc. 7th International Meshing Roundtable, 1998, pp. 505–513. +- **Use.** The boundary residual of each layer's map interpolated into the interior. + +**B8.4 Roca and Sarrate 2010, least-squares sweeping.** **Verdict: BUILD ON (inner nodes).** +- **Citation.** X. Roca, J. Sarrate. *An automatic and general least-squares projection procedure for sweep + meshing.* Engineering with Computers 26(4):391–406, 2010. DOI 10.1007/s00366-009-0172-z (checked, repository + record). Conference version: Proc. 15th IMR, 2006, pp. 487–506. +- **Use.** Each layer's interior nodes from the source mesh by the least-squares affine map of the source + ring onto the layer's ring (`hex::sweep::Affine`). Related: Roca, Sarrate, Huerta, *A new least-squares + approximation of affine mappings for sweep algorithms*, Eng. Comput. 26(3):327–337, 2010. + +**B8.5 Interval assignment: Tam and Armstrong 1993; Mitchell 1997/2021.** **Verdict: BUILD ON.** +- **Citations.** T. K. H. Tam, C. G. Armstrong. *Finite element mesh control by integer programming.* IJNME + 36:2581–2605, 1993. S. A. Mitchell. *High fidelity interval assignment.* Proc. 6th IMR, 1997 (IJCGA + 10(4):399–415, 2000). S. A. Mitchell. *Incremental interval assignment by integer linear algebra.* + International Meshing Roundtable 2021 (Zenodo 5559025); extended in Computer-Aided Design 158:103485, 2023. +- **Use.** Division counts of all curves together: side equalities of mapped faces and linking faces, + even (or multiple-of-four) loop sums, minimum loop sums. `hex::interval` uses incremental repair from the + rounded goals with increments only, chosen by squared relative error; it is bounded and reports an + infeasible system instead of looping. Mitchell's integer-linear-algebra solver (Hermite normal form plus + nullspace improvement) is the upgrade if the increment-only scheme proves too coarse. + +**B8.6 Whiteley, White, Benzley, Blacker 1996, submapping.** **Verdict: BUILD ON (vertex classification).** +- **Citation.** M. Whiteley, D. White, S. Benzley, T. Blacker. *Two and three-quarter dimensional meshing + facilitators.* Engineering with Computers 12:144–154, 1996. +- **Use.** Corners by interior angle (≤ 135° end, ≥ 225° reflex); cross-section regions are cut from + their reflex corners to the boundary, which decomposes rectilinear sections into rectangles (an explicit + form of submapping's virtual decomposition), then mapped. + +**B8.7 Gordon and Thiel 1982, trivariate transfinite interpolation.** **Verdict: BUILD ON (mapped blocks).** +- **Citation.** W. J. Gordon, L. C. Thiel. *Transfinite mappings and their application to grid generation.* + Applied Mathematics and Computation 10–11:171–233, 1982. +- **Use.** Interior nodes of a logical cube from its six face grids (`hex::mapped`). + +**B8.8 Support.** Generalized winding number for the inside test (A. Jacobson, L. Kavan, O. Sorkine-Hornung, +*Robust inside-outside segmentation using generalized winding numbers*, ACM TOG 32(4):33, 2013; solid angle +of a triangle: Van Oosterom and Strackee 1983); closest point on a triangle (C. Ericson, *Real-Time +Collision Detection*, 2005, §5.1.5). + ## Part C: Datasets and benchmarks | Dataset | Content | Licence | Usable for ModelCraft? | diff --git a/docs/parity.md b/docs/parity.md index 35e1edb..604c9b7 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -34,8 +34,8 @@ products are run, no documentation text is copied). Status: **done**, **started* | Mesh | Free tet meshing (tet4, tet10 with mid-side nodes on the geometry) | M2 | | started | | Mesh | Quality metrics, certified validity, verify per set | M1–M2 | yes | done (metrics, gate) | | Mesh | 2D quad meshing of faces | M5 | yes | started (`modelcraft_mesh::quad`; `mesh.planar`: 2D models of a polygon or a planar face, quad4/quad8/tri3/tri6, edges tagged by their neighbouring face) | -| Mesh | Mapped hex of six-sided regions | M5 | yes | planned | -| Mesh | Sweep hex (extrusions, revolutions) | M5 | yes | planned | +| Mesh | Mapped hex of six-sided regions | M5 | yes | done (`mesh.hex.mapped`) | +| Mesh | Sweep hex (extrusions, revolutions) | M5 | yes | done (`mesh.hex.sweep`: multi-level extrusions, revolutions, one-to-one sweeps) | | Mesh | Automatic decomposition into sweepable and mappable regions | M5 | yes | planned | | Mesh | Conforming interfaces, tet–hex transition with pyramids | M5 | yes | planned | | Mesh | Octree/grid-dual all-hex for any part | M6 | yes | planned |