explicit: symmetry about inclined planes (#242)

At an inclined plane's nodes the velocity along its normal is removed every
increment, from the components no prescribed velocity holds (planes meeting
at a node orthonormalized, so one pass is exact); a prescribed motion off the
plane is refused.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo

crates/explicit/src/integrate.rs

@@ -12,8 +12,10 @@
 //!    condition: T. J. R. Hughes, R. L. Taylor, J. L. Sackman, A. Curnier, W. Kanoknukulchai,
 //!    *A finite element method for a class of contact-impact problems*, CMAME 8(3):249–276,
 //!    1976), with Coulomb friction as a tangential impulse at most μ times the normal one (the
-//!    tangential velocity stops, or loses μ Δv_n when it slides), and prescribed degrees of
-//!    freedom take their velocity. The change a constraint makes is the impulse of its reaction,
+//!    tangential velocity stops, or loses μ Δv_n when it slides), prescribed degrees of
+//!    freedom take their velocity, and a symmetry plane not normal to an axis removes its nodes'
+//!    velocity along its normal (from their free components: the lumped mass is the same on a
+//!    node's three translations, so this is the mass-weighted projection onto the constraint). The change a constraint makes is the impulse of its reaction,
 //!    r^n = M Δv / Δt^n;
 //! 4. energies: the works of the internal, contact and external forces by the trapezoidal rule,
 //!    W += ½ Δu · (f^{n−1} + f^n) (exactly the strain energy for a linear elastic material),
@@ -77,7 +79,7 @@ struct Output {
 }
 
 /// Velocity constraints over an increment of length `dt` (rigid walls when dt > 0, then
-/// prescribed degrees of freedom); returns the kinetic energy the walls' friction took (`mass`
+/// prescribed degrees of freedom, then inclined symmetry planes); returns the kinetic energy the walls' friction took (`mass`
 /// per degree of freedom).
 fn constrain(sys: &System, u: &[f64], v: &mut [f64], dt: f64, mass: &[f64]) -> f64 {
     let mut dissipated = 0.0;
@@ -107,6 +109,9 @@ fn constrain(sys: &System, u: &[f64], v: &mut [f64], dt: f64, mass: &[f64]) -> f
             *c = *p;
         }
     }
+    for row in &sys.planes {
+        row.project(v);
+    }
     dissipated
 }
 

crates/explicit/src/system.rs

@@ -9,7 +9,8 @@ use fieldcraft_fem::plasticity::J2;
 use fieldcraft_fem::uniform::{Hex8, Tet4, modulus_bound};
 use fieldcraft_fem::{DofKind, DofLayout, DofSet, face};
 use fieldcraft_model::{
-    Elastic, ElementKind, ExplicitDynamic, Friction, HourglassKind, Integration, Load, Model, Resolved, SectionKind, Step, StressRate, Target,
+    BoundaryCondition, Elastic, ElementKind, ExplicitDynamic, Friction, HourglassKind, Integration, Load, Model, Resolved, SectionKind, Step,
+    StressRate, Target,
 };
 
 use crate::contact::Interface;
@@ -100,6 +101,36 @@ pub struct Wall {
     pub friction: Friction,
 }
 
+/// One velocity constraint of the symmetry planes not normal to an axis at a node: c · v = 0
+/// over the node's translations, `q` its part on the components no prescribed velocity holds,
+/// of unit length and perpendicular to the node's other rows' (Gram–Schmidt over the planes that
+/// meet at the node). [`SymmetryRow::project`] removes q (c · v) from the velocity: with the same
+/// lumped mass on a node's three translations this is the mass-weighted closest velocity that
+/// satisfies the planes and the prescribed values (the constrained update of a lumped-mass
+/// central difference scheme), exactly, in one pass over the rows.
+#[derive(Debug, Clone)]
+pub struct SymmetryRow {
+    pub node: usize,
+    pub q: [f64; 3],
+    pub c: [f64; 3],
+}
+
+impl SymmetryRow {
+    pub fn project(&self, v: &mut [f64]) {
+        let at = |k: usize| NODE_DOFS * self.node + k;
+        let along: f64 = (0..3).map(|k| self.c[k] * v.get(at(k)).copied().unwrap_or(0.0)).sum();
+        for k in 0..3 {
+            if let Some(x) = v.get_mut(at(k)) {
+                *x -= along * self.q[k];
+            }
+        }
+    }
+}
+
+/// A free part below this fraction of its row's size is round-off: the row repeats the node's
+/// earlier ones, or the prescribed velocities hold it.
+const REPEATED: f64 = 1e-10;
+
 /// The discrete system of an explicit dynamic step.
 ///
 /// Its elements are solids, so every node carries the three translations and the layout is
@@ -130,6 +161,8 @@ pub struct System {
     /// Prescribed velocity per degree of freedom (`None` = free).
     pub prescribed: Vec<Option<f64>>,
     pub walls: Vec<Wall>,
+    /// The symmetry planes not normal to an axis, per node.
+    pub planes: Vec<SymmetryRow>,
     /// The model's contacts ([`crate::contact`]).
     pub contacts: Vec<Interface>,
     /// Velocity per degree of freedom at the start (prescribed velocities included).
@@ -293,6 +326,7 @@ impl System {
             discrete_increment: None,
             prescribed: vec![None; dofs],
             walls: Vec::new(),
+            planes: Vec::new(),
             contacts: Vec::new(),
             initial_velocity: vec![0.0; dofs],
             external: vec![0.0; dofs],
@@ -541,6 +575,7 @@ impl System {
     /// Prescribed velocities from the step's boundary conditions: zero for held degrees of
     /// freedom, value / end time for a displacement reached over the step.
     fn prescribe(&mut self, model: &Model, resolved: &Resolved, step: &Step, end_time: f64) -> Result<()> {
+        let mut normals: std::collections::BTreeMap<usize, Vec<[f64; 3]>> = std::collections::BTreeMap::new();
         for name in &step.boundary_conditions {
             let bc = model.boundary_condition(name).ok_or_else(|| invalid(format!("boundary condition `{name}` is missing")))?;
             if let Target::ReferencePoint(rp) = bc.target() {
@@ -549,6 +584,14 @@ impl System {
             if bc.is_thermal() {
                 return Err(Error::NotSupported(format!("temperature boundary conditions (`{name}`) in explicit dynamics")));
             }
+            if let BoundaryCondition::Symmetry { plane, .. } = bc
+                && let Some(normal) = plane.inclined()
+            {
+                for node in resolved.target_nodes(bc.target()).iter().filter_map(|l| self.node(*l)) {
+                    normals.entry(node).or_default().push(normal);
+                }
+                continue;
+            }
             for &label in resolved.target_nodes(bc.target()) {
                 let Some(node) = self.node(label) else { continue };
                 for kind in self.layout.kinds(node).iter() {
@@ -571,6 +614,34 @@ impl System {
                 }
             }
         }
+        // Each node's planes: orthonormal free parts, or rows the prescribed velocities satisfy.
+        for (node, rows) in normals {
+            let p = [0, 1, 2].map(|k| self.prescribed.get(NODE_DOFS * node + k).copied().flatten());
+            let mut done: Vec<SymmetryRow> = Vec::new();
+            for n in rows {
+                let mut c = n;
+                for r in &done {
+                    let f: f64 = (0..3).map(|k| r.q[k] * if p[k].is_none() { c[k] } else { 0.0 }).sum();
+                    for k in 0..3 {
+                        c[k] -= f * r.c[k];
+                    }
+                }
+                let q = [0, 1, 2].map(|k| if p[k].is_none() { c[k] } else { 0.0 });
+                let length = q.iter().map(|x| x * x).sum::<f64>().sqrt();
+                if length > REPEATED {
+                    done.push(SymmetryRow { node, q: q.map(|x| x / length), c: c.map(|x| x / length) });
+                    continue;
+                }
+                // Nothing free is left: the prescribed velocities must satisfy the row.
+                let along: f64 = (0..3).map(|k| c[k] * p[k].unwrap_or(0.0)).sum();
+                let size: f64 = (0..3).map(|k| (n[k] * p[k].unwrap_or(0.0)).abs()).sum();
+                if along.abs() > 1e-9 * size {
+                    let label = self.labels.get(node).copied().unwrap_or_default();
+                    return Err(invalid(format!("node {label} moves off a symmetry plane (normal {n:?}) by its prescribed velocities")));
+                }
+            }
+            self.planes.extend(done);
+        }
         Ok(())
     }
 
@@ -729,12 +800,15 @@ impl System {
                 }
             }
         }
-        // A prescribed velocity holds from the start.
+        // A prescribed velocity holds from the start, and so do the symmetry planes.
         for (v, p) in self.initial_velocity.iter_mut().zip(&self.prescribed) {
             if let Some(p) = p {
                 *v = *p;
             }
         }
+        for row in &self.planes {
+            row.project(&mut self.initial_velocity);
+        }
         let size = self.xyz.iter().flatten().fold(0.0f64, |m, c| m.max(c.abs())).max(1.0);
         for w in &settings.rigid_walls {
             let len = w.normal.iter().map(|c| c * c).sum::<f64>().sqrt();

crates/explicit/tests/inclined_symmetry.rs

@@ -0,0 +1,164 @@
+//! Symmetry about inclined planes in explicit dynamics (#242): the velocity along the plane's
+//! normal is removed at its nodes every increment. A bar turned by a rotation R moves as R times
+//! the unturned bar on the axis-aligned planes; and a half (or quarter) of a bar whose mesh is
+//! mirror-symmetric moves as the whole at the nodes they share. One plane, and two planes meeting
+//! at an edge.
+
+use fieldcraft_explicit::{ExplicitResult, Limits, solve_explicit};
+use fieldcraft_model::*;
+
+/// A rotation by `angle` about the axis `a` (Rodrigues).
+fn rotation(a: [f64; 3], angle: f64) -> [[f64; 3]; 3] {
+    let l = (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt();
+    let a = a.map(|c| c / l);
+    let (s, c) = angle.sin_cos();
+    let k = [[0.0, -a[2], a[1]], [a[2], 0.0, -a[0]], [-a[1], a[0], 0.0]];
+    let mut r = [[0.0; 3]; 3];
+    for i in 0..3 {
+        for j in 0..3 {
+            r[i][j] = if i == j { c } else { 0.0 } + (1.0 - c) * a[i] * a[j] + s * k[i][j];
+        }
+    }
+    r
+}
+
+fn apply(r: &[[f64; 3]; 3], v: [f64; 3]) -> [f64; 3] {
+    [0, 1, 2].map(|i| (0..3).map(|j| r[i][j] * v[j]).sum())
+}
+
+fn norm(v: [f64; 3]) -> f64 {
+    v.iter().map(|c| c * c).sum::<f64>().sqrt()
+}
+
+/// A 6 × 1 × 1 bar of one-point hex8 cells along x, clamped at x = 0, under the body force
+/// `force` for 0.2 time units. `whole[a]` spans −1..1 across axis a (y or z) with two cells,
+/// otherwise 0..1 with one cell and a symmetry plane at 0. Everything is turned by `turn`.
+fn bar(whole: [bool; 2], force: [f64; 3], turn: &[[f64; 3]; 3]) -> (Model, ExplicitResult) {
+    let mut m = Model::default();
+    let cells = [6usize, if whole[0] { 2 } else { 1 }, if whole[1] { 2 } else { 1 }];
+    let start = [0.0, if whole[0] { -1.0 } else { 0.0 }, if whole[1] { -1.0 } else { 0.0 }];
+    let label = |i: usize, j: usize, k: usize| (1 + i + (cells[0] + 1) * (j + (cells[1] + 1) * k)) as u32;
+    for k in 0..=cells[2] {
+        for j in 0..=cells[1] {
+            for i in 0..=cells[0] {
+                m.set_node(label(i, j, k), apply(turn, [i as f64, start[1] + j as f64, start[2] + k as f64])).unwrap();
+            }
+        }
+    }
+    let corners = [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]];
+    let mut next = 1;
+    for k in 0..cells[2] {
+        for j in 0..cells[1] {
+            for i in 0..cells[0] {
+                m.set_element(next, ElementKind::Hex8, corners.iter().map(|c| label(i + c[0], j + c[1], k + c[2])).collect()).unwrap();
+                next += 1;
+            }
+        }
+    }
+    let all: Vec<u32> = m.elements().keys().copied().collect();
+    m.add_set("all", Set { elements: all, ..Set::default() }).unwrap();
+    let nodes = |f: &dyn Fn(usize, usize, usize) -> bool| -> Vec<u32> {
+        let mut out = Vec::new();
+        for k in 0..=cells[2] {
+            for j in 0..=cells[1] {
+                for i in 0..=cells[0] {
+                    if f(i, j, k) {
+                        out.push(label(i, j, k));
+                    }
+                }
+            }
+        }
+        out.sort_unstable();
+        out
+    };
+    m.add_set("root", Set { nodes: nodes(&|i, _, _| i == 0), ..Set::default() }).unwrap();
+    m.add_set("tip", Set { nodes: nodes(&|i, _, _| i == cells[0]), ..Set::default() }).unwrap();
+    m.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: 1.0e4, poisson: 0.3 }), density: Some(1.0), ..Material::default() })
+        .unwrap();
+    let formulation = Formulation { integration: Integration::Reduced, ..Formulation::default() };
+    m.add_section("bar", Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation }).unwrap();
+    let settings = ExplicitDynamic::new(0.2);
+    let mut step = Step::new("Pull", Procedure::ExplicitDynamic(settings));
+    m.add_boundary_condition("clamp", BoundaryCondition::Fixed { target: Target::Set("root".into()) }).unwrap();
+    step.boundary_conditions.insert("clamp".into());
+    for (a, name) in [(0usize, "y0"), (1, "z0")] {
+        if whole[a] {
+            continue;
+        }
+        m.add_set(name, Set { nodes: nodes(&|_, j, k| [j, k][a] == 0), ..Set::default() }).unwrap();
+        let mut normal = [0.0; 3];
+        normal[a + 1] = 1.0;
+        let plane = Symmetry::from_normal(apply(turn, normal)).unwrap();
+        m.add_boundary_condition(name, BoundaryCondition::Symmetry { target: Target::Set(name.into()), plane }).unwrap();
+        step.boundary_conditions.insert(name.into());
+    }
+    m.add_load("pull", Load::BodyForce { set: "all".into(), vector: apply(turn, force) }).unwrap();
+    step.loads.insert("pull".into());
+    m.add_step(step, None).unwrap();
+    let r = solve_explicit(&m, None, &Limits::default()).unwrap();
+    (m, r)
+}
+
+fn key(x: [f64; 3]) -> [i64; 3] {
+    x.map(|c| (c * 1e6).round() as i64)
+}
+
+/// The largest difference between `b`'s last-frame displacements and `turn` times `a`'s at the
+/// same (unturned) places, over the nodes of `b`, relative to `a`'s largest.
+fn difference(a: &(Model, ExplicitResult), b: &(Model, ExplicitResult), turn: &[[f64; 3]; 3]) -> f64 {
+    let ua = &a.1.frames.last().unwrap().displacement;
+    let at: std::collections::BTreeMap<[i64; 3], [f64; 3]> = a.1.nodes.iter().zip(ua).map(|(l, u)| (key(a.0.node(*l).unwrap().xyz), *u)).collect();
+    let size = ua.iter().fold(0.0f64, |m, u| m.max(norm(*u)));
+    let unturn = |x: [f64; 3]| [0, 1, 2].map(|i| (0..3).map(|j| turn[j][i] * x[j]).sum::<f64>());
+    let ub = &b.1.frames.last().unwrap().displacement;
+    b.1.nodes
+        .iter()
+        .zip(ub)
+        .map(|(l, u)| {
+            let want = apply(turn, at[&key(unturn(b.0.node(*l).unwrap().xyz))]);
+            norm([0, 1, 2].map(|k| u[k] - want[k])) / size
+        })
+        .fold(0.0, f64::max)
+}
+
+#[test]
+fn turned_planes_move_as_the_axis_aligned_ones_and_the_whole() {
+    let identity = rotation([0.0, 0.0, 1.0], 0.0);
+    let turn = rotation([1.0, 2.0, 3.0], 0.7);
+    // Bending about one plane; an axial pull and a bend about both, two planes meeting at an edge.
+    for (whole, force) in [([false, true], [0.0, 0.0, -1.0]), ([false, false], [4.0, 0.0, 0.0])] {
+        let axis = bar(whole, force, &identity);
+        let turned = bar(whole, force, &turn);
+        let full = bar([true, true], force, &identity);
+        let full_turned = bar([true, true], force, &turn);
+        // The part moves as the whole, turned or not, to round-off: the planes hold exactly.
+        let (mirrored, mirrored_turned) = (difference(&full, &axis, &identity), difference(&full_turned, &turned, &identity));
+        // Turning the model changes the stable increment in its last digits, which the central
+        // difference scheme carries over its increments: the turned part may differ from the
+        // unturned by as much as the turned whole does (its increments may end a step apart), no more.
+        let (rotated, own) = (difference(&axis, &turned, &turn), difference(&full, &full_turned, &turn));
+        let tip = axis.1.frames.last().unwrap().displacement.iter().fold(0.0f64, |m, u| m.max(norm(*u)));
+        eprintln!(
+            "planes {whole:?}: tip {tip:.4e}; part vs whole {mirrored:.1e} (turned {mirrored_turned:.1e}); turned vs axis-aligned {rotated:.1e} (the whole: {own:.1e})"
+        );
+        assert!(tip > 1e-4, "{tip}");
+        assert!(mirrored < 1e-12 && mirrored_turned < 1e-12, "{mirrored} {mirrored_turned}");
+        assert!(rotated < 1.01 * own + 1e-12 && own < 1e-3, "{rotated} {own}");
+        assert!(turned.1.completed);
+    }
+}
+
+/// A prescribed velocity that leaves no free component along the normal must move along the
+/// plane; one that does not is refused.
+#[test]
+fn a_prescribed_motion_off_the_plane_is_refused() {
+    let turn = rotation([0.0, 0.0, 1.0], 0.5);
+    let (mut m, _) = bar([false, true], [0.0, 0.0, -1.0], &turn);
+    let nodes = m.set("y0").unwrap().nodes.iter().copied().filter(|l| m.set("tip").unwrap().nodes.contains(l)).collect();
+    m.add_set("pushed", Set { nodes, ..Set::default() }).unwrap();
+    let push = BoundaryCondition::Displacement { target: Target::Set("pushed".into()), values: [Some(0.01), Some(0.01), None, None, None, None] };
+    m.add_boundary_condition("push", push).unwrap();
+    m.set_boundary_condition_active("Pull", "push", true).unwrap();
+    let e = solve_explicit(&m, None, &Limits::default()).unwrap_err().to_string();
+    assert!(e.contains("moves off a symmetry plane"), "{e}");
+}