diff --git a/crates/engine/src/cmd/forms.rs b/crates/engine/src/cmd/forms.rs index 5abbb28..da1f41e 100644 --- a/crates/engine/src/cmd/forms.rs +++ b/crates/engine/src/cmd/forms.rs @@ -1418,7 +1418,12 @@ pub const NODE_RELAX: &[Field] = &[ labels("nodes", "Nodes", LabelOf::Node).help("Instead of a sphere: these nodes.").hint("none"), labels("elements", "Elements", LabelOf::Element).help("Instead of a sphere: the nodes of these elements.").hint("none"), references("bodies", "Bodies", C::Body).help("Only the nodes of these bodies.").hint("every body"), - int("iterations", "Iterations", 1, 100).help("Sweeps over the nodes.").def("5"), + choice("mode", "Mode", &[("guarded", "Guarded (local)"), ("strong", "Strong (optimise the region)")]) + .help("Guarded: each move keeps the worst element around the node from getting worse. Strong: optimise every node of the region together toward scaled Jacobian 1, nodes on CAD edges held.") + .def("\"guarded\""), + range("below", "Below", D::None, 0.0, 1.0).help("Strong: relax the elements below this scaled Jacobian (and the layers around them).").hint("every element of the region"), + int("layers", "Layers", 0, 20).help("Strong: rings of elements added around the elements below the threshold.").def("2"), + int("iterations", "Iterations", 1, 2000).help("Sweeps over the nodes (guarded: 1 to 100, default 5; strong: default 200).").hint("5 (strong: 200)"), range("strength", "Strength", D::None, 0.0, 1.0).help("How far each node moves toward its relaxed position, 0 to 1.").def("1"), choice("falloff", "Falloff", &[("smooth", "Smooth"), ("linear", "Linear"), ("none", "None")]) .help("How the strength falls from the centre to the radius, so the edit blends into the mesh.") diff --git a/crates/engine/src/cmd/relax.rs b/crates/engine/src/cmd/relax.rs index af7beed..45ea34d 100644 --- a/crates/engine/src/cmd/relax.rs +++ b/crates/engine/src/cmd/relax.rs @@ -28,7 +28,7 @@ use crate::{CommandKind, CommandSpec, EngineError, Result, Session}; pub const RELAX: CommandSpec = CommandSpec { id: "mesh.node.relax", label: "Relax Nodes", - params: r#"{"center": [x, y, z], "radius": mm (the sphere of influence) | "around_element": label (a sphere round that element; radius default twice its longest edge) | "faces": [{"body": "Base", "face": "F3"}] (persistent B-rep face names, as model.set.create takes them: every node on those faces and the nodes of the elements touching them; pick faces in the view) | "nodes": [labels] | "elements": [labels] (their nodes), "bodies"?: ["Body1"] (only their nodes), "iterations"?: 5 (1..=100), "strength"?: 1 (0..1], "falloff"?: "smooth" (default: weight 1 at the centre to 0 at the radius) | "linear" | "none", "field"?: "Name" (the relaxation field to keep; default Relax1, Relax2, …), "reapply"?: false (apply the field again after each remesh of these bodies)} — moves the nodes to improve the worst elements around them (smart Laplacian, then optimization of the worst element; a move is kept only when the worst element around the node does not get worse, and never inverts one). Nodes on a B-rep face slide in its UV parameters and stay on it, nodes on a B-rep edge slide along it, nodes on vertices stay; nodes on partition cut faces move freely (shared by both cells); mid-side nodes follow, on the geometry. Reports the nodes moved, the largest displacement, the region's smallest scaled Jacobian, largest edge ratio and worst elements before and after, the field kept, and the listed faces no mesh is tagged with (`faces_without_mesh`; an error when none is). One undo step"#, + params: r#"{"center": [x, y, z], "radius": mm (the sphere of influence) | "around_element": label (a sphere round that element; radius default twice its longest edge) | "faces": [{"body": "Base", "face": "F3"}] (persistent B-rep face names, as model.set.create takes them: every node on those faces and the nodes of the elements touching them; pick faces in the view) | "nodes": [labels] | "elements": [labels] (their nodes), "bodies"?: ["Body1"] (only their nodes; a partitioned body's cells count as it), "iterations"?: 5 (1..=100), "strength"?: 1 (0..1], "falloff"?: "smooth" (default: weight 1 at the centre to 0 at the radius) | "linear" | "none", "field"?: "Name" (the relaxation field to keep; default Relax1, Relax2, …), "reapply"?: false (apply the field again after each remesh of these bodies), "mode"?: "guarded" (default) | "strong", "below"?: 0.7 (strong: the elements below this scaled Jacobian), "layers"?: 2 (strong: rings of elements added around them, 0..=20)} — moves the nodes to improve the worst elements around them (smart Laplacian, then optimization of the worst element; a move is kept only when the worst element around the node does not get worse, and never inverts one). `"mode": "strong"` optimises the whole region instead (`bodies` alone is the whole body; `iterations` are sweeps, default 200, 1..=2000): nodes on B-rep edges and vertices are held, nodes on faces slide on them, interior nodes are free, and the node positions minimise Σ (1/SJ − 1)² over the elements (a barrier: nothing inverts), so single elements may get worse on the way; the best state passed through is kept, and nothing moves unless it is valid and its smallest scaled Jacobian is not below the start. Its reply adds the mean scaled Jacobian and the objective, min and mean per sweep (`history`). Nodes on a B-rep face slide in its UV parameters and stay on it, nodes on a B-rep edge slide along it, nodes on vertices stay; nodes on partition cut faces move freely (shared by both cells); mid-side nodes follow, on the geometry. Reports the nodes moved, the largest displacement, the region's smallest scaled Jacobian, largest edge ratio and worst elements before and after, the field kept, and the listed faces no mesh is tagged with (`faces_without_mesh`; an error when none is). One undo step"#, form: Some(super::forms::NODE_RELAX), kind: CommandKind::Model, requires: &[Requirement::Mesh], @@ -281,8 +281,11 @@ fn design_geometry(s: &Session, nodes: &BTreeSet) -> Result fn body_nodes(model: &Model, p: &Value) -> Result, BTreeSet)>> { let Some(bodies) = super::opt::>(p, "bodies")? else { return Ok(None) }; let mut nodes = BTreeSet::new(); + // A body cut into cells by partitions is meshed as its cells (`Base/a`, `Base/b/a`, …). + let meshed: BTreeSet = model.classification().entities().iter().map(|r| r.body.clone()).collect(); for b in &bodies { - let elements = model.body_elements(b); + let cell = format!("{b}/"); + let elements: Vec = meshed.iter().filter(|m| *m == b || m.starts_with(&cell)).flat_map(|m| model.body_elements(m)).collect(); if elements.is_empty() { return Err(EngineError::Params(format!("body `{}` has no mesh", b.chars().take(80).collect::()))); } @@ -497,7 +500,156 @@ pub fn sphere(s: &Session, params: &Value) -> Option<([f64; 3], f64)> { Some((c, r)) } +/// Elements below this scaled Jacobian, with `layers` rings of elements around them: their +/// nodes (only of `only` when given). +fn below_nodes(model: &Model, below: f64, layers: usize, only: Option<&BTreeSet>) -> Result> { + let mut nodes = BTreeSet::new(); + for e in model.elements().values() { + if !fieldcraft_mesh::quality::is_measured(e.kind) || only.is_some_and(|o| !e.nodes.iter().all(|n| o.contains(n))) { + continue; + } + let x: Vec<[f64; 3]> = e.nodes.iter().filter_map(|n| model.node(*n).map(|p| p.xyz)).collect(); + let q = fieldcraft_mesh::quality::element_quality(e.kind, &x).ok().and_then(|q| { + q.get(fieldcraft_mesh::quality::Metric::ScaledJacobian).or_else(|| q.get(fieldcraft_mesh::quality::Metric::MinAngle).map(|a| a / 60.0)) + }); + if q.is_none_or(|q| q < below) { + nodes.extend(e.nodes.iter().copied()); + } + } + for _ in 0..layers { + let mut grown = nodes.clone(); + for e in model.elements().values() { + if e.nodes.iter().any(|n| nodes.contains(n)) { + grown.extend(e.nodes.iter().copied()); + } + } + nodes = grown; + } + if let Some(only) = only { + nodes.retain(|n| only.contains(n)); + } + Ok(nodes) +} + +/// The most layers `layers` may add. +const MAX_LAYERS: usize = 20; +/// Default sweeps of the strong mode. +const STRONG_SWEEPS: usize = 200; + +/// `mesh.node.relax` with `"mode": "strong"`: optimise the region globally ([`relax::optimize`]). +fn relax_strong(s: &mut Session, p: &Value) -> Result { + let only = body_nodes(&s.model, p)?; + let below = match p.get("below") { + None => None, + Some(_) => { + let b = super::optional_number(p, "below", 0.0)?; + if !(b > 0.0 && b <= 1.0) { + return Err(EngineError::param("below", ProblemKind::NotInRange { min: 0.0, max: 1.0 })); + } + Some(b) + } + }; + let layers = match p.get("layers") { + None => 2, + Some(v) => v + .as_u64() + .and_then(|l| usize::try_from(l).ok()) + .filter(|l| *l <= MAX_LAYERS) + .ok_or_else(|| EngineError::param("layers", ProblemKind::NotInRange { min: 0.0, max: MAX_LAYERS as f64 }))?, + }; + let listed = ["center", "radius", "around_element", "faces", "nodes", "elements"].iter().any(|k| p.get(*k).is_some()); + let mut nodes: BTreeSet = if listed { + weights(&s.model, p)?.0.into_keys().collect() + } else if let Some((_, n)) = &only { + n.clone() + } else if below.is_none() { + return Err(EngineError::Params( + "say which nodes to relax: give `bodies`, `below`, `center` and `radius`, `around_element`, `faces`, `nodes` or `elements`".into(), + )); + } else { + s.model.nodes().keys().copied().collect() + }; + if let Some(b) = below { + let low = below_nodes(&s.model, b, layers, only.as_ref().map(|(_, n)| n))?; + nodes.retain(|n| low.contains(n)); + } + let nodes: BTreeSet = relax::uniform_weights(&s.model, &nodes, 1.0)?.into_keys().collect(); + if nodes.is_empty() { + return Err(EngineError::Params("no mesh nodes in the region to relax".into())); + } + let sweeps = match p.get("iterations") { + None => STRONG_SWEEPS, + Some(v) => v + .as_u64() + .and_then(|i| usize::try_from(i).ok()) + .filter(|i| (1..=relax::MAX_SWEEPS).contains(i)) + .ok_or_else(|| EngineError::param("iterations", ProblemKind::NotInRange { min: 1.0, max: relax::MAX_SWEEPS as f64 }))?, + }; + let g = design_geometry(s, &touched(&s.model, &nodes.iter().map(|&n| (n, 1.0)).collect()))?; + let r = relax::optimize(&s.model, &nodes, sweeps, &g)?; + let field_name = match super::opt::(p, "field")? { + Some(n) => n, + None => (1..).map(|i| format!("Relax{i}")).find(|n| s.model.relaxation(n).is_none()).unwrap_or_else(|| "Relax".into()), + }; + let bodies = only.map(|(b, _)| b).unwrap_or_default(); + let reapply = super::opt::(p, "reapply")?.unwrap_or(false); + let as_relaxed = Relaxed { positions: r.positions.clone(), ..Relaxed::default() }; + if !r.positions.is_empty() { + let cells = relax::field_cells(&s.model, &as_relaxed); + s.model.set_relaxation(&field_name, Relaxation { bodies, cells, reapply })?; + } + for (n, x) in &r.positions { + s.model.set_node(*n, *x)?; + } + let summary = format!( + "strong: moved {} node{} (largest {:.3e}) over {} sweep{}; scaled Jacobian min {:.3} → {:.3}, mean {:.3} → {:.3}{}", + r.moved, + if r.moved == 1 { "" } else { "s" }, + r.max_displacement, + r.history.len().saturating_sub(1), + if r.history.len() == 2 { "" } else { "s" }, + r.before.min_quality, + r.after.min_quality, + r.mean_before, + r.mean_after, + if r.positions.is_empty() { "; nothing better found, nothing moved".to_string() } else { format!("; kept as field {field_name}") } + ); + let history: Vec = r + .history + .iter() + .map(|h| json!({"phase": match h.phase { relax::Phase::Smooth => "smooth", relax::Phase::MinMax => "min_max" }, "objective": h.objective, "min": h.min_quality, "mean": h.mean_quality, "invalid": h.invalid})) + .collect(); + let mut before = quality_json(&r.before); + let mut after = quality_json(&r.after); + if let (Some(b), Some(a)) = (before.as_object_mut(), after.as_object_mut()) { + b.insert("mean_scaled_jacobian".into(), json!(r.mean_before)); + a.insert("mean_scaled_jacobian".into(), json!(r.mean_after)); + } + Ok(json!({ + "mode": "strong", + "moved": r.moved, + "nodes_changed": r.positions.len(), + "max_displacement": r.max_displacement, + "iterations": r.history.len().saturating_sub(1), + "kept_sweep": r.kept, + "history": history, + "before": before, + "after": after, + "field": (!r.positions.is_empty()).then_some(field_name), + "faces_without_mesh": unmeshed_faces(&s.model, p), + "summary": summary, + })) +} + fn relax_nodes(s: &mut Session, p: &Value) -> Result { + match super::opt::(p, "mode")?.as_deref() { + None | Some("guarded") => {} + Some("strong") => return relax_strong(s, p), + Some(_) => return Err(EngineError::param("mode", ProblemKind::NotAChoice)), + } + if p.get("below").is_some() || p.get("layers").is_some() { + return Err(EngineError::Params("`below` and `layers` go with \"mode\": \"strong\"".into())); + } let (w, sphere) = weights(&s.model, p)?; if w.is_empty() { return Err(EngineError::Params("no mesh nodes in the region to relax".into())); diff --git a/crates/engine/src/cmd/relax/tests.rs b/crates/engine/src/cmd/relax/tests.rs index e5b21cb..4b41002 100644 --- a/crates/engine/src/cmd/relax/tests.rs +++ b/crates/engine/src/cmd/relax/tests.rs @@ -260,3 +260,100 @@ fn relaxing_fillet_faces_improves_them() { s.execute("mesh.node.relax", &json!({"elements": [e]})).unwrap(); assert!(s.execute("mesh.node.relax", &json!({"set": "All"})).is_err()); } + +/// Whether a node of the 20 mm box lies on two of its planes (an edge or a vertex of the box, +/// or a tangent line of its fillet). +fn on_box_edge(p: [f64; 3]) -> bool { + p.iter().filter(|c| c.abs() < 1e-6 || (**c - 20.0).abs() < 1e-6).count() >= 2 +} + +/// Strong mode on the whole filleted box: the smallest and the mean scaled Jacobian do not drop, +/// nothing inverts, the fillet's nodes stay on its cylinder, the plane faces' nodes on their +/// planes, and the nodes on the box's edges and the fillet's tangent lines do not move. The +/// reply carries the objective, min and mean per sweep. +#[test] +fn strong_relaxation_of_a_filleted_box_keeps_faces_and_edges() { + let (_d, mut s, _) = filleted_box(); + let before = s.model.clone(); + let r = s.execute("mesh.node.relax", &json!({"mode": "strong", "bodies": ["B"], "iterations": 40})).unwrap(); + let (b, a) = (&r["before"], &r["after"]); + assert!(a["min_scaled_jacobian"].as_f64().unwrap() > b["min_scaled_jacobian"].as_f64().unwrap(), "{r}"); + assert!(a["mean_scaled_jacobian"].as_f64().unwrap() > b["mean_scaled_jacobian"].as_f64().unwrap(), "{r}"); + assert_eq!(a["invalid"], 0, "{r}"); + assert_eq!(invalid(&mut s), 0); + let history = r["history"].as_array().unwrap(); + assert!(history.len() > 2 && history.iter().all(|h| h["objective"].is_number() && h["min"].is_number() && h["mean"].is_number()), "{r}"); + for (n, p) in s.model.nodes() { + let (old, new) = (before.node(*n).unwrap().xyz, p.xyz); + if on_box_edge(old) { + assert_eq!(old, new, "node {n} on an edge moved"); + } + if old[0] > 18.0 && old[2] > 18.0 && (fillet_radius(old) - 2.0).abs() < 1e-6 { + assert!((fillet_radius(new) - 2.0).abs() < 1e-6, "node {n} left the fillet"); + } + for k in 0..3 { + for plane in [0.0, 20.0] { + if (old[k] - plane).abs() < 1e-9 { + assert!((new[k] - plane).abs() < 1e-6, "node {n} left the plane {k} = {plane}"); + } + } + } + } +} + +/// Strong mode on a distorted hex plate: over the guarded mode's sphere it lifts the smallest +/// scaled Jacobian as far as the guarded mode (to the elements on the sphere's rim, which +/// neither moves) (with nodes on edges held, where the guarded mode slides +/// them); by quality threshold (the elements below 0.5 and two layers round them) it untangles +/// every element. Nothing inverts, nodes on the hole stay on its cylinder and the caps' on their +/// planes, nodes on the hole's edges do not move. +#[test] +fn strong_relaxation_untangles_a_distorted_plate() { + let mut s = plate(); + s.execute("mesh.hex.generate", &json!({"size": 2.0})).unwrap(); + perturb(&mut s, ON_HOLE, 6.0, 0.35); + let perturbed = s.model.clone(); + let sphere = json!({"center": ON_HOLE, "radius": 7.0, "iterations": 20}); + let guarded = s.execute("mesh.node.relax", &sphere).unwrap(); + let min = |r: &Value, at: &str| r[at]["min_scaled_jacobian"].as_f64().unwrap(); + for strong in [json!({"mode": "strong", "center": ON_HOLE, "radius": 7.0}), json!({"mode": "strong", "below": 0.5, "layers": 2})] { + s.model = perturbed.clone(); + let r = s.execute("mesh.node.relax", &strong).unwrap(); + assert!(min(&r, "after") > 0.5 && min(&r, "before") < 0.0, "{strong}: {r}"); + if strong.get("center").is_some() { + assert!(min(&r, "after") >= min(&guarded, "after") - 1e-9, "strong {r}\nguarded {guarded}"); + } + assert_eq!(r["after"]["invalid"], 0, "{r}"); + assert_eq!(invalid(&mut s), 0); + check_on_geometry(&perturbed, &s.model); + for (n, p) in s.model.nodes() { + let old = perturbed.node(*n).unwrap().xyz; + if on_faces(old).len() >= 2 { + assert_eq!(old, p.xyz, "node {n} on the hole's edge moved"); + } + } + } +} + +#[test] +fn strong_mode_parameters_are_checked() { + let mut s = plate(); + s.execute("mesh.hex.generate", &json!({"size": 4.0})).unwrap(); + for p in [ + json!({"mode": "strong"}), + json!({"mode": "wild", "nodes": [1]}), + json!({"mode": "strong", "below": 0}), + json!({"mode": "strong", "below": 1.5}), + json!({"mode": "strong", "below": 0.5, "layers": 21}), + json!({"mode": "strong", "below": 0.5, "layers": -1}), + json!({"mode": "strong", "bodies": ["nope"]}), + json!({"mode": "strong", "nodes": [1], "iterations": 0}), + json!({"mode": "strong", "nodes": [1], "iterations": 2001}), + json!({"below": 0.5}), + json!({"nodes": [1], "layers": 1}), + ] { + assert!(s.execute("mesh.node.relax", &p).is_err(), "{p}"); + } + // An undistorted mesh above the threshold: nothing to relax. + assert!(s.execute("mesh.node.relax", &json!({"mode": "strong", "below": 0.01})).is_err()); +} diff --git a/docs/mcp.md b/docs/mcp.md index 7e81d84..5404ebd 100644 --- a/docs/mcp.md +++ b/docs/mcp.md @@ -210,6 +210,11 @@ When `mesh.quality` shows a handful of spikes in an otherwise good mesh: vertex nodes stay; nodes on a partition cut face move freely (both cells share them). The reply gives the nodes moved, the largest displacement, and the region's smallest scaled Jacobian, largest edge ratio and worst elements before and after. One undo step. + For a whole region at once, `mesh.node.relax {"mode": "strong", "bodies": ["Base"]}` (or + `"below": 0.7, "layers": 2` for the elements below 0.7 and two rings round them) optimises + every node together toward scaled Jacobian 1, edge and vertex nodes held, single elements + allowed to get worse on the way; it keeps the best state, never one with a lower smallest + scaled Jacobian, and replies the objective, min and mean per sweep (`history`). 2. The relaxation is kept as a field of the model (`field`, default `Relax1`, …; saved with it). After a remesh, `mesh.node.relax.apply {"field": "Relax1"}` puts the relaxed shape on the new mesh (each node takes the displacement interpolated in the old cell it lies in, on its face, diff --git a/docs/meshing/literature.md b/docs/meshing/literature.md index af7a988..d50f1eb 100644 --- a/docs/meshing/literature.md +++ b/docs/meshing/literature.md @@ -970,6 +970,16 @@ input points), Steiner points kept half the local spacing from the boundary. mesh takes the displacement at each node's parametric coordinates in the old cell holding it (barycentric or trilinear shape functions), slides it onto its face or edge, undoes any move that would invert an element, and smooths once (`mesh.node.relax.apply`, or `reapply` after every remesh). +- **Strong mode** (`"mode": "strong"`, `fieldcraft-mesh::relax::optimize`). An optimisation of every corner node + of the region at once: interior nodes in space, face nodes in their face's parameters, edge and vertex nodes held. + It minimises Σ φ(SJ) over the elements by nonlinear Gauss–Seidel (per node: the neighbours' centroid, then a + central-difference gradient step with backtracking, each kept when it lowers the node's part of the sum), first + with φ = (1/SJ − 1)², then with φ = SJ⁻¹⁶ (a smooth stand-in for the smallest SJ); an invalid element costs + far more than any valid one (a barrier, and an untangling objective while elements are inverted; P. Knupp, + *Hexahedral and tetrahedral mesh untangling*, Eng. Comput. 17, 2001; L. Freitag, P. Plassmann, IJNME 49, 2000). + A term ½ ln²(V/V₀) keeps element sizes, and no element may exceed the region's largest edge ratio before: the + scaled Jacobian alone lets elements flatten against the boundary. The best state passed through is kept, and + only if it is valid and its smallest SJ is not below the start. ---