From 72ebd5bc7fa6bfd547dc8562ba601d5eadf61b9d Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Fri, 16 Oct 2026 13:36:42 -0400 Subject: [PATCH] mesh.cells: tet cells meet hex cells through pyramids With "transition": "pyramid" (the default; "tie" keeps the old join) a tet cell is meshed after the hex cells, with a pyramid on every hex face of its cuts and its tetrahedra on the pyramids' sides: one conforming mesh, no tie. Where it cannot be built the cell says why (no_pyramids) and is tied. The apex height is 0.4 of the quad's edge, so pyramids on faces meeting at a right angle do not cross. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo --- crates/engine/src/catalogue.rs | 3 +- crates/engine/src/cmd/cells.rs | 127 ++++++++++++---- crates/engine/src/cmd/cells/transition.rs | 172 ++++++++++++++++++++++ crates/engine/src/cmd/forms.rs | 3 + crates/engine/src/cmd/mesh.rs | 12 +- crates/mesh/src/transition.rs | 28 ++-- 6 files changed, 296 insertions(+), 49 deletions(-) create mode 100644 crates/engine/src/cmd/cells/transition.rs diff --git a/crates/engine/src/catalogue.rs b/crates/engine/src/catalogue.rs index 8e4deca..18a76f3 100644 --- a/crates/engine/src/catalogue.rs +++ b/crates/engine/src/catalogue.rs @@ -359,7 +359,7 @@ pub static CATALOGUE: &[Capability] = &[ "Partitioned Mesh", Mesh, "Structured Hex", - "One mesh of a partitioned part, a method per cell (hex, tet, or all hex as conforming blocks); cells joined across partitions", + "One mesh of a partitioned part, a method per cell (hex, tet, or all hex as conforming blocks); cells joined across partitions, tet cells to hex cells through conforming pyramids", &["mesh.hex.blocks"], ), ready( @@ -393,7 +393,6 @@ pub static CATALOGUE: &[Capability] = &[ "M5", "Find sweepable and mappable regions of 3D parts (extruded cross-sections are already cut into blocks by Hex Mesh)", ), - soon("mesh.transition", "Transition", Mesh, "Structured Hex", "M5", "Conforming tet–hex interfaces with pyramids"), ready( "mesh.hex.octree", "Octree Hex", diff --git a/crates/engine/src/cmd/cells.rs b/crates/engine/src/cmd/cells.rs index c9e6010..7b36b7e 100644 --- a/crates/engine/src/cmd/cells.rs +++ b/crates/engine/src/cmd/cells.rs @@ -35,14 +35,14 @@ use fieldcraft_model::{ElementFace, GeometryRef, Model, Set, Tie}; use serde_json::{Value, json}; use super::hex::{Input, add_cells, add_hex, hex_body}; -use super::mesh::{BodyTets, add_tets, mesh_body, quality_options}; +use super::mesh::{add_tets, mesh_body, quality_options}; use crate::requirement::Requirement; use crate::{CommandKind, CommandSpec, EngineError, Result, Session}; pub const CELLS: CommandSpec = CommandSpec { id: "mesh.cells", label: "Partitioned Mesh", - params: r#"{"methods": {"Body1/a": "tet", "Body1/b": "swept"} (per cell; default `default`), "default": "auto", "bodies": ["Body1/a", "Body1/b"] (default: every body of the open design), "size": 2.0 (element size, mm; default 1/25 of the bounding-box diagonal of all of them), "order": 1 | 2 (hex8 and tet4, or hex20 and tet10), "tolerance": 0.05 (chordal tolerance, mm)} — replaces the cells' meshes (as mesh.tet.generate; other bodies keep theirs) with one mesh of the cells, each with its method: "auto" (structured hex where a method fits, else tet), "hex" (any structured hex method), "swept", "extrude", "sweep", "revolve", "mapped" or "tet"; or "default": "blocks" (with no `methods`): every cell a mapped block, division counts assigned over all cells together, exactly mesh.hex.blocks with the same parameters. Cells on either side of a partition (faces `partition::a` and `:b`) are joined: where their element faces match node for node the nodes are merged (conforming); otherwise, as between a hex and a tet cell, they are tied: each cell's tied cut faces make one tie, `tie_` (`/` as `_`; sets `_main`, the facing cells' cut faces, and `_secondary`, the cell's own), hex cells main to tet cells and otherwise by name, so a node is a secondary node of one such tie only; the ties an earlier run made on these cells are replaced. The cells of one design body that may sweep (auto, hex, swept, extrude, revolve) are swept together first: one cross-section and one set of layers for all of them, so they share their nodes across every cut, T-junctions included, with no edge seeds (`swept_together` names the body; `not_swept_together` says why not, and those cells are meshed one by one). Reports each cell's method (and why auto fell back to tet), its `quality` (volume, smallest and mean scaled Jacobian, largest edge ratio), and each interface: `merged` with the nodes merged or shared, or `tied` (with `why` when both cells are of one design body: the tet cell and why auto fell back, or why the body's cells were not swept together)"#, + params: r#"{"methods": {"Body1/a": "tet", "Body1/b": "swept"} (per cell; default `default`), "default": "auto", "bodies": ["Body1/a", "Body1/b"] (default: every body of the open design), "size": 2.0 (element size, mm; default 1/25 of the bounding-box diagonal of all of them), "order": 1 | 2 (hex8 and tet4, or hex20 and tet10), "tolerance": 0.05 (chordal tolerance, mm), "transition": "pyramid" | "tie" (how a tet cell meets hex cells; default pyramid)} — replaces the cells' meshes (as mesh.tet.generate; other bodies keep theirs) with one mesh of the cells, each with its method: "auto" (structured hex where a method fits, else tet), "hex" (any structured hex method), "swept", "extrude", "sweep", "revolve", "mapped" or "tet"; or "default": "blocks" (with no `methods`): every cell a mapped block, division counts assigned over all cells together, exactly mesh.hex.blocks with the same parameters. Cells on either side of a partition (faces `partition::a` and `:b`) are joined: where their element faces match node for node the nodes are merged (conforming); between a tet cell and the hex cells across its partitions, with "transition": "pyramid" (linear meshes), the tet cell is meshed after them with a pyramid (pyramid5) on every hex element face of the cut, its tetrahedra sharing the pyramids' sides, so the mesh is conforming there (`joined`: `pyramids`, with their count); where that cannot be built (a cut the hex faces do not cover exactly, quadratic elements, a pyramid that would reach through a thin cell) the cell says why in `no_pyramids` and is tied; otherwise, as between a hex and a tet cell with "transition": "tie", they are tied: each cell's tied cut faces make one tie, `tie_` (`/` as `_`; sets `_main`, the facing cells' cut faces, and `_secondary`, the cell's own), hex cells main to tet cells and otherwise by name, so a node is a secondary node of one such tie only; the ties an earlier run made on these cells are replaced. The cells of one design body that may sweep (auto, hex, swept, extrude, revolve) are swept together first: one cross-section and one set of layers for all of them, so they share their nodes across every cut, T-junctions included, with no edge seeds (`swept_together` names the body; `not_swept_together` says why not, and those cells are meshed one by one). Reports each cell's method (and why auto fell back to tet), its `quality` (volume, smallest and mean scaled Jacobian, largest edge ratio), and each interface: `merged` with the nodes merged or shared, or `tied` (with `why` when both cells are of one design body: the tet cell and why auto fell back, or why the body's cells were not swept together)"#, form: Some(super::forms::CELLS), kind: CommandKind::Model, requires: &[Requirement::Design], @@ -80,10 +80,22 @@ impl Method { } } -/// A meshed cell. +mod transition; + +/// A meshed cell; tet cells are meshed after the hex cells are in the model, so they can meet +/// them through pyramids. enum Meshed<'a> { - Hex(Input<'a>, HexMesh), - Tet(BodyTets), + Hex(Box<(Input<'a>, HexMesh)>), + Tet, +} + +/// How a tet cell joins the hex cells across its partitions. +#[derive(Clone, Copy, PartialEq, Eq)] +enum Transition { + /// Pyramids on the hex faces, conforming (falls back to a tie with the reason). + Pyramid, + /// A tie (dual mortar), each mesh as its method made it. + Tie, } fn cells(s: &mut Session, p: &Value) -> Result { @@ -105,6 +117,11 @@ fn cells(s: &mut Session, p: &Value) -> Result { _ => return Err(EngineError::invalid("order", "must be 1 or 2")), }, }; + let transition = match p.get("transition").map(|v| v.as_str()) { + None | Some(Some("pyramid")) => Transition::Pyramid, + Some(Some("tie")) => Transition::Tie, + Some(_) => return Err(EngineError::invalid("transition", "must be \"pyramid\" or \"tie\"")), + }; let default = p.get("default").map(Method::parse).transpose()?.unwrap_or(Method::Auto); let methods: BTreeMap = match p.get("methods") { None => BTreeMap::new(), @@ -215,20 +232,19 @@ fn cells(s: &mut Session, p: &Value) -> Result { let mesh = hex_body(&hex_params, &b, m, quadratic, size)?; Ok((b, mesh)) }; - let tet = |s: &Session| mesh_body(s, &tet_params, input()?, &cdt, &tet_options, quality, true); let (cell, mut info) = match method { Method::Hex(m) => { let (b, mesh) = hex(m)?; let info = json!({"method": mesh.report.method.name()}); - (Meshed::Hex(b, mesh), info) + (Meshed::Hex(Box::new((b, mesh))), info) } - Method::Tet => (Meshed::Tet(tet(s)?), json!({"method": "tet"})), + Method::Tet => (Meshed::Tet, json!({"method": "tet"})), Method::Auto => match hex(HexMethod::Auto) { Ok((b, mesh)) => { let info = json!({"method": mesh.report.method.name()}); - (Meshed::Hex(b, mesh), info) + (Meshed::Hex(Box::new((b, mesh))), info) } - Err(e) => (Meshed::Tet(tet(s)?), json!({"method": "tet", "why": e.to_string()})), + Err(e) => (Meshed::Tet, json!({"method": "tet", "why": e.to_string()})), }, }; if let Some(o) = info.as_object_mut() { @@ -259,19 +275,50 @@ fn cells(s: &mut Session, p: &Value) -> Result { } interfaces.extend(swept_interfaces(&design, inputs, cm)); } - let mut k = 0; - for (i, _) in shapes.iter().enumerate() { - if joint.get(i).copied().flatten().is_some() { - continue; - } - let Some((cell, mut info)) = meshed.get_mut(k).map(|(c, v)| (c, std::mem::take(v))) else { break }; - k += 1; - let (kind, elements, nodes) = match &*cell { - Meshed::Hex(b, m) => { + // The cells meshed alone, in the order of `shapes`: the hex ones first, then the tet ones, + // which may meet the hex cells through pyramids. + let alone: Vec = (0..shapes.len()).filter(|&i| joint.get(i).copied().flatten().is_none()).collect(); + let mut hex_cells: Vec = swept.iter().flat_map(|(_, _, inputs, _)| inputs.iter().map(|b| b.name.clone())).collect(); + let mut tet_cells = Vec::new(); + for (&i, (cell, info)) in alone.iter().zip(meshed.iter_mut()) { + let mut info = std::mem::take(info); + match &*cell { + Meshed::Hex(hex) => { + let (b, m) = &**hex; add_hex(model, b, m)?; - (m.kind, m.element_count(), m.nodes.len()) + hex_cells.push(b.name.clone()); + if let Some(o) = info.as_object_mut() { + o.insert("kind".into(), json!(m.kind.name())); + o.insert("elements".into(), json!(m.element_count())); + o.insert("nodes".into(), json!(m.nodes.len())); + } + report.push((i, info)); } - Meshed::Tet(t) => { + Meshed::Tet => tet_cells.push((i, info)), + } + } + for (i, mut info) in tet_cells { + let Some(shape) = shapes.get(i) else { continue }; + let input = Input { seeds: seeds.clone(), ..super::hex::shape_input(shape, Some(tolerance))? }; + let pyramids = if transition == Transition::Tie { + Err("`transition` is \"tie\"".to_string()) + } else if quadratic { + Err("pyramid transitions are linear (pyramid5); order 2 is tied".to_string()) + } else { + transition::mesh_with_pyramids(&design, model, &tet_params, &input, &hex_cells, &cdt, &tet_options)? + }; + let (kind, elements, nodes) = match pyramids { + Ok(c) => { + interfaces.extend(c.interfaces.iter().cloned()); + let kind = c.tets.mesh.kind; + let (elements, nodes) = transition::add(model, c)?; + (kind, elements, nodes) + } + Err(why) => { + if let Some(o) = info.as_object_mut() { + o.insert("no_pyramids".into(), json!(why)); + } + let t = mesh_body(Some(&design), &tet_params, input, &cdt, &tet_options, quality, true)?; add_tets(model, &t.name, &t.face_names, &t.mesh)?; (t.mesh.kind, t.mesh.elements.len(), t.mesh.points.len()) } @@ -660,7 +707,7 @@ mod tests { fn hex_and_tet_cells_are_tied_and_matching_hex_cells_merged() { let (_dir, mut s) = cut_plate(); // A tet cell next to a hex cell: tied, the hex cell main. - let r = s.execute("mesh.cells", &json!({"methods": {"Body1/a": "tet", "Body1/b": "mapped"}, "size": 2.5})).unwrap(); + let r = s.execute("mesh.cells", &json!({"methods": {"Body1/a": "tet", "Body1/b": "mapped"}, "size": 2.5, "transition": "tie"})).unwrap(); assert_eq!(r["cells"][0]["kind"], "tet4", "{r}"); assert_eq!(r["cells"][1]["kind"], "hex8", "{r}"); let i = &r["interfaces"][0]; @@ -679,7 +726,7 @@ mod tests { let separate = r["cells"].as_array().unwrap().iter().map(|c| c["nodes"].as_u64().unwrap()).sum::(); assert_eq!(r["nodes"].as_u64().unwrap() + i["nodes"].as_u64().unwrap(), separate, "{r}"); // Rerunning replaces the tie it made rather than adding another. - s.execute("mesh.cells", &json!({"methods": {"Body1/a": "tet"}, "size": 2.5})).unwrap(); + s.execute("mesh.cells", &json!({"methods": {"Body1/a": "tet"}, "size": 2.5, "transition": "tie"})).unwrap(); assert_eq!(s.model.ties().len(), 1); } @@ -762,11 +809,10 @@ mod tests { } /// A 40 × 10 × 10 bar in uniaxial tension, cut in the middle: one half hex, the other tet, - /// tied. With ν = 0, symmetry conditions on three faces and the end load spread by consistent - /// weights, the stress is uniform, which both elements and the dual mortar tie carry exactly: - /// the end moves F L / (E A). - #[test] - fn a_tied_hex_tet_bar_in_tension_stretches_as_the_formula_says() { + /// joined by `transition` (`joined` in the reply). With ν = 0, symmetry conditions on three + /// faces and the end load spread by consistent weights, the stress is uniform, which the + /// elements, the pyramids and the dual mortar tie all carry exactly: the end moves F L / (E A). + fn hex_tet_bar_in_tension(transition: &str, joined: &str) -> Session { let commands = json!({"commands": [ {"command": "solid.box", "params": {"length": 40, "width": 10, "height": 10, "body_name": "Bar"}}, ]}); @@ -786,8 +832,8 @@ mod tests { point[axis] = mid; s.execute("geometry.datum", &json!({"name": "Mid", "by": {"type": "point_normal", "point": point, "normal": normal}})).unwrap(); s.execute("geometry.partition", &json!({"name": "Cut", "datum": "Mid"})).unwrap(); - let r = s.execute("mesh.cells", &json!({"methods": {"Bar/a": "mapped", "Bar/b": "tet"}, "size": 2.5})).unwrap(); - assert_eq!(r["interfaces"][0]["joined"], "tied", "{r}"); + let r = s.execute("mesh.cells", &json!({"methods": {"Bar/a": "mapped", "Bar/b": "tet"}, "size": 2.5, "transition": transition})).unwrap(); + assert_eq!(r["interfaces"][0]["joined"], joined, "{r}"); let m = &s.model; let on = |k: usize, v: f64| -> Vec { m.nodes().iter().filter(|(_, n)| (n.xyz[k] - v).abs() < 1e-9).map(|(l, _)| *l).collect() }; @@ -829,6 +875,27 @@ mod tests { let u = n["displacement"][axis].as_f64().unwrap(); assert!((u - expected).abs() < 1e-6 * expected, "{u} vs {expected}: {n}"); } + s + } + + #[test] + fn a_tied_hex_tet_bar_in_tension_stretches_as_the_formula_says() { + let s = hex_tet_bar_in_tension("tie", "tied"); + assert_eq!(s.model.ties().len(), 1); + } + + /// Through pyramids the bar is one conforming mesh: no tie, no free faces inside, every + /// element valid. + #[test] + fn a_hex_tet_bar_joined_by_pyramids_is_conforming_and_stretches_as_the_formula_says() { + let mut s = hex_tet_bar_in_tension("pyramid", "pyramids"); + assert!(s.model.ties().is_empty()); + let kinds: BTreeSet<&str> = s.model.elements().values().map(|e| e.kind.name()).collect(); + assert_eq!(kinds.into_iter().collect::>(), ["hex8", "pyramid5", "tet4"]); + let d = s.execute("mesh.check.duplicates", &json!({})).unwrap(); + assert_eq!(d["counts"], json!({"coincident": 0, "cracks": 0, "duplicate_elements": 0, "free_edges": 0}), "{d}"); + let q = s.execute("mesh.quality", &json!({})).unwrap(); + assert_eq!(q["invalid"]["count"], 0, "{q}"); } /// A box cut by planes into cells, `methods` per cell, a heat source in every cell and diff --git a/crates/engine/src/cmd/cells/transition.rs b/crates/engine/src/cmd/cells/transition.rs new file mode 100644 index 0000000..5fc5924 --- /dev/null +++ b/crates/engine/src/cmd/cells/transition.rs @@ -0,0 +1,172 @@ +//! Pyramid transitions between a tet cell and the hex cells across its partitions (#439). +//! +//! The tet cell is meshed after the hex cells: its cut faces that face hex cells are replaced by +//! a pyramid on every hex element face there ([`fieldcraft_mesh::transition`]), the tetrahedra +//! fill the rest and keep the pyramids' sides, and the pyramids join the hex nodes to the tet +//! nodes. The mesh is conforming across the cut, so no tie is needed there. + +use std::collections::{BTreeMap, BTreeSet}; + +use fieldcraft_mesh::cdt::{CdtOptions, tetrahedralize}; +use fieldcraft_mesh::surface::BoundarySurface; +use fieldcraft_mesh::tet::{QualityMesh, QualityOptions, quality_mesh}; +use fieldcraft_mesh::transition::{Transition, apex_points, pyramid_surface}; +use fieldcraft_model::{ElementKind, GeometryRef, Model}; +use serde_json::{Value, json}; + +use super::super::hex::Input; +use super::super::mesh::{BodyTets, add_tets, analytic_surfaces}; +use crate::{EngineError, Result}; + +/// A tet cell meshed against the pyramids on its hex neighbours' faces. +pub(super) struct Pyramids { + pub(super) tets: BodyTets, + transition: Transition, + /// Each pyramid's base: the hex nodes of its face, oriented into the tet cell. + bases: Vec<[u32; 4]>, + /// The interfaces it makes (reply entries). + pub(super) interfaces: Vec, +} + +/// The base name of a cut face (`partition:P:a` of `partition:P:a#1`): its partition and side. +fn cut_of(name: &str) -> Option<(&str, usize)> { + let base = name.split_once('#').map_or(name, |(b, _)| b); + let (partition, side) = base.strip_prefix("partition:")?.rsplit_once(':')?; + Some((partition, usize::from(side == "b"))) +} + +/// Mesh tet cell `input` with pyramids on the faces of the hex cells `hex` (already in `model`) +/// across its cut faces. `Ok(Err(why))` when there is no such face or the transition cannot be +/// built (the caller meshes the cell plainly and ties it); limits and a cancel are errors. +pub(super) fn mesh_with_pyramids( + design: &fieldcraft_geometry::Design, + model: &Model, + p: &Value, + input: &Input, + hex: &[String], + cdt: &CdtOptions, + options: &QualityOptions, +) -> Result> { + // The hex element faces across each of the cell's cut faces. + let sides = super::sides(model, hex); + let mut cut: Vec = Vec::new(); + let mut quads: Vec<[u32; 4]> = Vec::new(); + let mut interfaces = Vec::new(); + let mut seen: BTreeSet<(String, String)> = BTreeSet::new(); + for (i, name) in input.face_names.iter().enumerate() { + let Some((partition, side)) = cut_of(name) else { continue }; + let Some(across) = sides.get(partition).and_then(|s| s.get(1 - side)) else { continue }; + if across.is_empty() { + continue; + } + cut.push(i as u32); + for other in across { + if !seen.insert((partition.to_string(), other.body.clone())) { + continue; + } + for f in &other.faces { + let nodes = model.face_nodes(*f).unwrap_or_default(); + let Ok(q) = <[u32; 4]>::try_from(nodes.as_slice()) else { + return Ok(Err(format!("`{}` has quadratic or non-quadrilateral faces on the cut; pyramid transitions are linear", other.body))); + }; + quads.push(q); + } + interfaces.push((partition.to_string(), other.body.clone(), other.faces.len())); + } + } + if quads.is_empty() { + return Ok(Err("no hex cell faces it across a partition".into())); + } + // The hex faces must cover the cut faces exactly (a cut shared with a tet cell, or a hex + // cell larger than this one, is tied instead). + let xyz = |n: u32| model.node(n).map(|x| x.xyz).unwrap_or([f64::NAN; 3]); + let positions: Vec<[[f64; 3]; 4]> = quads.iter().map(|q| q.map(xyz)).collect(); + let quad_area: f64 = positions.iter().map(|q| 0.5 * norm(cross(sub(q[2], q[0]), sub(q[3], q[1])))).sum(); + let m = &input.mesh; + let cut_area: f64 = m + .triangles + .iter() + .zip(&m.tri_face) + .filter(|(_, f)| cut.contains(f)) + .map(|(t, _)| { + let c = t.map(|v| m.positions.get(v as usize).copied().unwrap_or([f64::NAN; 3])); + 0.5 * norm(cross(sub(c[1], c[0]), sub(c[2], c[0]))) + }) + .sum(); + if !((quad_area - cut_area).abs() <= 1e-3 * cut_area) { + return Ok(Err(format!("the hex cells' faces ({quad_area:.6} mm²) do not cover its cut faces ({cut_area:.6} mm²)"))); + } + let tag = u32::try_from(input.face_names.len()).map_err(|_| EngineError::Params("too many faces".into()))?; + let welding = super::super::mesh::surface_options(p, m)?; + let built = || -> fieldcraft_mesh::Result<(QualityMesh, Transition, usize, usize, f64)> { + let surface = BoundarySurface::new(&m.positions, &m.triangles, &m.tri_face, &welding)?; + let t = pyramid_surface(&surface, &cut, &positions, tag)?; + // Already welded: nothing more merges. + let exact = fieldcraft_mesh::surface::SurfaceOptions { weld_tolerance: 0.0, weld_distance: 0.0, ..welding }; + let region = BoundarySurface::new(&t.positions, &t.triangles, &t.tri_face, &exact)?; + let tets = tetrahedralize(®ion, cdt)?; + let mut analytic = design.body(&input.name).map(analytic_surfaces).unwrap_or_default(); + analytic.resize(tag as usize, None); + analytic.push(None); + let o = QualityOptions { seeds: input.seeds.clone(), frozen_faces: vec![tag], ..options.clone() }; + let mesh = quality_mesh(®ion, &tets, analytic, &o)?; + Ok((mesh, t, tets.steiner.len(), tets.stats.cavities, region.volume())) + }; + let (mesh, transition, steiner, cavities, volume) = match built() { + Ok(x) => x, + Err(fieldcraft_mesh::Error::Model(fieldcraft_model::Error::Cancelled)) => return Err(EngineError::Model(fieldcraft_model::Error::Cancelled)), + Err(e @ fieldcraft_mesh::Error::Limit(_)) => return Err(EngineError::Params(e.to_string())), + Err(e) => return Ok(Err(e.to_string())), + }; + if let Err(e) = apex_points(&mesh, &transition) { + return Ok(Err(e.to_string())); + } + let tets = BodyTets::new(input.name.clone(), input.face_names.clone(), mesh, steiner, cavities, volume); + let interfaces = interfaces + .into_iter() + .map(|(partition, other, faces)| json!({"partition": partition, "cells": [other, input.name], "joined": "pyramids", "pyramids": faces})) + .collect(); + Ok(Ok(Pyramids { tets, transition, bases: quads, interfaces })) +} + +/// Add the cell's tetrahedra and pyramids to `model`: the tetrahedra's nodes on the hex faces +/// are the hex nodes, the pyramids join them to the apexes. +pub(super) fn add(model: &mut Model, mut c: Pyramids) -> Result<(usize, usize)> { + let tag = c.tets.face_names.len() as u32; + let apexes = apex_points(&c.tets.mesh, &c.transition).map_err(|e| EngineError::Params(e.to_string()))?; + // The pyramid sides are inside the part: no geometry face. + c.tets.mesh.boundary.retain(|b| b.face != tag); + let (first_node, _) = add_tets(model, &c.tets.name, &c.tets.face_names, &c.tets.mesh)?; + // The tet mesh keeps the hex nodes' coordinates exactly: merge its copies into them. + let hex: BTreeMap<[u64; 3], u32> = c.bases.iter().flatten().filter_map(|&n| model.node(n).map(|x| (x.xyz.map(f64::to_bits), n))).collect(); + let mut merge = BTreeMap::new(); + for (i, x) in c.tets.mesh.points.iter().enumerate() { + if let Some(&h) = hex.get(&x.map(f64::to_bits)) { + merge.insert(first_node + i as u32, h); + } + } + let shared = merge.len(); + model.merge_nodes(&merge)?; + let (_, first) = crate::cmd::fresh_labels(model, 0, c.bases.len())?; + let body = model.geometry_tag(&GeometryRef::body(c.tets.name.as_str())?)?; + for (k, (base, apex)) in c.bases.iter().zip(apexes).enumerate() { + let label = first + k as u32; + let mut nodes = base.to_vec(); + nodes.push(first_node + apex); + model.set_element(label, ElementKind::Pyramid5, nodes)?; + model.tag_element(label, body)?; + } + Ok((c.tets.mesh.elements.len() + c.bases.len(), c.tets.mesh.points.len() - shared)) +} + +fn sub(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [a[0] - b[0], a[1] - b[1], a[2] - b[2]] +} + +fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]] +} + +fn norm(a: [f64; 3]) -> f64 { + (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt() +} diff --git a/crates/engine/src/cmd/forms.rs b/crates/engine/src/cmd/forms.rs index 60c18f7..83b77b4 100644 --- a/crates/engine/src/cmd/forms.rs +++ b/crates/engine/src/cmd/forms.rs @@ -1169,6 +1169,9 @@ pub const CELLS: &[Field] = &[ HEX_SIZE, HEX_ORDER, TOLERANCE, + choice("transition", "Hex–tet joint", &[("pyramid", "Pyramids (conforming)"), ("tie", "Tie")]) + .help("How a tet cell meets the hex cells across its partitions: a layer of pyramids on the hex faces, one conforming mesh (linear elements), or a tie between the two meshes.") + .def("\"pyramid\""), ]; pub const HEX_PROJECT: &[Field] = &[ diff --git a/crates/engine/src/cmd/mesh.rs b/crates/engine/src/cmd/mesh.rs index ef2adfa..091b8e9 100644 --- a/crates/engine/src/cmd/mesh.rs +++ b/crates/engine/src/cmd/mesh.rs @@ -255,10 +255,16 @@ pub(super) struct BodyTets { surface_volume: f64, } +impl BodyTets { + pub(super) fn new(name: String, face_names: Vec, mesh: QualityMesh, steiner: usize, cavities: usize, surface_volume: f64) -> Self { + BodyTets { name, face_names, mesh, steiner, cavities, surface_volume } + } +} + /// Mesh one body or imported surface; the model is not touched. With `several` bodies an error /// names the body. pub(super) fn mesh_body( - s: &Session, + design: Option<&fieldcraft_geometry::Design>, p: &Value, input: super::hex::Input, cdt: &CdtOptions, @@ -269,7 +275,7 @@ pub(super) fn mesh_body( let mesh = || -> Result { // The exact faces of a design body; an imported surface keeps to its triangles. let analytic = match p.get("surface") { - None => s.design.as_deref().and_then(|d| d.body(&input.name)).map(analytic_surfaces).unwrap_or_default(), + None => design.and_then(|d| d.body(&input.name)).map(analytic_surfaces).unwrap_or_default(), Some(_) => Vec::new(), }; fieldcraft_model::progress::stage(&format!("{}: boundary", input.name), None, ""); @@ -303,7 +309,7 @@ fn tet_generate(s: &mut Session, p: &Value) -> Result { let several = inputs.len() > 1; let mut meshed = Vec::with_capacity(inputs.len()); for input in inputs { - meshed.push(mesh_body(s, p, input, &cdt_options, &options, quality, several)?); + meshed.push(mesh_body(s.design.as_deref(), p, input, &cdt_options, &options, quality, several)?); } // Replace the bodies' meshes (replace.rs); other bodies and the definitions stay. diff --git a/crates/mesh/src/transition.rs b/crates/mesh/src/transition.rs index b1f776d..03c3a8d 100644 --- a/crates/mesh/src/transition.rs +++ b/crates/mesh/src/transition.rs @@ -4,7 +4,7 @@ //! A hexahedron's face is a quadrilateral and a tetrahedron's a triangle, so a tet mesh cannot //! share a hex mesh's faces. A pyramid can: its base is the quadrilateral and its four sides are //! triangles. [`pyramid_surface`] puts one pyramid on every interface quadrilateral (its apex -//! inside the tet region, about 0.6 of the quadrilateral's mean edge above its centre) and builds +//! inside the tet region, [`APEX_HEIGHT`] of the quadrilateral's mean edge above its centre) and builds //! the surface the tetrahedra fill: the tet region's own surface without the interface faces, //! closed by the pyramids' sides. Where the interface meets the tet region's other faces (the //! rim), those faces' triangles are re-stitched to the hex mesh's rim nodes: each rim vertex of @@ -35,8 +35,10 @@ use crate::tet::QualityMesh; use crate::tet::complex::face_toward; use crate::{Error, Result}; -/// The apex's height above the quadrilateral's centre, over the quadrilateral's mean edge. -pub const APEX_HEIGHT: f64 = 0.6; +/// The apex's height above the quadrilateral's centre, over the quadrilateral's mean edge. Below +/// one half, so the sides of two pyramids on faces that meet at a right angle (a box-shaped tet +/// cell's corner between hex cells) lean apart instead of crossing. +pub const APEX_HEIGHT: f64 = 0.4; /// The tet region's surface for a pyramid transition ([`pyramid_surface`]). #[derive(Debug, Clone)] @@ -294,10 +296,10 @@ mod tests { use crate::surface::SurfaceOptions; use crate::tet::{QualityOptions, quality_mesh}; - /// A box [0, 2] × [0, 2] × [0, 1] as twelve triangles (two per side, outward), sides tagged + /// A box [0, 2] × [0, 2] × [0, height] as twelve triangles (two per side, outward), sides tagged /// 0..6 (−x, +x, −y, +y, −z, +z). - fn slab() -> BoundarySurface { - let p: Vec<[f64; 3]> = (0..8).map(|i| [2.0 * (i & 1) as f64, 2.0 * ((i >> 1) & 1) as f64, ((i >> 2) & 1) as f64]).collect(); + fn slab(height: f64) -> BoundarySurface { + let p: Vec<[f64; 3]> = (0..8).map(|i| [2.0 * (i & 1) as f64, 2.0 * ((i >> 1) & 1) as f64, height * ((i >> 2) & 1) as f64]).collect(); let quads = [([0, 4, 6, 2], 0), ([1, 3, 7, 5], 1), ([0, 1, 5, 4], 2), ([2, 6, 7, 3], 3), ([0, 2, 3, 1], 4), ([4, 5, 7, 6], 5)]; let (mut t, mut f) = (Vec::new(), Vec::new()); for (q, tag) in quads { @@ -322,7 +324,7 @@ mod tests { #[test] fn pyramids_on_a_face_are_kept_by_the_tet_mesh() { - let s = slab(); + let s = slab(1.0); let quads = bottom_quads(); let t = pyramid_surface(&s, &[4], &quads, 6).unwrap(); assert_eq!(t.apexes.len(), 16); @@ -342,16 +344,14 @@ mod tests { #[test] fn a_pyramid_through_a_thin_region_is_refused() { - let s = slab(); - // 1 × 1 quads: apexes 0.6 up, inside the 1 mm slab; one 2 × 2 quad: its apex 1.2 up, outside. + // 1 × 1 quads: apexes 0.4 up, inside a 1 mm slab, through a 0.3 mm one. let quads: Vec<[[f64; 3]; 4]> = (0..4) .map(|k| ((k % 2) as f64, (k / 2) as f64)) .map(|(x, y)| [[x, y, 0.0], [x + 1.0, y, 0.0], [x + 1.0, y + 1.0, 0.0], [x, y + 1.0, 0.0]]) .collect(); - let t = pyramid_surface(&s, &[4], &quads, 6).unwrap(); - let big: Vec<[[f64; 3]; 4]> = vec![[[0.0, 0.0, 0.0], [2.0, 0.0, 0.0], [2.0, 2.0, 0.0], [0.0, 2.0, 0.0]]]; - let tall = pyramid_surface(&s, &[4], &big, 6).unwrap(); - assert!(BoundarySurface::new(&t.positions, &t.triangles, &t.tri_face, &SurfaceOptions::default()).is_ok()); - assert!(BoundarySurface::new(&tall.positions, &tall.triangles, &tall.tri_face, &SurfaceOptions::default()).is_err()); + for (height, ok) in [(1.0, true), (0.3, false)] { + let t = pyramid_surface(&slab(height), &[4], &quads, 6).unwrap(); + assert_eq!(BoundarySurface::new(&t.positions, &t.triangles, &t.tri_face, &SurfaceOptions::default()).is_ok(), ok, "{height}"); + } } }