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 <[email protected]>
Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo

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<Option<CornerCopy>> {
+) -> Result<std::result::Result<CornerCopy, Option<String>>> {
     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::<f64>().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<CornerCopy> {
     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::<f64>() / 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<u32, u32> = 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<u32> = 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<PlanCut> = 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::<f64>() / 8.0;
-    let low_blocks = sub_topology(&t, |b| centre(b) < x1);
     let middle_blocks = sub_topology(&t, |b| (x1..x2).contains(&centre(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() {

crates/engine/src/cmd/project.rs

@@ -544,11 +544,60 @@ fn hex_project(s: &mut Session, raw: &Value) -> Result<Value> {
         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<fieldcraft_mesh::project::ProjectedMesh> {
+    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<Value> {
+    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::<f64>(), 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;
 

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<FixedPatch>,
+    /// 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<P> {
+        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<P> {
+        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<usize> {
     let mut parent: Vec<usize> = (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<u32>> = 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<usize> = 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<u32, BTreeSet<usize>> = BTreeMap::new();
     let mut edge_on: BTreeMap<[u32; 2], BTreeSet<usize>> = 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,