diff --git a/crates/engine/src/cmd/project.rs b/crates/engine/src/cmd/project.rs index 7469200..879add9 100644 --- a/crates/engine/src/cmd/project.rs +++ b/crates/engine/src/cmd/project.rs @@ -19,7 +19,7 @@ use crate::{CommandKind, CommandSpec, EngineError, Result, Session}; pub const HEX_PROJECT: CommandSpec = CommandSpec { id: "mesh.hex.project", label: "Projection Blocking", - params: r#"{"body": "Body1", "size": 2.0 (element size, mm; default 1/25 of the body's diagonal), "template": "grid" (default) | "sphere" | "butterfly" | "blocks", "center": [x, y, z], "radius": r, "axis": [0, 0, 1], "levels": [z0, z1, …] (sphere and butterfly: default from the body's bounding box; butterfly levels are distances along the axis from the centre), "inset": 0.55 (inner square or cube corner at this fraction of the radius), "grid": {"planes": {"x": [...], "y": [...], "z": [...]} or "cells": [ni, nj, nk] (even planes over the bounding box, default [1, 1, 1]), "carve": [{"i": [a, b], "j": [c, d], "k": [e, f]}] (block index ranges, inclusive, removed)}, "blocks": {"vertices": [[x, y, z], ...], "blocks": [[v1..v8], ...] (0-based vertices, hex8 order)}, "operations": [{"op": "o_grid", "block": b, "open": [faces], "inset": 0.5 or [i1, i2, i3], "extrude": {"offset": [dx, dy, dz], "layers": 1, "faces": "all" (the open face's butterfly: a round branch) | "inner" (its inner face: a square stub)}} | {"op": "extrude", "faces": [[block, face], ...], "offset": [dx, dy, dz], "layers": 1} | {"op": "connect", "from": [block, face], "to": [block, face], "layers": 1} (blocks bridging two faces: closes a handle) | {"op": "remove", "blocks": [...]} (carve; later blocks renumber) | {"op": "merge", "distance": d} | {"op": "split", "block": b, "direction": 1|2|3 or "x"|"y"|"z", "at": [fractions]} (cuts the block and the whole sheet of blocks the cut runs through, so the topology stays conforming; fractions from the block's side-0 face, or from the low end along an axis; the block keeps its first slab, the others are appended)] (applied in order; blocks are 0-based, new blocks appended, the O-grid's inner block keeping the old index; a face is its hex8 number 1-6 or the direction it faces, "+x" … "-z"), "preview": true (return the topology without meshing), "faces": [{"block": b, "face": 1-6, "to": ["face names"]}] (targets set by hand; other outer block faces go to the smooth group of B-rep faces nearest them), "divisions": [{"block": b, "direction": 1|2|3, "count": n}], "smoothing": "elliptic" (default) | "elliptic-controlled" (Thomas–Middlecoff) | "laplace" | "none", "iterations": 100, "min_quality": 0.2 (smallest scaled Jacobian accepted), "merge": distance (join nodes of other bodies' meshes within it; default 1e-6 of the diagonal, 0 never), "tolerance": chordal tolerance of the tessellation} — replaces the body's mesh with an all-hex mesh of the block topology: outer block faces projected onto their B-rep faces, edges between faces with different targets onto the curve where those meet (tangent-plane iteration), vertices on three onto their point, nodes placed by transfinite interpolation and smoothed with the boundary held on the geometry. Elements are tagged with the body and boundary faces with their B-rep faces. Reports blocks, divisions, targets and quality before and after smoothing; a corner that flattens on a curved face is an error that says to use an O-grid or template there"#, + params: r#"{"body": "Body1", "size": 2.0 (element size, mm; default 1/25 of the body's diagonal), "template": "grid" (default) | "sphere" | "butterfly" | "blocks", "center": [x, y, z], "radius": r, "axis": [0, 0, 1], "levels": [z0, z1, …] (sphere and butterfly: default from the body's bounding box; butterfly levels are distances along the axis from the centre), "inset": 0.55 (inner square or cube corner at this fraction of the radius), "grid": {"planes": {"x": [...], "y": [...], "z": [...]} or "cells": [ni, nj, nk] (even planes over the bounding box, default [1, 1, 1]), "carve": [{"i": [a, b], "j": [c, d], "k": [e, f]}] (block index ranges, inclusive, removed)}, "blocks": {"vertices": [[x, y, z], ...], "blocks": [[v1..v8], ...] (0-based vertices, hex8 order)}, "operations": [{"op": "o_grid", "block": b, "open": [faces], "inset": 0.5 or [i1, i2, i3], "extrude": {"offset": [dx, dy, dz], "layers": 1, "faces": "all" (the open face's butterfly: a round branch) | "inner" (its inner face: a square stub)}} | {"op": "extrude", "faces": [[block, face], ...], "offset": [dx, dy, dz], "layers": 1} | {"op": "connect", "from": [block, face], "to": [block, face], "layers": 1} (blocks bridging two faces: closes a handle) | {"op": "remove", "blocks": [...]} (carve; later blocks renumber) | {"op": "merge", "distance": d} | {"op": "split", "block": b, "direction": 1|2|3 or "x"|"y"|"z", "at": [fractions]} (cuts the block and the whole sheet of blocks the cut runs through, so the topology stays conforming; fractions from the block's side-0 face, or from the low end along an axis; the block keeps its first slab, the others are appended)] (applied in order; blocks are 0-based, new blocks appended, the O-grid's inner block keeping the old index; a face is its hex8 number 1-6 or the direction it faces, "+x" … "-z"), "preview": true (return the topology without meshing), "faces": [{"block": b, "face": 1-6, "to": ["face names"]}] (targets set by hand; other outer block faces go to the smooth group of B-rep faces nearest them), "divisions": [{"block": b, "direction": 1|2|3, "count": n}], "smoothing": "elliptic" (default) | "elliptic-controlled" (Thomas–Middlecoff) | "laplace" | "none", "iterations": 100, "min_quality": 0.2 (smallest scaled Jacobian accepted), "max_arc_angle": 30 (degrees: the most a curved face may turn across one element; the divisions grow to keep it, 0 for no limit, at most 90), "merge": distance (join nodes of other bodies' meshes within it; default 1e-6 of the diagonal, 0 never), "tolerance": chordal tolerance of the tessellation} — replaces the body's mesh with an all-hex mesh of the block topology: outer block faces projected onto their B-rep faces, edges between faces with different targets onto the curve where those meet (tangent-plane iteration), vertices on three onto their point, nodes placed by transfinite interpolation and smoothed with the boundary held on the geometry. Elements are tagged with the body and boundary faces with their B-rep faces. Reports blocks, divisions, targets and quality before and after smoothing; a corner that flattens on a curved face is an error that says to use an O-grid or template there"#, form: Some(super::forms::HEX_PROJECT), kind: CommandKind::Model, requires: &[Requirement::Geometry], @@ -240,6 +240,7 @@ struct ProjectParams { smoothing: Option, iterations: Option, min_quality: Option, + max_arc_angle: Option, merge: Option, // Read from the raw parameters by the surface preparation (mesh.tet.generate's options). #[allow(dead_code)] @@ -509,7 +510,20 @@ fn hex_project(s: &mut Session, raw: &Value) -> Result { if !(-1.0..=1.0).contains(&min_quality) { return Err(EngineError::invalid("min_quality", "from -1 to 1")); } - let opts = ProjectOptions { size, divisions, faces, smoothing, iterations, max_elements: fieldcraft_model::MAX_ELEMENTS / 8, min_quality }; + let max_arc_angle = p.max_arc_angle.unwrap_or(fieldcraft_mesh::project::DEFAULT_MAX_ARC_ANGLE.to_degrees()); + if !(0.0..=90.0).contains(&max_arc_angle) { + return Err(EngineError::invalid("max_arc_angle", "from 0 to 90 degrees")); + } + let opts = ProjectOptions { + size, + max_arc_angle: max_arc_angle.to_radians(), + divisions, + faces, + smoothing, + iterations, + max_elements: fieldcraft_model::MAX_ELEMENTS / 8, + min_quality, + }; let merge = merge_distance(p.merge, diag)?; project_body(s, &p.body, b, &t, &opts, merge) } diff --git a/crates/mesh/src/project/mod.rs b/crates/mesh/src/project/mod.rs index b8c32e8..e380276 100644 --- a/crates/mesh/src/project/mod.rs +++ b/crates/mesh/src/project/mod.rs @@ -56,12 +56,22 @@ const VERTEX_SWEEPS: usize = 500; /// Rounds of arc-length redistribution of an edge's nodes. const EDGE_ROUNDS: usize = 4; +/// The default largest turn of the surface across one element: 30°. +pub const DEFAULT_MAX_ARC_ANGLE: f64 = std::f64::consts::PI / 6.0; + +/// Samples per side of an outer block face when measuring how much its target turns. +const CURVATURE_SAMPLES: usize = 32; + /// How to project. #[derive(Debug, Clone, PartialEq)] pub struct ProjectOptions { /// The element size: each class of parallel block edges gets `round(length / size)` /// divisions (at least one), from its longest projected edge. pub size: f64, + /// The largest angle (radians) the surface may turn through across one element on a curved + /// target face: a class takes enough divisions that its elements are no longer than the + /// face's radius of curvature times this (0 for no limit). + pub max_arc_angle: f64, /// Division counts set by hand: `(block, direction 0..3, count)`; the count applies to the /// whole class of parallel edges. pub divisions: Vec<(usize, usize, usize)>, @@ -80,6 +90,7 @@ impl Default for ProjectOptions { fn default() -> Self { ProjectOptions { size: 1.0, + max_arc_angle: DEFAULT_MAX_ARC_ANGLE, divisions: Vec::new(), faces: Vec::new(), smoothing: Smoothing::Elliptic, @@ -376,7 +387,62 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o } } } - let mut class_n: Vec = class_len.iter().map(|l| ((l / opts.size).round() as usize).clamp(1, crate::block::MAX_DIVISIONS)).collect(); + // Curved targets: across each outer block face, along each of its two directions, the + // largest normal curvature of the target in that direction (at a grid of points projected + // onto it, from the turn of the B-rep face's normal over a short step on that face) times + // the face's length there, over the largest turn an element may span, is the least count + // of that direction's class. + let probe = 1e-3 * opts.size; + let mut class_arc = vec![0usize; nclass]; + if opts.max_arc_angle > 0.0 { + for &(b, d, s) in &outer { + let Some(block) = t.blocks.get(b) else { continue }; + let Some(&tg) = face_target.get(&face_key(block, d, s)) else { continue }; + let Some(set) = targets.get(tg) else { continue }; + let c: [P; 8] = std::array::from_fn(|k| vpos.get(block[k] as usize).copied().unwrap_or_default()); + let (e1, e2) = other(d); + let m = CURVATURE_SAMPLES; + let mut grid: Vec> = Vec::with_capacity((m + 1) * (m + 1)); + for j in 0..=m { + for i in 0..=m { + let mut u = [0.0; 3]; + u[d] = s as f64; + u[e1] = i as f64 / m as f64; + u[e2] = j as f64 / m as f64; + grid.push(target::on_target(g, set, trilinear(&c, u))); + } + } + let at = |i: usize, j: usize| grid.get(j * (m + 1) + i).copied().flatten(); + for (dir, axis) in [(0, e1), (1, e2)] { + let (mut kappa, mut length) = (0.0f64, 0.0f64); + for line in 0..=m { + let mut run = 0.0; + for k in 0..m { + let (p, q) = if dir == 0 { (at(k, line), at(k + 1, line)) } else { (at(line, k), at(line, k + 1)) }; + let (Some((pa, na, fa)), Some((pb, _, _))) = (p, q) else { continue }; + let step = dist(pa, pb); + run += step; + if step > tol { + let x = add(pa, scale(crate::vec3::sub(pb, pa), probe / step)); + let Some(qx) = g.project(fa, x) else { continue }; + let Some(nx) = g.normal(fa, qx) else { continue }; + let h = dist(pa, qx); + if h > 0.5 * probe { + kappa = kappa.max(crate::vec3::dot(na, nx).clamp(-1.0, 1.0).acos() / h); + } + } + } + length = length.max(run); + } + let need = (kappa * length / opts.max_arc_angle).ceil(); + if let Some(slot) = classes.get(b).and_then(|cl| cl.get(axis)).and_then(|&k| class_arc.get_mut(k)) { + *slot = (*slot).max(need.min(crate::block::MAX_DIVISIONS as f64) as usize); + } + } + } + } + let mut class_n: Vec = + class_len.iter().zip(&class_arc).map(|(l, &arc)| ((l / opts.size).round() as usize).max(arc).clamp(1, crate::block::MAX_DIVISIONS)).collect(); for &(b, d, n) in &opts.divisions { let cls = classes.get(b).and_then(|c| c.get(d)).ok_or_else(|| Error::Invalid(format!("no direction {d} of block {b}")))?; if !(1..=crate::block::MAX_DIVISIONS).contains(&n) {