Mesh: merge coincident nodes, keeping tags and sets on the survivor

Model::merge_nodes replaces merged nodes by their survivors in elements (node
order kept, so face and edge tags hold) and in set node lists, moves a vertex
tag to an untagged survivor, and refuses a merge that would collapse an
element. fieldcraft_mesh::join::coincident_nodes / merge_coincident find the
nodes within a distance of a lower-labelled survivor, deterministically.
Meshing cell by cell and mirroring use them; node merging commands will too.

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo

crates/mesh/src/join.rs

@@ -1,7 +1,8 @@
 //! Joining new elements to the mesh already in a model: which new points coincide with existing
 //! nodes ([`weld`]), and whether the new elements then meet the mesh face to face
-//! ([`check_conforming`]). Mesh generators that add to a model (`mesh.block`) use both; node
-//! merging and orphan-mesh import will too.
+//! ([`check_conforming`]). Mesh generators that add to a model (`mesh.block`) use both. Nodes that
+//! coincide are merged by [`coincident_nodes`] and [`merge_coincident`] (`mesh.cells`,
+//! `mesh.mirror`); orphan-mesh import will use them too.
 
 use std::collections::{BTreeMap, BTreeSet, HashMap};
 
@@ -91,6 +92,56 @@ pub fn weld(model: &Model, points: &[P], tol: f64) -> Result<Vec<Option<u32>>> {
         .collect()
 }
 
+/// Most nodes [`coincident_nodes`] looks at.
+pub const MAX_MERGE_NODES: usize = 50_000_000;
+
+/// The nodes among `candidates` (every node when `None`) that lie within `tol` of a survivor,
+/// each mapped to it. Nodes are taken in ascending label order; one with no survivor within reach
+/// becomes a survivor (the nearest wins, the lowest label on a tie). Deterministic: the labels
+/// decide, never the order of a hash map.
+pub fn coincident_nodes(model: &Model, candidates: Option<&BTreeSet<u32>>, tol: f64) -> Result<BTreeMap<u32, u32>> {
+    if !(tol.is_finite() && tol > 0.0) {
+        return Err(Error::Invalid(format!("the merge distance must be a positive number, not {tol}")));
+    }
+    let nodes: Vec<(u32, P)> = match candidates {
+        Some(c) => c.iter().filter_map(|&l| model.node(l).map(|n| (l, n.xyz))).collect(),
+        None => model.nodes().iter().map(|(&l, n)| (l, n.xyz)).collect(),
+    };
+    if nodes.len() > MAX_MERGE_NODES {
+        return Err(Error::Invalid(format!("{} nodes to merge; at most {MAX_MERGE_NODES}", nodes.len())));
+    }
+    let mut grid: Grid<(u32, P)> = Grid::new(tol);
+    let mut map = BTreeMap::new();
+    for (label, p) in nodes {
+        // The survivor of the nearest earlier node within reach (the lowest label on a tie).
+        let near = grid
+            .near(p)
+            .map(|(l, q)| (dot(sub(q, p), sub(q, p)), l))
+            .filter(|(d2, _)| *d2 <= tol * tol)
+            .min_by(|a, b| a.0.total_cmp(&b.0).then(a.1.cmp(&b.1)))
+            .map(|(_, l)| map.get(&l).copied().unwrap_or(l));
+        match near {
+            Some(survivor) => {
+                map.insert(label, survivor);
+            }
+            None => grid.insert(p, p, (label, p)),
+        }
+    }
+    Ok(map)
+}
+
+/// Merge the nodes among `candidates` (every node when `None`) that coincide within `tol`
+/// ([`coincident_nodes`]) into their survivors ([`Model::merge_nodes`]: elements, set node lists
+/// and vertex tags follow the survivor); refused, with nothing changed, when that would collapse
+/// an element. Returns the merged nodes and their survivors.
+pub fn merge_coincident(model: &mut Model, candidates: Option<&BTreeSet<u32>>, tol: f64) -> Result<BTreeMap<u32, u32>> {
+    let map = coincident_nodes(model, candidates, tol)?;
+    if !map.is_empty() {
+        model.merge_nodes(&map).map_err(|e| Error::Invalid(e.to_string()))?;
+    }
+    Ok(map)
+}
+
 /// A face or edge by its sorted node labels.
 fn key(nodes: impl Iterator<Item = u32>) -> Vec<u32> {
     let mut k: Vec<u32> = nodes.collect();
@@ -196,3 +247,47 @@ pub fn check_conforming(model: &Model, new: &[(ElementKind, Vec<u32>)], joined:
     }
     Ok(())
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// Two unit hex8 cubes side by side, each with its own nodes: the four on x = 1 coincide.
+    fn two_cubes(gap: f64) -> Model {
+        let mut m = Model::default();
+        let mut label = 0;
+        for (e, x0) in [(1u32, 0.0), (2, 1.0 + gap)] {
+            let first = label + 1;
+            // Hex8 order: the bottom face counter-clockwise, then the top.
+            for [x, y, z] in [[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [1., 1., 1.], [0., 1., 1.]] {
+                label += 1;
+                m.set_node(label, [x0 + x, y, z]).unwrap();
+            }
+            m.set_element(e, ElementKind::Hex8, (first..first + 8).collect()).unwrap();
+        }
+        m
+    }
+
+    #[test]
+    fn coincident_nodes_merge_into_the_lower_labels() {
+        let mut m = two_cubes(0.0);
+        let map = merge_coincident(&mut m, None, 1e-9).unwrap();
+        assert_eq!(map.len(), 4);
+        assert!(map.iter().all(|(from, to)| from > to && *from > 8 && *to <= 8));
+        assert_eq!(m.nodes().len(), 12);
+        // A gap wider than the distance merges nothing; candidates limit the merge.
+        let mut apart = two_cubes(1e-3);
+        assert!(merge_coincident(&mut apart, None, 1e-6).unwrap().is_empty());
+        let mut some = two_cubes(0.0);
+        let only: BTreeSet<u32> = [2, 10].into();
+        assert_eq!(merge_coincident(&mut some, Some(&only), 1e-9).unwrap().len(), 0, "2 and 10 are not the same corner");
+        // A distance that would put two corners of one element together is refused.
+        let mut whole = two_cubes(0.0);
+        let before = whole.clone();
+        assert!(merge_coincident(&mut whole, None, 1.5).is_err());
+        assert_eq!(whole, before);
+        for bad in [0.0, -1.0, f64::NAN, f64::INFINITY] {
+            assert!(coincident_nodes(&m, None, bad).is_err());
+        }
+    }
+}

crates/model/src/geometry.rs

@@ -244,6 +244,11 @@ impl Classification {
         self.elements.remove(&label);
     }
 
+    /// Tag a node with a tag already interned (the caller checked both).
+    pub(crate) fn tag_node_unchecked(&mut self, label: u32, tag: u32) {
+        self.nodes.insert(label, tag);
+    }
+
     pub(crate) fn forget_node(&mut self, label: u32) {
         self.nodes.remove(&label);
     }

crates/model/src/mesh.rs

@@ -341,6 +341,61 @@ impl Model {
         Ok(())
     }
 
+    /// Merge nodes: each node `from` in `map` is replaced by the node `map[from]` (its survivor)
+    /// in every element and in the node lists of sets, then removed. The survivor keeps its
+    /// coordinates; it takes the merged node's geometry vertex tag when it has none. Survivors
+    /// must exist and not be merged themselves. Refused, with nothing changed, when an element
+    /// would hold one node twice (the merge would collapse it).
+    pub fn merge_nodes(&mut self, map: &BTreeMap<u32, u32>) -> Result<()> {
+        for (&from, &to) in map {
+            if from == to {
+                return Err(Error::Invalid(format!("node {from} cannot be merged into itself")));
+            }
+            if !self.nodes.contains_key(&from) {
+                return Err(Error::NotFound(format!("node {from}")));
+            }
+            if !self.nodes.contains_key(&to) || map.contains_key(&to) {
+                return Err(Error::Invalid(format!("node {to}, the survivor of node {from}, must be a node that is not merged")));
+            }
+        }
+        let target = |n: u32| map.get(&n).copied().unwrap_or(n);
+        let mut renamed: Vec<(u32, Vec<u32>)> = Vec::new();
+        for (&label, e) in &self.elements {
+            if !e.nodes.iter().any(|n| map.contains_key(n)) {
+                continue;
+            }
+            let nodes: Vec<u32> = e.nodes.iter().map(|&n| target(n)).collect();
+            if let Some(n) = nodes.iter().enumerate().find_map(|(i, n)| nodes.get(..i).filter(|b| b.contains(n)).map(|_| n)) {
+                return Err(Error::Invalid(format!("merging would collapse element {label}: node {n} would appear twice")));
+            }
+            renamed.push((label, nodes));
+        }
+        for (label, nodes) in renamed {
+            // The node order is kept, so the element's face and edge tags still hold.
+            if let Some(e) = self.elements.get_mut(&label) {
+                e.nodes = nodes;
+            }
+        }
+        for set in self.sets.values_mut() {
+            if set.nodes.iter().any(|n| map.contains_key(n)) {
+                let mut nodes: Vec<u32> = set.nodes.iter().map(|&n| target(n)).collect();
+                nodes.sort_unstable();
+                nodes.dedup();
+                set.nodes = nodes;
+            }
+        }
+        for (&from, &to) in map {
+            if let Some(tag) = self.classification.nodes().get(&from).copied()
+                && !self.classification.nodes().contains_key(&to)
+            {
+                self.classification.tag_node_unchecked(to, tag);
+            }
+            self.classification.forget_node(from);
+            self.nodes.remove(&from);
+        }
+        Ok(())
+    }
+
     /// Element counts by kind, in [`ElementKind::ALL`] order, kinds with none left out.
     pub fn element_counts(&self) -> Vec<(ElementKind, usize)> {
         ElementKind::ALL.into_iter().map(|k| (k, self.elements.values().filter(|e| e.kind == k).count())).filter(|(_, n)| *n > 0).collect()
@@ -470,4 +525,45 @@ mod tests {
         }
         assert_eq!(ElementKind::from_name("C3D8"), None);
     }
+
+    #[test]
+    fn merging_nodes_rewrites_elements_sets_and_tags_and_refuses_a_collapse() {
+        use crate::{GeometryRef, Set};
+        let mut m = Model::default();
+        // Two tet4 sharing a face whose nodes are duplicated (4, 5, 6 copy 1, 2, 3).
+        for (l, p) in [
+            (1, [0., 0., 0.]),
+            (2, [1., 0., 0.]),
+            (3, [0., 1., 0.]),
+            (7, [0., 0., 1.]),
+            (4, [0., 0., 0.]),
+            (5, [1., 0., 0.]),
+            (6, [0., 1., 0.]),
+            (8, [0., 0., -1.]),
+        ] {
+            m.set_node(l, p).unwrap();
+        }
+        m.set_element(1, ElementKind::Tet4, vec![1, 2, 3, 7]).unwrap();
+        m.set_element(2, ElementKind::Tet4, vec![4, 6, 5, 8]).unwrap();
+        let face = m.geometry_tag(&GeometryRef::face("B", "bottom").unwrap()).unwrap();
+        m.tag_element_face(2, 1, face).unwrap();
+        let vertex = m.geometry_tag(&GeometryRef::new("B", crate::geometry::Entity::Vertex("v".into())).unwrap()).unwrap();
+        m.tag_node(4, vertex).unwrap();
+        m.add_set("pins", Set { nodes: vec![1, 4, 5], ..Set::default() }).unwrap();
+
+        let collapse: BTreeMap<u32, u32> = [(8, 7), (7, 1)].into();
+        assert!(m.merge_nodes(&collapse).is_err());
+        assert!(m.merge_nodes(&[(1, 2)].into()).is_err(), "collapses element 1");
+        assert!(m.merge_nodes(&[(9, 1)].into()).is_err());
+        let before = m.clone();
+        assert!(m.merge_nodes(&[(4, 1), (5, 2), (6, 6)].into()).is_err());
+        assert_eq!(m, before, "a refused merge changes nothing");
+
+        m.merge_nodes(&[(4, 1), (5, 2), (6, 3)].into()).unwrap();
+        assert_eq!(m.nodes().len(), 5);
+        assert_eq!(m.element(2).unwrap().nodes, vec![1, 3, 2, 8]);
+        assert_eq!(m.set("pins").unwrap().nodes, vec![1, 2]);
+        assert_eq!(m.classification().nodes().get(&1), Some(&vertex));
+        assert_eq!(m.classification().elements().get(&2).and_then(|t| t.faces.get(&1)), Some(&face));
+    }
 }