From 755576de8a922d0ad0faac2f002129fc122da459 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Thu, 15 Oct 2026 12:57:18 -0400 Subject: [PATCH] mesh.cells: sweep the cells of a body together, conforming with no seeds (#288) hex_mesh_swept_cells lays out one cross-section for several cells: every cell's walls (its cut faces included) are imprinted together, every cell's levels are levels of all, one interval assignment and one set of layers, and each region of each level interval is filled by the cell it lies in. Cells meet node for node across every cut, T-junctions included. mesh.cells sweeps the cells of each design body together first and meshes cell by cell only what does not sweep together (saying why); the reply gives every cell's volume, scaled Jacobian and edge ratio. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo --- crates/engine/src/cmd/cells.rs | 230 ++++++++++++++++++++++++-- crates/engine/src/cmd/hex.rs | 40 +++-- crates/geometry/src/design.rs | 10 ++ crates/mesh/src/hex/cells.rs | 113 +++++++++++++ crates/mesh/src/hex/layered.rs | 290 +++++++++++++++++++++------------ crates/mesh/src/hex/mod.rs | 33 ++-- 6 files changed, 580 insertions(+), 136 deletions(-) diff --git a/crates/engine/src/cmd/cells.rs b/crates/engine/src/cmd/cells.rs index 00a8128..7ae7373 100644 --- a/crates/engine/src/cmd/cells.rs +++ b/crates/engine/src/cmd/cells.rs @@ -2,6 +2,14 @@ //! typically made by `geometry.partition`) with its own method, all into one mesh, and the cells //! joined across the partitions that divide them. //! +//! **Sweeping together (#288).** The cells of one body of the design whose method may sweep +//! (`auto`, `hex`, `swept`, `extrude`, `revolve`) are first meshed together +//! ([`hex_mesh_swept_cells`]): one cross-section with every cell's cut faces imprinted, one +//! interval assignment and one set of layers, so cells meet node for node across every cut — a +//! cut that ends on another included — with no edge seeds, and their nodes are shared. When they +//! do not sweep together in one frame (the reply says why, per cell, as `not_swept_together`), +//! each cell is meshed alone as below. +//! //! **Joining.** Two cells on either side of a partition `P` meet on faces named //! `partition:P:a` and `partition:P:b`. Where the element faces of the two sides match node for //! node (the same divisions and kind on both sides, e.g. a mesh and its mirror image) the @@ -17,12 +25,14 @@ use std::collections::{BTreeMap, BTreeSet}; use fieldcraft_mesh::cdt::CdtOptions; -use fieldcraft_mesh::hex::{HexMesh, HexMethod}; +use fieldcraft_mesh::hex::{Cell, CellsMesh, HexMesh, HexMethod, HexOptions, hex_mesh_swept_cells}; +use fieldcraft_mesh::size::UniformSize; +use fieldcraft_mesh::surface::BoundarySurface; use fieldcraft_model::geometry::Entity; use fieldcraft_model::{ElementFace, GeometryRef, Model, Set, Tie}; use serde_json::{Value, json}; -use super::hex::{Input, add_hex, hex_body}; +use super::hex::{Input, add_cells, add_hex, hex_body}; use super::mesh::{BodyTets, add_tets, mesh_body, quality_options}; use crate::requirement::Requirement; use crate::{CommandKind, CommandSpec, EngineError, Result, Session}; @@ -30,7 +40,7 @@ use crate::{CommandKind, CommandSpec, EngineError, Result, Session}; pub const CELLS: CommandSpec = CommandSpec { id: "mesh.cells", label: "Partitioned Mesh", - params: r#"{"methods": {"Body1/a": "tet", "Body1/b": "swept"} (per cell; default `default`), "default": "auto", "bodies": ["Body1/a", "Body1/b"] (default: every body of the open design), "size": 2.0 (element size, mm; default 1/25 of the bounding-box diagonal of all of them), "order": 1 | 2 (hex8 and tet4, or hex20 and tet10), "tolerance": 0.05 (chordal tolerance, mm)} — replaces the cells' meshes (as mesh.tet.generate; other bodies keep theirs) with one mesh of the cells, each with its method: "auto" (structured hex where a method fits, else tet), "hex" (any structured hex method), "swept", "extrude", "sweep", "revolve", "mapped" or "tet"; or "default": "blocks" (with no `methods`): every cell a mapped block, division counts assigned over all cells together, exactly mesh.hex.blocks with the same parameters. Cells on either side of a partition (faces `partition::a` and `:b`) are joined: where their element faces match node for node the nodes are merged (conforming); otherwise, as between a hex and a tet cell, a set on each side and a tie between them are made (`tie__`, sets `_main` and `_secondary`: the side with fewer nodes is the main surface), replacing ties of those names from an earlier run. Reports each cell's method (and why auto fell back to tet), and each interface: `merged` with the nodes merged, or `tied`"#, + params: r#"{"methods": {"Body1/a": "tet", "Body1/b": "swept"} (per cell; default `default`), "default": "auto", "bodies": ["Body1/a", "Body1/b"] (default: every body of the open design), "size": 2.0 (element size, mm; default 1/25 of the bounding-box diagonal of all of them), "order": 1 | 2 (hex8 and tet4, or hex20 and tet10), "tolerance": 0.05 (chordal tolerance, mm)} — replaces the cells' meshes (as mesh.tet.generate; other bodies keep theirs) with one mesh of the cells, each with its method: "auto" (structured hex where a method fits, else tet), "hex" (any structured hex method), "swept", "extrude", "sweep", "revolve", "mapped" or "tet"; or "default": "blocks" (with no `methods`): every cell a mapped block, division counts assigned over all cells together, exactly mesh.hex.blocks with the same parameters. Cells on either side of a partition (faces `partition::a` and `:b`) are joined: where their element faces match node for node the nodes are merged (conforming); otherwise, as between a hex and a tet cell, a set on each side and a tie between them are made (`tie__`, sets `_main` and `_secondary`: the side with fewer nodes is the main surface), replacing ties of those names from an earlier run. The cells of one design body that may sweep (auto, hex, swept, extrude, revolve) are swept together first: one cross-section and one set of layers for all of them, so they share their nodes across every cut, T-junctions included, with no edge seeds (`swept_together` names the body; `not_swept_together` says why not, and those cells are meshed one by one). Reports each cell's method (and why auto fell back to tet), its `quality` (volume, smallest and mean scaled Jacobian, largest edge ratio), and each interface: `merged` with the nodes merged or shared, or `tied`"#, form: Some(super::forms::CELLS), kind: CommandKind::Model, requires: &[Requirement::Design], @@ -148,9 +158,55 @@ fn cells(s: &mut Session, p: &Value) -> Result { let cdt = CdtOptions { max_steiner: s.limits.max_steiner, max_tets: tet_options.max_tets }; let seeds = super::seeds::edge_seeds(s)?; + let method_of = |n: &str| methods.get(n).copied().unwrap_or(default); + // The cells of one body that may sweep are swept together first (#288): one cross-section + // and one set of layers for all of them, so they meet node for node across their cuts. Cells + // that do not sweep together are meshed one by one, as below. + let mut groups: BTreeMap<(String, HexMethod), Vec> = BTreeMap::new(); + for (i, shape) in shapes.iter().enumerate() { + let sweep = match method_of(&shape.name) { + Method::Auto | Method::Hex(HexMethod::Auto | HexMethod::Swept) => Some(HexMethod::Swept), + Method::Hex(m @ (HexMethod::Extrude | HexMethod::Revolve)) => Some(m), + Method::Hex(_) | Method::Tet => None, + }; + if let (Some(m), Some(source)) = (sweep, design.source_of(&shape.name)) { + groups.entry((source.to_string(), m)).or_default().push(i); + } + } + let mut joint: Vec> = vec![None; shapes.len()]; + let mut not_together: BTreeMap = BTreeMap::new(); + let mut swept: Vec<(String, Vec, Vec, CellsMesh)> = Vec::new(); + for ((source, m), members) in groups { + if members.len() < 2 { + continue; + } + let inputs: Vec = members + .iter() + .filter_map(|&i| shapes.get(i)) + .map(|shape| Ok(Input { seeds: seeds.clone(), ..super::hex::shape_input(shape, Some(tolerance))? })) + .collect::>()?; + match swept_together(p, &inputs, m, quadratic, size, tolerance, &seeds)? { + Ok(cm) => { + for &i in &members { + if let Some(slot) = joint.get_mut(i) { + *slot = Some(swept.len()); + } + } + swept.push((source, members, inputs, cm)); + } + Err(why) => { + for &i in &members { + not_together.insert(i, why.clone()); + } + } + } + } let mut meshed: Vec<(Meshed, Value)> = Vec::with_capacity(shapes.len()); - for shape in shapes { - let method = methods.get(&shape.name).copied().unwrap_or(default); + for (i, shape) in shapes.iter().enumerate() { + if joint.get(i).copied().flatten().is_some() { + continue; + } + let method = method_of(&shape.name); let input = || Ok::<_, EngineError>(Input { seeds: seeds.clone(), ..super::hex::shape_input(shape, Some(tolerance))? }); let hex = |m: HexMethod| -> Result<(Input<'static>, HexMesh)> { let b = input()?; @@ -175,6 +231,9 @@ fn cells(s: &mut Session, p: &Value) -> Result { }; if let Some(o) = info.as_object_mut() { o.insert("body".into(), json!(shape.name)); + if let Some(why) = not_together.get(&i) { + o.insert("not_swept_together".into(), json!(why)); + } } meshed.push((cell, info)); } @@ -183,9 +242,29 @@ fn cells(s: &mut Session, p: &Value) -> Result { // and the definitions stay. let model = &mut s.model; let cleared = super::replace::clear(model, names.iter().map(String::as_str)); - let mut report = Vec::with_capacity(meshed.len()); - for (cell, mut info) in meshed { - let (kind, elements, nodes) = match &cell { + let mut report: Vec<(usize, Value)> = Vec::with_capacity(shapes.len()); + let mut interfaces = Vec::new(); + for (source, members, inputs, cm) in &swept { + add_cells(model, inputs, cm)?; + let m = &cm.mesh; + for (c, (&i, b)) in members.iter().zip(inputs).enumerate() { + let elements: Vec = (0..m.element_count()).filter(|&e| cm.element_cell.get(e).is_some_and(|&x| x as usize == c)).collect(); + let nodes: BTreeSet = elements.iter().filter_map(|&e| m.element(e)).flatten().copied().collect(); + report.push(( + i, + json!({"body": b.name, "method": m.report.method.name(), "swept_together": source, "description": m.report.description, "kind": m.kind.name(), "elements": elements.len(), "nodes": nodes.len()}), + )); + } + interfaces.extend(swept_interfaces(&design, inputs, cm)); + } + let mut k = 0; + for (i, _) in shapes.iter().enumerate() { + if joint.get(i).copied().flatten().is_some() { + continue; + } + let Some((cell, mut info)) = meshed.get_mut(k).map(|(c, v)| (c, std::mem::take(v))) else { break }; + k += 1; + let (kind, elements, nodes) = match &*cell { Meshed::Hex(b, m) => { add_hex(model, b, m)?; (m.kind, m.element_count(), m.nodes.len()) @@ -200,13 +279,118 @@ fn cells(s: &mut Session, p: &Value) -> Result { o.insert("elements".into(), json!(elements)); o.insert("nodes".into(), json!(nodes)); } - report.push(info); + report.push((i, info)); } - let interfaces = join(model, &names, COINCIDENT * diagonal)?; + report.sort_by_key(|r| r.0); + let report: Vec = report + .into_iter() + .map(|(_, mut info)| { + let quality = info.get("body").and_then(Value::as_str).map(|b| cell_quality(model, b)); + if let (Some(o), Some(q)) = (info.as_object_mut(), quality) { + o.insert("quality".into(), q); + } + info + }) + .collect(); + interfaces.extend(join(model, &names, COINCIDENT * diagonal)?); let (nodes, elements) = (model.nodes().len(), model.elements().len()); Ok(cleared.report(s, json!({"cells": report, "interfaces": interfaces, "nodes": nodes, "elements": elements, "size": size}))) } +/// Sweep the cells `inputs` together ([`hex_mesh_swept_cells`]) with `method` (extrude, +/// revolve, or swept: either). `Ok(Err(why))` when they do not sweep together; limits and a +/// cancel are errors. +fn swept_together( + p: &Value, + inputs: &[Input], + method: HexMethod, + quadratic: bool, + size: Option, + tolerance: f64, + seeds: &[fieldcraft_mesh::seeds::EdgeSeed], +) -> Result> { + let h = size.ok_or_else(|| EngineError::Params("no element size".into()))?; + let mut surfaces = Vec::with_capacity(inputs.len()); + for b in inputs { + let options = super::mesh::surface_options(p, &b.mesh)?; + match BoundarySurface::new(&b.mesh.positions, &b.mesh.triangles, &b.mesh.tri_face, &options) { + Ok(surface) => surfaces.push(surface), + Err(e) => return Ok(Err(format!("body `{}`: {e}", b.name))), + } + } + let cells: Vec = surfaces.iter().zip(inputs).map(|(surface, b)| Cell { name: &b.name, surface, faces: &b.surfaces }).collect(); + let options = HexOptions { + method, + quadratic, + chordal_tolerance: tolerance, + max_elements: fieldcraft_model::MAX_ELEMENTS / 8, + seeds: seeds.to_vec(), + ..HexOptions::default() + }; + match hex_mesh_swept_cells(&cells, &UniformSize(h), &options) { + Ok(cm) => Ok(Ok(cm)), + Err(fieldcraft_mesh::Error::Model(fieldcraft_model::Error::Cancelled)) => Err(EngineError::Model(fieldcraft_model::Error::Cancelled)), + Err(e @ fieldcraft_mesh::Error::Limit(_)) => Err(EngineError::Params(e.to_string())), + Err(e) => Ok(Err(e.to_string())), + } +} + +/// The interfaces of cells swept together: each pair of cells whose elements share faces, with +/// the partition between them and the nodes they share (one mesh: merged by construction). +fn swept_interfaces(design: &fieldcraft_geometry::Design, inputs: &[Input], cm: &CellsMesh) -> Vec { + let m = &cm.mesh; + let mut owner: BTreeMap<[u32; 4], u32> = BTreeMap::new(); + let mut pairs: BTreeMap<(u32, u32), BTreeSet> = BTreeMap::new(); + for e in 0..m.element_count() { + let (Some(h), Some(&c)) = (m.element(e), cm.element_cell.get(e)) else { continue }; + for local in fieldcraft_model::ElementKind::Hex8.faces() { + let mut key = [0u32; 4]; + for (slot, &i) in key.iter_mut().zip(local.iter()) { + *slot = h.get(i).copied().unwrap_or(0); + } + key.sort_unstable(); + match owner.insert(key, c) { + Some(o) if o != c => pairs.entry((o.min(c), o.max(c))).or_default().extend(key), + _ => {} + } + } + } + // The partitions whose cut faces a cell has. + let cuts = |c: u32| -> BTreeSet { + let name = inputs.get(c as usize).map_or("", |b| b.name.as_str()); + let faces = design.body(name).map(|b| b.faces().iter().map(|f| f.name.clone()).collect::>()).unwrap_or_default(); + faces.iter().filter_map(|f| f.strip_prefix("partition:")?.rsplit_once(':').map(|(p, _)| p.to_string())).collect() + }; + pairs + .into_iter() + .map(|((a, b), nodes)| { + let (na, nb) = (inputs.get(a as usize).map_or("", |x| x.name.as_str()), inputs.get(b as usize).map_or("", |x| x.name.as_str())); + let partition = cuts(a).intersection(&cuts(b)).next().cloned().unwrap_or_default(); + json!({"partition": partition, "cells": [na, nb], "joined": "merged", "nodes": nodes.len(), "swept_together": true}) + }) + .collect() +} + +/// A cell's mesh in numbers: its volume, smallest and mean scaled Jacobian and largest edge +/// ratio. +fn cell_quality(model: &Model, body: &str) -> Value { + let (mut volume, mut min, mut sum, mut n, mut ratio) = (0.0, f64::INFINITY, 0.0, 0usize, 0.0f64); + for label in model.body_elements(body) { + let Some(e) = model.element(label) else { continue }; + let xyz: Vec<[f64; 3]> = e.nodes.iter().filter_map(|&l| model.node(l).map(|x| x.xyz)).collect(); + let Ok(q) = fieldcraft_mesh::quality::element_quality(e.kind, &xyz) else { continue }; + volume += q.size; + ratio = ratio.max(q.get(fieldcraft_mesh::quality::Metric::EdgeRatio).unwrap_or(0.0)); + if let Some(sj) = q.get(fieldcraft_mesh::quality::Metric::ScaledJacobian) { + min = min.min(sj); + sum += sj; + n += 1; + } + } + let min = if n > 0 { Some(min) } else { None }; + json!({"volume": volume, "min_scaled_jacobian": min, "mean_scaled_jacobian": if n > 0 { Some(sum / n as f64) } else { None }, "max_edge_ratio": ratio}) +} + /// One side of an interface: a cell's element faces on one side of a partition. struct Side { body: String, @@ -491,4 +675,30 @@ mod tests { assert!((u - expected).abs() < 1e-6 * expected, "{u} vs {expected}: {n}"); } } + + /// Cells cut from one box by planes, one of them (`y = 10`) cutting only the cell left of + /// `x = 10`, so a cut ends on another (a T-junction): with no seeds the cells sweep together + /// and every interface merges node for node (#288), and the mesh fills the box. + #[test] + fn cells_of_a_body_sweep_together_and_merge_without_seeds() { + let (_dir, mut s) = super::super::datum::tests::scripted(json!([ + {"command": "solid.box", "params": {"length": 30, "width": 20, "height": 10, "corner": [0, 0, 0], "body_name": "Box"}}, + ])); + super::super::datum::tests::cut_at(&mut s, &[(0, 10.0), (0, 21.0)]); + s.execute("geometry.datum", &json!({"name": "Y", "by": {"type": "offset", "plane": "zx", "offset": 7.5}})).unwrap(); + s.execute("geometry.partition", &json!({"name": "PY", "datum": "Y", "bodies": ["Box/a"]})).unwrap(); + let r = s.execute("mesh.cells", &json!({"default": "auto", "size": 2.0})).unwrap(); + let cells = r["cells"].as_array().unwrap(); + assert_eq!(cells.len(), 4, "{r}"); + assert!(cells.iter().all(|c| c["swept_together"] == "Box" && c["kind"] == "hex8"), "{r}"); + let interfaces = r["interfaces"].as_array().unwrap(); + assert_eq!(interfaces.len(), 4, "{r}"); + assert!(interfaces.iter().all(|i| i["joined"] == "merged"), "{r}"); + assert!(s.model.ties().is_empty()); + let q = s.execute("mesh.quality", &json!({})).unwrap(); + let volume = q["volume"].as_f64().unwrap(); + assert!((volume - 6000.0).abs() < 1e-6, "{q}"); + let d = s.execute("mesh.check.duplicates", &json!({})).unwrap(); + assert_eq!((d["counts"]["coincident"].as_u64(), d["counts"]["cracks"].as_u64()), (Some(0), Some(0)), "{d}"); + } } diff --git a/crates/engine/src/cmd/hex.rs b/crates/engine/src/cmd/hex.rs index f508b21..3f2695e 100644 --- a/crates/engine/src/cmd/hex.rs +++ b/crates/engine/src/cmd/hex.rs @@ -228,22 +228,7 @@ pub(super) fn hex_blocks(s: &mut Session, p: &Value) -> Result { // Replace the cells' meshes (replace.rs); other bodies and the definitions stay. let model = &mut s.model; let cleared = super::replace::clear(model, inputs.iter().map(|b| b.name.as_str())); - let (first_node, first_element) = super::fresh_labels(model, m.nodes.len(), m.element_count())?; - for (i, xyz) in m.nodes.iter().enumerate() { - model.set_node(first_node + i as u32, *xyz)?; - } - let mut body_tags = Vec::with_capacity(inputs.len()); - for b in &inputs { - body_tags.push(model.geometry_tag(&GeometryRef::body(b.name.as_str())?)?); - } - for e in 0..m.element_count() { - let cell = cm.element_cell.get(e).map_or(0, |&c| c as usize); - let nodes: Vec = m.element(e).unwrap_or(&[]).iter().map(|&i| first_node + i).collect(); - let l = first_element + e as u32; - model.set_element(l, m.kind, nodes)?; - model.tag_element(l, body_tags.get(cell).copied().ok_or_else(|| EngineError::Params("element without a cell".into()))?)?; - } - tag_faces(model, &inputs, |e| cm.element_cell.get(e as usize).map_or(0, |&c| c as usize), &m.boundary, first_element)?; + add_cells(model, &inputs, &cm)?; let blocks: Vec = inputs.iter().zip(&cm.blocks).map(|(b, d)| json!({"body": b.name, "divisions": d, "elements": d.iter().product::()})).collect(); let reply = json!({ @@ -355,6 +340,29 @@ pub(super) fn hex_body(p: &Value, b: &Input, method: HexMethod, quadratic: bool, hex_mesh(&surface, &b.surfaces, &UniformSize(h), &options).map_err(|e| EngineError::of_body(&b.name, e)) } +/// Add a mesh of several cells (`inputs[c]` is cell `c`) after the model's largest labels, its +/// nodes shared between the cells, tagging each element with its cell's body and the boundary +/// element faces with the cell's B-rep faces. +pub(super) fn add_cells(model: &mut Model, inputs: &[Input], cm: &fieldcraft_mesh::hex::CellsMesh) -> Result<()> { + let m = &cm.mesh; + let (first_node, first_element) = super::fresh_labels(model, m.nodes.len(), m.element_count())?; + for (i, xyz) in m.nodes.iter().enumerate() { + model.set_node(first_node + i as u32, *xyz)?; + } + let mut body_tags = Vec::with_capacity(inputs.len()); + for b in inputs { + body_tags.push(model.geometry_tag(&GeometryRef::body(b.name.as_str())?)?); + } + for e in 0..m.element_count() { + let cell = cm.element_cell.get(e).map_or(0, |&c| c as usize); + let nodes: Vec = m.element(e).unwrap_or(&[]).iter().map(|&i| first_node + i).collect(); + let l = first_element + e as u32; + model.set_element(l, m.kind, nodes)?; + model.tag_element(l, body_tags.get(cell).copied().ok_or_else(|| EngineError::Params("element without a cell".into()))?)?; + } + tag_faces(model, inputs, |e| cm.element_cell.get(e as usize).map_or(0, |&c| c as usize), &m.boundary, first_element) +} + /// Add one body's hex mesh after the model's largest labels, tagging the elements with the body /// and the boundary element faces with its B-rep faces; returns the first node and element /// labels. diff --git a/crates/geometry/src/design.rs b/crates/geometry/src/design.rs index 6e27361..493e4cc 100644 --- a/crates/geometry/src/design.rs +++ b/crates/geometry/src/design.rs @@ -76,6 +76,16 @@ impl Design { pub fn body(&self, name: &str) -> Option<&Shape> { self.bodies.iter().find(|b| b.name == name) } + + /// The body of the design as opened that `name` is, or that partitions cut it from + /// (`Base` for `Base/a/b1`). + pub fn source_of(&self, name: &str) -> Option<&str> { + self.source + .iter() + .map(|b| b.name.as_str()) + .filter(|s| name == *s || name.strip_prefix(s).is_some_and(|rest| rest.starts_with('/'))) + .max_by_key(|s| s.len()) + } } /// Open a design from its file bytes (a `.solidcraft` container or a JSON design) and rebuild it. diff --git a/crates/mesh/src/hex/cells.rs b/crates/mesh/src/hex/cells.rs index ac845f7..44bad3e 100644 --- a/crates/mesh/src/hex/cells.rs +++ b/crates/mesh/src/hex/cells.rs @@ -130,6 +130,119 @@ pub fn hex_mesh_cells(cells: &[Cell], field: &dyn SizeField, opts: &HexOptions) Ok(CellsMesh { mesh, element_cell, blocks, interfaces: shared }) } +/// Mesh `cells` (bodies that together make a part, meeting face to face) into one conforming +/// all-hex mesh by sweeping them together, along one direction or about one axis (`opts.method`: +/// [`HexMethod::Extrude`], [`HexMethod::Revolve`], or [`HexMethod::Auto`] / [`HexMethod::Swept`] +/// for the first that fits). +/// +/// This is many-to-many sweeping ([`super::layered`]) over the cells at once: the walls of every +/// cell (its cut faces included) are imprinted into one cross-section and the levels of every cell +/// are levels of all, so one interval assignment gives every curve and every level interval one +/// count, and each region of each interval is filled by the cell it lies in. Where two cells meet +/// they share the curves and layers there, so their elements meet node for node with no seeds; +/// a cut that ends on another cut (a T-junction) is a curve ending on a curve like any other. +/// Each element belongs to one cell and boundary faces are tagged with that cell's faces. Cells +/// that do not all sweep in one frame are [`Error::Unsupported`] with the first reason. +pub fn hex_mesh_swept_cells(cells: &[Cell], field: &dyn SizeField, opts: &HexOptions) -> Result { + let methods: &[HexMethod] = match opts.method { + HexMethod::Auto | HexMethod::Swept => &[HexMethod::Extrude, HexMethod::Revolve], + HexMethod::Extrude => &[HexMethod::Extrude], + HexMethod::Revolve => &[HexMethod::Revolve], + m => return Err(Error::Invalid(format!("cells are swept together by extrusion or revolution; method `{}` does not apply", m.name()))), + }; + 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 cells.is_empty() || cells.len() > MAX_CELLS { + return Err(Error::Invalid(format!("{} cells; 1 to {MAX_CELLS} are meshed together", cells.len()))); + } + if let Some(c) = cells.iter().find(|c| c.surface.shells().len() != 1) { + return Err(Error::Unsupported(format!("cell `{}` has {} shells; a cell is one solid without voids", c.name, c.surface.shells().len()))); + } + let surfaces: Vec<&BoundarySurface> = cells.iter().map(|c| c.surface).collect(); + super::check_estimate(&surfaces, field, opts)?; + let breps: Vec = cells.iter().map(|c| Brep::new(c.surface, c.faces)).collect::>()?; + let trees: Vec = cells.iter().map(|c| TriangleTree::new(c.surface)).collect::>()?; + // One classification tolerance for all, from their joint extent. + let (mut lo, mut hi) = ([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]); + for p in cells.iter().flat_map(|c| c.surface.points()) { + for k in 0..3 { + lo[k] = lo[k].min(p[k]); + hi[k] = hi[k].max(p[k]); + } + } + let tol = super::RELATIVE_TOLERANCE * dist(lo, hi); + if !(tol > 0.0 && tol.is_finite()) { + return Err(Error::Invalid("the cells have no extent".into())); + } + let ctxs: Vec = breps.iter().zip(&trees).map(|(brep, tree)| Ctx::with_tolerance(brep, tree, field, opts, tol)).collect::>()?; + let names: Vec<&str> = cells.iter().map(|c| c.name).collect(); + let mut reasons: Vec = Vec::new(); + for &method in methods { + let setups = match super::candidate_setups(&ctxs, method) { + Ok(s) => s, + Err(why) => { + reasons.push(format!("{}: {why}", method.name())); + continue; + } + }; + let mut last = None; + for setup in setups { + fieldcraft_model::progress::check()?; + let prepared = match super::layered::prepare(&ctxs, setup) { + Ok(p) => p, + Err(e @ Error::Limit(_)) => return Err(e), + Err(e) => { + last = Some(e); + continue; + } + }; + // As a single body's sweep: on a face mesher's failure, again with finer curves. + for refine in 0..=super::MAX_REFINE { + fieldcraft_model::progress::check()?; + let made = super::layered::mesh(&ctxs, &prepared, refine) + .and_then(|(raw, of)| Ok((finish(&Cells { ctxs: &ctxs, of: &of, names: &names }, raw, method)?, of))); + match made { + Ok((mesh, element_cell)) => { + let interfaces = shared_faces(&mesh, &element_cell); + return Ok(CellsMesh { mesh, element_cell, blocks: Vec::new(), interfaces }); + } + Err(e @ Error::Limit(_)) => return Err(e), + Err(e @ Error::Failed(_)) => last = Some(e), + Err(e) => { + last = Some(e); + break; + } + } + } + } + if let Some(e) = last { + reasons.push(format!("{}: {}", method.name(), super::plain(&e))); + } + } + Err(Error::Unsupported(format!("the cells do not sweep together ({})", reasons.join("; ")))) +} + +/// The pairs of cells that share element faces in `mesh`, counted once each. +fn shared_faces(mesh: &HexMesh, element_cell: &[u32]) -> usize { + let mut owner: std::collections::HashMap<[u32; 4], u32> = std::collections::HashMap::new(); + let mut pairs = std::collections::BTreeSet::new(); + for e in 0..mesh.element_count() { + let (Some(h), Some(&c)) = (mesh.element(e), element_cell.get(e)) else { continue }; + for k in 0..6 { + let mut key = super::hex_face(k).map(|i| h.get(i).copied().unwrap_or(0)); + key.sort_unstable(); + match owner.insert(key, c) { + Some(o) if o != c => { + pairs.insert((o.min(c), o.max(c))); + } + _ => {} + } + } + } + pairs.len() +} + /// Vertices, edges and faces shared by the cells. struct Topology { /// Shared vertices: where they are. diff --git a/crates/mesh/src/hex/layered.rs b/crates/mesh/src/hex/layered.rs index f43be6d..a43b72b 100644 --- a/crates/mesh/src/hex/layered.rs +++ b/crates/mesh/src/hex/layered.rs @@ -221,8 +221,8 @@ fn role(brep: &Brep, face: &Face, frame: &Frame, tol: f64) -> std::result::Resul } /// Candidate extrusion directions: the normals of planar faces and the axes of cylinders, most -/// common first. -pub(crate) fn translations(brep: &Brep) -> Vec { +/// common first, over the faces of every body given. +pub(crate) fn translations(breps: &[&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 }; @@ -236,7 +236,7 @@ pub(crate) fn translations(brep: &Brep) -> Vec { None => dirs.push((d, weight)), } }; - for f in &brep.faces { + for f in breps.iter().flat_map(|b| &b.faces) { if let Some((_, n)) = f.plane { push(n, f.area); } else if let Some(Surface::Cylinder { axis, .. }) = f.surface { @@ -244,76 +244,56 @@ pub(crate) fn translations(brep: &Brep) -> Vec { } } dirs.sort_by(|a, b| b.1.total_cmp(&a.1)); - let o = brep.point(0); + let Some(o) = breps.first().map(|b| b.point(0)) else { return Vec::new() }; 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 { +/// Candidate axes of revolution: the axes of the analytic surfaces of revolution of every body +/// given. +pub(crate) fn rotations(breps: &[&Brep]) -> Vec { let mut axes: Vec<(P, P)> = Vec::new(); - for f in &brep.faces { + let diagonal = breps.iter().map(|b| b.diagonal).fold(0.0, f64::max); + for f in breps.iter().flat_map(|b| &b.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) { + if !axes.iter().any(|(p, b)| parallel(*b, a) && line_distance(o, *p, *b) <= RELATIVE_AXIS * 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. +/// The frame's classification of the bodies swept together: face roles (per body) and levels. pub(crate) struct Setup { frame: Frame, - roles: Vec, + roles: Vec>, /// Sorted levels (heights or angles in `[0, 2π)`). levels: Vec, /// Sweep coordinates closer than this are one level. h_tol: f64, } -/// 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); +/// Classify faces and edges and find the levels of the bodies (one context each, sharing one +/// tolerance) swept together; an `Unsupported` error says why a body does not sweep in this +/// frame, naming it when there are several. +pub(crate) fn setup(ctxs: &[Ctx], frame: &Frame) -> Result { + let tol = ctxs.first().map_or(0.0, |c| c.tol); + let r_max = ctxs.iter().flat_map(|c| c.brep.surface.points()).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); - } + let (mut roles, mut heights) = (Vec::with_capacity(ctxs.len()), Vec::new()); + for (c, ctx) in ctxs.iter().enumerate() { + let named = |e: Error| match e { + Error::Unsupported(m) if ctxs.len() > 1 => Error::Unsupported(format!("cell {c}: {m}")), + e => e, + }; + roles.push(classify(ctx, frame, tol, &mut heights).map_err(named)?); } heights.sort_by(f64::total_cmp); let mut levels: Vec = Vec::new(); @@ -336,6 +316,42 @@ pub(crate) fn setup(ctx: &Ctx, frame: &Frame) -> Result { Ok(Setup { frame: *frame, roles, levels, h_tol }) } +/// One body's face roles in the frame; the heights of its edges across the sweep and of its caps +/// are added to `heights`. +fn classify(ctx: &Ctx, frame: &Frame, tol: f64, heights: &mut Vec) -> Result> { + let brep = ctx.brep; + 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"))); + } + } + 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); + } + } + Ok(roles) +} + impl Setup { /// The number of levels the frame sweeps through (for ordering candidates). pub(crate) fn level_count(&self) -> usize { @@ -356,50 +372,74 @@ impl SizeField for PlaneField<'_> { } } -/// The cross-section of a frame the body sweeps in: the arrangement (cut into blocks), -/// the level intervals and which of them each region fills. It does not depend on the element -/// size, so refinement retries reuse it. +/// The cross-section of a frame the bodies sweep in: the arrangement (cut into blocks), +/// the level intervals and which body fills each region in each of them. 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>, + filled: Vec>>, used: Vec, + /// The first arrangement tag of each body's faces (curve tags number the faces of all bodies + /// one after the other). + offsets: Vec, } -/// Lay out the cross-section of `setup`'s frame. -pub(crate) fn prepare(ctx: &Ctx, setup: Setup) -> Result { - let brep = ctx.brep; +impl Prepared { + /// The body and face of an arrangement tag. + fn face(&self, ctxs: &[Ctx], tag: usize) -> Option<(usize, usize)> { + let c = self.offsets.iter().rposition(|&o| o <= tag)?; + let f = tag - self.offsets.get(c)?; + ctxs.get(c).filter(|ctx| f < ctx.brep.faces.len()).map(|_| (c, f)) + } +} + +/// Lay out the cross-section of `setup`'s frame, for the bodies of `ctxs` (those `setup` was +/// made for) swept together: the walls of every body are imprinted on each other, so where two +/// bodies meet on a wall they get the same curves, and each region of each level interval is +/// filled by the body it lies in. +pub(crate) fn prepare(ctxs: &[Ctx], setup: Setup) -> Result { + let ctx = ctxs.first().ok_or_else(|| Error::Failed("no body to sweep".into()))?; let tol = ctx.tol; let frame = &setup.frame; - let (roles, levels) = (&setup.roles, &setup.levels); - // Lower outlines of the walls, and every vertex. + let levels = &setup.levels; + // Lower outlines of the walls, and every vertex, of every body. 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) { + let mut vertices: Vec = Vec::new(); + let mut offsets = Vec::with_capacity(ctxs.len()); + let mut offset = 0; + for (c, body) in ctxs.iter().enumerate() { + let brep = body.brep; + let roles = setup.roles.get(c).map(Vec::as_slice).unwrap_or(&[]); + offsets.push(offset); + 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; } - // 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)); + 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(v) = brep.third_vertex(e, f) else { continue }; + let (qc, hc) = frame.coords(brep.point(v)); + if frame.dh(hc, h0) > 0.0 && frame.span(qc, frame.dh(hc, h0)) > tol { + polylines.push((pts.iter().map(|x| x.0).collect(), offset + f)); + } } } + vertices.extend( + 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), + ); + offset += brep.faces.len(); } - 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 { @@ -407,7 +447,7 @@ pub(crate) fn prepare(ctx: &Ctx, setup: Setup) -> Result { }; 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 mut inside = Occupants::new(ctxs, tol); let mut arr = arrangement::build(&vertices, &polylines, tol, coincide)?; let (mut filled, mut used) = occupancy(&arr, frame, &intervals, &mut inside)?; // Cut the filled regions into blocks the structured schemes fit, and lay the arrangement out @@ -415,20 +455,21 @@ pub(crate) fn prepare(ctx: &Ctx, setup: Setup) -> Result { let h_ref = levels.first().copied().unwrap_or(0.0); let plane = PlaneField { field: ctx.field, frame: *frame, h: h_ref }; let settings = decompose::Settings { field: &plane, tol, opts: &ctx.opts.quad }; - let locked = seeded_curves(ctx, frame, h_ref, &intervals, setup.h_tol, &arr); + let locked = seeded_curves(ctxs, frame, h_ref, &intervals, setup.h_tol, &arr); let cut_lines = decompose::cuts(&arr, &used, &locked, &settings)?; if !cut_lines.is_empty() { polylines.extend(cut_lines.into_iter().map(|c| (c, decompose::CUT))); arr = arrangement::build(&vertices, &polylines, tol, coincide)?; (filled, used) = occupancy(&arr, frame, &intervals, &mut inside)?; } - Ok(Prepared { setup, arr, intervals, filled, used }) + Ok(Prepared { setup, arr, intervals, filled, used, offsets }) } -/// 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; +/// Mesh the prepared cross-section, with its curves `2^refine` times finer than the size asks; +/// with the body (index into `ctxs`) of each hexahedron. +pub(crate) fn mesh(ctxs: &[Ctx], prepared: &Prepared, refine: u32) -> Result<(Raw, Vec)> { + let ctx = ctxs.first().ok_or_else(|| Error::Failed("no body to sweep".into()))?; + 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. @@ -446,7 +487,7 @@ pub(crate) fn mesh(ctx: &Ctx, prepared: &Prepared, refine: u32) -> Result { ia.var(goal * fine, 1); } let schemes = section::schemes(arr, used, &mut ia); - let seeds = Seeds::of(ctx, frame, h_ref, intervals, setup.h_tol, &(0..arr.curves.len()).map(curve3d).collect::>())?; + let seeds = Seeds::of(ctxs, frame, h_ref, intervals, setup.h_tol, &(0..arr.curves.len()).map(curve3d).collect::>())?; for (c, seed) in seeds.curves.iter().enumerate() { if let Some((seed, _)) = seed { ia.seed(c, seed); @@ -493,7 +534,12 @@ pub(crate) fn mesh(ctx: &Ctx, prepared: &Prepared, refine: u32) -> Result { let wall = arr .curves .get(c) - .and_then(|curve| curve.tags.iter().find_map(|&f| brep.faces.get(f).and_then(|x| x.surface))) + .and_then(|curve| { + curve.tags.iter().find_map(|&t| { + let (c, f) = prepared.face(ctxs, t)?; + ctxs.get(c)?.brep.faces.get(f)?.surface + }) + }) .filter(|s| !matches!(s, Surface::Plane { .. })); Ok(pts3.iter().map(|p| frame.coords(wall.and_then(|s| s.project(*p)).unwrap_or(*p)).0).collect()) }; @@ -523,7 +569,7 @@ pub(crate) fn mesh(ctx: &Ctx, prepared: &Prepared, refine: u32) -> Result { 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 { + if f.is_some() { count = count.saturating_add(quads.len().saturating_mul(layers.get(k).copied().unwrap_or(0))); } } @@ -543,25 +589,26 @@ pub(crate) fn mesh(ctx: &Ctx, prepared: &Prepared, refine: u32) -> Result { }) }; let mut hexes: Vec<[u32; 8]> = Vec::with_capacity(count); + let mut cells: Vec = 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 Some(cell) = *f else { 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]]); + cells.push(cell); } } } } 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 bodies = if ctxs.len() > 1 { format!(" ({} cells together)", ctxs.len()) } else { String::new() }; let description = format!( - "{}: {} level{}, {} region{} ({} structured, {} free), {} layer{}", + "{}{bodies}: {} level{}, {} region{} ({} structured, {} free), {} layer{}", frame.describe(), levels.len(), if levels.len() == 1 { "" } else { "s" }, @@ -572,7 +619,7 @@ pub(crate) fn mesh(ctx: &Ctx, prepared: &Prepared, refine: u32) -> Result { total_layers.saturating_sub(usize::from(!frame.periodic())), if total_layers == 2 { "" } else { "s" }, ); - Ok(Raw { nodes, hexes, description }) + Ok((Raw { nodes, hexes, description }, cells)) } /// The edge seeds as a layered mesh takes them: on the cross-section's curves (seeds of edges @@ -633,10 +680,19 @@ fn split_over(nodes: &[f64], bounds: &[f64]) -> Option>> { } impl<'s> Seeds<'s> { - fn of(ctx: &'s Ctx, frame: &Frame, h_ref: f64, intervals: &[(f64, f64)], h_tol: f64, curves: &[Vec

]) -> Result> { + fn of(ctxs: &'s [Ctx], frame: &Frame, h_ref: f64, intervals: &[(f64, f64)], h_tol: f64, curves: &[Vec

]) -> Result> { + let mut out = Seeds { curves: vec![None; curves.len()], layers: vec![None; intervals.len()] }; + for ctx in ctxs { + out.add(ctx, frame, h_ref, intervals, h_tol, curves)?; + } + Ok(out) + } + + /// Add the seeds on the edges of one body. + fn add(&mut self, ctx: &'s Ctx, frame: &Frame, h_ref: f64, intervals: &[(f64, f64)], h_tol: f64, curves: &[Vec

]) -> Result<()> { let brep = ctx.brep; let tol = ctx.tol; - let mut out = Seeds { curves: vec![None; curves.len()], layers: vec![None; intervals.len()] }; + let out = self; let stol = super::seed_tolerance(ctx.opts, tol); let conflict = |a: &EdgeSeed, b: &EdgeSeed, why: &str| { Error::Invalid(format!( @@ -710,26 +766,60 @@ impl<'s> Seeds<'s> { *slot = Some(LayerSeed { seed, fractions }); } } - Ok(out) + Ok(()) } } -/// 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)> { +/// The bodies swept together, for inside tests: each body's test and bounding box. +struct Occupants<'a> { + bodies: Vec<(Inside<'a>, P, P)>, +} + +impl<'a> Occupants<'a> { + fn new(ctxs: &[Ctx<'a>], tol: f64) -> Self { + let bodies = ctxs + .iter() + .map(|c| { + let (mut lo, mut hi) = ([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]); + for p in c.brep.surface.points() { + for k in 0..3 { + lo[k] = lo[k].min(p[k] - tol); + hi[k] = hi[k].max(p[k] + tol); + } + } + (Inside::new(c.brep.surface), lo, hi) + }) + .collect(); + Occupants { bodies } + } + + /// The first body `x` lies in. + fn body_at(&mut self, x: P) -> Result> { + for (i, (inside, lo, hi)) in self.bodies.iter_mut().enumerate() { + if (0..3).all(|k| x[k] >= lo[k] && x[k] <= hi[k]) && inside.contains(x)? { + return Ok(u32::try_from(i).ok()); + } + } + Ok(None) + } +} + +/// Which body fills each region of the arrangement in each level interval (inside tests at the +/// region's centre, mid-interval), and which regions any body fills. +fn occupancy(arr: &Arrangement, frame: &Frame, intervals: &[(f64, f64)], inside: &mut Occupants) -> 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()); + 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)))?); + row.push(inside.body_at(frame.point(c, 0.5 * (lo + hi)))?); } filled.push(row); } - let used: Vec = filled.iter().map(|r| r.iter().any(|&x| x)).collect(); + let used: Vec = filled.iter().map(|r| r.iter().any(Option::is_some)).collect(); if !used.iter().any(|&u| u) { return Err(Error::Unsupported("no region of the cross-section lies inside the body".into())); } @@ -737,10 +827,10 @@ fn occupancy(arr: &Arrangement, frame: &Frame, intervals: &[(f64, f64)], inside: } /// Which curves of `arr` carry an edge seed (a cut may not end inside them). -fn seeded_curves(ctx: &Ctx, frame: &Frame, h_ref: f64, intervals: &[(f64, f64)], h_tol: f64, arr: &Arrangement) -> Vec { +fn seeded_curves(ctxs: &[Ctx], frame: &Frame, h_ref: f64, intervals: &[(f64, f64)], h_tol: f64, arr: &Arrangement) -> Vec { let curves: Vec> = arr.curves.iter().map(|cv| cv.points.iter().map(|&p| frame.point(arr.points.get(p).copied().unwrap_or([0.0; 2]), h_ref)).collect()).collect(); - match Seeds::of(ctx, frame, h_ref, intervals, h_tol, &curves) { + match Seeds::of(ctxs, frame, h_ref, intervals, h_tol, &curves) { Ok(seeds) => seeds.curves.iter().map(Option::is_some).collect(), Err(_) => vec![false; arr.curves.len()], } diff --git a/crates/mesh/src/hex/mod.rs b/crates/mesh/src/hex/mod.rs index c061781..a9e5e2f 100644 --- a/crates/mesh/src/hex/mod.rs +++ b/crates/mesh/src/hex/mod.rs @@ -22,7 +22,9 @@ //! **Partitioned parts.** [`hex_mesh_cells`] ([`cells`]) meshes a part split into cells that //! meet face to face (SolidCraft split bodies) as one conforming mesh: every cell a mapped block, //! interval assignment over all cells together, each shared face meshed once and its nodes -//! merged by identity. +//! merged by identity. [`hex_mesh_swept_cells`] sweeps such cells together instead: one +//! cross-section imprinted with every cell's walls, one interval assignment and one set of +//! layers, so cells that share a sweep conform across their cut faces. //! //! **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 @@ -51,7 +53,7 @@ mod tests; use std::collections::HashMap; -pub use cells::{Cell, CellsMesh, GROWTH, MAX_CELLS, hex_mesh_cells}; +pub use cells::{Cell, CellsMesh, GROWTH, MAX_CELLS, hex_mesh_cells, hex_mesh_swept_cells}; pub use geom::Surface; use crate::seeds::EdgeSeed; @@ -197,7 +199,17 @@ pub(crate) struct Ctx<'a> { impl<'a> Ctx<'a> { pub(crate) fn new(brep: &'a Brep<'a>, tree: &'a TriangleTree<'a>, field: &'a dyn SizeField, opts: &'a HexOptions) -> Result> { - let tol = RELATIVE_TOLERANCE * brep.diagonal; + Ctx::with_tolerance(brep, tree, field, opts, RELATIVE_TOLERANCE * brep.diagonal) + } + + /// A context with the classification tolerance `tol` (bodies meshed together share one). + pub(crate) fn with_tolerance( + brep: &'a Brep<'a>, + tree: &'a TriangleTree<'a>, + field: &'a dyn SizeField, + opts: &'a HexOptions, + tol: f64, + ) -> Result> { let curves: Vec> = (0..brep.edges.len()).map(|e| brep.edge_points(e, true)).collect(); let edge_seeds = crate::seeds::match_curves(&opts.seeds, &curves, seed_tolerance(opts, tol))?; Ok(Ctx { brep, tree, field, opts, tol, edge_seeds }) @@ -256,7 +268,7 @@ pub(crate) struct Raw { const MAX_REFINE: u32 = 2; /// Classification tolerance relative to the bounding-box diagonal. -const RELATIVE_TOLERANCE: f64 = 1e-6; +pub(crate) 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); @@ -290,7 +302,7 @@ pub fn hex_mesh(surface: &BoundarySurface, faces: &[Option], field: &dy let mut last = None; match method { HexMethod::Extrude | HexMethod::Revolve => { - let setups = match candidate_setups(&ctx, method) { + let setups = match candidate_setups(std::slice::from_ref(&ctx), method) { Ok(s) => s, Err(why) => { reasons.push(format!("{}: {why}", method.name())); @@ -300,7 +312,7 @@ pub fn hex_mesh(surface: &BoundarySurface, faces: &[Option], field: &dy for setup in setups { fieldcraft_model::progress::check()?; // The cross-section is laid out once per frame; refinement retries reuse it. - let prepared = match layered::prepare(&ctx, setup) { + let prepared = match layered::prepare(std::slice::from_ref(&ctx), setup) { Ok(p) => p, Err(e @ Error::Limit(_)) => return Err(e), Err(e) => { @@ -308,7 +320,7 @@ pub fn hex_mesh(surface: &BoundarySurface, faces: &[Option], field: &dy continue; } }; - match with_refinement(&ctx, method, |refine| layered::mesh(&ctx, &prepared, refine)) { + match with_refinement(&ctx, method, |refine| layered::mesh(std::slice::from_ref(&ctx), &prepared, refine).map(|m| m.0)) { Ok(mesh) => return Ok(mesh), Err(e @ Error::Limit(_)) => return Err(e), Err(e) => last = Some(e), @@ -351,12 +363,13 @@ fn check_estimate(surfaces: &[&BoundarySurface], field: &dyn SizeField, opts: &H } /// 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) }; +pub(crate) fn candidate_setups(ctxs: &[Ctx], method: HexMethod) -> std::result::Result, String> { + let breps: Vec<&Brep> = ctxs.iter().map(|c| c.brep).collect(); + let frames = if method == HexMethod::Extrude { layered::translations(&breps) } else { layered::rotations(&breps) }; 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) { + match layered::setup(ctxs, &f) { Ok(s) => fitting.push((s.level_count(), k, s)), Err(e) => { first_reason.get_or_insert_with(|| plain(&e));