diff --git a/crates/explicit/src/integrate.rs b/crates/explicit/src/integrate.rs index 7e8a866..7a9fefc 100644 --- a/crates/explicit/src/integrate.rs +++ b/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 } diff --git a/crates/explicit/src/system.rs b/crates/explicit/src/system.rs index ac2947c..fa30156 100644 --- a/crates/explicit/src/system.rs +++ b/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>, pub walls: Vec, + /// The symmetry planes not normal to an axis, per node. + pub planes: Vec, /// The model's contacts ([`crate::contact`]). pub contacts: Vec, /// 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> = 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 = 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::().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::().sqrt(); diff --git a/crates/explicit/tests/inclined_symmetry.rs b/crates/explicit/tests/inclined_symmetry.rs new file mode 100644 index 0000000..ef8938e --- /dev/null +++ b/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::().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 = 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 { + 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::()); + 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}"); +}