From db0f0567a0a57b9d339c7c76cfda7affd13e0790 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Fri, 16 Oct 2026 10:58:02 -0400 Subject: [PATCH] mesh: the corner copy takes the corner from the whole cell's mesh, sheets at the cuts held planar Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo --- crates/engine/src/cmd/plan.rs | 100 ++++++++++++++++++++++++------- crates/engine/src/cmd/project.rs | 49 +++++++++++++++ crates/mesh/src/project/mod.rs | 73 +++++++++++++++++++++- 3 files changed, 196 insertions(+), 26 deletions(-) diff --git a/crates/engine/src/cmd/plan.rs b/crates/engine/src/cmd/plan.rs index 93f5533..329e75c 100644 --- a/crates/engine/src/cmd/plan.rs +++ b/crates/engine/src/cmd/plan.rs @@ -33,7 +33,7 @@ use std::collections::BTreeSet; use fieldcraft_geometry::{Curve, Design, FaceSurface, FeatureNode, Shape}; -use fieldcraft_mesh::project::Topology; +use fieldcraft_mesh::project::{ProjectOptions, Topology}; use fieldcraft_model::geometry::Entity; use fieldcraft_model::{CutRole, Datum, MeshPlan, Partition, PlanCut, PlanFeature, PlanMesher, PlanRing, PlanUnmapped, RingKind}; use serde_json::{Value, json}; @@ -1339,13 +1339,14 @@ fn apply_symmetric(s: &mut Session, name: &str, plan: &MeshPlan, corners: bool) undo(s)?; return Ok(None); } - let copied = if corners { copy_corner(s, name, plan, &cuts, &cell, &cell_plan, tolerance)? } else { None }; + let copied = if corners { copy_corner(s, name, plan, &cuts, &cell, &cell_plan, tolerance)? } else { Err(None) }; let (r, mut meshed, corner) = match copied { - Some(c) => (c.reply, c.cells, Some(c.summary)), - None => { + Ok(c) => (c.reply, c.cells, Some(c.summary)), + Err(why) => { + let corner = why.map(|w| json!({"skipped": w})); let cell_name = MeshPlan::default_name(&cell); s.model.set_mesh_plan(&cell_name, cell_plan.clone())?; - (apply_projection(s, &cell_name, &cell_plan, tolerance)?, vec![cell.clone()], None) + (apply_projection(s, &cell_name, &cell_plan, tolerance)?, vec![cell.clone()], corner) } }; // The cells now (the corner copy cuts the body further). @@ -1489,18 +1490,18 @@ fn copy_corner( cell: &str, cell_plan: &MeshPlan, tolerance: f64, -) -> Result> { +) -> Result>> { let design = s.design.as_deref().ok_or_else(|| EngineError::Params("open a design first".into()))?; - let (Some(shape), Some(d)) = (design.body(cell), cell_plan.direction) else { return Ok(None) }; + let (Some(shape), Some(d)) = (design.body(cell), cell_plan.direction) else { return Ok(Err(None)) }; let layout = match layout_of(design, cell_plan)? { Some(Ok(l)) => l, - _ => return Ok(None), + _ => return Ok(Err(None)), }; let m = shape.measure()?; let (lo, hi) = (m.bbox_min, m.bbox_max); let diagonal = (0..3).map(|k| (hi[k] - lo[k]).powi(2)).sum::().sqrt(); let eps = 1e-6 * diagonal.max(1.0); - let Some(mirror) = Mirror::of(shape)? else { return Ok(None) }; + let Some(mirror) = Mirror::of(shape)? else { return Ok(Err(None)) }; let mut found = None; for a in in_plane_axes(d).into_iter().filter(|a| AXES.iter().any(|x| parallel(*x, *a))) { if sym.iter().any(|c| parallel(c.normal, a)) { @@ -1513,14 +1514,14 @@ fn copy_corner( break; } } - let Some((a, at, c)) = found else { return Ok(None) }; + let Some((a, at, c)) = found else { return Ok(Err(None)) }; let before = s.snapshot(); - match corner_meshes(s, name, plan, sym, &layout.topology, a, at, c, tolerance) { - Ok(copy) => Ok(Some(copy)), + match corner_meshes(s, name, plan, sym, cell, &layout.topology, a, at, c, tolerance) { + Ok(copy) => Ok(Ok(copy)), // A corner copy that does not mesh leaves the cell to be meshed whole. - Err(_) => { + Err(e) => { s.restore(before); - Ok(None) + Ok(Err(Some(e.to_string()))) } } } @@ -1532,6 +1533,7 @@ fn corner_meshes( name: &str, plan: &MeshPlan, sym: &[PlanCut], + cell: &str, t: &Topology, a: V3, at: [f64; 2], @@ -1539,6 +1541,58 @@ fn corner_meshes( tolerance: f64, ) -> Result { let [x1, x2] = at; + // The whole cell meshed first, its layout split at both cuts: the low corner is taken from + // this mesh, so it is as good as the whole cell's. + let mut t = t.clone(); + split_at(&mut t, a, x1)?; + split_at(&mut t, a, x2)?; + let opts = fieldcraft_mesh::project::ProjectOptions { size: plan.size, ..Default::default() }; + let whole_body = super::project::open_body(s, cell, Some(tolerance), &json!({}))?; + let merge = super::project::merge_distance(None, whole_body.diagonal())?; + // The sheets at the cuts held planar, so the corner's cut face is the plane. + let sheets = ProjectOptions { planes: at.iter().map(|&x| (a.map(|k| k * x), a)).collect(), ..opts.clone() }; + let whole = super::project::project_only(&whole_body, cell, &t, &sheets)?; + let whole_min = whole.report.min_scaled_jacobian; + let along = |p: &[f64; 3]| dot(*p, a); + let hex_centre = |h: &[u32; 8]| h.iter().map(|&v| whole.nodes.get(v as usize).map_or(0.0, along)).sum::() / 8.0; + // The low corner's hexes, over their own nodes; the nodes on the cut sheet put on the plane. + let mut index: std::collections::BTreeMap = std::collections::BTreeMap::new(); + let mut corner_nodes: Vec<[f64; 3]> = Vec::new(); + let mut corner_hexes: Vec<[u32; 8]> = Vec::new(); + let mut beyond: std::collections::BTreeSet = std::collections::BTreeSet::new(); + for h in &whole.hexes { + if hex_centre(h) >= x1 { + beyond.extend(h.iter().copied()); + } + } + for h in whole.hexes.iter().filter(|h| hex_centre(h) < x1) { + corner_hexes.push(h.map(|v| { + *index.entry(v).or_insert_with(|| { + let mut p = whole.nodes.get(v as usize).copied().unwrap_or_default(); + if beyond.contains(&v) { + let off = x1 - along(&p); + p = std::array::from_fn(|k| p[k] + off * a[k]); + } + corner_nodes.push(p); + corner_nodes.len() as u32 - 1 + }) + })); + } + let corner_sj = corner_hexes + .iter() + .map(|h| { + let x = h.map(|v| corner_nodes.get(v as usize).copied().unwrap_or_default()); + fieldcraft_mesh::quality::element_quality(fieldcraft_model::ElementKind::Hex8, &x) + .ok() + .and_then(|q| q.get(fieldcraft_mesh::quality::Metric::ScaledJacobian)) + .unwrap_or(-1.0) + }) + .fold(f64::INFINITY, f64::min); + if corner_hexes.is_empty() || corner_sj < whole_min - 1e-9 { + return Err(EngineError::Params(format!( + "the corner taken from the whole cell's mesh has a smallest scaled Jacobian of {corner_sj}, below the whole cell's {whole_min}" + ))); + } let mut cuts: Vec = sym.to_vec(); for x in at { cuts.push(PlanCut { @@ -1573,16 +1627,10 @@ fn corner_meshes( } }; let (low, middle, high) = (pick(&|x| x < x1)?, pick(&|x| x > x1 && x < x2)?, pick(&|x| x > x2)?); - let mut t = t.clone(); - split_at(&mut t, a, x1)?; - split_at(&mut t, a, x2)?; let centre = |b: &[u32; 8]| b.iter().map(|&v| t.vertices.get(v as usize).map_or(0.0, |p| dot(*p, a))).sum::() / 8.0; - let low_blocks = sub_topology(&t, |b| centre(b) < x1); let middle_blocks = sub_topology(&t, |b| (x1..x2).contains(¢re(b))); - let opts = fieldcraft_mesh::project::ProjectOptions { size: plan.size, ..Default::default() }; let body = super::project::open_body(s, &low, Some(tolerance), &json!({}))?; - let merge = super::project::merge_distance(None, body.diagonal())?; - let corner = super::project::project_body(s, &low, body, &low_blocks, &opts, merge)?; + super::project::add_hexes(s, &body, &low, &corner_nodes, &corner_hexes, merge)?; (super::mirror::MIRROR.run)(s, &json!({"plane": {"point": a.map(|k| k * c), "normal": a}, "bodies": [low]}))?; let body = super::project::open_body(s, &middle, Some(tolerance), &json!({}))?; let (fixed, fixed_faces) = cut_face_patches(s, &body, &middle, plan.size)?; @@ -1591,7 +1639,13 @@ fn corner_meshes( } let opts = fieldcraft_mesh::project::ProjectOptions { fixed, ..opts }; let mut mesh = super::project::project_body(s, &middle, body, &middle_blocks, &opts, merge)?; - let corner_sj = corner["min_scaled_jacobian"].as_f64(); + let middle_sj = mesh["min_scaled_jacobian"].as_f64().unwrap_or(-1.0); + if middle_sj < whole_min - 1e-3 { + return Err(EngineError::Params(format!( + "the middle meshed round the corners has a smallest scaled Jacobian of {middle_sj}, below the whole cell's {whole_min}" + ))); + } + let corner_sj = Some(corner_sj); let summary = json!({ "axis": a, "cuts": at, @@ -1599,7 +1653,7 @@ fn corner_meshes( "image": high, "middle": middle, "fixed": fixed_faces, - "min_scaled_jacobian": {"corner": corner_sj, "middle": mesh["min_scaled_jacobian"]}, + "min_scaled_jacobian": {"whole_cell": whole_min, "corner": corner_sj, "middle": mesh["min_scaled_jacobian"]}, }); let min = [corner_sj, mesh["min_scaled_jacobian"].as_f64()].into_iter().flatten().reduce(f64::min); if let Some(o) = mesh.as_object_mut() { diff --git a/crates/engine/src/cmd/project.rs b/crates/engine/src/cmd/project.rs index 7aef6c9..59abcb5 100644 --- a/crates/engine/src/cmd/project.rs +++ b/crates/engine/src/cmd/project.rs @@ -544,11 +544,60 @@ fn hex_project(s: &mut Session, raw: &Value) -> Result { max_elements: fieldcraft_model::MAX_ELEMENTS / 8, min_quality, fixed, + planes: Vec::new(), }; let merge = merge_distance(p.merge, diag)?; project_body(s, &p.body, b, &t, &opts, merge) } +/// The projection mesh of topology `t` over body `b`, not added to the model. +pub(crate) fn project_only(b: &Body, name: &str, t: &Topology, opts: &ProjectOptions) -> Result { + project_mesh(t, &b.surface, &b.target, opts).map_err(|e| EngineError::of_body(name, e)) +} + +/// Replace body `name`'s mesh (opened as `b`) with `nodes` and `hexes` made for it elsewhere +/// (a part of another body's mesh): its outer faces (those no other of the hexes shares) tagged +/// with the B-rep face of `b` nearest their middle. Nodes within `merge` of other meshes' join +/// them. +pub(crate) fn add_hexes(s: &mut Session, b: &Body, name: &str, nodes: &[[f64; 3]], hexes: &[[u32; 8]], merge: f64) -> Result { + use fieldcraft_mesh::project::Geometry; + let faces = ElementKind::Hex8.faces(); + let mut count: BTreeMap<[u32; 4], usize> = BTreeMap::new(); + let quad = |h: &[u32; 8], f: &[usize]| -> [u32; 4] { [0, 1, 2, 3].map(|i| f.get(i).and_then(|&l| h.get(l)).copied().unwrap_or(0)) }; + for h in hexes { + for f in faces { + let mut k = quad(h, f); + k.sort_unstable(); + *count.entry(k).or_insert(0) += 1; + } + } + let face_count = b.input.face_names.len() as u32; + let mut boundary = Vec::new(); + for (e, h) in hexes.iter().enumerate() { + for (fi, f) in faces.iter().enumerate() { + let q = quad(h, f); + let mut k = q; + k.sort_unstable(); + if count.get(&k) != Some(&1) { + continue; + } + let mid = q.iter().fold([0.0; 3], |a, &v| { + let p = nodes.get(v as usize).copied().unwrap_or_default(); + [a[0] + 0.25 * p[0], a[1] + 0.25 * p[1], a[2] + 0.25 * p[2]] + }); + let nearest = (0..face_count) + .filter_map(|f| b.target.project(f, mid).map(|p| ((0..3).map(|i| (p[i] - mid[i]).powi(2)).sum::(), f))) + .min_by(|x, y| x.0.total_cmp(&y.0)); + if let Some((_, brep_face)) = nearest { + boundary.push(BoundaryFace { element: e as u32, face: fi as u8 + 1, brep_face }); + } + } + } + let _ = name; + let (reply, _, _) = add_mesh(s, &b.input, nodes, hexes, &boundary, merge)?; + Ok(reply) +} + /// Fraction of the element size within which a node of another mesh lies on a face. const ON_FACE: f64 = 0.05; diff --git a/crates/mesh/src/project/mod.rs b/crates/mesh/src/project/mod.rs index d0fc6b6..0fb5c49 100644 --- a/crates/mesh/src/project/mod.rs +++ b/crates/mesh/src/project/mod.rs @@ -94,6 +94,10 @@ pub struct ProjectOptions { /// Already-meshed patches of the boundary taken as given ([`fixed`]): the block faces on /// them take their divisions and nodes from the patch, held through smoothing. pub fixed: Vec, + /// Planes `(point, unit normal)` inner block faces are held on: a sheet of block faces whose + /// corners lie on one (within a tenth of the size) keeps its nodes on it through smoothing, so + /// the mesh can be cut there into two meshes with planar faces. + pub planes: Vec<(P, P)>, } impl Default for ProjectOptions { @@ -108,6 +112,7 @@ impl Default for ProjectOptions { max_elements: 2_000_000, min_quality: DEFAULT_MIN_QUALITY, fixed: Vec::new(), + planes: Vec::new(), } } } @@ -201,6 +206,34 @@ fn claim_narrow_faces(boundary: &mut [BoundaryFace], centres: &[(P, f64, usize)] } } +/// Most planes inner faces are held on. +pub const MAX_PLANES: usize = 64; + +/// The geometry with planes added as faces `first`, `first + 1`, …. +struct Planed<'a> { + g: &'a dyn Geometry, + first: u32, + planes: &'a [(P, P)], +} + +impl Geometry for Planed<'_> { + fn project(&self, face: u32, p: P) -> Option

{ + match face.checked_sub(self.first) { + Some(k) => { + let (o, n) = self.planes.get(k as usize)?; + Some(crate::vec3::sub(p, scale(*n, crate::vec3::dot(crate::vec3::sub(p, *o), *n)))) + } + None => self.g.project(face, p), + } + } + fn normal(&self, face: u32, p: P) -> Option

{ + match face.checked_sub(self.first) { + Some(k) => self.planes.get(k as usize).map(|x| x.1), + None => self.g.normal(face, p), + } + } +} + /// The face groups joined across smooth B-rep edges. fn smooth_groups(surface: &BoundarySurface, faces: usize) -> Vec { let mut parent: Vec = (0..faces).collect(); @@ -302,6 +335,13 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o let groups = smooth_groups(surface, face_count); let diag = t.diagonal().max(opts.size); let tol = 1e-9 * diag; + if opts.planes.len() > MAX_PLANES + || opts.planes.iter().any(|(p, n)| !(p.iter().chain(n).all(|x| x.is_finite()) && (crate::vec3::norm(*n) - 1.0).abs() < 1e-6)) + { + return Err(Error::Invalid(format!("at most {MAX_PLANES} planes, each a finite point and a unit normal"))); + } + let planed = Planed { g, first: face_count as u32, planes: &opts.planes }; + let g: &dyn Geometry = &planed; // Targets of the outer block faces. let explicit: BTreeMap<(usize, u8), Vec> = opts.faces.iter().map(|(b, f, s)| ((*b, *f), s.clone())).collect(); @@ -370,10 +410,31 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o face_target.insert(key, target_of(set)); } } + // Inner block faces on a plane: a target of their own. + let mut sheets: Vec<(usize, usize, usize)> = Vec::new(); + let mut sheet_targets: BTreeSet = BTreeSet::new(); + for (b, block) in t.blocks.iter().enumerate() { + for (d, s) in FACE_SIDES { + let key = face_key(block, d, s); + if uses.get(&key) != Some(&2) || face_target.contains_key(&key) { + continue; + } + let on = opts.planes.iter().position(|(p, n)| { + key.iter() + .all(|&v| t.vertices.get(v as usize).is_some_and(|x| crate::vec3::dot(crate::vec3::sub(*x, *p), *n).abs() <= 0.1 * opts.size)) + }); + if let Some(k) = on { + let tg = target_of(vec![face_count as u32 + k as u32]); + face_target.insert(key, tg); + sheet_targets.insert(tg); + sheets.push((b, d, s)); + } + } + } // Targets of the vertices and edges on the outside. let mut vertex_on: BTreeMap> = BTreeMap::new(); let mut edge_on: BTreeMap<[u32; 2], BTreeSet> = BTreeMap::new(); - for &(b, d, s) in &outer { + for &(b, d, s) in outer.iter().chain(&sheets) { let block = t.blocks.get(b).copied().unwrap_or([0; 8]); let tg = face_target.get(&face_key(&block, d, s)).copied().unwrap_or(0); let cyc = face_cycle(&block, d, s); @@ -690,6 +751,7 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o // The outer quads (outward, hex8 face order) and the target each lies on. let hex8_faces = ElementKind::Hex8.faces(); let mut quads: Vec<([u32; 4], u32, u8, usize)> = Vec::new(); + let mut sheet_quads: Vec<[u32; 4]> = Vec::new(); for (e, (hex, (b, cell))) in lat.hexes.iter().zip(&lat.hex_cell).enumerate() { let block = t.blocks.get(*b as usize).copied().unwrap_or([0; 8]); let nd = lat.divisions.get(*b as usize).copied().unwrap_or([1; 3]); @@ -700,10 +762,15 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o let Some(&tg) = face_target.get(&face_key(&block, d, s)) else { continue }; let local = hex8_faces.get(fi).copied().unwrap_or(&[]); let q = [0, 1, 2, 3].map(|i| local.get(i).and_then(|&l| hex.get(l)).copied().unwrap_or(0)); - quads.push((q, e as u32, fi as u8 + 1, tg)); + if sheet_targets.contains(&tg) { + sheet_quads.push(q); + } else { + quads.push((q, e as u32, fi as u8 + 1, tg)); + } } } - let boundary_quads: Vec<[u32; 4]> = quads.iter().map(|q| q.0).collect(); + // Nodes on a sheet slide along it among their neighbours on it, as boundary nodes do. + let boundary_quads: Vec<[u32; 4]> = quads.iter().map(|q| q.0).chain(sheet_quads).collect(); let smoothing = smooth::smooth( smooth::Mesh { nodes: &mut nodes, hexes: &lat.hexes, constraint: &constraint, held: &held, boundary: &boundary_quads }, g,