forked from 1b-agent/fieldcraft-upstream
engine, ui: mesh.node.relax and relaxation fields, Relax Nodes tool with sphere and preview (#398)
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo
This commit is contained in:
jyvikoandClaude Opus 5.5 authored and Irakli Kourtis committed
1 parent
75f9e5ae6e
commit
94fdbf0fca
19 files changed
+1038
-9
No files matched your search
@@ -1264,6 +1264,29 @@ pub const NODE_MERGE: &[Field] = &[
|
||||
|
||||
pub const CHECK_DUPLICATES: &[Field] = &[MERGE_TOLERANCE];
|
||||
|
||||
pub const NODE_RELAX: &[Field] = &[
|
||||
xyz("center", "Centre", D::Length).help("The centre of the sphere of influence; click the mesh to place it.").hint("click the mesh"),
|
||||
positive("radius", "Radius", D::Length).help("Nodes closer to the centre than this move; Alt + scroll over the view changes it."),
|
||||
int("around_element", "Around element", 1, 99_999_999)
|
||||
.help("Instead of a centre: a sphere round this element (radius default twice its longest edge).")
|
||||
.hint("none"),
|
||||
labels("nodes", "Nodes", LabelOf::Node).help("Instead of a sphere: these nodes.").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"),
|
||||
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.")
|
||||
.def("\"smooth\""),
|
||||
text("field", "Keep as").help("The relaxation field to keep, to apply again after a remesh.").hint("Relax1, Relax2, …"),
|
||||
flag("reapply", "Reapply after remesh").help("Apply the field again by itself after each remesh.").def("false"),
|
||||
];
|
||||
|
||||
pub const NODE_RELAX_APPLY: &[Field] = &[
|
||||
text("field", "Field").req().help("The relaxation field kept by Relax Nodes."),
|
||||
references("bodies", "Bodies", C::Body).help("Only the nodes of these bodies.").hint("the field's bodies"),
|
||||
range("polish", "Polish", D::None, 0.0, 1.0).help("Strength of one guarded smoothing sweep after the field is applied; 0 for none.").def("0.3"),
|
||||
];
|
||||
|
||||
pub const RENUMBER: &[Field] = &[
|
||||
choice("what", "Renumber", &[("both", "Nodes and elements"), ("nodes", "Nodes"), ("elements", "Elements")])
|
||||
.help("Which labels change.")
|
||||
|
||||
@@ -18,10 +18,12 @@ pub use datum::edit_params as datum_edit_params;
|
||||
pub use datum::preview as datum_preview;
|
||||
pub use mesh::analytic_surfaces;
|
||||
pub use nodes::merge_preview as node_merge_preview;
|
||||
pub use relax::{Preview as RelaxPreview, preview as relax_preview, sphere as relax_sphere};
|
||||
mod mirror;
|
||||
mod model;
|
||||
pub(crate) mod nodes;
|
||||
pub(crate) mod octree;
|
||||
pub(crate) mod relax;
|
||||
pub mod replace;
|
||||
mod results;
|
||||
pub(crate) mod seeds;
|
||||
@@ -86,6 +88,8 @@ pub static ALL: &[CommandSpec] = &[
|
||||
nodes::MERGE,
|
||||
nodes::DUPLICATES,
|
||||
nodes::RENUMBER,
|
||||
relax::RELAX,
|
||||
relax::APPLY,
|
||||
tags::TAG,
|
||||
tags::TAG_CLEAR,
|
||||
model::UNITS,
|
||||
|
||||
@@ -0,0 +1,540 @@
|
||||
//! Node relaxation (#398): `mesh.node.relax` moves the nodes inside a sphere of influence (or
|
||||
//! listed nodes, elements or a set) to fix the worst elements there, nodes on the geometry
|
||||
//! staying on it ([`fieldcraft_mesh::relax`]), and keeps what it did as a relaxation field of the
|
||||
//! model; `mesh.node.relax.apply` does a field again on the current mesh (after a remesh), and a
|
||||
//! field marked `reapply` is applied by itself after every remesh of its bodies.
|
||||
//!
|
||||
//! **Where a node may go** comes from the mesh's geometry tags: a node on one B-rep face slides
|
||||
//! in that face's parameters, a node where two faces meet slides along their edge, a node on a
|
||||
//! geometry vertex or on three or more faces stays. Faces of the same surface (the two halves a
|
||||
//! kernel splits a cylinder into, or one face as cut into cells) count once. Partition cut faces
|
||||
//! (`partition:<name>:a` and `:b`) are not geometry: a node on one moves freely, and the cells on
|
||||
//! both sides share it, so the mesh stays conforming across the cut. A node on the skin whose
|
||||
//! geometry is unknown (an imported mesh) stays.
|
||||
|
||||
use std::collections::{BTreeMap, BTreeSet};
|
||||
|
||||
use fieldcraft_mesh::relax::{self, Constraints, Falloff, Freedom, RegionQuality, Relaxed};
|
||||
use fieldcraft_mesh::surface::{BoundarySurface, SurfaceOptions};
|
||||
use fieldcraft_mesh::tet::geometry::{AnalyticSurface, Tessellation};
|
||||
use fieldcraft_model::geometry::Entity;
|
||||
use fieldcraft_model::{Model, Relaxation};
|
||||
use serde_json::{Value, json};
|
||||
|
||||
use crate::requirement::Requirement;
|
||||
use crate::schema::ProblemKind;
|
||||
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) | "nodes": [labels] | "elements": [labels] | "set": "Name", "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, and the field kept. One undo step"#,
|
||||
form: Some(super::forms::NODE_RELAX),
|
||||
kind: CommandKind::Model,
|
||||
requires: &[Requirement::Mesh],
|
||||
run: relax_nodes,
|
||||
};
|
||||
|
||||
pub const APPLY: CommandSpec = CommandSpec {
|
||||
id: "mesh.node.relax.apply",
|
||||
label: "Apply Relaxation",
|
||||
params: r#"{"field": "Name", "bodies"?: ["Body1"] (default the field's), "polish"?: 0.3 (strength of one guarded smoothing sweep after; 0 for none)} — applies a relaxation field kept by mesh.node.relax to the current mesh (typically after a remesh): each node inside the relaxed region takes the displacement interpolated from the field at its parametric coordinates in the old cells, on its face or edge as mesh.node.relax moves it; a move that would invert an element is undone; then one guarded smoothing sweep. Reports the nodes moved and the largest displacement"#,
|
||||
form: Some(super::forms::NODE_RELAX_APPLY),
|
||||
kind: CommandKind::Model,
|
||||
requires: &[Requirement::Mesh],
|
||||
run: apply,
|
||||
};
|
||||
|
||||
/// The commands after which `reapply` fields are applied again.
|
||||
const REMESH: &[&str] = &[
|
||||
"mesh.tet.generate",
|
||||
"mesh.hex.generate",
|
||||
"mesh.hex.sweep",
|
||||
"mesh.hex.mapped",
|
||||
"mesh.hex.blocks",
|
||||
"mesh.hex.octree",
|
||||
"mesh.cells",
|
||||
"mesh.block",
|
||||
];
|
||||
|
||||
const DEFAULT_ITERATIONS: usize = 5;
|
||||
|
||||
/// The region's quality as JSON.
|
||||
fn quality_json(q: &RegionQuality) -> Value {
|
||||
json!({
|
||||
"elements": q.elements,
|
||||
"min_scaled_jacobian": q.min_quality,
|
||||
"max_edge_ratio": if q.max_edge_ratio.is_finite() { json!(q.max_edge_ratio) } else { Value::Null },
|
||||
"invalid": q.invalid,
|
||||
"worst": q.worst.iter().map(|(e, v)| json!({"element": e, "quality": v})).collect::<Vec<_>>(),
|
||||
})
|
||||
}
|
||||
|
||||
/// Faces a node lies on (geometry references), from the element faces, shells and vertex tags.
|
||||
struct NodeTags {
|
||||
faces: BTreeMap<u32, BTreeSet<(String, String)>>,
|
||||
vertex: BTreeSet<u32>,
|
||||
/// Nodes on a partition cut face.
|
||||
cut: BTreeSet<u32>,
|
||||
}
|
||||
|
||||
fn node_tags(model: &Model, nodes: &BTreeSet<u32>) -> NodeTags {
|
||||
let c = model.classification();
|
||||
let mut t = NodeTags { faces: BTreeMap::new(), vertex: BTreeSet::new(), cut: BTreeSet::new() };
|
||||
let note = |n: u32, tag: u32, t: &mut NodeTags| {
|
||||
let Some(g) = c.entity(tag) else { return };
|
||||
if let Entity::Face(name) = &g.entity {
|
||||
if name.starts_with("partition:") {
|
||||
t.cut.insert(n);
|
||||
} else {
|
||||
t.faces.entry(n).or_default().insert((g.body.clone(), name.clone()));
|
||||
}
|
||||
}
|
||||
};
|
||||
for (&label, tags) in c.elements() {
|
||||
let Some(e) = model.element(label) else { continue };
|
||||
if !e.nodes.iter().any(|n| nodes.contains(n)) {
|
||||
continue;
|
||||
}
|
||||
if e.kind.dimension() == 2
|
||||
&& let Some(r) = tags.region
|
||||
{
|
||||
for &n in &e.nodes {
|
||||
note(n, r, &mut t);
|
||||
}
|
||||
}
|
||||
for (&face, &tag) in &tags.faces {
|
||||
let Some(local) = e.kind.faces().get(usize::from(face).wrapping_sub(1)) else { continue };
|
||||
for &i in local.iter() {
|
||||
if let Some(&n) = e.nodes.get(i) {
|
||||
note(n, tag, &mut t);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
for (&n, &tag) in c.nodes() {
|
||||
if nodes.contains(&n) && c.entity(tag).is_some_and(|g| matches!(g.entity, Entity::Vertex(_))) {
|
||||
t.vertex.insert(n);
|
||||
}
|
||||
}
|
||||
t
|
||||
}
|
||||
|
||||
/// Whether two exact surfaces are the same surface.
|
||||
fn same_surface(a: &AnalyticSurface, b: &AnalyticSurface, tol: f64) -> bool {
|
||||
use AnalyticSurface as S;
|
||||
let sub = |x: [f64; 3], y: [f64; 3]| [x[0] - y[0], x[1] - y[1], x[2] - y[2]];
|
||||
let dot = |x: [f64; 3], y: [f64; 3]| x[0] * y[0] + x[1] * y[1] + x[2] * y[2];
|
||||
let len = |x: [f64; 3]| dot(x, x).sqrt();
|
||||
let unit = |x: [f64; 3]| {
|
||||
let l = len(x);
|
||||
if l > 0.0 { [x[0] / l, x[1] / l, x[2] / l] } else { x }
|
||||
};
|
||||
let parallel = |x: [f64; 3], y: [f64; 3]| (dot(unit(x), unit(y)).abs() - 1.0).abs() < 1e-9;
|
||||
let on_axis = |p: [f64; 3], o: [f64; 3], a: [f64; 3]| {
|
||||
let d = sub(p, o);
|
||||
let a = unit(a);
|
||||
let along = dot(d, a);
|
||||
len([d[0] - a[0] * along, d[1] - a[1] * along, d[2] - a[2] * along]) <= tol
|
||||
};
|
||||
match (*a, *b) {
|
||||
(S::Plane { normal: n1, point: p1 }, S::Plane { normal: n2, point: p2 }) => parallel(n1, n2) && dot(sub(p2, p1), unit(n1)).abs() <= tol,
|
||||
(S::Cylinder { origin: o1, axis: a1, radius: r1 }, S::Cylinder { origin: o2, axis: a2, radius: r2 }) => {
|
||||
(r1 - r2).abs() <= tol && parallel(a1, a2) && on_axis(o2, o1, a1)
|
||||
}
|
||||
(S::Sphere { center: c1, radius: r1 }, S::Sphere { center: c2, radius: r2 }) => len(sub(c1, c2)) <= tol && (r1 - r2).abs() <= tol,
|
||||
(S::Cone { apex: p1, axis: a1, half_angle: h1 }, S::Cone { apex: p2, axis: a2, half_angle: h2 }) => {
|
||||
len(sub(p1, p2)) <= tol && dot(unit(a1), unit(a2)) > 1.0 - 1e-9 && (h1 - h2).abs() <= 1e-9
|
||||
}
|
||||
(S::Torus { center: c1, axis: a1, major: m1, minor: n1 }, S::Torus { center: c2, axis: a2, major: m2, minor: n2 }) => {
|
||||
len(sub(c1, c2)) <= tol && parallel(a1, a2) && (m1 - m2).abs() <= tol && (n1 - n2).abs() <= tol
|
||||
}
|
||||
_ => false,
|
||||
}
|
||||
}
|
||||
|
||||
/// The design's geometry as relaxation constraints.
|
||||
struct DesignGeometry {
|
||||
/// Per body: its tessellation with exact surfaces.
|
||||
bodies: Vec<Tessellation>,
|
||||
/// Face ids: (body, B-rep face).
|
||||
faces: Vec<(usize, u32)>,
|
||||
/// Curve ids: (body, the two B-rep faces).
|
||||
curves: Vec<(usize, [u32; 2])>,
|
||||
freedom: BTreeMap<u32, Freedom>,
|
||||
}
|
||||
|
||||
impl Constraints for DesignGeometry {
|
||||
fn freedom(&self, node: u32) -> Option<Freedom> {
|
||||
self.freedom.get(&node).copied()
|
||||
}
|
||||
fn slide_face(&self, face: u32, p: [f64; 3], d: [f64; 3]) -> Option<[f64; 3]> {
|
||||
let &(b, f) = self.faces.get(face as usize)?;
|
||||
self.bodies.get(b)?.slide_face(f, p, d)
|
||||
}
|
||||
fn slide_curve(&self, curve: u32, p: [f64; 3], d: [f64; 3]) -> Option<[f64; 3]> {
|
||||
let &(b, f) = self.curves.get(curve as usize)?;
|
||||
self.bodies.get(b)?.slide_curve(f, p, d)
|
||||
}
|
||||
}
|
||||
|
||||
/// The constraints of `nodes` from the mesh's tags and the open design.
|
||||
fn design_geometry(s: &Session, nodes: &BTreeSet<u32>) -> Result<DesignGeometry> {
|
||||
let tags = node_tags(&s.model, nodes);
|
||||
let mut g = DesignGeometry { bodies: Vec::new(), faces: Vec::new(), curves: Vec::new(), freedom: BTreeMap::new() };
|
||||
// Body name → (index into `g.bodies`, its exact surfaces), built when first needed.
|
||||
let mut loaded: BTreeMap<String, Option<(usize, Vec<Option<AnalyticSurface>>)>> = BTreeMap::new();
|
||||
let mut load = |name: &str, g: &mut DesignGeometry| -> Result<Option<(usize, Vec<Option<AnalyticSurface>>)>> {
|
||||
if let Some(l) = loaded.get(name) {
|
||||
return Ok(l.clone());
|
||||
}
|
||||
let shape = s.design.as_deref().and_then(|d| d.body(name));
|
||||
let entry = match shape {
|
||||
None => None,
|
||||
Some(shape) => {
|
||||
let tol = (shape.size()? * 1e-3).max(1e-6);
|
||||
let mesh = shape.tessellate(tol)?;
|
||||
let surface = BoundarySurface::new(
|
||||
&mesh.positions,
|
||||
&mesh.triangles,
|
||||
&mesh.tri_face,
|
||||
&SurfaceOptions { weld_distance: mesh.weld_distance, ..SurfaceOptions::default() },
|
||||
)?;
|
||||
let analytic = super::mesh::analytic_surfaces(shape);
|
||||
g.bodies.push(Tessellation::new(&surface, analytic.clone()));
|
||||
Some((g.bodies.len() - 1, analytic))
|
||||
}
|
||||
};
|
||||
loaded.insert(name.to_string(), entry.clone());
|
||||
Ok(entry)
|
||||
};
|
||||
let tol = 1e-6 * s.model.nodes().values().fold(0.0f64, |m, n| m.max(n.xyz.iter().fold(0.0f64, |a, x| a.max(x.abs())))).max(1.0);
|
||||
let mut face_ids: BTreeMap<(usize, u32), u32> = BTreeMap::new();
|
||||
let mut curve_ids: BTreeMap<(usize, [u32; 2]), u32> = BTreeMap::new();
|
||||
for &n in nodes {
|
||||
if tags.vertex.contains(&n) {
|
||||
g.freedom.insert(n, Freedom::Fixed);
|
||||
continue;
|
||||
}
|
||||
let Some(refs) = tags.faces.get(&n) else {
|
||||
if tags.cut.contains(&n) {
|
||||
g.freedom.insert(n, Freedom::Free);
|
||||
}
|
||||
continue;
|
||||
};
|
||||
// Resolve each face; group those on one surface.
|
||||
let mut groups: Vec<Vec<(usize, u32, Option<AnalyticSurface>)>> = Vec::new();
|
||||
let mut unknown = false;
|
||||
for (body, face) in refs {
|
||||
let Some((b, analytic)) = load(body, &mut g)? else {
|
||||
unknown = true;
|
||||
continue;
|
||||
};
|
||||
let Some(f) = s.design.as_deref().and_then(|d| d.body(body)).and_then(|sh| sh.face_index(face)) else {
|
||||
unknown = true;
|
||||
continue;
|
||||
};
|
||||
let surface = analytic.get(f as usize).copied().flatten();
|
||||
let group = groups.iter_mut().find(|grp| {
|
||||
grp.iter().any(|(gb, gf, gs)| (*gb == b && *gf == f) || matches!((gs, surface), (Some(x), Some(y)) if same_surface(x, &y, tol)))
|
||||
});
|
||||
match group {
|
||||
Some(grp) => grp.push((b, f, surface)),
|
||||
None => groups.push(vec![(b, f, surface)]),
|
||||
}
|
||||
}
|
||||
let freedom = match groups.as_slice() {
|
||||
_ if unknown => Freedom::Fixed,
|
||||
[one] => {
|
||||
let Some(&(b, f, _)) = one.first() else { continue };
|
||||
let next = g.faces.len() as u32;
|
||||
let id = *face_ids.entry((b, f)).or_insert(next);
|
||||
if id == next {
|
||||
g.faces.push((b, f));
|
||||
}
|
||||
Freedom::Face(id)
|
||||
}
|
||||
[x, y] => {
|
||||
// A body with a face of each surface carries their edge.
|
||||
let pair = x.iter().find_map(|&(bx, fx, _)| y.iter().find(|(by, _, _)| *by == bx).map(|&(_, fy, _)| (bx, [fx.min(fy), fx.max(fy)])));
|
||||
match pair {
|
||||
Some(key) => {
|
||||
let next = g.curves.len() as u32;
|
||||
let id = *curve_ids.entry(key).or_insert(next);
|
||||
if id == next {
|
||||
g.curves.push(key);
|
||||
}
|
||||
Freedom::Curve(id)
|
||||
}
|
||||
None => Freedom::Fixed,
|
||||
}
|
||||
}
|
||||
_ => Freedom::Fixed,
|
||||
};
|
||||
g.freedom.insert(n, freedom);
|
||||
}
|
||||
Ok(g)
|
||||
}
|
||||
|
||||
/// The nodes of `bodies`' elements, when given.
|
||||
fn body_nodes(model: &Model, p: &Value) -> Result<Option<(Vec<String>, BTreeSet<u32>)>> {
|
||||
let Some(bodies) = super::opt::<Vec<String>>(p, "bodies")? else { return Ok(None) };
|
||||
let mut nodes = BTreeSet::new();
|
||||
for b in &bodies {
|
||||
let elements = model.body_elements(b);
|
||||
if elements.is_empty() {
|
||||
return Err(EngineError::Params(format!("body `{}` has no mesh", b.chars().take(80).collect::<String>())));
|
||||
}
|
||||
for e in elements {
|
||||
if let Some(el) = model.element(e) {
|
||||
nodes.extend(el.nodes.iter().copied());
|
||||
}
|
||||
}
|
||||
}
|
||||
Ok(Some((bodies, nodes)))
|
||||
}
|
||||
|
||||
fn point(p: &Value, key: &str) -> Result<Option<[f64; 3]>> {
|
||||
match super::opt::<[f64; 3]>(p, key)? {
|
||||
Some(c) if c.iter().all(|x| x.is_finite()) => Ok(Some(c)),
|
||||
Some(_) => Err(EngineError::param(key, ProblemKind::NotANumber)),
|
||||
None => Ok(None),
|
||||
}
|
||||
}
|
||||
|
||||
/// The sphere around element `label`: its centroid, and twice its longest corner edge.
|
||||
fn around_element(model: &Model, label: u32) -> Result<([f64; 3], f64)> {
|
||||
let e = model.element(label).ok_or_else(|| EngineError::Params(format!("no element {label}")))?;
|
||||
let x: Vec<[f64; 3]> = e.nodes.iter().take(e.kind.corner_count()).filter_map(|n| model.node(*n).map(|p| p.xyz)).collect();
|
||||
if x.is_empty() {
|
||||
return Err(EngineError::Params(format!("element {label} has no nodes")));
|
||||
}
|
||||
let c = x.iter().fold([0.0; 3], |s, p| [s[0] + p[0], s[1] + p[1], s[2] + p[2]]).map(|v| v / x.len() as f64);
|
||||
let mut longest = 0.0f64;
|
||||
for edge in e.kind.edges() {
|
||||
let at = |i: usize| edge.get(i).and_then(|&k| x.get(k)).copied();
|
||||
if let (Some(a), Some(b)) = (at(0), at(1)) {
|
||||
longest = longest.max(((a[0] - b[0]).powi(2) + (a[1] - b[1]).powi(2) + (a[2] - b[2]).powi(2)).sqrt());
|
||||
}
|
||||
}
|
||||
Ok((c, 2.0 * longest))
|
||||
}
|
||||
|
||||
/// The nodes to move and their weights, and the sphere when there is one.
|
||||
pub(crate) fn weights(model: &Model, p: &Value) -> Result<(BTreeMap<u32, f64>, Option<([f64; 3], f64)>)> {
|
||||
let strength = super::optional_number(p, "strength", 1.0)?;
|
||||
if !(strength > 0.0 && strength <= 1.0) {
|
||||
return Err(EngineError::param("strength", ProblemKind::NotInRange { min: 0.0, max: 1.0 }));
|
||||
}
|
||||
let falloff = match super::opt::<String>(p, "falloff")? {
|
||||
None => Falloff::Smooth,
|
||||
Some(f) => Falloff::from_name(&f).ok_or_else(|| EngineError::param("falloff", ProblemKind::NotAChoice))?,
|
||||
};
|
||||
let only = body_nodes(model, p)?.map(|(_, n)| n);
|
||||
let around = match super::opt::<u32>(p, "around_element")? {
|
||||
Some(e) => Some(around_element(model, e)?),
|
||||
None => None,
|
||||
};
|
||||
let center = point(p, "center")?.or(around.map(|a| a.0));
|
||||
let radius = match p.get("radius") {
|
||||
Some(_) => Some(super::optional_number(p, "radius", 0.0)?)
|
||||
.filter(|r| *r > 0.0 && r.is_finite())
|
||||
.ok_or_else(|| EngineError::param("radius", ProblemKind::NotPositive))
|
||||
.map(Some)?,
|
||||
None => around.map(|a| a.1),
|
||||
};
|
||||
let listed = ["nodes", "elements", "set"].iter().filter(|k| p.get(**k).is_some()).count();
|
||||
if let (Some(c), Some(r)) = (center, radius) {
|
||||
if listed > 0 {
|
||||
return Err(EngineError::Params(
|
||||
"give a sphere (`center` and `radius`, or `around_element`) or `nodes`, `elements` or `set`, not both".into(),
|
||||
));
|
||||
}
|
||||
return Ok((relax::sphere_weights(model, c, r, strength, falloff, only.as_ref())?, Some((c, r))));
|
||||
}
|
||||
if center.is_some() != radius.is_some() {
|
||||
return Err(EngineError::Params("a sphere of influence needs both `center` and `radius`".into()));
|
||||
}
|
||||
let mut nodes: BTreeSet<u32> = BTreeSet::new();
|
||||
if p.get("nodes").is_some() {
|
||||
nodes.extend(super::labels(p, "nodes", fieldcraft_model::MAX_NODES)?);
|
||||
}
|
||||
let mut elements: Vec<u32> = Vec::new();
|
||||
if p.get("elements").is_some() {
|
||||
elements.extend(super::labels(p, "elements", fieldcraft_model::MAX_ELEMENTS)?);
|
||||
}
|
||||
if let Some(name) = super::opt::<String>(p, "set")? {
|
||||
let set = model.resolve_set(&name)?;
|
||||
nodes.extend(set.nodes.iter().copied());
|
||||
elements.extend(set.elements.iter().copied());
|
||||
elements.extend(set.faces.iter().map(|f| f.element));
|
||||
}
|
||||
for e in elements {
|
||||
let el = model.element(e).ok_or_else(|| EngineError::Params(format!("no element {e}")))?;
|
||||
nodes.extend(el.nodes.iter().copied());
|
||||
}
|
||||
if listed == 0 {
|
||||
return Err(EngineError::Params("say which nodes to relax: `center` and `radius`, `around_element`, `nodes`, `elements` or `set`".into()));
|
||||
}
|
||||
if let Some(only) = &only {
|
||||
nodes.retain(|n| only.contains(n));
|
||||
}
|
||||
Ok((relax::uniform_weights(model, &nodes, strength)?, None))
|
||||
}
|
||||
|
||||
fn iterations(p: &Value) -> Result<usize> {
|
||||
match p.get("iterations") {
|
||||
None => Ok(DEFAULT_ITERATIONS),
|
||||
Some(v) => v
|
||||
.as_u64()
|
||||
.and_then(|i| usize::try_from(i).ok())
|
||||
.filter(|i| (1..=relax::MAX_ITERATIONS).contains(i))
|
||||
.ok_or_else(|| EngineError::param("iterations", ProblemKind::NotInRange { min: 1.0, max: relax::MAX_ITERATIONS as f64 })),
|
||||
}
|
||||
}
|
||||
|
||||
/// The nodes the relaxation may touch: the weighted ones and their elements' nodes.
|
||||
fn touched(model: &Model, weights: &BTreeMap<u32, f64>) -> BTreeSet<u32> {
|
||||
let mut out: BTreeSet<u32> = weights.keys().copied().collect();
|
||||
for e in model.elements().values() {
|
||||
if e.nodes.iter().any(|n| weights.contains_key(n)) {
|
||||
out.extend(e.nodes.iter().copied());
|
||||
}
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
/// What `mesh.node.relax` with `params` would do, without doing it: the open form's preview.
|
||||
pub struct Preview {
|
||||
pub sphere: Option<([f64; 3], f64)>,
|
||||
pub relaxed: Relaxed,
|
||||
}
|
||||
|
||||
pub fn preview(s: &Session, params: &Value) -> Option<Preview> {
|
||||
let (w, sphere) = weights(&s.model, params).ok()?;
|
||||
let g = design_geometry(s, &touched(&s.model, &w)).ok()?;
|
||||
let relaxed = relax::relax(&s.model, &w, iterations(params).ok()?, &g).ok()?;
|
||||
Some(Preview { sphere, relaxed })
|
||||
}
|
||||
|
||||
/// The sphere `params` describe, for the form's gizmo before the rest is valid.
|
||||
pub fn sphere(s: &Session, params: &Value) -> Option<([f64; 3], f64)> {
|
||||
let around = params.get("around_element").and_then(Value::as_u64).and_then(|e| around_element(&s.model, u32::try_from(e).ok()?).ok());
|
||||
let c = point(params, "center").ok().flatten().or(around.map(|a| a.0))?;
|
||||
let r = params.get("radius").and_then(Value::as_f64).filter(|r| *r > 0.0 && r.is_finite()).or(around.map(|a| a.1))?;
|
||||
Some((c, r))
|
||||
}
|
||||
|
||||
fn relax_nodes(s: &mut Session, p: &Value) -> Result<Value> {
|
||||
let (w, sphere) = weights(&s.model, p)?;
|
||||
if w.is_empty() {
|
||||
return Err(EngineError::Params("no mesh nodes in the region to relax".into()));
|
||||
}
|
||||
let iterations = iterations(p)?;
|
||||
let g = design_geometry(s, &touched(&s.model, &w))?;
|
||||
let r = relax::relax(&s.model, &w, iterations, &g)?;
|
||||
let field_name = match super::opt::<String>(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 = match body_nodes(&s.model, p)? {
|
||||
Some((b, _)) => b,
|
||||
None => Vec::new(),
|
||||
};
|
||||
let reapply = super::opt::<bool>(p, "reapply")?.unwrap_or(false);
|
||||
if !r.positions.is_empty() {
|
||||
let cells = relax::field_cells(&s.model, &r);
|
||||
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!(
|
||||
"moved {} node{} (largest {:.3e}); scaled Jacobian min {:.3} → {:.3}{}",
|
||||
r.moved,
|
||||
if r.moved == 1 { "" } else { "s" },
|
||||
r.max_displacement,
|
||||
r.before.min_quality,
|
||||
r.after.min_quality,
|
||||
if r.positions.is_empty() { String::new() } else { format!("; kept as field {field_name}") }
|
||||
);
|
||||
Ok(json!({
|
||||
"moved": r.moved,
|
||||
"nodes_changed": r.positions.len(),
|
||||
"max_displacement": r.max_displacement,
|
||||
"iterations": r.iterations,
|
||||
"sphere": sphere.map(|(c, r)| json!({"center": c, "radius": r})),
|
||||
"before": quality_json(&r.before),
|
||||
"after": quality_json(&r.after),
|
||||
"field": (!r.positions.is_empty()).then_some(field_name),
|
||||
"summary": summary,
|
||||
}))
|
||||
}
|
||||
|
||||
/// Apply field `name` to the nodes of `bodies` (or of the whole mesh).
|
||||
fn apply_field(s: &mut Session, name: &str, bodies: &[String], polish: f64) -> Result<Relaxed> {
|
||||
let field = s
|
||||
.model
|
||||
.relaxation(name)
|
||||
.cloned()
|
||||
.ok_or_else(|| EngineError::Params(format!("no relaxation field `{}`", name.chars().take(80).collect::<String>())))?;
|
||||
let mut nodes = BTreeSet::new();
|
||||
let bodies = if bodies.is_empty() { &field.bodies[..] } else { bodies };
|
||||
if bodies.is_empty() {
|
||||
nodes.extend(s.model.nodes().keys().copied());
|
||||
} else {
|
||||
for b in bodies {
|
||||
for e in s.model.body_elements(b) {
|
||||
if let Some(el) = s.model.element(e) {
|
||||
nodes.extend(el.nodes.iter().copied());
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Only nodes the field reaches need constraints.
|
||||
nodes.retain(|n| s.model.node(*n).is_some_and(|p| relax::field_displacement(&field.cells, p.xyz).is_some()));
|
||||
let g = design_geometry(s, &nodes)?;
|
||||
let r = relax::apply_field(&s.model, &field.cells, Some(&nodes), &g, polish)?;
|
||||
for (n, x) in &r.positions {
|
||||
s.model.set_node(*n, *x)?;
|
||||
}
|
||||
Ok(r)
|
||||
}
|
||||
|
||||
fn apply(s: &mut Session, p: &Value) -> Result<Value> {
|
||||
let name: String = super::req(p, "field")?;
|
||||
let bodies = super::opt::<Vec<String>>(p, "bodies")?.unwrap_or_default();
|
||||
let polish = super::optional_number(p, "polish", 0.3)?;
|
||||
if !(0.0..=1.0).contains(&polish) {
|
||||
return Err(EngineError::param("polish", ProblemKind::NotInRange { min: 0.0, max: 1.0 }));
|
||||
}
|
||||
let r = apply_field(s, &name, &bodies, polish)?;
|
||||
Ok(json!({
|
||||
"moved": r.moved,
|
||||
"nodes_changed": r.positions.len(),
|
||||
"max_displacement": r.max_displacement,
|
||||
"summary": format!("applied field {name}: moved {} node{} (largest {:.3e})", r.moved, if r.moved == 1 { "" } else { "s" }, r.max_displacement),
|
||||
}))
|
||||
}
|
||||
|
||||
/// After a remesh (`id`), apply again every field marked `reapply`; its reply lists them.
|
||||
pub(crate) fn after_remesh(s: &mut Session, id: &str, mut reply: Value) -> Result<Value> {
|
||||
if !REMESH.contains(&id) {
|
||||
return Ok(reply);
|
||||
}
|
||||
let fields: Vec<String> = s.model.relaxations().iter().filter(|(_, f)| f.reapply).map(|(n, _)| n.clone()).collect();
|
||||
let mut applied = Vec::new();
|
||||
for name in fields {
|
||||
let r = apply_field(s, &name, &[], 0.3)?;
|
||||
applied.push(json!({"field": name, "moved": r.moved}));
|
||||
}
|
||||
if !applied.is_empty()
|
||||
&& let Some(o) = reply.as_object_mut()
|
||||
{
|
||||
o.insert("relaxations_applied".into(), Value::Array(applied));
|
||||
}
|
||||
Ok(reply)
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests;
|
||||
@@ -0,0 +1,194 @@
|
||||
use super::*;
|
||||
|
||||
/// The plate (40 × 30 × 10) with its Ø10 hole round (15, 15): a point on the hole's wall.
|
||||
const ON_HOLE: [f64; 3] = [20.0, 15.0, 5.0];
|
||||
const HOLE: ([f64; 2], f64) = ([15.0, 15.0], 5.0);
|
||||
|
||||
fn plate() -> Session {
|
||||
super::super::datum::tests::plate().0
|
||||
}
|
||||
|
||||
fn hole_radius(p: [f64; 3]) -> f64 {
|
||||
((p[0] - HOLE.0[0]).powi(2) + (p[1] - HOLE.0[1]).powi(2)).sqrt()
|
||||
}
|
||||
|
||||
/// Which analytic face a node of the plate lies on, if any.
|
||||
fn on_faces(p: [f64; 3]) -> Vec<&'static str> {
|
||||
let mut f = Vec::new();
|
||||
if (hole_radius(p) - HOLE.1).abs() < 1e-6 {
|
||||
f.push("hole");
|
||||
}
|
||||
if p[2].abs() < 1e-6 || (p[2] - 10.0).abs() < 1e-6 {
|
||||
f.push("cap");
|
||||
}
|
||||
f
|
||||
}
|
||||
|
||||
/// Spikes near the hole: interior nodes pushed about; nodes on the hole turned round it and
|
||||
/// moved along its axis (still on the cylinder), nodes on the caps moved in their plane.
|
||||
fn perturb(s: &mut Session, centre: [f64; 3], radius: f64, amount: f64) {
|
||||
let nodes: Vec<(u32, [f64; 3])> = s.model.nodes().iter().map(|(n, p)| (*n, p.xyz)).collect();
|
||||
for (k, (n, p)) in nodes.into_iter().enumerate() {
|
||||
let d2: f64 = (0..3).map(|i| (p[i] - centre[i]).powi(2)).sum();
|
||||
if d2 > radius * radius {
|
||||
continue;
|
||||
}
|
||||
let wiggle = |m: usize| amount * ((k * m % 7) as f64 - 3.0) / 3.0;
|
||||
let faces = on_faces(p);
|
||||
let q = match faces.as_slice() {
|
||||
[] => [p[0] + wiggle(3), p[1] + wiggle(5), p[2] + wiggle(11)],
|
||||
["hole"] => {
|
||||
let a = wiggle(3) / HOLE.1;
|
||||
let (x, y) = (p[0] - HOLE.0[0], p[1] - HOLE.0[1]);
|
||||
[HOLE.0[0] + x * a.cos() - y * a.sin(), HOLE.0[1] + x * a.sin() + y * a.cos(), p[2] + wiggle(5)]
|
||||
}
|
||||
["cap"] if hole_radius(p) > HOLE.1 + 1e-6 => [p[0] + wiggle(3), p[1] + wiggle(5), p[2]],
|
||||
_ => continue,
|
||||
};
|
||||
s.model.set_node(n, q).unwrap();
|
||||
}
|
||||
}
|
||||
|
||||
fn check_on_geometry(before: &Model, after: &Model) {
|
||||
for (n, p) in after.nodes() {
|
||||
let old = before.node(*n).unwrap().xyz;
|
||||
for f in on_faces(old) {
|
||||
match f {
|
||||
"hole" => assert!((hole_radius(p.xyz) - HOLE.1).abs() < 1e-6, "node {n} left the hole: r = {}", hole_radius(p.xyz)),
|
||||
_ => assert!((p.xyz[2] - old[2]).abs() < 1e-6, "node {n} left its cap"),
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fn invalid(s: &mut Session) -> u64 {
|
||||
s.execute("mesh.quality", &json!({})).unwrap()["invalid"]["count"].as_u64().unwrap_or(0)
|
||||
}
|
||||
|
||||
/// The issue's hex test on a design: spikes near the hole relax, every surface node stays on its
|
||||
/// face (the hole's cylinder exactly), quality rises, nothing inverts, nothing outside the sphere
|
||||
/// moves; one undo step; the field is kept and saved with the model.
|
||||
#[test]
|
||||
fn hex_spikes_near_a_hole_relax_on_the_geometry() {
|
||||
let mut s = plate();
|
||||
s.execute("mesh.hex.generate", &json!({"size": 2.0})).unwrap();
|
||||
perturb(&mut s, ON_HOLE, 3.0, 0.35);
|
||||
let perturbed = s.model.clone();
|
||||
let r = s.execute("mesh.node.relax", &json!({"center": ON_HOLE, "radius": 4.0, "iterations": 10})).unwrap();
|
||||
let (before, after) = (r["before"]["min_scaled_jacobian"].as_f64().unwrap(), r["after"]["min_scaled_jacobian"].as_f64().unwrap());
|
||||
assert!(after > before + 0.05, "{r}");
|
||||
assert_eq!(r["after"]["invalid"], 0, "{r}");
|
||||
assert!(r["moved"].as_u64().unwrap() > 0 && r["field"] == "Relax1", "{r}");
|
||||
check_on_geometry(&perturbed, &s.model);
|
||||
for (n, p) in s.model.nodes() {
|
||||
let old = perturbed.node(*n).unwrap().xyz;
|
||||
if old != p.xyz {
|
||||
let d2: f64 = (0..3).map(|i| (old[i] - ON_HOLE[i]).powi(2)).sum();
|
||||
assert!(d2 < 16.0, "node {n} outside the sphere moved");
|
||||
}
|
||||
}
|
||||
assert_eq!(invalid(&mut s), 0);
|
||||
// The field is saved with the model.
|
||||
let text = serde_json::to_string(&s.model).unwrap();
|
||||
let back: Model = serde_json::from_str(&text).unwrap();
|
||||
assert_eq!(back.relaxation("Relax1"), s.model.relaxation("Relax1"));
|
||||
s.execute("edit.undo", &json!({})).unwrap();
|
||||
assert_eq!(s.model, perturbed);
|
||||
}
|
||||
|
||||
/// The same on tet10 with curved mid-side nodes on the hole: they stay on the cylinder.
|
||||
#[test]
|
||||
fn tet10_relaxes_with_curved_mid_side_nodes_on_the_hole() {
|
||||
let mut s = plate();
|
||||
s.execute("mesh.tet.generate", &json!({"size": 2.5, "order": "quadratic"})).unwrap();
|
||||
perturb(&mut s, ON_HOLE, 3.0, 0.2);
|
||||
// Put the mid-side nodes back on their edges' middles, on the hole where the edge is.
|
||||
let mids: Vec<(u32, u32, u32)> = s
|
||||
.model
|
||||
.elements()
|
||||
.values()
|
||||
.flat_map(|e| e.kind.edges().iter().map(|ed| (e.nodes[ed[0]], e.nodes[ed[1]], e.nodes[ed[2]])).collect::<Vec<_>>())
|
||||
.collect();
|
||||
let original = s.model.clone();
|
||||
for (a, b, m) in mids {
|
||||
let (pa, pb, pm) = (original.node(a).unwrap().xyz, original.node(b).unwrap().xyz, original.node(m).unwrap().xyz);
|
||||
let mut q = [(pa[0] + pb[0]) / 2.0, (pa[1] + pb[1]) / 2.0, (pa[2] + pb[2]) / 2.0];
|
||||
if on_faces(pm).contains(&"hole") {
|
||||
let r = hole_radius(q);
|
||||
q = [HOLE.0[0] + (q[0] - HOLE.0[0]) * HOLE.1 / r, HOLE.0[1] + (q[1] - HOLE.0[1]) * HOLE.1 / r, q[2]];
|
||||
}
|
||||
s.model.set_node(m, q).unwrap();
|
||||
}
|
||||
let perturbed = s.model.clone();
|
||||
let r = s.execute("mesh.node.relax", &json!({"center": ON_HOLE, "radius": 4.0, "iterations": 8})).unwrap();
|
||||
assert!(r["after"]["min_scaled_jacobian"].as_f64().unwrap() >= r["before"]["min_scaled_jacobian"].as_f64().unwrap(), "{r}");
|
||||
assert_eq!(r["after"]["invalid"], 0, "{r}");
|
||||
check_on_geometry(&perturbed, &s.model);
|
||||
}
|
||||
|
||||
/// A partition cut is not geometry: nodes on it move freely, shared by the cells on both sides.
|
||||
#[test]
|
||||
fn nodes_on_a_partition_cut_move_freely() {
|
||||
let mut s = plate();
|
||||
s.execute("geometry.datum", &json!({"name": "Cut", "by": {"type": "offset", "plane": "yz", "offset": 27}})).unwrap();
|
||||
s.execute("geometry.partition", &json!({"name": "P", "datum": "Cut"})).unwrap();
|
||||
s.execute("mesh.cells", &json!({"size": 2.0, "default": "hex"})).unwrap();
|
||||
let centre = [27.0, 15.0, 5.0];
|
||||
perturb(&mut s, centre, 2.5, 0.3);
|
||||
let before = s.model.clone();
|
||||
s.execute("mesh.node.relax", &json!({"center": centre, "radius": 3.5, "iterations": 10})).unwrap();
|
||||
let moved_off = s
|
||||
.model
|
||||
.nodes()
|
||||
.iter()
|
||||
.filter(|(n, p)| {
|
||||
let old = before.node(**n).unwrap().xyz;
|
||||
(old[0] - 27.0).abs() < 1e-6 && on_faces(old).is_empty() && (p.xyz[0] - 27.0).abs() > 1e-6
|
||||
})
|
||||
.count();
|
||||
assert!(moved_off > 0, "nodes on the cut are free");
|
||||
assert_eq!(invalid(&mut s), 0);
|
||||
}
|
||||
|
||||
/// A field marked `reapply` is applied again after a remesh; `mesh.node.relax.apply` applies it
|
||||
/// on demand; around the worst element works.
|
||||
#[test]
|
||||
fn a_kept_field_is_applied_again_after_a_remesh() {
|
||||
let mut s = plate();
|
||||
s.execute("mesh.hex.generate", &json!({"size": 3.0})).unwrap();
|
||||
perturb(&mut s, ON_HOLE, 4.0, 0.3);
|
||||
let r = s.execute("mesh.node.relax", &json!({"center": ON_HOLE, "radius": 5.0, "field": "Spike", "reapply": true})).unwrap();
|
||||
assert_eq!(r["field"], "Spike");
|
||||
let r = s.execute("mesh.hex.generate", &json!({"size": 1.5})).unwrap();
|
||||
assert_eq!(r["relaxations_applied"][0]["field"], "Spike", "{r}");
|
||||
assert_eq!(invalid(&mut s), 0);
|
||||
let r = s.execute("mesh.node.relax.apply", &json!({"field": "Spike"})).unwrap();
|
||||
assert!(r["summary"].as_str().unwrap().contains("Spike"));
|
||||
let worst = s.execute("mesh.quality", &json!({})).unwrap()["metrics"]["scaled_jacobian"]["worst"].as_u64().unwrap();
|
||||
let r = s.execute("mesh.node.relax", &json!({"around_element": worst})).unwrap();
|
||||
assert!(r["after"]["min_scaled_jacobian"].as_f64() >= r["before"]["min_scaled_jacobian"].as_f64(), "{r}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn hostile_parameters_never_panic() {
|
||||
let mut s = plate();
|
||||
s.execute("mesh.hex.generate", &json!({"size": 4.0})).unwrap();
|
||||
for p in [
|
||||
json!({}),
|
||||
json!({"center": [0, 0, 0]}),
|
||||
json!({"center": [1e300, 0, 0], "radius": 1}),
|
||||
json!({"center": [0, 0, 0], "radius": -1}),
|
||||
json!({"center": [0, 0, 0], "radius": 1, "nodes": [1]}),
|
||||
json!({"nodes": [1], "strength": 2}),
|
||||
json!({"nodes": [1], "iterations": 0}),
|
||||
json!({"nodes": [1], "iterations": 1000}),
|
||||
json!({"nodes": [1], "falloff": "cubic"}),
|
||||
json!({"around_element": 99999999}),
|
||||
json!({"set": "nope"}),
|
||||
json!({"nodes": [1], "bodies": ["nope"]}),
|
||||
] {
|
||||
assert!(s.execute("mesh.node.relax", &p).is_err(), "{p}");
|
||||
}
|
||||
assert!(s.execute("mesh.node.relax.apply", &json!({"field": "nope"})).is_err());
|
||||
assert!(s.execute("mesh.node.relax.apply", &json!({})).is_err());
|
||||
}
|
||||
Reference in new issue
Block a user