diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs index 010be25..0a3893b 100644 --- a/crates/decks/src/abaqus.rs +++ b/crates/decks/src/abaqus.rs @@ -148,7 +148,7 @@ use std::fmt::Write as _; use fieldcraft_model::{ BeamProfile, BoundaryCondition, ContactSurfaces, CouplingKind, Elastic, Element, ElementKind, Hardening, Hyperelastic, InitialPenetration, - Integration, Load, Model, OutputVariable, Plastic, Procedure, Resolved, Section, SectionKind, Step, Symmetry, Target, + Integration, Load, Model, OutputVariable, Plastic, Procedure, Resolved, Section, SectionKind, Step, Target, }; use crate::{Error, Result, unsupported}; @@ -193,6 +193,7 @@ pub fn write_resolved(model: &Model, r: &Resolved) -> Result { d.couplings(model, r); d.contacts(model, r)?; d.ties(model)?; + d.equations(model, r)?; d.steps(model, r)?; Ok(d.out) } @@ -370,9 +371,36 @@ fn check_expressible(model: &Model) -> Result<()> { bad.insert(format!("coupling `{name}`: an Abaqus distributing coupling always couples translations 1-3")); } } + for name in inclined_symmetries(model) { + if !model.steps().iter().all(|s| s.boundary_conditions.contains(name)) { + bad.insert(format!( + "boundary condition `{name}`: a symmetry plane not normal to X, Y or Z active in some steps only (it is written as `*EQUATION`, which is model data and holds in every step)" + )); + } + } unsupported(FORMAT, bad) } +/// The symmetry conditions on inclined planes some step uses: linear equations rather than held +/// degrees of freedom. +fn inclined_symmetries(model: &Model) -> BTreeSet<&String> { + model + .boundary_conditions() + .iter() + .filter(|(n, bc)| { + matches!(bc, BoundaryCondition::Symmetry { plane, .. } if plane.axis().is_none()) && !model.boundary_condition_users(n).is_empty() + }) + .map(|(n, _)| n) + .collect() +} + +/// The boundary conditions of `names` that `*BOUNDARY` writes (inclined symmetry planes are +/// `*EQUATION`s). +fn held_conditions(model: &Model, names: &BTreeSet) -> BTreeSet { + let inclined = inclined_symmetries(model); + names.iter().filter(|n| !inclined.contains(n)).cloned().collect() +} + /// The property of a group of discrete elements. enum DiscreteProperty { Mass(f64), @@ -1076,10 +1104,12 @@ impl Deck { match bc { BoundaryCondition::Fixed { .. } => self.line(&format!("{target}, ENCASTRE")), BoundaryCondition::Symmetry { plane, .. } => { - let s = match plane { - Symmetry::X => "XSYMM", - Symmetry::Y => "YSYMM", - Symmetry::Z => "ZSYMM", + let s = match plane.axis() { + Some(0) => "XSYMM", + Some(1) => "YSYMM", + Some(2) => "ZSYMM", + // An inclined plane is an `*EQUATION` (model data), never here. + _ => continue, }; self.line(&format!("{target}, {s}")); } @@ -1097,8 +1127,74 @@ impl Deck { Ok(()) } + /// The inclined symmetry planes: one `*EQUATION` per node and held direction (n · u = 0, and + /// the rotations about the plane's axes at shell and beam nodes), whose first term Abaqus + /// eliminates; it is chosen among the degrees of freedom no boundary condition, coupling + /// or tie constrains ([`crate::symmetry_equations`]). + fn equations(&mut self, model: &Model, r: &Resolved) -> Result<()> { + let inclined = inclined_symmetries(model); + if inclined.is_empty() { + return Ok(()); + } + let steps = model.steps(); + // Per node and degree of freedom: the value it is held at in each step that holds it. + let mut held: BTreeMap<(u32, u8), Vec> = BTreeMap::new(); + for step in steps { + for name in held_conditions(model, &step.boundary_conditions) { + let Some(bc) = model.boundary_condition(&name) else { continue }; + for &node in r.target_nodes(bc.target()) { + for dof in bc.dofs().list() { + let v = if matches!(step.procedure, Procedure::Frequency { .. } | Procedure::Buckle { .. }) { + 0.0 + } else { + bc.value(dof).unwrap_or(0.0) + }; + held.entry((node, dof)).or_default().push(v); + } + } + } + } + let mut dependent: BTreeSet<(u32, u8)> = BTreeSet::new(); + for c in model.couplings().values().filter(|c| c.kind == CouplingKind::Kinematic) { + for &node in r.sets.get(&c.set).map(|s| s.nodes.as_slice()).unwrap_or_default() { + dependent.extend(c.dofs.list().into_iter().map(|d| (node, d))); + } + } + for t in model.ties().values() { + for &node in r.sets.get(&t.secondary).map(|s| s.nodes.as_slice()).unwrap_or_default() { + dependent.extend((1..=6).map(|d| (node, d))); + } + } + let state = |node: u32, dof: u8| { + if dependent.contains(&(node, dof)) { + return crate::Held::Sometimes; + } + match held.get(&(node, dof)) { + None => crate::Held::Free, + Some(v) if v.len() == steps.len() && v.iter().all(|x| *x == v[0]) => crate::Held::Always(v[0]), + Some(_) => crate::Held::Sometimes, + } + }; + let mut problems = BTreeSet::new(); + let equations = crate::symmetry_equations(model, r, inclined, &model.rotating_nodes(), state, &mut problems); + unsupported(FORMAT, problems)?; + self.comment("Symmetry planes not normal to X, Y or Z: per node, no motion along the normal (and no rotation about the plane's axes)"); + for e in &equations { + self.line("*EQUATION"); + self.line(&e.terms.len().to_string()); + let mut fields = Vec::new(); + for (dof, c) in &e.terms { + fields.push(format!("{}, {dof}", e.node)); + fields.push(number(*c)?); + } + self.line(&fields.join(", ")); + } + Ok(()) + } + /// Boundary conditions of a step against the state it starts from. fn step_boundaries(&mut self, model: &Model, r: &Resolved, base: &BTreeSet, active: &BTreeSet) -> Result<()> { + let (base, active) = (&held_conditions(model, base), &held_conditions(model, active)); if active == base { self.comment("boundary conditions: as in the base state"); Ok(()) @@ -1188,7 +1284,8 @@ impl Deck { BTreeSet::new() } else { let in_every_step = |n: &String| steps.iter().all(|s| s.boundary_conditions.contains(n)); - model.boundary_conditions().iter().filter(|(n, bc)| holds_at_zero(bc) && in_every_step(n)).map(|(n, _)| n.clone()).collect() + let all = model.boundary_conditions().iter().filter(|(n, bc)| holds_at_zero(bc) && in_every_step(n)).map(|(n, _)| n.clone()).collect(); + held_conditions(model, &all) }; if !initial.is_empty() { self.comment("Initial state: boundary conditions active in every step"); diff --git a/crates/decks/src/lib.rs b/crates/decks/src/lib.rs index f157f7b..faad399 100644 --- a/crates/decks/src/lib.rs +++ b/crates/decks/src/lib.rs @@ -16,7 +16,7 @@ pub mod nastran; use std::collections::{BTreeMap, BTreeSet}; use fieldcraft_model::nodal::edge_pressure; -use fieldcraft_model::{Diagnostic, Load, Model, Resolved, ResolvedDirection}; +use fieldcraft_model::{BoundaryCondition, Diagnostic, Load, Model, Resolved, ResolvedDirection}; #[derive(Debug, Clone, PartialEq, thiserror::Error)] pub enum Error { @@ -176,6 +176,120 @@ pub fn edge_pressure_loads<'a>(model: &Model, r: &Resolved, active: impl IntoIte Ok(out) } +/// How a degree of freedom is constrained apart from the inclined symmetry planes, for choosing +/// the dependent one of a plane's equation (formats eliminate it, so it may carry no other +/// constraint). +#[derive(Debug, Clone, Copy, PartialEq)] +pub(crate) enum Held { + Free, + /// Held at this value wherever the equation applies: it cannot depend on the others, and an + /// equation of held ones alone is redundant when the values satisfy it. + Always(f64), + /// Held in some steps only, or made dependent by a coupling or tie: it cannot be dependent, + /// nor can the equation be dropped. + Sometimes, +} + +/// A linear equation Σ cⱼ u(node, dofⱼ) = 0 of an inclined symmetry plane, dependent degree of +/// freedom first. +#[derive(Debug, Clone, PartialEq)] +pub(crate) struct Equation { + pub node: u32, + pub terms: Vec<(u8, f64)>, +} + +/// A coefficient below this fraction of its row's largest is round-off of the elimination. +const NEGLIGIBLE: f64 = 1e-12; + +/// The equations of the inclined symmetry conditions `names` (others are skipped), node by node +/// in label order: n · u = 0 and, at a node in `rotating`, t₁ · θ = 0 and t₂ · θ = 0 for the +/// plane's in-plane axes ([`fieldcraft_model::Symmetry`]). The rows of one node (several planes +/// meet at an edge) are reduced against each other, so each has a dependent degree of freedom +/// of its own: the largest coefficient among those `held` says are free. A row of held degrees +/// of freedom only is left out when their values satisfy it; otherwise it is a problem, as is a +/// row with none free that cannot be left out. +pub(crate) fn symmetry_equations<'a>( + model: &Model, + r: &Resolved, + names: impl IntoIterator, + rotating: &BTreeSet, + held: impl Fn(u32, u8) -> Held, + problems: &mut BTreeSet, +) -> Vec { + let mut rows: BTreeMap> = BTreeMap::new(); + let mut seen = BTreeSet::new(); + for name in names { + if !seen.insert(name.as_str()) { + continue; + } + let Some(BoundaryCondition::Symmetry { target, plane }) = model.boundary_condition(name) else { continue }; + let (Some(n), Some([t1, t2])) = (plane.inclined(), plane.in_plane()) else { continue }; + for &node in r.target_nodes(target) { + let at = rows.entry(node).or_default(); + at.push((name, [n[0], n[1], n[2], 0.0, 0.0, 0.0])); + if rotating.contains(&node) || target.reference_point().is_some() { + at.push((name, [0.0, 0.0, 0.0, t1[0], t1[1], t1[2]])); + at.push((name, [0.0, 0.0, 0.0, t2[0], t2[1], t2[2]])); + } + } + } + let mut out = Vec::new(); + for (node, rows) in rows { + // The equations of this node so far: (dependent index, coefficients). + let mut done: Vec<(usize, [f64; 6])> = Vec::new(); + for (name, mut c) in rows { + let scale = c.iter().fold(0.0f64, |m, v| m.max(v.abs())); + for (d, e) in &done { + let f = c[*d] / e[*d]; + for k in 0..6 { + c[k] -= f * e[k]; + } + c[*d] = 0.0; + } + for v in &mut c { + if v.abs() <= NEGLIGIBLE * scale { + *v = 0.0; + } + } + if c.iter().all(|v| *v == 0.0) { + continue; + } + let dof = |k: usize| k as u8 + 1; + let pivot = (0..6).filter(|&k| c[k] != 0.0 && held(node, dof(k)) == Held::Free).fold(None, |best: Option, k| match best { + Some(b) if c[b].abs() >= c[k].abs() => Some(b), + _ => Some(k), + }); + let Some(pivot) = pivot else { + let mut residual = 0.0; + let mut size = 0.0; + let mut always = true; + for k in (0..6).filter(|&k| c[k] != 0.0) { + match held(node, dof(k)) { + Held::Always(v) => { + residual += c[k] * v; + size += (c[k] * v).abs(); + } + _ => always = false, + } + } + if !always { + problems.insert(format!( + "boundary condition `{name}`: at node {node} every degree of freedom the inclined symmetry plane relates is held in some steps only or follows a coupling or tie, so none can be eliminated by its equation" + )); + } else if residual.abs() > 1e-9 * size { + problems.insert(format!("boundary condition `{name}`: node {node} is held off the inclined symmetry plane")); + } + continue; + }; + let mut terms = vec![(dof(pivot), c[pivot])]; + terms.extend((0..6).filter(|&k| k != pivot && c[k] != 0.0).map(|k| (dof(k), c[k]))); + out.push(Equation { node, terms }); + done.push((pivot, c)); + } + } + out +} + /// The global degree of freedom (1–6) a spring or dashpot direction acts along, when it is one /// (a vector or line along an axis counts: the sign does not matter to a spring). pub(crate) fn axis_dof(d: ResolvedDirection) -> Option { diff --git a/crates/decks/src/lsdyna/mod.rs b/crates/decks/src/lsdyna/mod.rs index 929da99..6c8dc72 100644 --- a/crates/decks/src/lsdyna/mod.rs +++ b/crates/decks/src/lsdyna/mod.rs @@ -40,7 +40,8 @@ //! LS-DYNA's defaults; a step without nonlinear geometry and hyperelastic materials in nonlinear //! steps are refused); boundary conditions //! as `*BOUNDARY_SPC_*_ID` (named) and `*BOUNDARY_PRESCRIBED_MOTION_*_ID` for non-zero -//! displacements; forces, moments and tractions as `*LOAD_NODE_SET`/`_POINT` from the +//! displacements, a symmetry plane not normal to X, Y or Z as `*BOUNDARY_SPC_*_ID` in a +//! `*DEFINE_COORDINATE_SYSTEM` whose x axis is the plane's normal; forces, moments and tractions as `*LOAD_NODE_SET`/`_POINT` from the //! resolved nodal values; pressure as `*LOAD_SEGMENT_SET`; gravity as `*LOAD_BODY_X/Y/Z` //! (with `*LOAD_BODY_PARTS` when it acts on a set); output as `*DATABASE_BINARY_D3PLOT`, //! `*DATABASE_SPCFORC` and `*DATABASE_EXTENT_BINARY`. @@ -1048,6 +1049,24 @@ fn boundary_condition( (r.reference_point_nodes.get(p).copied().ok_or_else(|| broken(format!("boundary condition `{name}`: reference point `{p}`")))?, "NODE") } }; + // A symmetry plane not normal to an axis: held in a coordinate system whose x axis is its + // normal and whose xy plane holds its first in-plane axis (translation x, rotations y and z). + if let BoundaryCondition::Symmetry { plane, .. } = bc + && let (Some(n), Some([t1, _])) = (plane.inclined(), plane.in_plane()) + { + *constraint_id = constraint_id.saturating_add(1); + let cid = *constraint_id; + d.line("*DEFINE_COORDINATE_SYSTEM"); + d.header(&["cid", "xo", "yo", "zo", "xl", "yl", "zl", "cidl"]); + d.card([Field::from(cid), Field::from(0.0), Field::from(0.0), Field::from(0.0), Field::from(n[0]), Field::from(n[1]), Field::from(n[2])]); + d.header(&["xp", "yp", "zp"]); + d.card([Field::from(t1[0]), Field::from(t1[1]), Field::from(t1[2])]); + d.line(&format!("*BOUNDARY_SPC_{suffix}_ID")); + id_card(d, constraint_id, name); + d.header(&["nsid", "cid", "dofx", "dofy", "dofz", "dofrx", "dofry", "dofrz"]); + d.card([Field::from(entity), Field::from(cid)].into_iter().chain([1, 0, 0, 0, 1, 1].map(|v: i32| Field::from(v)))); + return Ok(()); + } let dofs = bc.dofs(); // Held degrees of freedom (in a frequency step every constrained one is held). let held: Vec = dofs.list().into_iter().filter(|dof| !is_static || bc.value(*dof).is_none_or(|v| v == 0.0)).collect(); diff --git a/crates/decks/src/nastran/mod.rs b/crates/decks/src/nastran/mod.rs index f807129..c910167 100644 --- a/crates/decks/src/nastran/mod.rs +++ b/crates/decks/src/nastran/mod.rs @@ -22,6 +22,9 @@ //! the basic X axis through a `CORD2R`; //! - boundary conditions: per step, one constraint set of `SPC1` (held degrees of freedom) and //! `SPC` (prescribed values); in a modes step every constrained degree of freedom is held; +//! symmetry about a plane not normal to X, Y or Z: one `MPC` set per step, an equation per node +//! and held direction (n · u = 0, and t · θ = 0 about the plane's two axes at shell and beam +//! grids), its dependent component one no `SPC` or rigid element constrains; //! - forces, moments and tractions: `FORCE` and `MOMENT` at the nodal values the model resolves, //! so every writer and the solver share one distribution; pressure: `PLOAD4` per element face; //! gravity: `GRAV`; a step with several loads combines them with `LOAD`; @@ -241,6 +244,8 @@ type Constraints = BTreeMap; 6]>; struct Subcase<'m> { step: &'m Step, spc: Option, + /// The equations of its inclined symmetry planes. + mpc: Option, load: Option, method: Option, } @@ -269,6 +274,8 @@ struct Deck<'m> { dependent: BTreeMap, /// Constraint sets by id. spc_sets: Vec<(u32, Constraints)>, + /// Multipoint constraint sets (inclined symmetry planes) by id. + mpc_sets: Vec<(u32, Vec)>, /// Load name → load set id. load_sids: BTreeMap<&'m str, u32>, /// Combinations: (id, load set ids). @@ -310,7 +317,7 @@ impl<'m> Deck<'m> { model, r, sol, - subcases: steps.into_iter().map(|step| Subcase { step, spc: None, load: None, method: None }).collect(), + subcases: steps.into_iter().map(|step| Subcase { step, spc: None, mpc: None, load: None, method: None }).collect(), problems, next_sid: 1, materials: BTreeMap::new(), @@ -321,6 +328,7 @@ impl<'m> Deck<'m> { discrete: Vec::new(), dependent: BTreeMap::new(), spc_sets: Vec::new(), + mpc_sets: Vec::new(), load_sids: BTreeMap::new(), combinations: Vec::new(), } @@ -587,6 +595,7 @@ impl<'m> Deck<'m> { /// One constraint set per distinct step constraint, checked against coupled nodes. fn plan_constraints(&mut self) -> Result<()> { let (model, r) = (self.model, self.r); + let rotating = model.rotating_nodes(); // Conflicts, one problem per pair of definitions: the nodes they meet on. let mut clashes: BTreeMap<(&str, &str, &str), BTreeSet> = BTreeMap::new(); let mut on_dependent: BTreeMap<(&str, &str), BTreeSet> = BTreeMap::new(); @@ -598,6 +607,10 @@ impl<'m> Deck<'m> { for name in &step.boundary_conditions { let Some(bc) = model.boundary_condition(name) else { continue }; let dofs = bc.dofs(); + // Inclined symmetry planes hold no global degree of freedom: they are equations. + if dofs.is_empty() { + continue; + } // In a modes step every constrained degree of freedom is held. let values: Vec<(u8, f64)> = dofs.list().into_iter().filter_map(|d| Some((d, if modes { 0.0 } else { bc.value(d)? }))).collect(); for &node in r.target_nodes(bc.target()) { @@ -622,6 +635,19 @@ impl<'m> Deck<'m> { } } } + let dependent = &self.dependent; + let held = |node: u32, dof: u8| match ( + dependent.get(&node), + set.get(&node).and_then(|v| v.get(usize::from(dof).checked_sub(1)?).copied().flatten()), + ) { + (Some((d, _)), _) if d.contains(dof) => crate::Held::Sometimes, + (_, Some(v)) => crate::Held::Always(v), + _ => crate::Held::Free, + }; + let equations = crate::symmetry_equations(model, r, &step.boundary_conditions, &rotating, held, &mut self.problems); + if !equations.is_empty() { + sub.mpc = Some(intern(&mut self.mpc_sets, &mut self.next_sid, equations)?); + } if !set.is_empty() { sub.spc = Some(intern(&mut self.spc_sets, &mut self.next_sid, set)?); } @@ -707,7 +733,7 @@ impl<'m> Deck<'m> { let _ = writeln!(out, "SUBCASE {}", i + 1); let title: String = format!(" SUBTITLE = {}", sub.step.name).chars().take(72).collect(); let _ = writeln!(out, "{title}"); - for (key, sid) in [("SPC", sub.spc), ("LOAD", sub.load), ("METHOD", sub.method)] { + for (key, sid) in [("SPC", sub.spc), ("MPC", sub.mpc), ("LOAD", sub.load), ("METHOD", sub.method)] { if let Some(sid) = sid { let _ = writeln!(out, " {key} = {sid}"); } @@ -715,6 +741,7 @@ impl<'m> Deck<'m> { for v in &sub.step.output { let request = match v { OutputVariable::Displacement => "DISPLACEMENT", + OutputVariable::ReactionForce if sub.mpc.is_some() => "SPCFORCES = ALL\n MPCFORCES", OutputVariable::ReactionForce => "SPCFORCES", OutputVariable::Stress => "STRESS", OutputVariable::Strain => "STRAIN", @@ -892,6 +919,23 @@ impl<'m> Deck<'m> { Card::new("SPC").int(*sid).int(node).int(component(Dofs::from_bits(mask))).real(value).write(out)?; } } + for (sid, equations) in &self.mpc_sets { + for sub in self.subcases.iter().filter(|s| s.mpc == Some(*sid)) { + comment(out, &format!("step {}: inclined symmetry planes, MPC {sid} (the first grid component is dependent)", sub.step.name)); + } + for e in equations { + // SID, then (G, C, A) triples: two on the first line (its last field blank), two + // per continuation after a blank first field. + let mut card = Card::new("MPC").int(*sid); + for (i, (dof, c)) in e.terms.iter().enumerate() { + if i % 2 == 0 && i > 0 { + card = card.blank().blank(); + } + card = card.int(e.node).int(i64::from(*dof)).real(*c); + } + card.write(out)?; + } + } Ok(()) } diff --git a/crates/decks/tests/inclined_symmetry.rs b/crates/decks/tests/inclined_symmetry.rs new file mode 100644 index 0000000..20b2266 --- /dev/null +++ b/crates/decks/tests/inclined_symmetry.rs @@ -0,0 +1,193 @@ +//! Symmetry about a plane not normal to X, Y or Z in the three decks (#242): Abaqus `*EQUATION`s, +//! Nastran `MPC` sets and an LS-DYNA `*BOUNDARY_SPC` in a local coordinate system, on a small +//! shell plate of our own whose diagonal x = y is the plane. + +use fieldcraft_decks::{Error, abaqus, lsdyna, nastran}; +use fieldcraft_model::{ + BoundaryCondition, Elastic, ElementFace, ElementKind, Formulation, Load, Material, Model, Procedure, Section, SectionKind, Set, Step, Symmetry, + Target, +}; + +/// A 2 × 2 quad4 shell plate in the XY plane (nodes 1–9, x and y in 0, 1, 2), clamped on x = 0, +/// the inclined plane (normal (1, −1, 0)) on the diagonal nodes 1, 5 and 9, under pressure; +/// steps `Load` (static) and `Modes`, `in_modes` saying whether the plane is active in both. +fn plate(in_modes: bool) -> Model { + let mut m = Model::default(); + for j in 0..3u32 { + for i in 0..3u32 { + m.set_node(1 + i + 3 * j, [f64::from(i), f64::from(j), 0.0]).unwrap(); + } + } + for (e, n) in [(1, [1, 2, 5, 4]), (2, [2, 3, 6, 5]), (3, [4, 5, 8, 7]), (4, [5, 6, 9, 8])] { + m.set_element(e, ElementKind::Quad4, n.to_vec()).unwrap(); + } + m.add_material( + "Steel", + Material { elastic: Some(Elastic::Isotropic { young: 210e3, poisson: 0.3 }), density: Some(7.85e-9), ..Material::default() }, + ) + .unwrap(); + m.add_set("all", Set { elements: vec![1, 2, 3, 4], ..Set::default() }).unwrap(); + let kind = SectionKind::Shell { thickness: 0.1, offset: 0.0, points: 5 }; + m.add_section("plate", Section { set: "all".into(), material: "Steel".into(), kind, formulation: Formulation::default() }).unwrap(); + m.add_set("edge", Set { nodes: vec![1, 4, 7], ..Set::default() }).unwrap(); + m.add_set("diagonal", Set { nodes: vec![1, 5, 9], ..Set::default() }).unwrap(); + let faces = (1..=4).map(|element| ElementFace { element, face: 1 }).collect(); + m.add_set("top", Set { faces, ..Set::default() }).unwrap(); + m.add_boundary_condition("clamp", BoundaryCondition::Fixed { target: Target::Set("edge".into()) }).unwrap(); + let plane = Symmetry::from_normal([1.0, -1.0, 0.0]).unwrap(); + m.add_boundary_condition("sym", BoundaryCondition::Symmetry { target: Target::Set("diagonal".into()), plane }).unwrap(); + m.add_load("q", Load::Pressure { set: "top".into(), magnitude: 0.01 }).unwrap(); + let mut load = Step::new("Load", Procedure::Static {}); + load.boundary_conditions = ["clamp", "sym"].map(String::from).into(); + load.loads = ["q".to_string()].into(); + m.add_step(load, None).unwrap(); + let mut modes = Step::new("Modes", Procedure::Frequency { modes: 3, min_hz: None, max_hz: None }); + modes.boundary_conditions = if in_modes { ["clamp", "sym"].map(String::from).into() } else { ["clamp".to_string()].into() }; + m.add_step(modes, None).unwrap(); + m +} + +/// A Nastran real: an exponent may come without `E`. +fn real(s: &str) -> f64 { + let split = s.char_indices().skip(1).find(|(i, c)| (*c == '+' || *c == '-') && !s[..*i].ends_with(['E', 'e'])); + match split { + Some((i, _)) => format!("{}E{}", &s[..i], &s[i..]).parse().unwrap(), + None => s.parse().unwrap(), + } +} + +fn dot(a: [f64; 3], b: [f64; 3]) -> f64 { + (0..3).map(|k| a[k] * b[k]).sum() +} + +/// Every equation is one of the plane's at a free node: its terms on one node, its coefficients +/// along the normal (translations) or along an in-plane axis (rotations), unit length; nodes 5 +/// and 9 get three each (n · u, t₁ · θ, t₂ · θ); node 1 is clamped, so its equations are left out; +/// no dependent degree of freedom repeats. +fn check(equations: &[(u32, Vec<(u8, f64)>)]) { + let n = [1.0 / 2f64.sqrt(), -1.0 / 2f64.sqrt(), 0.0]; + let mut dependents = std::collections::BTreeSet::new(); + let mut per_node = std::collections::BTreeMap::::new(); + for (node, terms) in equations { + assert!(dependents.insert((*node, terms[0].0)), "{node} {terms:?}"); + *per_node.entry(*node).or_default() += 1; + let mut v = [0.0; 3]; + let rotation = terms[0].0 > 3; + for (dof, c) in terms { + assert_eq!(*dof > 3, rotation, "{terms:?}"); + v[usize::from((dof - 1) % 3)] = *c; + } + assert!((dot(v, v) - 1.0).abs() < 1e-12, "{terms:?}"); + if rotation { + assert!(dot(v, n).abs() < 1e-12, "{terms:?}"); + } else { + assert!((dot(v, n).abs() - 1.0).abs() < 1e-12, "{terms:?}"); + } + } + assert_eq!(per_node, [(5, 3), (9, 3)].into()); +} + +#[test] +fn abaqus_writes_equations_in_the_model_data() { + let deck = abaqus::write(&plate(true)).unwrap(); + let lines: Vec<&str> = deck.lines().collect(); + let mut equations = Vec::new(); + for (i, l) in lines.iter().enumerate() { + if *l == "*EQUATION" { + let count: usize = lines[i + 1].trim().parse().unwrap(); + let f: Vec<&str> = lines[i + 2].split(',').map(str::trim).collect(); + assert_eq!(f.len(), 3 * count, "{}", lines[i + 2]); + let node: u32 = f[0].parse().unwrap(); + let terms = f.chunks(3).map(|t| (t[1].parse().unwrap(), t[2].parse().unwrap())).collect(); + assert!(f.chunks(3).all(|t| t[0] == f[0])); + equations.push((node, terms)); + } + } + check(&equations); + // Model data, before the first step; the plane is no `*BOUNDARY` line. + let first_step = deck.find("*STEP").unwrap(); + assert!(deck.rfind("*EQUATION").unwrap() < first_step); + assert!(!deck.contains("SYMM"), "{deck}"); + + let e = abaqus::write(&plate(false)).unwrap_err(); + assert!(matches!(&e, Error::Unsupported { what, .. } if what.contains("boundary condition `sym`") && what.contains("some steps only")), "{e}"); +} + +#[test] +fn nastran_writes_one_mpc_set_per_step() { + // A static deck and a modes deck (one solution per run). + for (steps, sol) in [("Load", 101), ("Modes", 103)] { + let m = fieldcraft_decks::with_steps(&plate(true), &[steps.to_string()]).unwrap(); + let deck = nastran::write(&m).unwrap(); + assert!(deck.contains(&format!("SOL {sol}"))); + let sid = deck.lines().find_map(|l| l.trim().strip_prefix("MPC = ")).unwrap().to_string(); + // MPC entries: SID, then (grid, component, coefficient) triples; the continuation after + // two triples starts with a blank field. Large-field lines are 16 columns per field. + let mut equations = Vec::new(); + let mut fields: Vec = Vec::new(); + let flush = |fields: &mut Vec, equations: &mut Vec<(u32, Vec<(u8, f64)>)>| { + if fields.is_empty() { + return; + } + assert_eq!(fields[0], sid); + let data: Vec<&String> = fields[1..].iter().enumerate().filter(|(i, _)| !matches!(i % 8, 6 | 7)).map(|(_, f)| f).collect(); + let triples: Vec<(u32, u8, f64)> = data + .chunks(3) + .filter(|t| t.len() == 3 && !t[0].is_empty()) + .map(|t| (t[0].parse().unwrap(), t[1].parse().unwrap(), real(t[2]))) + .collect(); + assert!(triples.iter().all(|t| t.0 == triples[0].0)); + equations.push((triples[0].0, triples.iter().map(|t| (t.1, t.2)).collect())); + fields.clear(); + }; + for line in deck.lines() { + let (head, rest) = line.split_at(8.min(line.len())); + let row: Vec = rest.as_bytes().chunks(16).map(|c| String::from_utf8_lossy(c).trim().to_string()).collect(); + let mut row = row; + row.resize(4, String::new()); + match head.trim() { + "MPC*" => { + flush(&mut fields, &mut equations); + fields.extend(row); + } + "*" if !fields.is_empty() => fields.extend(row), + _ => flush(&mut fields, &mut equations), + } + } + flush(&mut fields, &mut equations); + check(&equations); + } +} + +#[test] +fn lsdyna_holds_the_plane_in_a_local_system() { + let deck = lsdyna::write(&plate(true)).unwrap().text; + let block = |keyword: &str| -> Vec> { + let at = deck.find(keyword).unwrap(); + deck[at..] + .lines() + .skip(1) + .filter(|l| !l.starts_with('$')) + .take_while(|l| !l.starts_with('*')) + .map(|l| l.as_bytes().chunks(20).filter_map(|c| String::from_utf8_lossy(c).trim().parse().ok()).collect()) + .collect() + }; + let system = block("*DEFINE_COORDINATE_SYSTEM"); + let (cid, x, p) = (system[0][0], [system[0][4], system[0][5], system[0][6]], [system[1][0], system[1][1], system[1][2]]); + assert!((x[0] + x[1]).abs() < 1e-12 && x[2] == 0.0 && (dot(x, x) - 1.0).abs() < 1e-12, "{x:?}"); + assert!(dot(x, p).abs() < 1e-12, "{p:?}"); + // The SPC after the system: its id card, then nsid, cid and the held local directions. + let at = deck.find("*DEFINE_COORDINATE_SYSTEM").unwrap(); + let spc = deck[at..].find("*BOUNDARY_SPC_SET_ID").unwrap() + at; + let card: Vec = deck[spc..] + .lines() + .skip(2) + .find(|l| !l.starts_with('$')) + .unwrap() + .as_bytes() + .chunks(20) + .map(|c| String::from_utf8_lossy(c).trim().parse().unwrap()) + .collect(); + assert_eq!(card[1], cid); + assert_eq!(&card[2..], [1.0, 0.0, 0.0, 0.0, 1.0, 1.0]); +} diff --git a/crates/model/src/bcs.rs b/crates/model/src/bcs.rs index 2d9f3a5..8bca0db 100644 --- a/crates/model/src/bcs.rs +++ b/crates/model/src/bcs.rs @@ -110,25 +110,101 @@ impl std::fmt::Display for Target { } } -/// The plane of a symmetry boundary condition, named by its normal axis. -#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] +/// The plane of a symmetry boundary condition, named by its normal: an axis (`"x"`), or any +/// direction (`{"normal": [1, 1, 0]}`, not necessarily of unit length). +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] #[serde(rename_all = "snake_case")] pub enum Symmetry { /// The plane normal to X: no X translation, no rotation about Y or Z. X, Y, Z, + /// A plane normal to a direction: no translation along it, no rotation about the two axes + /// in the plane (n · u = 0, and θ ∥ n at nodes with rotations). + Normal([f64; 3]), } +/// A normal whose component along an axis is within this of 1 (unit length) is that axis. +const ALONG_AXIS: f64 = 1e-12; + impl Symmetry { - /// The fixed degrees of freedom. - pub fn dofs(self) -> Dofs { + /// The plane normal to `normal`: [`Symmetry::X`], `Y` or `Z` when it lies along an axis, + /// otherwise [`Symmetry::Normal`] with the unit normal; `None` for a zero or non-finite one. + pub fn from_normal(normal: [f64; 3]) -> Option { + let unit = unit(normal)?; + Some(match (0..3).find(|&k| (unit[k].abs() - 1.0).abs() < ALONG_AXIS) { + Some(0) => Symmetry::X, + Some(1) => Symmetry::Y, + Some(2) => Symmetry::Z, + _ => Symmetry::Normal(unit), + }) + } + + /// The unit normal (`None` for a zero or non-finite [`Symmetry::Normal`]). + pub fn unit_normal(self) -> Option<[f64; 3]> { match self { - Symmetry::X => Dofs(0b110_001), - Symmetry::Y => Dofs(0b101_010), - Symmetry::Z => Dofs(0b011_100), + Symmetry::X => Some([1.0, 0.0, 0.0]), + Symmetry::Y => Some([0.0, 1.0, 0.0]), + Symmetry::Z => Some([0.0, 0.0, 1.0]), + Symmetry::Normal(n) => unit(n), } } + + /// The axis (0–2) the normal lies along, so the condition holds global degrees of freedom; + /// `None` for an inclined plane, which is a linear constraint instead. + pub fn axis(self) -> Option { + let n = self.unit_normal()?; + (0..3).find(|&k| (n[k].abs() - 1.0).abs() < ALONG_AXIS) + } + + /// The held global degrees of freedom: those of an axis-aligned plane; none for an inclined + /// one, which holds combinations of them ([`Symmetry::unit_normal`]). + pub fn dofs(self) -> Dofs { + match self.axis() { + Some(0) => Dofs(0b110_001), + Some(1) => Dofs(0b101_010), + Some(2) => Dofs(0b011_100), + _ => Dofs::none(), + } + } + + /// The inclined plane's unit normal; `None` for one normal to an axis. + pub fn inclined(self) -> Option<[f64; 3]> { + if self.axis().is_some() { None } else { self.unit_normal() } + } + + /// Two unit vectors in the plane, perpendicular to each other and to the unit normal n: + /// t₁ = n × e / |n × e| with e the axis n leans on least, and t₂ = n × t₁ (`None` for a + /// zero or non-finite normal). + pub fn in_plane(self) -> Option<[[f64; 3]; 2]> { + let n = self.unit_normal()?; + let least = (0..3).fold(0, |m, k| if n[k].abs() < n[m].abs() { k } else { m }); + let mut e = [0.0; 3]; + e[least] = 1.0; + let cross = |a: [f64; 3], b: [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]]; + let t1 = unit(cross(n, e))?; + Some([t1, cross(n, t1)]) + } + + fn validate(self) -> Result<()> { + if let Symmetry::Normal(n) = self + && unit(n).is_none() + { + return Err(Error::Invalid(format!("the symmetry plane's normal must be a finite, non-zero vector, got {n:?}"))); + } + Ok(()) + } +} + +/// `v` scaled to unit length; `None` for a zero or non-finite vector. +fn unit(v: [f64; 3]) -> Option<[f64; 3]> { + let scale = v.iter().fold(0.0f64, |m, c| m.max(c.abs())); + if !(scale > 0.0 && v.iter().all(|c| c.is_finite())) { + return None; + } + let w = v.map(|c| c / scale); + let length = w.iter().map(|c| c * c).sum::().sqrt(); + Some(w.map(|c| c / length)) } /// A boundary condition. @@ -198,6 +274,9 @@ impl BoundaryCondition { } pub fn validate(&self) -> Result<()> { + if let BoundaryCondition::Symmetry { plane, .. } = self { + plane.validate()?; + } if let BoundaryCondition::Temperature { target, value } = self { finite("temperature", *value)?; if let Target::ReferencePoint(r) = target { @@ -285,4 +364,35 @@ mod tests { let parsed: BoundaryCondition = serde_json::from_value(json.clone()).unwrap(); assert_eq!(serde_json::to_value(&parsed).unwrap(), json); } + + #[test] + fn inclined_symmetry_planes() { + let tilt = Symmetry::from_normal([1.0, 1.0, 0.0]).unwrap(); + let h = 1.0 / 2f64.sqrt(); + assert_eq!(tilt, Symmetry::Normal([h, h, 0.0])); + assert_eq!((tilt.axis(), tilt.dofs(), tilt.inclined()), (None, Dofs::none(), Some([h, h, 0.0]))); + assert_eq!(Symmetry::from_normal([0.0, -2.0, 0.0]), Some(Symmetry::Y)); + assert_eq!(Symmetry::Normal([0.0, 0.0, 3.0]).dofs(), Symmetry::Z.dofs()); + assert_eq!(Symmetry::Normal([0.0, 0.0, 3.0]).inclined(), None); + for bad in [[0.0; 3], [f64::NAN, 1.0, 0.0], [f64::INFINITY, 0.0, 0.0]] { + assert_eq!(Symmetry::from_normal(bad), None); + let bc = BoundaryCondition::Symmetry { target: Target::Set("s".into()), plane: Symmetry::Normal(bad) }; + assert!(bc.validate().is_err(), "{bad:?}"); + } + // The tiniest and the largest finite normals still give a unit one. + assert_eq!(Symmetry::from_normal([1e-300, 1e-300, 0.0]), Some(Symmetry::Normal([h, h, 0.0]))); + assert_eq!(Symmetry::from_normal([f64::MAX, f64::MAX, 0.0]), Some(Symmetry::Normal([h, h, 0.0]))); + for plane in [tilt, Symmetry::X, Symmetry::Normal([0.3, -2.0, 7.0])] { + let n = plane.unit_normal().unwrap(); + let [a, b] = plane.in_plane().unwrap(); + let dot = |u: [f64; 3], v: [f64; 3]| (0..3).map(|k| u[k] * v[k]).sum::(); + for (u, v, want) in [(a, a, 1.0), (b, b, 1.0), (a, b, 0.0), (a, n, 0.0), (b, n, 0.0)] { + assert!((dot(u, v) - want).abs() < 1e-15, "{plane:?}"); + } + } + let json = serde_json::json!({"type": "symmetry", "target": {"set": "s"}, "plane": {"normal": [1.0, 1.0, 0.0]}}); + let parsed: BoundaryCondition = serde_json::from_value(json.clone()).unwrap(); + assert_eq!(parsed, BoundaryCondition::Symmetry { target: Target::Set("s".into()), plane: Symmetry::Normal([1.0, 1.0, 0.0]) }); + assert_eq!(serde_json::to_value(&parsed).unwrap(), json); + } } diff --git a/crates/model/src/check.rs b/crates/model/src/check.rs index 0a6d284..6728def 100644 --- a/crates/model/src/check.rs +++ b/crates/model/src/check.rs @@ -176,6 +176,19 @@ impl Model { .collect() } + /// The mesh nodes with rotational degrees of freedom under the model's own sections: those of + /// shells and beams, not of plane stress, plane strain or axisymmetric elements. + pub fn rotating_nodes(&self) -> BTreeSet { + let planar: BTreeSet = self + .sections + .values() + .filter(|s| s.kind.extent().is_some()) + .filter_map(|s| self.resolve_set_partial(&s.set)) + .flat_map(|r| r.elements) + .collect(); + self.rotational_nodes(|l| planar.contains(&l)) + } + fn analyze(&self) -> (Resolved, Vec) { let mut d = Diagnostics::default(); let mut resolved = Resolved::default(); diff --git a/crates/render/src/glyphs.rs b/crates/render/src/glyphs.rs index 8e88c53..815bbd0 100644 --- a/crates/render/src/glyphs.rs +++ b/crates/render/src/glyphs.rs @@ -354,19 +354,6 @@ fn scan(items: &[T]) -> impl Iterator { items.iter().step_by(items.len().div_ceil(MAX_SCAN).max(1)) } -/// The nodes that can rotate: shell and beam nodes (not those of 2D sections: plane stress, -/// plane strain, axisymmetric), as the model check counts them. -fn rotating_nodes(model: &Model) -> BTreeSet { - let planar: BTreeSet = model - .sections() - .values() - .filter(|s| s.kind.extent().is_some()) - .filter_map(|s| model.resolve_set_partial(&s.set)) - .flat_map(|r| r.elements) - .collect(); - model.rotational_nodes(|l| planar.contains(&l)) -} - /// The corner node labels of an element face, in the face's order. fn face_corners(model: &Model, f: ElementFace) -> Option> { let kind = model.element(f.element)?.kind; @@ -589,7 +576,7 @@ impl Builder<'_> { } let (spots, total) = self.spots(bc.target()); if dofs & 0b111_000 != 0 && self.rotating.is_none() { - self.rotating = Some(rotating_nodes(self.model)); + self.rotating = Some(self.model.rotating_nodes()); } let none = BTreeSet::new(); let rotating = self.rotating.as_ref().unwrap_or(&none); diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index 0212e3d..4b121f5 100644 --- a/crates/solver/src/analysis.rs +++ b/crates/solver/src/analysis.rs @@ -12,7 +12,7 @@ use serde::Serialize; use fieldcraft_parallel as parallel; use crate::cholesky::{Factor, FactorError, UpperColumns}; -use crate::constraints::{Budget, Reduction, coupling_rows}; +use crate::constraints::{Budget, Reduction, coupling_rows, symmetry_rows}; use crate::eigen::{Operator, Wanted, lanczos}; use crate::ldlt::{Kind, Symmetric}; use crate::ordering::nested_dissection; @@ -342,6 +342,13 @@ pub(crate) struct Setup<'a> { pub(crate) fn reduction(s: &Setup, values: bool, order: &[usize], budget: Budget, warnings: &mut Vec) -> Result { let mut rows = coupling_rows(s.model, s.resolved, s.mesh)?; rows.extend(crate::tie::rows(s.model, s.resolved, s.mesh, warnings)?); + // Reference point directions exist when a coupling, tie or discrete element uses them. + let used: std::collections::BTreeSet = rows.iter().flat_map(|r| r.terms.iter().map(|t| t.0)).chain(s.mesh.discrete.dofs()).collect(); + let node_dofs = s.mesh.layout.node_dofs(); + // Boundary conditions first, as the elimination takes them. + let mut symmetry = symmetry_rows(s.model, s.resolved, s.mesh, s.step, |d| d < node_dofs || used.contains(&d))?; + symmetry.append(&mut rows); + let rows = symmetry; Reduction::new( &s.mesh.layout, s.mesh.discrete.dofs(), diff --git a/crates/solver/src/constraints.rs b/crates/solver/src/constraints.rs index d30e861..0986e61 100644 --- a/crates/solver/src/constraints.rs +++ b/crates/solver/src/constraints.rs @@ -53,7 +53,7 @@ use std::collections::BTreeMap; use fieldcraft_fem::dense::inverse3; use fieldcraft_fem::{DofKind, DofLayout, DofSet}; -use fieldcraft_model::{CouplingKind, Model, Resolved}; +use fieldcraft_model::{BoundaryCondition, CouplingKind, Model, Resolved, Step}; use crate::sparse::Csr; use crate::system::Mesh; @@ -74,6 +74,9 @@ pub struct Row<'m> { #[derive(Debug, Clone, Copy, PartialEq, Eq)] pub enum Source { Coupling, + /// A symmetry boundary condition on an inclined plane: n · u = 0, and θ ∥ n at nodes with + /// rotations. + Symmetry, /// A tie: its secondary node follows a weighted average of main nodes. Tie, } @@ -82,6 +85,7 @@ impl Source { fn what(self) -> &'static str { match self { Source::Coupling => "coupling", + Source::Symmetry => "boundary condition", Source::Tie => "tie", } } @@ -224,6 +228,64 @@ pub fn coupling_rows<'m>(model: &'m Model, resolved: &Resolved, mesh: &Mesh) -> Ok(rows) } +/// The constraint rows of a step's symmetry conditions on inclined planes (those on a plane +/// normal to an axis hold global degrees of freedom and are prescribed values instead), in the +/// step's order, nodes by label. At every target node the translation along the unit normal n +/// is zero, n · u = 0, eliminating the translation of the largest |nₖ|; at a node with +/// rotations (shells, beams, a reference point whose rotations exist) the rotation about each of +/// two in-plane axes t₁, t₂ is zero too, tᵢ · θ = 0: a mirror image turns only about the normal. +/// `exists` says whether a reference point direction exists (mesh node ones always do). +pub fn symmetry_rows<'m>(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &'m Step, exists: impl Fn(usize) -> bool) -> Result>> { + let mut rows = Vec::new(); + for name in &step.boundary_conditions { + let Some(BoundaryCondition::Symmetry { target, plane }) = model.boundary_condition(name) else { continue }; + let (Some(n), Some([t1, t2])) = (plane.inclined(), plane.in_plane()) else { continue }; + let rp = target.reference_point(); + if rp.is_some_and(|r| mesh.point(r).is_none()) { + // Reported by `prescribed` as holding nothing, like any condition on such a point. + continue; + } + for &label in resolved.target_nodes(target) { + let slot = match rp { + Some(r) => mesh.point(r).map(|k| mesh.layout.point_slot(k)), + None => mesh.layout.slot(label), + }; + let Some(slot) = slot else { continue }; + let kinds = mesh.layout.kinds(slot); + let dof = |kind: Option| kind.and_then(|k| kinds.contains(k).then(|| mesh.layout.dof(slot, k)).flatten()).filter(|d| exists(*d)); + let row = |axis: [f64; 3], rotation: bool| -> Result>> { + let mut terms = Vec::with_capacity(3); + let mut missing = false; + for (k, c) in axis.into_iter().enumerate().filter(|(_, c)| *c != 0.0) { + match dof(if rotation { DofKind::rotation(k) } else { DofKind::translation(k) }) { + Some(d) => terms.push((d, c)), + None => missing = true, + } + } + if missing { + // No rotation exists at a solid node: only its translation is held. + if rotation && terms.is_empty() { + return Ok(None); + } + return Err(Error::Invalid(format!( + "boundary condition `{name}`: node {label} cannot move along the symmetry plane's normal {n:?} (it has no {} along some axis the normal leans on); in a 2D model the plane must contain the Z axis", + if rotation { "rotation" } else { "translation" } + ))); + } + let pivot = terms.iter().fold(None, |best: Option<(usize, f64)>, t| match best { + Some(b) if b.1.abs() >= t.1.abs() => Some(b), + _ => Some(*t), + }); + Ok(pivot.map(|(pivot, _)| Row { terms, pivot, coupling: name, source: Source::Symmetry })) + }; + rows.extend(row(n, false)?); + rows.extend(row(t1, true)?); + rows.extend(row(t2, true)?); + } + } + Ok(rows) +} + /// The degrees of freedom that exist: every mesh node's, and the reference point directions /// some constraint row or `extra` (the discrete elements') involves. fn active_dofs(layout: &DofLayout, extra: impl Iterator, rows: &[Row]) -> Vec { @@ -346,6 +408,10 @@ pub struct Reduction { prescribed_motion: Vec, /// Whether any degree of freedom depends on others (T is not a plain selection). has_dependents: bool, + /// The degrees of freedom an inclined symmetry plane supports (those of its rows that are not + /// prescribed): their share of K u − f is the plane's reaction on them, n λ, and is reported + /// where it acts rather than carried to prescribed degrees of freedom. + supports: Vec, budget: Budget, } @@ -447,7 +513,7 @@ impl Reduction { for row in rows { let scale = row.terms.iter().fold(0.0f64, |m, t| m.max(t.1.abs())); if !(scale > 0.0 && scale.is_finite()) { - return Err(Error::Invalid(format!("coupling `{}`: a constraint with no finite coefficient", row.coupling))); + return Err(Error::Invalid(format!("{} `{}`: a constraint with no finite coefficient", row.source.what(), row.coupling))); } // The row in current free and prescribed degrees of freedom. let mut expanded = Vec::with_capacity(row.terms.len()); @@ -497,7 +563,8 @@ impl Reduction { let size: f64 = terms.iter().map(|t| value(t).abs()).sum(); if residual.abs() > 1e-9 * size.max(f64::MIN_POSITIVE) && residual.abs() > 0.0 { return Err(Error::Invalid(format!( - "coupling `{}` contradicts the prescribed values (boundary conditions or other couplings hold its nodes in a way it cannot follow)", + "{} `{}` contradicts the prescribed values (boundary conditions or couplings hold its nodes in a way it cannot follow)", + row.source.what(), row.coupling ))); } @@ -572,7 +639,17 @@ impl Reduction { } let prescribed_motion = (0..n).map(|d| fixed.row(d).map(|(p, c)| c * values.get(p).copied().unwrap_or(0.0)).sum()).collect(); let columns = free.transpose(eq_dof.len()); - Ok(Reduction { eq_dof, active, free, columns, fixed, prescribed_motion, has_dependents: !exprs.is_empty(), budget }) + let mut supports = vec![false; n]; + for row in rows.iter().filter(|r| r.source == Source::Symmetry) { + for &(dof, _) in &row.terms { + if state.get(dof) != Some(&State::Fixed) + && let Some(x) = supports.get_mut(dof) + { + *x = true; + } + } + } + Ok(Reduction { eq_dof, active, free, columns, fixed, prescribed_motion, has_dependents: !exprs.is_empty(), supports, budget }) } /// Free equations. @@ -737,10 +814,17 @@ impl Reduction { } /// T_pᵀ g for a vector g over every degree of freedom (the reactions of the out-of-balance - /// forces g = f_int − f_ext), zero where nothing is prescribed. + /// forces g = f_int − f_ext), zero where nothing is prescribed; at the degrees of freedom an + /// inclined symmetry plane supports, g itself. pub fn reactions_of(&self, g: &[f64]) -> Vec { let mut out = vec![0.0; self.prescribed_motion.len()]; for (dof, v) in g.iter().enumerate() { + if self.supports.get(dof).copied().unwrap_or(false) { + if let Some(x) = out.get_mut(dof) { + *x += v; + } + continue; + } for (p, c) in self.fixed.row(dof) { if let Some(x) = out.get_mut(p) { *x += c * v; @@ -750,24 +834,23 @@ impl Reduction { out } - /// The reactions at the prescribed degrees of freedom, T_pᵀ (K u − f), over every degree of - /// freedom (zero where nothing is prescribed). + /// The reactions at the prescribed degrees of freedom, T_pᵀ (K u − f), and at those an + /// inclined symmetry plane supports ([`Reduction::reactions_of`]), over every degree of + /// freedom (zero elsewhere). pub fn reactions(&self, k: &Csr, u: &[f64], f: &[f64]) -> Vec { - let mut out = vec![0.0; self.prescribed_motion.len()]; - for (dof, fv) in f.iter().enumerate() { - if self.fixed.is_empty_row(dof) { - continue; - } - let (cols, vals) = k.row(dof); - let ku: f64 = cols.iter().zip(vals).map(|(j, v)| v * u.get(*j).copied().unwrap_or(0.0)).sum(); - let r = ku - fv; - for (p, c) in self.fixed.row(dof) { - if let Some(x) = out.get_mut(p) { - *x += c * r; + let residual: Vec = f + .iter() + .enumerate() + .map(|(dof, fv)| { + if self.fixed.is_empty_row(dof) && !self.supports.get(dof).copied().unwrap_or(false) { + return 0.0; } - } - } - out + let (cols, vals) = k.row(dof); + let ku: f64 = cols.iter().zip(vals).map(|(j, v)| v * u.get(*j).copied().unwrap_or(0.0)).sum(); + ku - fv + }) + .collect(); + self.reactions_of(&residual) } } diff --git a/crates/solver/tests/inclined_symmetry.rs b/crates/solver/tests/inclined_symmetry.rs new file mode 100644 index 0000000..14a2d0b --- /dev/null +++ b/crates/solver/tests/inclined_symmetry.rs @@ -0,0 +1,289 @@ +//! Symmetry about an inclined plane (#242): n · u = 0 on the plane, and zero rotation about the +//! two axes in it at shell nodes, eliminated as linear constraints. +//! +//! The checks need no reference program: a finite element model turned by a rotation R solves +//! to R times the displacements of the unturned one (the elements are frame-indifferent), so an +//! inclined plane is checked against the axis-aligned plane it is the image of; and half of a +//! part whose mesh is mirror-symmetric about the plane, with the condition and half the load, +//! solves exactly as the whole does at the nodes they share. + +mod common; + +use std::collections::BTreeMap; + +use common::{full, grid, iso, near, reduced}; +use fieldcraft_model::{ + BoundaryCondition, Elastic, ElementFace, ElementKind, Formulation, Load, Material, Model, Procedure, Section, SectionKind, Set, Step, Symmetry, + Target, +}; +use fieldcraft_solver::{FrequencyOptions, Limits, StaticResult, solve_frequency, solve_static}; + +/// A rotation by `angle` about the unit 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 mut r = [[0.0; 3]; 3]; + for i in 0..3 { + for j in 0..3 { + let k = [[0.0, -a[2], a[1]], [a[2], 0.0, -a[0]], [-a[1], a[0], 0.0]][i][j]; + r[i][j] = if i == j { c } else { 0.0 } + (1.0 - c) * a[i] * a[j] + s * k; + } + } + 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() +} + +fn key(x: [f64; 3]) -> [i64; 3] { + x.map(|c| (c * 1e6).round() as i64) +} + +/// The displacements of a result by node position (rounded). +fn by_position(model: &Model, r: &StaticResult) -> BTreeMap<[i64; 3], ([f64; 3], f64)> { + r.nodes.iter().enumerate().map(|(i, l)| (key(model.nodes()[l].xyz), (r.displacement[i], r.von_mises[i]))).collect() +} + +/// A cantilever L × 2b × h along x (y from −b to b when `whole`, 0 to b otherwise) clamped at +/// x = 0, under a downward traction on its end; the half has the symmetry condition on y = 0. +/// Everything is turned by `turn` (nodes, the traction and the plane's normal). +fn cantilever(kind: ElementKind, formulation: Formulation, whole: bool, turn: &[[f64; 3]; 3]) -> (Model, StaticResult) { + let (l, b, h) = (10.0, 1.0, 1.0); + let n = if whole { [8, 4, 2] } else { [8, 2, 2] }; + let mut g = grid(n, kind, |p| apply(turn, [l * p[0], if whole { b * (2.0 * p[1] - 1.0) } else { b * p[1] }, h * p[2]])); + g.define(iso(210e3, 0.3), None, formulation, Procedure::Static {}); + g.fixed("root", |p| near(p[0], 0.0)); + g.face_set("tip", |p| near(p[0], 1.0)); + g.load("push", Load::Traction { set: "tip".into(), vector: apply(turn, [0.0, 0.0, -10.0]) }); + if !whole { + let normal = apply(turn, [0.0, 1.0, 0.0]); + g.symmetry("mid", Symmetry::from_normal(normal).unwrap(), |p| near(p[1], 0.0)); + } + let r = solve_static(&g.model, None, &Limits::default()).unwrap(); + (g.model, r) +} + +/// The turned half cantilever with its inclined plane solves to the turned displacements of the +/// unturned half (and, for the mirror-symmetric hexahedral meshes, of the whole), at every node; +/// the von Mises stress is the same; the plane's reactions balance the load. +#[test] +fn a_turned_half_solid_matches_the_axis_aligned_half_and_the_whole() { + let turn = rotation([1.0, 2.0, 3.0], 0.7); + let identity = rotation([0.0, 0.0, 1.0], 0.0); + for (kind, formulation, mirrored) in [ + (ElementKind::Hex8, full(), true), + (ElementKind::Hex20, full(), true), + (ElementKind::Hex20, reduced(), true), + (ElementKind::Tet10, full(), false), + (ElementKind::Wedge6, full(), false), + ] { + let (m_axis, axis) = cantilever(kind, formulation, false, &identity); + let (m_turned, turned) = cantilever(kind, formulation, false, &turn); + let size = axis.displacement.iter().fold(0.0f64, |m, u| m.max(norm(*u))); + let stress = axis.von_mises.iter().fold(0.0f64, |m, s| m.max(*s)); + let reference = by_position(&m_axis, &axis); + let mut worst = (0.0f64, 0.0f64); + let unturn = |x: [f64; 3]| [0, 1, 2].map(|i| (0..3).map(|j| turn[j][i] * x[j]).sum::()); + for (i, l) in turned.nodes.iter().enumerate() { + let (u, s) = reference[&key(unturn(m_turned.nodes()[l].xyz))]; + let want = apply(&turn, u); + let got = turned.displacement[i]; + worst.0 = worst.0.max(norm([0, 1, 2].map(|k| got[k] - want[k])) / size); + worst.1 = worst.1.max((turned.von_mises[i] - s).abs() / stress); + } + eprintln!( + "{kind:?} {formulation:?}: tip deflection {size:.6e}, turned vs axis-aligned: displacement {:.1e}, von Mises {:.1e}", + worst.0, worst.1 + ); + assert!(worst.0 < 1e-9 && worst.1 < 1e-8, "{kind:?}: {worst:?}"); + for r in [&axis, &turned] { + let imbalance = (0..3).map(|c| (r.applied_total[c] + r.reaction_total[c]).abs()).fold(0.0, f64::max); + assert!(imbalance < 1e-8 * norm(r.applied_total), "{kind:?}: {:?} vs {:?}", r.applied_total, r.reaction_total); + assert!(!r.warnings.iter().any(|w| w.contains("reactions")), "{:?}", r.warnings); + } + if mirrored { + let (m_whole, whole) = cantilever(kind, formulation, true, &identity); + let all = by_position(&m_whole, &whole); + let mut worst = 0.0f64; + for (x, (u, _)) in &reference { + let w = all[x].0; + worst = worst.max(norm([0, 1, 2].map(|k| u[k] - w[k])) / size); + } + eprintln!("{kind:?}: half with the plane vs whole: displacement {worst:.1e}"); + assert!(worst < 1e-9, "{kind:?}: {worst}"); + } + } +} + +/// A shell model: one material and step; `n` × `n` cells of the quarter of a simply supported +/// square plate (side 1), two triangles per cell split along x = y, so the mesh is mirror- +/// symmetric about the diagonal plane x = y; `eighth` keeps the cells below the diagonal and +/// puts the inclined symmetry plane (normal (1, −1, 0)) on it. +fn plate(n: usize, eighth: bool, procedure: Procedure) -> (Model, BTreeMap) { + let (a, t) = (1.0, 0.01); + let mut m = Model::default(); + m.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: 1e7, poisson: 0.3 }), density: Some(1.0), ..Material::default() }) + .unwrap(); + m.add_step(Step::new("Step", procedure), None).unwrap(); + let label = |i: usize, j: usize| (1 + i + (n + 1) * j) as u32; + let mut xyz = BTreeMap::new(); + for j in 0..=n { + for i in 0..=n { + if eighth && j > i { + continue; + } + let p = [i as f64 / n as f64 * a / 2.0, j as f64 / n as f64 * a / 2.0, 0.0]; + m.set_node(label(i, j), p).unwrap(); + xyz.insert(label(i, j), p); + } + } + let mut elements = Vec::new(); + for j in 0..n { + for i in 0..n { + let (p, q, r, s) = (label(i, j), label(i + 1, j), label(i + 1, j + 1), label(i, j + 1)); + for (tri, below) in [([p, q, r], true), ([p, r, s], false)] { + if eighth && (j > i || (i == j && !below)) { + continue; + } + let e = elements.len() as u32 + 1; + m.set_element(e, ElementKind::Tri3, tri.to_vec()).unwrap(); + elements.push(e); + } + } + } + let faces = elements.iter().map(|e| ElementFace { element: *e, face: 1 }).collect(); + m.add_set("plate", Set { elements, ..Set::default() }).unwrap(); + m.add_section( + "plate", + Section { + set: "plate".into(), + material: "Mat".into(), + kind: SectionKind::Shell { thickness: t, offset: 0.0, points: 5 }, + formulation: Formulation::default(), + }, + ) + .unwrap(); + let on = |x: f64, c: f64| (x - c).abs() < 1e-9; + let free = None; + let hold = |m: &mut Model, name: &str, bc: BoundaryCondition, f: &dyn Fn([f64; 3]) -> bool| { + let nodes = xyz.iter().filter(|(_, x)| f(**x)).map(|(l, _)| *l).collect(); + m.add_set(name, Set { nodes, ..Set::default() }).unwrap(); + m.add_boundary_condition(name, bc).unwrap(); + m.set_boundary_condition_active("Step", name, true).unwrap(); + }; + let held = |values| move |name: &str| BoundaryCondition::Displacement { target: Target::Set(name.into()), values }; + // Hard simple support on x = 0 and y = 0; symmetry on x = a/2 and y = a/2. + let x0 = held([free, free, Some(0.0), Some(0.0), free, free]); + let y0 = held([free, free, Some(0.0), free, Some(0.0), free]); + let sym = |name: &str, plane| BoundaryCondition::Symmetry { target: Target::Set(name.into()), plane }; + if !eighth { + hold(&mut m, "x0", x0("x0"), &|p| on(p[0], 0.0)); + hold(&mut m, "ys", sym("ys", Symmetry::Y), &|p| on(p[1], a / 2.0)); + } + hold(&mut m, "y0", y0("y0"), &|p| on(p[1], 0.0)); + hold(&mut m, "xs", sym("xs", Symmetry::X), &|p| on(p[0], a / 2.0)); + if eighth { + hold(&mut m, "diagonal", sym("diagonal", Symmetry::from_normal([1.0, -1.0, 0.0]).unwrap()), &|p| on(p[0], p[1])); + } + if matches!(m.steps()[0].procedure, Procedure::Static {}) { + m.add_set("top", Set { faces, ..Set::default() }).unwrap(); + m.add_load("q", Load::Pressure { set: "top".into(), magnitude: 1.0 }).unwrap(); + m.set_load_active("Step", "q", true).unwrap(); + } + (m, xyz) +} + +/// An eighth of the simply supported plate, cut on the diagonal with the inclined plane, solves +/// as the quarter (deflections and rotations at every shared node, to round-off; the first +/// frequency, whose mode is symmetric about the diagonal), and both match Navier's series. +#[test] +fn an_eighth_plate_on_the_diagonal_plane_matches_the_quarter() { + let n = 16; + let (mq, quarter) = plate(n, false, Procedure::Static {}); + let (me, eighth) = plate(n, true, Procedure::Static {}); + let rq = solve_static(&mq, None, &Limits::default()).unwrap(); + let re = solve_static(&me, None, &Limits::default()).unwrap(); + let at = |r: &StaticResult, l: u32| r.nodes.binary_search(&l).unwrap(); + let size = rq.displacement.iter().fold(0.0f64, |m, u| m.max(norm(*u))); + let turn = rq.rotation.iter().fold(0.0f64, |m, u| m.max(norm(*u))); + let (mut du, mut dr) = (0.0f64, 0.0f64); + for l in eighth.keys() { + assert!(quarter.contains_key(l)); + let (a, b) = (at(&rq, *l), at(&re, *l)); + du = du.max(norm([0, 1, 2].map(|k| rq.displacement[a][k] - re.displacement[b][k])) / size); + dr = dr.max(norm([0, 1, 2].map(|k| rq.rotation[a][k] - re.rotation[b][k])) / turn); + } + // The diagonal turns only about its normal: θx + θy = 0 and θz = 0 there. + for (l, x) in &eighth { + if (x[0] - x[1]).abs() < 1e-9 { + let r = re.rotation[at(&re, *l)]; + assert!((r[0] + r[1]).abs() < 1e-12 * turn && r[2].abs() < 1e-12 * turn, "node {l}: {r:?}"); + } + } + let centre = *eighth.iter().find(|(_, x)| (x[0] - 0.5).abs() < 1e-9 && (x[1] - 0.5).abs() < 1e-9).unwrap().0; + let w = -re.displacement[at(&re, centre)][2]; + let d = 1e7 * 0.01f64.powi(3) / (12.0 * (1.0 - 0.09)); + let navier = navier(1.0, 1.0, d); + eprintln!("eighth vs quarter plate: displacement {du:.1e}, rotation {dr:.1e}; centre deflection {w:.6e} vs Navier {navier:.6e}"); + assert!(du < 1e-9 && dr < 1e-9, "{du} {dr}"); + assert!((w - navier).abs() < 0.01 * navier, "{w} vs {navier}"); + let imbalance = (0..3).map(|c| (re.applied_total[c] + re.reaction_total[c]).abs()).fold(0.0, f64::max); + assert!(imbalance < 1e-9 * norm(re.applied_total), "{:?} vs {:?}", re.applied_total, re.reaction_total); + assert!(!re.warnings.iter().any(|w| w.contains("reactions")), "{:?}", re.warnings); + // An eighth carries an eighth of the plate's load (the quarter a quarter). + assert!((2.0 * re.applied_total[2] - rq.applied_total[2]).abs() < 1e-12 * rq.applied_total[2].abs()); + + let frequency = |eighth| { + let (m, _) = plate(n, eighth, Procedure::Frequency { modes: 1, min_hz: None, max_hz: None }); + solve_frequency(&m, None, &FrequencyOptions::default(), &Limits::default()).unwrap().modes[0].frequency + }; + let (fq, fe) = (frequency(false), frequency(true)); + eprintln!("first frequency: quarter {fq:.9}, eighth {fe:.9}"); + assert!((fq - fe).abs() < 1e-8 * fq, "{fq} vs {fe}"); +} + +/// Navier's series for the centre deflection of a simply supported square plate under uniform +/// pressure (S. Timoshenko, S. Woinowsky-Krieger, *Theory of Plates and Shells*, 2nd ed., +/// McGraw-Hill 1959, §30). +fn navier(q: f64, a: f64, d: f64) -> f64 { + let mut s = 0.0; + for m in (1..400).step_by(2) { + for n in (1..400).step_by(2) { + let (mf, nf) = (f64::from(m), f64::from(n)); + let sign = if ((m + n) / 2 - 1) % 2 == 0 { 1.0 } else { -1.0 }; + s += sign / (mf * nf * (mf * mf + nf * nf).powi(2)); + } + } + 16.0 * q * a.powi(4) / (std::f64::consts::PI.powi(6) * d) * s +} + +/// A node that is also held, or listed in two planes, is no conflict; a prescribed motion off +/// the plane is. +#[test] +fn conflicts() { + // The clamp's nodes on y = 0 are also on the plane: held there, the plane's rows are redundant. + let (_, r) = cantilever(ElementKind::Hex8, full(), false, &rotation([0.0, 0.0, 1.0], 0.3)); + assert!(!r.warnings.iter().any(|w| w.contains("reactions")), "{:?}", r.warnings); + let mut g = grid([4, 1, 1], ElementKind::Hex8, |p| [4.0 * p[0], p[1], p[2]]); + g.define(iso(1e3, 0.3), None, full(), Procedure::Static {}); + g.fixed("root", |p| near(p[0], 0.0)); + g.symmetry("a", Symmetry::Normal([1.0, 1.0, 0.0]), |p| near(p[1], 0.0) && near(p[0], 1.0)); + g.symmetry("b", Symmetry::Normal([2.0, 2.0, 0.0]), |p| near(p[1], 0.0) && near(p[0], 1.0)); + g.face_set("tip", |p| near(p[0], 1.0)); + g.load("push", Load::Traction { set: "tip".into(), vector: [0.0, 0.0, -1.0] }); + let r = solve_static(&g.model, None, &Limits::default()).unwrap(); + let u = r.displacement[r.nodes.binary_search(&g.nodes_where(|p| near(p[1], 0.0) && near(p[0], 1.0))[0]).unwrap()]; + assert!((u[0] + u[1]).abs() < 1e-12 * norm(u), "{u:?}"); + // A prescribed motion the plane cannot follow is a contradiction, reported by name. + let nodes = g.nodes_where(|p| near(p[1], 0.0) && near(p[0], 1.0)); + g.model.add_set("pulled", Set { nodes, ..Set::default() }).unwrap(); + g.bc("pull", BoundaryCondition::Displacement { target: Target::Set("pulled".into()), values: [Some(0.1), Some(0.1), None, None, None, None] }); + let e = solve_static(&g.model, None, &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("boundary condition `") && e.contains("contradicts"), "{e}"); +}